Skip to content

STRIDE: banded LU for the Δ′ boundary value problem #351

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
matt-pharr merged 2 commits into from
Oct 23, 2026
Merged

STRIDE: banded LU for the Δ′ boundary value problem #351

Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
6c17d27
STRIDE - ENHANCEMENT - dprime_banded option: banded LU solve of the D…
StuartBenjamin Oct 22, 2026
c125d36
STRIDE - ENHANCEMENT - Always solve the Delta' BVP banded; drop the d…
StuartBenjamin Oct 23, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
52 changes: 31 additions & 21 deletions stride/ode.F
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1046,11 +1046,11 @@ SUBROUTINE ode_calc_delta_prime
REAL(r8), DIMENSION(:,:), ALLOCATABLE :: atol

INTEGER :: nMat, nSp, i, j, k, jsing, ksing, s2, m2
INTEGER :: iChg, info, ipert0, dRow
REAL(r8) :: rcond, mnorm, ZLANGE
INTEGER :: iChg, info, ipert0, dRow, kl, ku
REAL(r8) :: rcond, mnorm, ZLANGB
INTEGER, DIMENSION(:), ALLOCATABLE :: isp, jsp
COMPLEX(r8), DIMENSION(:), ALLOCATABLE :: asp, b, x
COMPLEX(r8), DIMENSION(:,:), ALLOCATABLE :: M, MInv, TempMat
COMPLEX(r8), DIMENSION(:,:), ALLOCATABLE :: MInv, TempMat
INTEGER, DIMENSION(:), ALLOCATABLE :: ipiv, ipivTemp
COMPLEX(r8), DIMENSION(:), ALLOCATABLE :: uwork, cwork,
$ cworkTemp
Expand All @@ -1077,6 +1077,9 @@ SUBROUTINE ode_calc_delta_prime
$ s2 * (m2-2) + !each side of each sing
$ (msing-1) * m2*m2*2 + !intervals sandwiched b/w sings
$ s2 !driven big solution coefficients (2sides * msing)
!rows come in blocks of 2*m2 per singular surface, each ending
!with that surface's two big solution rows, so the matrix is
!banded with kl = ku = 3*mpert-1

ALLOCATE(b(nMat),x(nMat),isp(nSp),jsp(nSp),asp(nSp))
ALLOCATE(TempMat(m2,m2),ipivTemp(m2),rworkTemp(2*m2),
Expand Down Expand Up @@ -1108,7 +1111,7 @@ SUBROUTINE ode_calc_delta_prime
DO i = 1,m2
DO j = 1,m2
k=k+1
isp(k) = i + (2*m2-2) * (jsing-1)
isp(k) = i + (2*m2) * (jsing-1)
jsp(k) = j + mpert + (2*m2) * (jsing-1)
asp(k) = s%uShootL(i,j)
ENDDO
Expand All @@ -1120,7 +1123,7 @@ SUBROUTINE ode_calc_delta_prime
DO i = 1,m2
IF (i/=ipert0 .AND. i/=ipert0+mpert) THEN
k=k+1
isp(k) = i + iChg + m2 + (2*m2-2) * (jsing-1)
isp(k) = i + iChg + m2 + (2*m2) * (jsing-1)
jsp(k) = i + mpert + (2*m2) * (jsing-1)
asp(k) = 1.0

Expand All @@ -1138,7 +1141,7 @@ SUBROUTINE ode_calc_delta_prime
DO i = 1,m2
DO j = 1,m2
k=k+1
isp(k) = i + (2*m2-2) * jsing
isp(k) = i + (2*m2) * jsing
jsp(k) = j + mpert + m2 + (2*m2)*(jsing-1)
asp(k) = -s%uShootR(i,j)
ENDDO
Expand All @@ -1150,7 +1153,7 @@ SUBROUTINE ode_calc_delta_prime
DO i = 1,m2
DO j = 1,mpert
k=k+1
isp(k) = i + (2*m2-2) * msing
isp(k) = i + (2*m2) * msing
jsp(k) = j + mpert + (2*m2) * msing
asp(k) = uEdge(i,j+mpert)
ENDDO
Expand All @@ -1163,12 +1166,12 @@ SUBROUTINE ode_calc_delta_prime
ipert0 = s%sing_col
!left side
k=k+1
isp(k) = isp(k-1)+1
isp(k) = 2*m2*jsing-1
jsp(k) = ipert0 + mpert + (2*m2)*(jsing-1)
asp(k) = 1.0
!right side
k=k+1
isp(k) = isp(k-1)+1
isp(k) = 2*m2*jsing
jsp(k) = ipert0 + mpert + m2 + (2*m2)*(jsing-1)
asp(k) = 1.0
ENDDO
Expand All @@ -1181,16 +1184,18 @@ SUBROUTINE ode_calc_delta_prime
c-----------------------------------------------------------------------
c check condition number of BVP matrix.
c-----------------------------------------------------------------------
ALLOCATE(M(nMat,nMat),MInv(nMat,nMat))
kl = MAXVAL(isp-jsp)
ku = MAXVAL(jsp-isp)
ALLOCATE(MInv(2*kl+ku+1,nMat))
ALLOCATE(ipiv(nMat),uwork(nMat),cwork(2*nMat),rwork(2*nMat))
M=0.0
MInv=0.0
DO k = 1,nSp
M(isp(k),jsp(k)) = asp(k)
MInv(kl+ku+1+isp(k)-jsp(k),jsp(k)) = asp(k)
ENDDO
MInv = M
CALL ZGETRF(nMat,nMat,MInv,nMat,ipiv,info)
mnorm = ZLANGE('O',nMat,nMat,MInv,nMat,rwork)
CALL ZGECON('O',nMat,MInv,nMat,mnorm,rcond,cwork,rwork,info)
mnorm = ZLANGB('O',nMat,kl,ku,MInv(kl+1,1),2*kl+ku+1,rwork)
CALL ZGBTRF(nMat,nMat,kl,ku,MInv,2*kl+ku+1,ipiv,info)
CALL ZGBCON('O',nMat,kl,ku,MInv,2*kl+ku+1,ipiv,mnorm,rcond,
$ cwork,rwork,info)
print *,"Delta' sparse matrix condition number = ",rcond
c-----------------------------------------------------------------------
c solve BVP matrix without sparse matrix.
Expand All @@ -1201,15 +1206,16 @@ SUBROUTINE ode_calc_delta_prime
b = 0.0
IF (i == 1) THEN
!Left side driving term
b(nMat-s2+2*jsing-1) = 1.0
b(2*m2*jsing-1) = 1.0
dRow = 2*jsing-1
ELSE
!Right side driving term
b(nMat-s2+2*jsing) = 1.0
b(2*m2*jsing) = 1.0
dRow = 2*jsing
ENDIF
x = b
CALL ZGETRS('N',nMat,1,MInv,nMat,ipiv,x,nMat,info)
CALL ZGBTRS('N',nMat,kl,ku,1,MInv,2*kl+ku+1,ipiv,x,nMat,
$ info)
IF (info /= 0) THEN
print *,"info was not zero! info=",info
ENDIF
Expand All @@ -1226,10 +1232,14 @@ SUBROUTINE ode_calc_delta_prime
delta_prime_mat(dRow,2*ksing) = x(k+ipert0+mpert)
ENDDO
c-----------------------------------------------------------------------
c error-check non-sparse matrix solution.
c error-check the solution with the sparse matrix.
c-----------------------------------------------------------------------
!Error-check big solution elements by backward solve.
b = ABS(MATMUL(M,x))
b = 0.0
DO k = 1,nSp
b(isp(k)) = b(isp(k)) + asp(k)*x(jsp(k))
ENDDO
b = ABS(b)
IF (ABS(1-SUM(ABS(b))) > big_soln_err_tol) THEN
print *,"Big sol'n error, surface #",jsing,"side",i
print *,"=",ABS(1-SUM(ABS(b)))
Expand Down
Loading