From 552005f910506b70c7640dabd913485dc9f4579e Mon Sep 17 00:00:00 2001 From: Antoine Marie Date: Wed, 14 Jun 2023 16:53:16 +0200 Subject: [PATCH] dgemm instead of matmul for nVV matrix multiplication in sort_ppRPA.f90 --- src/RPA/sort_ppRPA.f90 | 23 ++++++++++++++++------- 1 file changed, 16 insertions(+), 7 deletions(-) diff --git a/src/RPA/sort_ppRPA.f90 b/src/RPA/sort_ppRPA.f90 index 4028c7f..f24743f 100644 --- a/src/RPA/sort_ppRPA.f90 +++ b/src/RPA/sort_ppRPA.f90 @@ -23,7 +23,9 @@ subroutine sort_ppRPA(nOO,nVV,Omega,Z,Omega1,X1,Y1,Omega2,X2,Y2) double precision,allocatable :: S1(:,:) double precision,allocatable :: S2(:,:) double precision,allocatable :: O1(:,:) - double precision,allocatable :: O2(:,:) + double precision,allocatable :: O2(:,:) + double precision,allocatable :: tmp1(:,:) + double precision,allocatable :: tmp2(:,:) integer,allocatable :: order1(:) integer,allocatable :: order2(:) @@ -86,7 +88,7 @@ subroutine sort_ppRPA(nOO,nVV,Omega,Z,Omega1,X1,Y1,Omega2,X2,Y2) end if - end do + end do if(minval(Omega1(:)) < 0d0 .or. ab /= nVV) call print_warning('You may have instabilities in pp-RPA!!') if(maxval(Omega2(:)) > 0d0 .or. ij /= nOO) call print_warning('You may have instabilities in pp-RPA!!') @@ -112,7 +114,7 @@ subroutine sort_ppRPA(nOO,nVV,Omega,Z,Omega1,X1,Y1,Omega2,X2,Y2) call quick_sort(Omega2(:),order2(:),nOO) call set_order(Z2(:,:),order2(:),nOO+nVV,nOO) - end if + end if ! Orthogonalize eigenvectors @@ -196,16 +198,23 @@ subroutine sort_ppRPA(nOO,nVV,Omega,Z,Omega1,X1,Y1,Omega2,X2,Y2) ! ij_start = ij ! end if -! end do - +! end do + allocate(S1(nVV,nVV),S2(nOO,nOO),O1(nVV,nVV),O2(nOO,nOO)) - S1 = + matmul(transpose(Z1),matmul(M,Z1)) + allocate(tmp1(nOO+nVV,nVV),tmp2(nOO+nVV,nOO)) + + call dgemm ('N', 'N', nOO+nVV, nVV, nOO+nVV, 1d0, M, nOO+nVV, Z1, nOO+nVV, 0d0, tmp1, nOO+nVV) + call dgemm ('T', 'N', nVV , nVV, nOO+nVV, 1d0, Z1, nOO+nVV, tmp1, nOO+nVV, 0d0, S1, nVV) + S2 = - matmul(transpose(Z2),matmul(M,Z2)) if(nVV > 0) call orthogonalization_matrix(1,nVV,S1,O1) if(nOO > 0) call orthogonalization_matrix(1,nOO,S2,O2) - Z1 = matmul(Z1,O1) + + write (*,*) 'OK SO FAR' + + call dgemm ('N', 'N', nVV, nVV, nVV, 1d0, Z1, nVV, O1, nVV, 0d0, Z1, nVV) Z2 = matmul(Z2,O2) ! Define submatrices X1, Y1, X2, & Y2