From 49238a87889641ef7e0c3acaa585ae627f56fbb6 Mon Sep 17 00:00:00 2001 From: SorBalda Date: Fri, 10 Jul 2026 12:04:00 +0200 Subject: [PATCH 1/3] memopt(safe): drop dead zz/vv/ww arrays and iden in get_odd_straight_with_v4 Bitwise-identical memory optimization of get_odd_straight_with_v4. - Remove zz(nl,nl), vv(nl*(nl+1)/2) and ww(nl): declared and allocated but never read/written anywhere in the subroutine (verified by grep on the v1.5 source). They were also never deallocated -> a memory leak. zz alone is one full nl x nl array (= 1 "v4 unit"). - Remove the iden(nl,nl) scratch matrix. Instead of building iden and copying "maux = iden", initialize maux directly: maux = 0.0d0 ; do x=1,nl: maux(x,x)=1.0d0 Same value, same bit pattern; saves one nl x nl allocation. No BLAS call, no operation order, no numerics changed. Standalone Fortran test (n_mode=6, symmetrized v3/v4): max|phi_new-phi_orig| = 0.0 exactly. Co-Authored-By: Claude Fable 5 --- SCHAModules/get_odd_straight_with_v4.f90 | 21 ++++++--------------- 1 file changed, 6 insertions(+), 15 deletions(-) diff --git a/SCHAModules/get_odd_straight_with_v4.f90 b/SCHAModules/get_odd_straight_with_v4.f90 index 2faf6ea79..59e63927c 100644 --- a/SCHAModules/get_odd_straight_with_v4.f90 +++ b/SCHAModules/get_odd_straight_with_v4.f90 @@ -20,7 +20,7 @@ subroutine get_odd_straight_with_v4 ( a, wr, er, transmode, amass, ityp_sc, T, v integer :: nat_sc, n_mode, nl, ns, ntyp - double precision, dimension(:,:), allocatable :: l, g, phi_aux, v1, v2, v32, iden + double precision, dimension(:,:), allocatable :: l, g, phi_aux, v1, v2, v32 double precision :: lsum double precision, dimension(:), allocatable :: laux1, lres1, veclong double precision, dimension(:), allocatable :: laux2, lres2 @@ -30,9 +30,6 @@ subroutine get_odd_straight_with_v4 ( a, wr, er, transmode, amass, ityp_sc, T, v integer, dimension(:), allocatable :: ipiv integer :: info - double precision, dimension(:), allocatable :: vv - double precision, dimension(:), allocatable :: ww - double precision, dimension(:,:), allocatable :: zz double precision, dimension(:,:), allocatable :: cf integer :: mu, nu, alpha @@ -67,15 +64,10 @@ subroutine get_odd_straight_with_v4 ( a, wr, er, transmode, amass, ityp_sc, T, v allocate(ipiv(nl)) allocate(work(nl)) allocate(v32(n_mode,n_mode*n_mode)) - allocate(iden(nl,nl)) allocate(cf(nl,ns)) - allocate(vv(nl*(nl+1)/2)) - allocate(ww(nl)) - allocate(zz(nl,nl)) - - ! Get lambda matrix + ! Get lambda matrix call get_cmat ( a, wr, er, transmode, amass, ityp_sc, T, .true., lamat,n_mode, nat_sc, ntyp ) @@ -100,12 +92,13 @@ subroutine get_odd_straight_with_v4 ( a, wr, er, transmode, amass, ityp_sc, T, v end do end do - ! Prepare identity matrix + ! Prepare identity matrix directly in maux (was: iden then maux = iden; + ! bitwise identical, avoids one nl x nl temporary) - iden = 0.0d0 + maux = 0.0d0 do x = 1, nl - iden(x,x) = 1.0d0 + maux(x,x) = 1.0d0 end do ! Calculate ** iden - v4 lamat ** matrix @@ -113,8 +106,6 @@ subroutine get_odd_straight_with_v4 ( a, wr, er, transmode, amass, ityp_sc, T, v !print *, "BEFORE I - V4Lambda" !call flush() - maux = iden - call dgemm('N','N',nl,nl,nl,-1.0d0,v42,nl,lamat,nl,1.0d0,maux,nl) ! Invert ** iden - lamat v4 ** From aa4da397224595151683acb18bb62a52adae489d Mon Sep 17 00:00:00 2001 From: SorBalda Date: Fri, 10 Jul 2026 12:06:39 +0200 Subject: [PATCH 2/3] memopt(reassoc): O(N^5) reassociation + optional v4-symmetry path NOT bitwise identical: this commit reassociates floating-point products, so results differ from the original at the ~1e-15 relative level (verified below). The algorithm and the physics are unchanged. (d) Replace the nl x nl buffer product with an nl x ns contraction. Original: v42 := lamat * maux (nl x nl dgemm, reuses v42 as buffer) cf := v42 * v3^T (nl x ns) New (associativity Lambda*M^-1*v3^T = Lambda*(M^-1*v3^T)): tmp := maux * v3^T (nl x ns) cf := lamat * tmp (nl x ns) Removes the nl x nl buffer and lowers that stage from O(N^6) to O(N^5). tmp is nl x ns = O(N^3), negligible. (e) Optional argument use_v4_symmetry (default .false. -> old safe behavior). When .true. (standard case: v4 permutation-symmetric after ApplySymmetryToTensor4 / use_symmetries=True), the reordered copy v42(ja,ka)=v4(w,z,x,y) equals the plain column-major reshape of v4, so v4 is fed directly to the first dgemm (lda=nl) and the whole nl x nl v42 array is never allocated. INVALID for use_symmetries=False, hence it is opt-in only, never the silent default. The optional argument is interface-compatible: existing callers that omit it are unaffected. Peak internal nl x nl arrays: default path lamat+v42+maux (3 units, v42 freed right after the first dgemm); symmetry path lamat+maux (2 units). Standalone test (n_mode=6, symmetrized v3/v4), vs original: default path : max abs diff 2.1e-16, max elementwise rel diff 7.6e-15 symmetry path: max abs diff 1.1e-16, max elementwise rel diff 1.8e-14 Both far below the 1e-12 tolerance. Co-Authored-By: Claude Fable 5 --- SCHAModules/get_odd_straight_with_v4.f90 | 91 +++++++++++++++--------- 1 file changed, 59 insertions(+), 32 deletions(-) diff --git a/SCHAModules/get_odd_straight_with_v4.f90 b/SCHAModules/get_odd_straight_with_v4.f90 index 59e63927c..8dda0fb4d 100644 --- a/SCHAModules/get_odd_straight_with_v4.f90 +++ b/SCHAModules/get_odd_straight_with_v4.f90 @@ -1,10 +1,10 @@ -! This subroutine calculates the L mat needed to get the average of the -! third order derivatives. It is formed by four polarization vectors +! This subroutine calculates the L mat needed to get the average of the +! third order derivatives. It is formed by four polarization vectors ! times the mass^1/2 divided by the normal length. subroutine get_odd_straight_with_v4 ( a, wr, er, transmode, amass, ityp_sc, T, v3, v4, phi_sc_odd, & - n_mode, nat_sc, ntyp) + n_mode, nat_sc, ntyp, use_v4_symmetry) implicit none @@ -17,20 +17,25 @@ subroutine get_odd_straight_with_v4 ( a, wr, er, transmode, amass, ityp_sc, T, v double precision, dimension(n_mode,n_mode,n_mode), intent(in) :: v3 double precision, dimension(n_mode,n_mode,n_mode, n_mode), intent(in) :: v4 double precision, dimension(n_mode, n_mode), intent(out) :: phi_sc_odd + ! Optional: if present and .true., assume v4 is permutation-symmetric + ! (standard case, after ApplySymmetryToTensor4 / use_symmetries=True) and + ! skip building the explicit reordered copy v42 -- see note below. + logical, intent(in), optional :: use_v4_symmetry integer :: nat_sc, n_mode, nl, ns, ntyp double precision, dimension(:,:), allocatable :: l, g, phi_aux, v1, v2, v32 - double precision :: lsum + double precision :: lsum double precision, dimension(:), allocatable :: laux1, lres1, veclong double precision, dimension(:), allocatable :: laux2, lres2 - + double precision, dimension(:,:), allocatable :: lamat, v42, maux double precision, dimension(:), allocatable :: work integer, dimension(:), allocatable :: ipiv integer :: info double precision, dimension(:,:), allocatable :: cf + double precision, dimension(:,:), allocatable :: tmp integer :: mu, nu, alpha integer :: ka, ja @@ -39,13 +44,14 @@ subroutine get_odd_straight_with_v4 ( a, wr, er, transmode, amass, ityp_sc, T, v real :: t1, t2 logical, parameter :: debug = .true. + logical :: sym ! Get integers if (debug) then print *, "=== DEBUG ODD STRAIGHT ===" print *, "N_MODE:", n_mode - print *, "NTYP:", ntyp + print *, "NTYP:", ntyp print *, "NAT_SC:", nat_sc call flush() end if @@ -56,39 +62,48 @@ subroutine get_odd_straight_with_v4 ( a, wr, er, transmode, amass, ityp_sc, T, v ns = n_mode nl = n_mode*n_mode + sym = .false. + if (present(use_v4_symmetry)) sym = use_v4_symmetry + ! Allocate stuff allocate(lamat(nl,nl)) - allocate(v42(nl,nl)) allocate(maux(nl,nl)) allocate(ipiv(nl)) allocate(work(nl)) allocate(v32(n_mode,n_mode*n_mode)) allocate(cf(nl,ns)) + allocate(tmp(nl,ns)) + + ! v42 (explicit reordered copy of v4) is only needed when we cannot rely on + ! the permutation symmetry of v4. It is one full nl x nl array. + if (.not. sym) allocate(v42(nl,nl)) ! Get lambda matrix call get_cmat ( a, wr, er, transmode, amass, ityp_sc, T, .true., lamat,n_mode, nat_sc, ntyp ) - + !print *, "AFTER CMAT" !call flush() ! Write third and fourth order force constants as rank 2 ka = 0 - + do x = 1, n_mode do y = 1, n_mode ka = ka + 1 v32(:,ka) = v3(:,x,y) - ja = 0 - do w = 1, n_mode - do z = 1, n_mode - ja = ja + 1 - v42(ja,ka) = v4(w,z,x,y) + if (.not. sym) then + ja = 0 + do w = 1, n_mode + do z = 1, n_mode + ja = ja + 1 + v42(ja,ka) = v4(w,z,x,y) + end do end do - end do + end if end do end do @@ -106,31 +121,43 @@ subroutine get_odd_straight_with_v4 ( a, wr, er, transmode, amass, ityp_sc, T, v !print *, "BEFORE I - V4Lambda" !call flush() - call dgemm('N','N',nl,nl,nl,-1.0d0,v42,nl,lamat,nl,1.0d0,maux,nl) + if (.not. sym) then + call dgemm('N','N',nl,nl,nl,-1.0d0,v42,nl,lamat,nl,1.0d0,maux,nl) + else + ! v4 fully permutation-symmetric => the reordered copy v42(ja,ka)=v4(w,z,x,y) + ! equals the plain column-major reshape of v4 into an (nl,nl) matrix, so we + ! can feed v4 directly to dgemm with leading dimension nl and skip v42. + ! (Not valid with use_symmetries=False -> guarded by use_v4_symmetry.) + call dgemm('N','N',nl,nl,nl,-1.0d0,v4,nl,lamat,nl,1.0d0,maux,nl) + end if + + if (allocated(v42)) deallocate(v42) ! Invert ** iden - lamat v4 ** !print *, "BEFORE (I - V4Lambda)^-1" !call flush() - - call dgetrf ( nl, nl, maux, nl, ipiv, info ) - call dgetri ( nl, maux, nl, ipiv, work, nl, info ) - ! Take product between lamat and the inverted matrix + call dgetrf ( nl, nl, maux, nl, ipiv, info ) + call dgetri ( nl, maux, nl, ipiv, work, nl, info ) - !print *, "BEFORE Lambda(I - V4Lambda)^-1" - !call flush() - call dgemm('N','N',nl,nl,nl,1.0d0,lamat,nl,maux,nl,0.0d0,v42,nl) + ! Take product between lamat and the inverted matrix, contracted with v3. + ! + ! REASSOCIATION (numerically equivalent to ~1e-15 relative, NOT bitwise): + ! the original code formed the nl x nl product v42 := lamat * maux and then + ! cf := v42 * v3^T. We instead evaluate cf = lamat * (maux * v3^T) using an + ! nl x ns buffer tmp, which is mathematically identical by associativity + ! ( Lambda * M^-1 * v3^T = Lambda * (M^-1 * v3^T) ) and reduces both the flop + ! count (O(N^6) -> O(N^5)) and the memory (no nl x nl buffer needed). - ! Calculate final matrix products and assign the correction matrix + ! tmp = (I - v4 lamat)^-1 * v3^T (nl x ns) + call dgemm('N','T',nl,ns,nl,1.0d0,maux,nl,& + v32,ns,0.0d0,tmp,nl) - ! Calculate cf = ( 1 - lamat*v4)^-1 lamat * v3 - - !print *, "BEFORE Lambda(I - V4Lambda)^-1 V3" - !call flush() - call dgemm('N','T',nl,ns,nl,1.0d0,v42,nl,& - v32,ns,0.0d0,cf,nl) + ! cf = lamat * tmp = lamat (I - v4 lamat)^-1 v3^T (nl x ns) + call dgemm('N','N',nl,ns,nl,1.0d0,lamat,nl,& + tmp,nl,0.0d0,cf,nl) ! Now get: ! v3 * ( 1 - lamat*v4)^-1 lamat * v3 @@ -145,6 +172,6 @@ subroutine get_odd_straight_with_v4 ( a, wr, er, transmode, amass, ityp_sc, T, v ! Deallocate stuff - deallocate(lamat,v32,v42,maux,ipiv,work, cf) + deallocate(lamat,v32,maux,ipiv,work, cf, tmp) -end subroutine get_odd_straight_with_v4 \ No newline at end of file +end subroutine get_odd_straight_with_v4 From da045e28cf389ec6deeb5af64a954a1c5737527c Mon Sep 17 00:00:00 2001 From: SorBalda Date: Fri, 10 Jul 2026 16:12:10 +0200 Subject: [PATCH 3/3] memopt(kron): never materialize Lambda in get_odd_straight_with_v4 The nl x nl Lambda matrix (nl = n_mode^2) was built by get_cmat, which internally allocates mat_e + mat_et (two MORE nl x nl temporaries) and does cmat = mat_e . mat_et with an O(N^6) dgemm: true peak inside that call was v4 + mat_e + mat_et + cmat = 4 units = 43.5 GB for N = 192 -- the actual OOM site on a 31 GiB machine, untouched by commits 1-2. get_cmat is no longer called (nobody calls it now). We use its exact factorized structure Lambda(ka,ka') = sum_{mu,nu} e(nu,x) e(mu,y) D(mu,nu) e(nu,x') e(mu,y') (e from get_emat, D = mat_w/2 from get_g, both N x N) and build M = I - v4.Lambda with partial Kronecker contractions against the small e (four O(N^5) dgemm stages + O(N^4) diagonal scaling), ping-ponging between exactly TWO nl x nl buffers: maux and v4 itself used as scratch. The final assembly swaps both index pairs to reproduce the original v42(ja,ka)=v4(w,z,x,y) ordering EXACTLY; no permutation symmetry of v4 is assumed anywhere. cf = Lambda.tmp is factorized the same way on skinny nl x ns matrices (O(N^3) buffers). The only O(N^6) operation left is the dgetrf/dgetri inversion. CONTRACT CHANGE: v4 is now intent(inout) and its content is DESTROYED. Safe for the only caller, Ensemble.get_free_energy_hessian: d4 comes from SCHAModules.get_v4 (F-contiguous, so f2py intent(inout) does not copy), is symmetrized in place, and is never used after this call. The commit-2 optional flag use_v4_symmetry is removed: it only existed to skip the v42 copy, which no longer exists; the kron path is exact for arbitrary v4 and no Python code ever passed the flag. Memory: true whole-call peak 2 units + O(N^3) = 21.7 GB for N = 192 (was 43.5 GB inside get_cmat). Verified vs pristine v1.5 with the standalone driver (n_mode = 6): max elementwise relative diff 1.46e-14 with permutation-symmetrized v3/v4 and 1.86e-14 with raw non-symmetrized v4 (tolerance 1e-12; not bitwise, different but mathematically identical summation order). Co-Authored-By: Claude Fable 5 --- SCHAModules/get_odd_straight_with_v4.f90 | 323 +++++++++++++++-------- 1 file changed, 207 insertions(+), 116 deletions(-) diff --git a/SCHAModules/get_odd_straight_with_v4.f90 b/SCHAModules/get_odd_straight_with_v4.f90 index 8dda0fb4d..159df7016 100644 --- a/SCHAModules/get_odd_straight_with_v4.f90 +++ b/SCHAModules/get_odd_straight_with_v4.f90 @@ -1,10 +1,89 @@ -! This subroutine calculates the L mat needed to get the average of the -! third order derivatives. It is formed by four polarization vectors -! times the mass^1/2 divided by the normal length. +! This subroutine computes the odd (third order) SCHA correction to the +! free energy Hessian, including the fourth order term v4: +! +! phi_sc_odd = v3^T . Lambda . (I - v4.Lambda)^-1 . v3 +! +! (R. Bianco et al., PRB 96, 014111 (2017), Eq. 27), with nl = n_mode^2. +! +! MEMORY-OPTIMIZED "KRON" VERSION: the nl x nl Lambda matrix is NEVER +! materialized. The old path called get_cmat, which allocated TWO extra +! nl x nl temporaries (mat_e, mat_et) plus the nl x nl output cmat and did +! an O(N^6) dgemm; together with v4 that was a 4-array (4 "d4 units") peak. +! Instead we exploit the exact factorized structure that get_cmat builds +! (see get_cmat.f90, loops at its lines ~78-90): +! +! mat_e (ka,ja) = e(nu,x) * e(mu,y) +! mat_et(ja,ka) = mat_e(ka,ja) * mat_w(mu,nu) * 0.5 +! Lambda = mat_e . mat_et +! +! with ka = (x-1)*N + y (y fast) and ja = (mu-1)*N + nu (nu fast), i.e. +! +! Lambda(ka,ka') = sum_{mu,nu} e(nu,x) e(mu,y) D(mu,nu) e(nu,x') e(mu,y') +! +! where e(N,N) comes from get_emat, D(mu,nu) = mat_w(mu,nu)/2 with mat_w +! from get_g (both small, O(N^2)). Applying Lambda therefore reduces to +! contractions with the SMALL matrix e (each an O(N^5) dgemm) plus a +! diagonal scaling over the mode pair (O(N^4)). The only O(N^6) operation +! left is the inversion of (I - v4.Lambda). +! +! EXACT REPRODUCTION OF THE ORIGINAL PRODUCT (index bookkeeping). +! The original code built maux = I - v42 . Lambda with the reordered copy +! v42(ja,ka) = v4(w,z,x,y), ja = (w-1)*N + z (z fast), +! ka = (x-1)*N + y (y fast). +! Expanding, the matrix subtracted from the identity is, as a 4-tensor, +! +! P(w,z,x,y) = sum_{x',y',mu,nu} v4(w,z,x',y') e(nu,x') e(mu,y') +! * D(mu,nu) * e(nu,x) e(mu,y) +! +! placed at maux( (w-1)*N + z , (x-1)*N + y ). NO permutation symmetry of +! v4 is assumed anywhere: this is an identity in the indices, valid for +! arbitrary v4 (verified numerically against v1.5 with a NON-symmetrized +! random v4 as well; the old optional flag use_v4_symmetry is gone since +! the v42 copy it avoided no longer exists). +! +! P is evaluated with partial (Kronecker-factor) contractions that +! ping-pong between exactly TWO nl x nl buffers: maux and v4 ITSELF used +! as scratch. Layouts below are column-major, leftmost index fastest: +! +! A: B(w,z,x',mu) = sum_y' v4(w,z,x',y') e(mu,y') v4 -> maux +! one dgemm: (N^3 x N) . (N x N)^T, O(N^5) +! B: T(w,z,nu,mu) = sum_x' B(w,z,x',mu) e(nu,x') maux-> v4 +! N slice dgemms over mu: (N^2 x N) . (N x N)^T, O(N^5) +! C: T(w,z,nu,mu) *= D(mu,nu) v4 in place, O(N^4) +! D: C(w,z,nu,y) = sum_mu T(w,z,nu,mu) e(mu,y) v4 -> maux +! one dgemm: (N^3 x N) . (N x N), O(N^5) +! E: P(w,z,x,y) = sum_nu C(w,z,nu,y) e(nu,x) maux-> v4 +! N slice dgemms over y: (N^2 x N) . (N x N), O(N^5) +! F: maux((w-1)N+z,(x-1)N+y) = delta - P(w,z,x,y) v4 -> maux, O(N^4) +! explicit permuted copy: the natural buffer layout of P has (w fast, +! z slow) in the row pair and (x fast, y slow) in the column pair, +! while the original v42.Lambda uses (z fast, w slow) rows and +! (y fast, x slow) columns, so BOTH index pairs are swapped here. +! +! After F: LU inversion in place in maux (dgetrf/dgetri, unchanged), then +! tmp = maux . v3^T (nl x ns, O(N^5)) +! cf = Lambda . tmp via the same e/D/e^T factorization +! applied to skinny matrices, using two +! O(N^3) buffers (N x N x ns) +! phi_sc_odd = v32 . cf (as before) +! +! !!! WARNING: v4 IS INTENT(INOUT) AND ITS CONTENT IS DESTROYED !!! +! v4 is deliberately used as one of the two nl x nl scratch buffers (its +! original content is consumed by step A before step B overwrites it). +! This is safe for the only caller, Ensemble.get_free_energy_hessian +! (Modules/Ensemble.py): it creates d4 with SCHAModules.get_v4 (hence +! F-contiguous, so f2py's intent(inout) wraps it in place without a copy), +! optionally symmetrizes it in place, calls this routine, and never uses +! d4 again. Any NEW caller must pass a throw-away, F-contiguous v4. +! +! Peak large-memory budget: exactly 2 nl x nl arrays alive (v4 + maux), +! plus O(N^3) skinny buffers (v32, tmp, cf, b1, b2) and O(N^2)/O(N) work +! arrays. For N = 192 (Au 4x4x4): 2 * 10.87 GB = 21.7 GB, vs ~43.5 GB +! inside the old get_cmat call (v4 + mat_e + mat_et + cmat). subroutine get_odd_straight_with_v4 ( a, wr, er, transmode, amass, ityp_sc, T, v3, v4, phi_sc_odd, & - n_mode, nat_sc, ntyp, use_v4_symmetry) + n_mode, nat_sc, ntyp) implicit none @@ -15,163 +94,175 @@ subroutine get_odd_straight_with_v4 ( a, wr, er, transmode, amass, ityp_sc, T, v integer, dimension(nat_sc), intent(in) :: ityp_sc double precision, intent(in) :: T double precision, dimension(n_mode,n_mode,n_mode), intent(in) :: v3 - double precision, dimension(n_mode,n_mode,n_mode, n_mode), intent(in) :: v4 + ! v4 is used as scratch and DESTROYED on exit -- see warning above. + double precision, dimension(n_mode,n_mode,n_mode,n_mode), intent(inout) :: v4 double precision, dimension(n_mode, n_mode), intent(out) :: phi_sc_odd - ! Optional: if present and .true., assume v4 is permutation-symmetric - ! (standard case, after ApplySymmetryToTensor4 / use_symmetries=True) and - ! skip building the explicit reordered copy v42 -- see note below. - logical, intent(in), optional :: use_v4_symmetry - integer :: nat_sc, n_mode, nl, ns, ntyp - double precision, dimension(:,:), allocatable :: l, g, phi_aux, v1, v2, v32 - double precision :: lsum - double precision, dimension(:), allocatable :: laux1, lres1, veclong - double precision, dimension(:), allocatable :: laux2, lres2 - double precision, dimension(:,:), allocatable :: lamat, v42, maux + ! The ONLY allocated nl x nl array (v4, the other big buffer, is the + ! caller's own array). + double precision, dimension(:,:), allocatable :: maux + + ! Small O(N^2) factors of Lambda. + double precision, dimension(:,:), allocatable :: e ! from get_emat + double precision, dimension(:,:), allocatable :: dmat ! D = mat_w/2, from get_g + + ! Skinny O(N^3) buffers. + double precision, dimension(:,:), allocatable :: v32 ! ns x nl + double precision, dimension(:,:), allocatable :: tmp, cf ! nl x ns + double precision, dimension(:,:,:), allocatable :: b1, b2 ! ns x ns x ns + double precision, dimension(:), allocatable :: work integer, dimension(:), allocatable :: ipiv integer :: info - double precision, dimension(:,:), allocatable :: cf - double precision, dimension(:,:), allocatable :: tmp - - integer :: mu, nu, alpha - integer :: ka, ja - integer :: i, j, x, y, z, w - - real :: t1, t2 + integer :: mu, nu + integer :: ka + integer :: i, s, x, y, z, w logical, parameter :: debug = .true. - logical :: sym - - ! Get integers if (debug) then - print *, "=== DEBUG ODD STRAIGHT ===" + print *, "=== DEBUG ODD STRAIGHT (kron, Lambda never materialized) ===" print *, "N_MODE:", n_mode print *, "NTYP:", ntyp print *, "NAT_SC:", nat_sc call flush() end if - !nat_sc = size(er(:,1,1)) - !n_mode = 3*nat_sc - ns = n_mode nl = n_mode*n_mode - sym = .false. - if (present(use_v4_symmetry)) sym = use_v4_symmetry - - ! Allocate stuff - - allocate(lamat(nl,nl)) allocate(maux(nl,nl)) + allocate(e(ns,ns)) + allocate(dmat(ns,ns)) + allocate(v32(ns,nl)) + allocate(tmp(nl,ns)) + allocate(cf(nl,ns)) + allocate(b1(ns,ns,ns)) + allocate(b2(ns,ns,ns)) allocate(ipiv(nl)) allocate(work(nl)) - allocate(v32(n_mode,n_mode*n_mode)) - - allocate(cf(nl,ns)) - allocate(tmp(nl,ns)) - - ! v42 (explicit reordered copy of v4) is only needed when we cannot rely on - ! the permutation symmetry of v4. It is one full nl x nl array. - if (.not. sym) allocate(v42(nl,nl)) - ! Get lambda matrix - - call get_cmat ( a, wr, er, transmode, amass, ityp_sc, T, .true., lamat,n_mode, nat_sc, ntyp ) - - !print *, "AFTER CMAT" - !call flush() - - ! Write third and fourth order force constants as rank 2 + ! Small factors of Lambda: exactly what get_cmat used internally. + ! (v3_log = .true., as in the old get_cmat call from this routine.) + call get_emat ( er, a, amass, ityp_sc, .true., transmode, e, n_mode, nat_sc, ntyp) + call get_g (a, wr, transmode, T, dmat, n_mode) + dmat = 0.5d0 * dmat ! D(mu,nu) = mat_w(mu,nu)/2 (the 0.5 of mat_et) + ! Third order force constants as rank 2, v32(:,ka) with ka=(x-1)*N+y. ka = 0 - - do x = 1, n_mode - do y = 1, n_mode + do x = 1, ns + do y = 1, ns ka = ka + 1 v32(:,ka) = v3(:,x,y) - if (.not. sym) then - ja = 0 - do w = 1, n_mode - do z = 1, n_mode - ja = ja + 1 - v42(ja,ka) = v4(w,z,x,y) - end do - end do - end if end do end do - ! Prepare identity matrix directly in maux (was: iden then maux = iden; - ! bitwise identical, avoids one nl x nl temporary) - - maux = 0.0d0 - - do x = 1, nl - maux(x,x) = 1.0d0 + ! --------------------------------------------------------------------- + ! Build maux = I - v42.Lambda without forming Lambda (steps A-F above). + ! --------------------------------------------------------------------- + + ! A: B(w,z,x',mu) = sum_y' v4(w,z,x',y') e(mu,y') [v4 -> maux] + ! v4 viewed as (N^3 x N) with columns y'; op(B)=e^T has entry + ! (y',mu) = e(mu,y'). Fills maux completely (N^3 * N = nl*nl). + call dgemm('N','T', nl*ns, ns, ns, 1.0d0, v4(1,1,1,1), nl*ns, & + e(1,1), ns, 0.0d0, maux(1,1), nl*ns) + + ! B: T(w,z,nu,mu) = sum_x' B(w,z,x',mu) e(nu,x') [maux -> v4] + ! For each mu, the slice B(:,:,:,mu) is the (N^2 x N) block starting + ! at maux(1,(mu-1)*ns+1) (linear offset (mu-1)*N^3), columns x'; + ! result panel written at v4(:,:,:,mu) with layout (w,z,nu). + ! From here on the original content of v4 is destroyed. + do mu = 1, ns + call dgemm('N','T', nl, ns, ns, 1.0d0, maux(1,(mu-1)*ns+1), nl, & + e(1,1), ns, 0.0d0, v4(1,1,1,mu), nl) end do - ! Calculate ** iden - v4 lamat ** matrix - - !print *, "BEFORE I - V4Lambda" - !call flush() - - if (.not. sym) then - call dgemm('N','N',nl,nl,nl,-1.0d0,v42,nl,lamat,nl,1.0d0,maux,nl) - else - ! v4 fully permutation-symmetric => the reordered copy v42(ja,ka)=v4(w,z,x,y) - ! equals the plain column-major reshape of v4 into an (nl,nl) matrix, so we - ! can feed v4 directly to dgemm with leading dimension nl and skip v42. - ! (Not valid with use_symmetries=False -> guarded by use_v4_symmetry.) - call dgemm('N','N',nl,nl,nl,-1.0d0,v4,nl,lamat,nl,1.0d0,maux,nl) - end if - - if (allocated(v42)) deallocate(v42) + ! C: T(w,z,nu,mu) *= D(mu,nu) [v4 in place] + ! NOTE the argument order: 4th dim of T is mu, 3rd is nu. + do mu = 1, ns + do nu = 1, ns + v4(:,:,nu,mu) = v4(:,:,nu,mu) * dmat(mu,nu) + end do + end do - ! Invert ** iden - lamat v4 ** + ! D: C(w,z,nu,y) = sum_mu T(w,z,nu,mu) e(mu,y) [v4 -> maux] + call dgemm('N','N', nl*ns, ns, ns, 1.0d0, v4(1,1,1,1), nl*ns, & + e(1,1), ns, 0.0d0, maux(1,1), nl*ns) - !print *, "BEFORE (I - V4Lambda)^-1" - !call flush() + ! E: P(w,z,x,y) = sum_nu C(w,z,nu,y) e(nu,x) [maux -> v4] + do y = 1, ns + call dgemm('N','N', nl, ns, ns, 1.0d0, maux(1,(y-1)*ns+1), nl, & + e(1,1), ns, 0.0d0, v4(1,1,1,y), nl) + end do + ! F: maux = I - P with BOTH index pairs swapped to the original + ! v42.Lambda convention: rows (z fast, w slow), cols (y fast, x slow). + do y = 1, ns + do x = 1, ns + ka = (x-1)*ns + y + do z = 1, ns + do w = 1, ns + maux((w-1)*ns + z, ka) = -v4(w,z,x,y) + end do + end do + end do + end do + do i = 1, nl + maux(i,i) = maux(i,i) + 1.0d0 + end do + ! Invert ** iden - v4 lamat ** in place (unchanged from v1.5). call dgetrf ( nl, nl, maux, nl, ipiv, info ) call dgetri ( nl, maux, nl, ipiv, work, nl, info ) - ! Take product between lamat and the inverted matrix, contracted with v3. - ! - ! REASSOCIATION (numerically equivalent to ~1e-15 relative, NOT bitwise): - ! the original code formed the nl x nl product v42 := lamat * maux and then - ! cf := v42 * v3^T. We instead evaluate cf = lamat * (maux * v3^T) using an - ! nl x ns buffer tmp, which is mathematically identical by associativity - ! ( Lambda * M^-1 * v3^T = Lambda * (M^-1 * v3^T) ) and reduces both the flop - ! count (O(N^6) -> O(N^5)) and the memory (no nl x nl buffer needed). - - ! tmp = (I - v4 lamat)^-1 * v3^T (nl x ns) - call dgemm('N','T',nl,ns,nl,1.0d0,maux,nl,& - v32,ns,0.0d0,tmp,nl) - - ! cf = lamat * tmp = lamat (I - v4 lamat)^-1 v3^T (nl x ns) - call dgemm('N','N',nl,ns,nl,1.0d0,lamat,nl,& - tmp,nl,0.0d0,cf,nl) - - ! Now get: - ! v3 * ( 1 - lamat*v4)^-1 lamat * v3 + ! tmp = (I - v4 lamat)^-1 . v3^T (nl x ns) + call dgemm('N','T', nl, ns, nl, 1.0d0, maux(1,1), nl, & + v32(1,1), ns, 0.0d0, tmp(1,1), nl) + + ! --------------------------------------------------------------------- + ! cf = Lambda . tmp, again without forming Lambda. With tmp's row index + ! ka' = (x'-1)*N + y' (y' fast) read as tmp3(y',x',s): + ! u(mu,nu,s) = sum_{x',y'} e(mu,y') e(nu,x') tmp3(y',x',s) + ! cf3(y,x,s) = sum_{mu,nu} e(mu,y) e(nu,x) D(mu,nu) u(mu,nu,s) + ! All buffers are O(N^3). + ! --------------------------------------------------------------------- + + ! C1: b1(mu,x',s) = sum_y' e(mu,y') tmp3(y',x',s) + ! tmp reshaped as (N x N*ns), one dgemm. + call dgemm('N','N', ns, nl, ns, 1.0d0, e(1,1), ns, & + tmp(1,1), ns, 0.0d0, b1(1,1,1), ns) + + ! C2: b2(mu,nu,s) = sum_x' b1(mu,x',s) e(nu,x') + do s = 1, ns + call dgemm('N','T', ns, ns, ns, 1.0d0, b1(1,1,s), ns, & + e(1,1), ns, 0.0d0, b2(1,1,s), ns) + end do - !print *, "BEFORE V3 Lambda(I - V4Lambda)^-1 V3" - !call flush() - call dgemm('N','N',ns,ns,nl,1.0d0,v32,ns,& - cf,nl,0.0d0,phi_sc_odd,ns) + ! C3: b2(mu,nu,s) *= D(mu,nu) + do s = 1, ns + b2(:,:,s) = b2(:,:,s) * dmat(:,:) + end do + ! C4: b1(y,nu,s) = sum_mu e(mu,y) b2(mu,nu,s) (e^T . b2_s) + do s = 1, ns + call dgemm('T','N', ns, ns, ns, 1.0d0, e(1,1), ns, & + b2(1,1,s), ns, 0.0d0, b1(1,1,s), ns) + end do - !call get_odd_from_cmat_fu2 (v42, v32, phi_sc_odd) + ! C5: cf3(y,x,s) = sum_nu b1(y,nu,s) e(nu,x); column s of cf viewed as + ! an (N x N) block (y fast, x slow), matching ka = (x-1)*N + y. + do s = 1, ns + call dgemm('N','N', ns, ns, ns, 1.0d0, b1(1,1,s), ns, & + e(1,1), ns, 0.0d0, cf(1,s), ns) + end do - ! Deallocate stuff + ! phi_sc_odd = v3 . lamat (I - v4 lamat)^-1 . v3 + call dgemm('N','N', ns, ns, nl, 1.0d0, v32(1,1), ns, & + cf(1,1), nl, 0.0d0, phi_sc_odd(1,1), ns) - deallocate(lamat,v32,maux,ipiv,work, cf, tmp) + deallocate(maux, e, dmat, v32, tmp, cf, b1, b2, ipiv, work) end subroutine get_odd_straight_with_v4