From 76c785e22d06fd809e5f54e332a7d2b5d41d0340 Mon Sep 17 00:00:00 2001 From: Matias Date: Fri, 17 Nov 2023 15:10:33 -0300 Subject: [PATCH] bug on PropagateRhorealHBlacs solved and routine cleaned --- src/dftbp/timedep/timeprop.F90 | 174 +++------------------------------ 1 file changed, 12 insertions(+), 162 deletions(-) diff --git a/src/dftbp/timedep/timeprop.F90 b/src/dftbp/timedep/timeprop.F90 index 1c24a96bf6..6d686b2e0b 100644 --- a/src/dftbp/timedep/timeprop.F90 +++ b/src/dftbp/timedep/timeprop.F90 @@ -2344,7 +2344,7 @@ subroutine propagateRhoRealHBlacs(this, rhoOld, rho, H1, Sinv, step) !> Time step in atomic units real(dp), intent(in) :: step - real(dp), allocatable :: T1R(:,:), T2R(:,:), T3R(:,:),T4R(:,:), T5R(:,:) + real(dp), allocatable :: T1R(:,:), T2R(:,:), T3R(:,:),T4R(:,:) integer :: nLocalCols, nLocalRows, i,j @@ -2355,166 +2355,53 @@ subroutine propagateRhoRealHBlacs(this, rhoOld, rho, H1, Sinv, step) allocate(T2R(nLocalRows,nLocalCols)) allocate(T3R(nLocalRows,nLocalCols)) allocate(T4R(nLocalRows,nLocalCols)) - allocate(T5R(nLocalRows,nLocalCols)) ! The code below takes into account that Sinv and H1 are real, this is twice as fast as the ! original above (propageteRho) ! get the real part of Sinv and H1 + T2R(:,:) = real(H1) T1R(:,:) = real(rho) - !call gemm(T3R,T1R,T2R) - call pblasfx_psymm(T2R, this%denseDesc%blacsOrbSqr, T1R, this%denseDesc%blacsOrbSqr,& & T3R, this%denseDesc%blacsOrbSqr, side="R") - -! =============== - ! write(*, *) "(H1) pblasfx_psymm" - ! do i =1,2 - ! do j=1,2 - ! write(*, *) i,j,H1(i,j) - ! enddo - ! enddo - - ! call gemm(T4R,T1R,T2R) - - ! write(*, *) "Real(Rho * H1 ) gemm" - ! do i =1,2 - ! do j=1,2 - ! write(*, *) i,j,T4R(i,j) - ! enddo - ! enddo -! ================= + ! T3R= Re[Rho]*Re[H] T2R(:,:) = T3R T1R(:,:) = real(Sinv) - ! Call ORIGINAL call pblasfx_psymm(T1R, this%denseDesc%blacsOrbSqr, T2R, this%denseDesc%blacsOrbSqr,& & T3R, this%denseDesc%blacsOrbSqr, side="R") - - ! Call ALTERNATIVO - ! call pblasfx_psymm(T2R, this%denseDesc%blacsOrbSqr, T1R, this%denseDesc%blacsOrbSqr,& - ! & T3R, this%denseDesc%blacsOrbSqr) - -! call gemm(T3R,T2R,T1R) - - ! T3R real(rho)HSinv - -! ================= + ! T3R= Re[Rho]*Re[H]*Re[Sinv] - ! write(*, *) "Real(Rho * H1 * Sinv) pblasfx_psymm" - ! do i =1,2 - ! do j=1,2 - ! write(*, *) i,j,T3R(i,j) - ! enddo - ! enddo - - ! call gemm(T5R,T2R,T1R) - - ! write(*, *) "Real(Rho * H1 * Sinv) gemm" - ! do i =1,2 - ! do j=1,2 - ! write(*, *) i,j,T5R(i,j) - ! enddo - ! enddo -! ================= - T1R(:,:) = aimag(rho) T2R(:,:) = real(H1) -! call gemm(T4R,T1R,T2R) call pblasfx_psymm(T2R, this%denseDesc%blacsOrbSqr, T1R, this%denseDesc%blacsOrbSqr,& & T4R, this%denseDesc%blacsOrbSqr, side="R") + ! T4R= Im[Rho]*Im[H] T2R(:,:) = T4R T1R(:,:) = real(Sinv) -! call gemm(T4R,T2R,T1R) - - ! write(*, *) "Imag(Rho * H1 ) pblasfx_psymm" - ! do i =1,2 - ! do j=1,2 - ! write(*, *) i,j,T4R(i,j) - ! enddo - ! enddo - - ! call original: - call pblasfx_psymm(T1R, this%denseDesc%blacsOrbSqr, T2R, this%denseDesc%blacsOrbSqr,& - & T4R, this%denseDesc%blacsOrbSqr, side="R") - - !Call alternativo - ! call pblasfx_psymm(T2R, this%denseDesc%blacsOrbSqr, T1R, this%denseDesc%blacsOrbSqr,& - ! & T4R, this%denseDesc%blacsOrbSqr) - -! ================= - ! write(*, *) "Imag(Rho * H1 * Sinv) pblasfx_psymm" - ! do i =1,2 - ! do j=1,2 - ! write(*, *) i,j,T4R(i,j) - ! enddo - ! enddo - - ! call gemm(T5R,T2R,T1R) - - ! write(*, *) "Imag(Rho * H1 * Sinv) gemm" - ! do i =1,2 - ! do j=1,2 - ! write(*, *) i,j,T5R(i,j) - ! enddo - ! enddo -! ================= - - ! T4R= imag(rho)HSinv + call pblasfx_psymm(T1R, this%denseDesc%blacsOrbSqr, T2R, this%denseDesc%blacsOrbSqr,& + & T4R, this%denseDesc%blacsOrbSqr, side="R") + ! T4R= Im[Rho]*Im[H]*Im[Sinv] call pblasfx_ptran(T3R, this%denseDesc%blacsOrbSqr, T1R, this%denseDesc%blacsOrbSqr) - ! T1R = transpose(T3R) call pblasfx_ptran(T4R, this%denseDesc%blacsOrbSqr, T2R, this%denseDesc%blacsOrbSqr) - ! T2R = transpose(T4R) -! ================= - ! write(*, *) "Real(Rho * H1 * Sinv)TRANS " - ! do i =1,2 - ! do j=1,2 - ! write(*, *) i,j,T1R(i,j) - ! enddo - ! enddo - - - ! write(*, *) "iamg(Rho * H1 * Sinv)TRANS " - ! do i =1,2 - ! do j=1,2 - ! write(*, *) i,j,T2R(i,j) - ! enddo - ! enddo -! ================= - -! ================= -! (rhoOld) - ! write(*, *) "(Rho_old)" - ! do i =1,2 - ! do j=1,2 - ! write(*, *) i,j,rhoOld(i,j) - ! enddo - ! enddo -! ================= - + ! T1R = Transpose(T3R) + ! T4R = Transpose(T2R) ! build the commutator combining the real and imaginary parts of the previous result - ! NEWWWW : !$OMP WORKSHARE rhoOld(:,:) = rhoOld + cmplx(0, step, dp) * (T3R + imag * T4R)& & + cmplx(0, -step, dp) * (T1R - imag * T2R) !$OMP END WORKSHARE - ! !$OMP WORKSHARE - ! rhoOld(:,:) = rhoOld + cmplx(0, -step, dp) * (T3R + imag * T4R)& - ! & + cmplx(0, step, dp) * (T1R - imag * T2R) - ! !$OMP END WORKSHARE - - deallocate(T1R) deallocate(T2R) @@ -2545,7 +2432,7 @@ subroutine propagateRhoRealH(this, rhoOld, rho, H1, Sinv, step) !> Time step in atomic units real(dp), intent(in) :: step - real(dp), allocatable :: T1R(:,:), T2R(:,:), T3R(:,:),T4R(:,:), T5R(:,:), T6R(:,:) + real(dp), allocatable :: T1R(:,:), T2R(:,:), T3R(:,:),T4R(:,:) integer :: i, j @@ -2553,8 +2440,6 @@ subroutine propagateRhoRealH(this, rhoOld, rho, H1, Sinv, step) allocate(T2R(this%nOrbs,this%nOrbs)) allocate(T3R(this%nOrbs,this%nOrbs)) allocate(T4R(this%nOrbs,this%nOrbs)) - allocate(T5R(this%nOrbs,this%nOrbs)) - ! The code below takes into account that Sinv and H1 are real, this is twice as fast as the ! original above (propageteRho) @@ -2562,14 +2447,7 @@ subroutine propagateRhoRealH(this, rhoOld, rho, H1, Sinv, step) ! get the real part of Sinv and H1 T1R(:,:) = real(H1) T2R(:,:) = real(Sinv) - - write(*, *) "H1" - do i =1,2 - do j=1,2 - write(*, *) i,j,H1(i,j) - enddo - enddo - + call gemm(T3R,T2R,T1R) T2R(:,:) = T3R @@ -2580,34 +2458,6 @@ subroutine propagateRhoRealH(this, rhoOld, rho, H1, Sinv, step) T1R(:,:) = aimag(rho) call gemm(T4R,T2R,T1R) - ! ================= - T6R = transpose(T3R) - T5R = transpose(T4R) - - - write(*, *) "TRANS Real(Sinv * H1 * Rho) " - do i =1,2 - do j=1,2 - - write(*, *) i,j,T6R(i,j) - enddo - enddo - - write(*, *) "TRANS iamg(Sinv * H1 * Rho) " - do i =1,2 - do j=1,2 - write(*, *) i,j,T5R(i,j) - enddo - enddo -! ================= - - write(*, *) "(Rho_old)" - do i =1,2 - do j=1,2 - write(*, *) i,j,rhoOld(i,j) - enddo - enddo - ! build the commutator combining the real and imaginary parts of the previous result !$OMP WORKSHARE rhoOld(:,:) = rhoOld + cmplx(0, -step, dp) * (T3R + imag * T4R)&