From 6c17d27247b38ac8e07dbd8ebcc5824fd34b0a77 Mon Sep 17 00:00:00 2001 From: StuartBenjamin Date: Thu, 22 Oct 2026 12:44:14 -0400 Subject: [PATCH 1/2] STRIDE - ENHANCEMENT - dprime_banded option: banded LU solve of the Delta' BVP The Delta' BVP matrix is block-banded except for the 2*msing big solution rows, which the dense solve puts last. With dprime_banded=t each surface keeps its two big solution rows in its block, so the matrix is banded with kl = ku = 3*mpert-1 and is solved with ZGBTRF/ZGBTRS (residual check from the sparse entries). Default f keeps the dense solve unchanged. Co-Authored-By: Claude Opus 5.5 --- input/stride.in | 1 + stride/ode.F | 76 +++++++++++++++++++++++++++++++++++-------------- stride/stride.F | 2 +- 3 files changed, 57 insertions(+), 22 deletions(-) diff --git a/input/stride.in b/input/stride.in index ff979efd7..fdbb3507f 100755 --- a/input/stride.in +++ b/input/stride.in @@ -61,6 +61,7 @@ calc_delta_prime=t ! When false, sing_order=2 is better for speed, (with less accuracy). calc_dp_with_vac=t ! When vac_flag=t, this will use vacuum edge conditions to calculate delta prime. big_soln_err_tol=1e-7 ! Threshold error for correct solution to the delta prime bvp. + dprime_banded=f ! Solve the delta prime bvp with a banded LU instead of a dense one: same result, much faster and smaller with many rational surfaces. integrate_riccati=f ! Use Riccati integration for dW. Cannot calculate delta_prime. riccati_bounce=t ! Toggle integration mode to control growth of matrix diff --git a/stride/ode.F b/stride/ode.F index 73bd2df80..5e5dca911 100644 --- a/stride/ode.F +++ b/stride/ode.F @@ -60,6 +60,7 @@ MODULE stride_ode_mod LOGICAL :: calc_delta_prime LOGICAL :: kill_big_soln_for_ideal_dW LOGICAL :: calc_dp_with_vac + LOGICAL :: dprime_banded=.FALSE. TYPE(sing_calculator), DIMENSION(:), TARGET, ALLOCATABLE :: scalc REAL(r8) :: axisPsi, outerPsi @@ -1046,8 +1047,8 @@ 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, rs, rb, kl, ku + REAL(r8) :: rcond, mnorm, ZLANGE, ZLANGB INTEGER, DIMENSION(:), ALLOCATABLE :: isp, jsp COMPLEX(r8), DIMENSION(:), ALLOCATABLE :: asp, b, x COMPLEX(r8), DIMENSION(:,:), ALLOCATABLE :: M, MInv, TempMat @@ -1077,6 +1078,11 @@ 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 per singular surface: the banded solve keeps each + !surface's two big solution rows in its block (the dense + !solve puts them last), so kl = ku = 3*mpert-1 + rs = 2*m2-2 + IF (dprime_banded) rs = 2*m2 ALLOCATE(b(nMat),x(nMat),isp(nSp),jsp(nSp),asp(nSp)) ALLOCATE(TempMat(m2,m2),ipivTemp(m2),rworkTemp(2*m2), @@ -1108,7 +1114,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 + rs * (jsing-1) jsp(k) = j + mpert + (2*m2) * (jsing-1) asp(k) = s%uShootL(i,j) ENDDO @@ -1120,7 +1126,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 + rs * (jsing-1) jsp(k) = i + mpert + (2*m2) * (jsing-1) asp(k) = 1.0 @@ -1138,7 +1144,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 + rs * jsing jsp(k) = j + mpert + m2 + (2*m2)*(jsing-1) asp(k) = -s%uShootR(i,j) ENDDO @@ -1150,7 +1156,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 + rs * msing jsp(k) = j + mpert + (2*m2) * msing asp(k) = uEdge(i,j+mpert) ENDDO @@ -1161,14 +1167,16 @@ SUBROUTINE ode_calc_delta_prime DO jsing = 1,msing s => scalc(jsing) ipert0 = s%sing_col + rb = nMat-s2+2*(jsing-1) + IF (dprime_banded) rb = rs*jsing-2 !left side k=k+1 - isp(k) = isp(k-1)+1 + isp(k) = rb+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) = rb+2 jsp(k) = ipert0 + mpert + m2 + (2*m2)*(jsing-1) asp(k) = 1.0 ENDDO @@ -1181,16 +1189,30 @@ SUBROUTINE ode_calc_delta_prime c----------------------------------------------------------------------- c check condition number of BVP matrix. c----------------------------------------------------------------------- - ALLOCATE(M(nMat,nMat),MInv(nMat,nMat)) ALLOCATE(ipiv(nMat),uwork(nMat),cwork(2*nMat),rwork(2*nMat)) - M=0.0 - DO k = 1,nSp - M(isp(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) + IF (dprime_banded) THEN + kl = MAXVAL(isp-jsp) + ku = MAXVAL(jsp-isp) + ALLOCATE(MInv(2*kl+ku+1,nMat)) + MInv=0.0 + DO k = 1,nSp + MInv(kl+ku+1+isp(k)-jsp(k),jsp(k)) = asp(k) + ENDDO + 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) + ELSE + ALLOCATE(M(nMat,nMat),MInv(nMat,nMat)) + M=0.0 + DO k = 1,nSp + M(isp(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) + ENDIF print *,"Delta' sparse matrix condition number = ",rcond c----------------------------------------------------------------------- c solve BVP matrix without sparse matrix. @@ -1201,15 +1223,19 @@ SUBROUTINE ode_calc_delta_prime b = 0.0 IF (i == 1) THEN !Left side driving term - b(nMat-s2+2*jsing-1) = 1.0 dRow = 2*jsing-1 ELSE !Right side driving term - b(nMat-s2+2*jsing) = 1.0 dRow = 2*jsing ENDIF + b(isp(nSp-s2+dRow)) = 1.0 !its big solution row x = b - CALL ZGETRS('N',nMat,1,MInv,nMat,ipiv,x,nMat,info) + IF (dprime_banded) THEN + CALL ZGBTRS('N',nMat,kl,ku,1,MInv,2*kl+ku+1,ipiv,x, + $ nMat,info) + ELSE + CALL ZGETRS('N',nMat,1,MInv,nMat,ipiv,x,nMat,info) + ENDIF IF (info /= 0) THEN print *,"info was not zero! info=",info ENDIF @@ -1229,7 +1255,15 @@ SUBROUTINE ode_calc_delta_prime c error-check non-sparse matrix solution. c----------------------------------------------------------------------- !Error-check big solution elements by backward solve. - b = ABS(MATMUL(M,x)) + IF (dprime_banded) THEN + b = 0.0 + DO k = 1,nSp + b(isp(k)) = b(isp(k)) + asp(k)*x(jsp(k)) + ENDDO + b = ABS(b) + ELSE + b = ABS(MATMUL(M,x)) + ENDIF 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))) diff --git a/stride/stride.F b/stride/stride.F index 585058ce7..76bc43e61 100644 --- a/stride/stride.F +++ b/stride/stride.F @@ -55,7 +55,7 @@ PROGRAM stride NAMELIST/stride_params/grid_packing,asymp_at_sing, $ integrate_riccati,calc_delta_prime,calc_dp_with_vac, $ axis_mid_pt_skew, - $ big_soln_err_tol, kill_big_soln_for_ideal_dW, + $ big_soln_err_tol, kill_big_soln_for_ideal_dW, dprime_banded, $ ric_dt,ric_tol,riccati_bounce,verbose_riccati_output, $ riccati_match_hamiltonian_evals,verbose_performance_output, $ fourfit_metric_parallel,vac_parallel,nIntervalsTot,nThreads From c125d36e847eddae893d52834bdb108b9aa04331 Mon Sep 17 00:00:00 2001 From: StuartBenjamin Date: Thu, 22 Oct 2026 22:50:33 -0400 Subject: [PATCH 2/2] STRIDE - ENHANCEMENT - Always solve the Delta' BVP banded; drop the dprime_banded option The banded LU gives bit-identical Delta' to the dense one, so it replaces it: the dprime_banded namelist key, the dense path and the stride.in line go. Each singular surface's two big solution rows sit at the end of its block of 2*m2 rows, the matrix is factored and solved with ZGBTRF/ZGBCON/ZGBTRS (kl = ku = 3*mpert-1), and the residual check uses the sparse entries. The printed condition number now uses the norm of the matrix itself, as ZGBCON expects (the dense path took the norm after factorization). Co-Authored-By: Claude Opus 5.5 --- input/stride.in | 1 - stride/ode.F | 90 ++++++++++++++++++------------------------------- stride/stride.F | 2 +- 3 files changed, 34 insertions(+), 59 deletions(-) diff --git a/input/stride.in b/input/stride.in index fdbb3507f..ff979efd7 100755 --- a/input/stride.in +++ b/input/stride.in @@ -61,7 +61,6 @@ calc_delta_prime=t ! When false, sing_order=2 is better for speed, (with less accuracy). calc_dp_with_vac=t ! When vac_flag=t, this will use vacuum edge conditions to calculate delta prime. big_soln_err_tol=1e-7 ! Threshold error for correct solution to the delta prime bvp. - dprime_banded=f ! Solve the delta prime bvp with a banded LU instead of a dense one: same result, much faster and smaller with many rational surfaces. integrate_riccati=f ! Use Riccati integration for dW. Cannot calculate delta_prime. riccati_bounce=t ! Toggle integration mode to control growth of matrix diff --git a/stride/ode.F b/stride/ode.F index 5e5dca911..111a82357 100644 --- a/stride/ode.F +++ b/stride/ode.F @@ -60,7 +60,6 @@ MODULE stride_ode_mod LOGICAL :: calc_delta_prime LOGICAL :: kill_big_soln_for_ideal_dW LOGICAL :: calc_dp_with_vac - LOGICAL :: dprime_banded=.FALSE. TYPE(sing_calculator), DIMENSION(:), TARGET, ALLOCATABLE :: scalc REAL(r8) :: axisPsi, outerPsi @@ -1047,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, rs, rb, kl, ku - REAL(r8) :: rcond, mnorm, ZLANGE, ZLANGB + 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 @@ -1078,11 +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 per singular surface: the banded solve keeps each - !surface's two big solution rows in its block (the dense - !solve puts them last), so kl = ku = 3*mpert-1 - rs = 2*m2-2 - IF (dprime_banded) rs = 2*m2 + !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), @@ -1114,7 +1111,7 @@ SUBROUTINE ode_calc_delta_prime DO i = 1,m2 DO j = 1,m2 k=k+1 - isp(k) = i + rs * (jsing-1) + isp(k) = i + (2*m2) * (jsing-1) jsp(k) = j + mpert + (2*m2) * (jsing-1) asp(k) = s%uShootL(i,j) ENDDO @@ -1126,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 + rs * (jsing-1) + isp(k) = i + iChg + m2 + (2*m2) * (jsing-1) jsp(k) = i + mpert + (2*m2) * (jsing-1) asp(k) = 1.0 @@ -1144,7 +1141,7 @@ SUBROUTINE ode_calc_delta_prime DO i = 1,m2 DO j = 1,m2 k=k+1 - isp(k) = i + rs * jsing + isp(k) = i + (2*m2) * jsing jsp(k) = j + mpert + m2 + (2*m2)*(jsing-1) asp(k) = -s%uShootR(i,j) ENDDO @@ -1156,7 +1153,7 @@ SUBROUTINE ode_calc_delta_prime DO i = 1,m2 DO j = 1,mpert k=k+1 - isp(k) = i + rs * msing + isp(k) = i + (2*m2) * msing jsp(k) = j + mpert + (2*m2) * msing asp(k) = uEdge(i,j+mpert) ENDDO @@ -1167,16 +1164,14 @@ SUBROUTINE ode_calc_delta_prime DO jsing = 1,msing s => scalc(jsing) ipert0 = s%sing_col - rb = nMat-s2+2*(jsing-1) - IF (dprime_banded) rb = rs*jsing-2 !left side k=k+1 - isp(k) = rb+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) = rb+2 + isp(k) = 2*m2*jsing jsp(k) = ipert0 + mpert + m2 + (2*m2)*(jsing-1) asp(k) = 1.0 ENDDO @@ -1189,30 +1184,18 @@ SUBROUTINE ode_calc_delta_prime c----------------------------------------------------------------------- c check condition number of BVP matrix. c----------------------------------------------------------------------- + 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)) - IF (dprime_banded) THEN - kl = MAXVAL(isp-jsp) - ku = MAXVAL(jsp-isp) - ALLOCATE(MInv(2*kl+ku+1,nMat)) - MInv=0.0 - DO k = 1,nSp - MInv(kl+ku+1+isp(k)-jsp(k),jsp(k)) = asp(k) - ENDDO - 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) - ELSE - ALLOCATE(M(nMat,nMat),MInv(nMat,nMat)) - M=0.0 - DO k = 1,nSp - M(isp(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) - ENDIF + MInv=0.0 + DO k = 1,nSp + MInv(kl+ku+1+isp(k)-jsp(k),jsp(k)) = asp(k) + ENDDO + 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. @@ -1223,19 +1206,16 @@ SUBROUTINE ode_calc_delta_prime b = 0.0 IF (i == 1) THEN !Left side driving term + b(2*m2*jsing-1) = 1.0 dRow = 2*jsing-1 ELSE !Right side driving term + b(2*m2*jsing) = 1.0 dRow = 2*jsing ENDIF - b(isp(nSp-s2+dRow)) = 1.0 !its big solution row x = b - IF (dprime_banded) THEN - CALL ZGBTRS('N',nMat,kl,ku,1,MInv,2*kl+ku+1,ipiv,x, - $ nMat,info) - ELSE - CALL ZGETRS('N',nMat,1,MInv,nMat,ipiv,x,nMat,info) - ENDIF + 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 @@ -1252,18 +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. - IF (dprime_banded) THEN - b = 0.0 - DO k = 1,nSp - b(isp(k)) = b(isp(k)) + asp(k)*x(jsp(k)) - ENDDO - b = ABS(b) - ELSE - b = ABS(MATMUL(M,x)) - ENDIF + 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))) diff --git a/stride/stride.F b/stride/stride.F index 76bc43e61..585058ce7 100644 --- a/stride/stride.F +++ b/stride/stride.F @@ -55,7 +55,7 @@ PROGRAM stride NAMELIST/stride_params/grid_packing,asymp_at_sing, $ integrate_riccati,calc_delta_prime,calc_dp_with_vac, $ axis_mid_pt_skew, - $ big_soln_err_tol, kill_big_soln_for_ideal_dW, dprime_banded, + $ big_soln_err_tol, kill_big_soln_for_ideal_dW, $ ric_dt,ric_tol,riccati_bounce,verbose_riccati_output, $ riccati_match_hamiltonian_evals,verbose_performance_output, $ fourfit_metric_parallel,vac_parallel,nIntervalsTot,nThreads