LCPQ/quack
4
1
Fork 0
mirror of https://github.com/pfloos/quack synced 2024-11-07 22:53:59 +01:00
Code Releases Activity

dgemm instead of matmul for nVV matrix multiplication in sort_ppRPA.f90

Browse Source
This commit is contained in:
Antoine Marie 2023-06-14 16:53:16 +02:00
parent 3941db73e3
commit 552005f910
1 changed files with 16 additions and 7 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

23
src/RPA/sort_ppRPA.f90
View File

@ -23,7 +23,9 @@ subroutine sort_ppRPA(nOO,nVV,Omega,Z,Omega1,X1,Y1,Omega2,X2,Y2)
double precision,allocatable :: S1(:,:) double precision,allocatable :: S1(:,:)
double precision,allocatable :: S2(:,:) double precision,allocatable :: S2(:,:)
double precision,allocatable :: O1(:,:) 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 :: order1(:)
integer,allocatable :: order2(:) integer,allocatable :: order2(:)
@ -86,7 +88,7 @@ subroutine sort_ppRPA(nOO,nVV,Omega,Z,Omega1,X1,Y1,Omega2,X2,Y2)
end if end if
end do end do
if(minval(Omega1(:)) < 0d0 .or. ab /= nVV) call print_warning('You may have instabilities in pp-RPA!!') 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!!') 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 quick_sort(Omega2(:),order2(:),nOO)
call set_order(Z2(:,:),order2(:),nOO+nVV,nOO) call set_order(Z2(:,:),order2(:),nOO+nVV,nOO)
end if end if
! Orthogonalize eigenvectors ! Orthogonalize eigenvectors
@ -196,16 +198,23 @@ subroutine sort_ppRPA(nOO,nVV,Omega,Z,Omega1,X1,Y1,Omega2,X2,Y2)
! ij_start = ij ! ij_start = ij
! end if ! end if
! end do ! end do
allocate(S1(nVV,nVV),S2(nOO,nOO),O1(nVV,nVV),O2(nOO,nOO)) 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)) S2 = - matmul(transpose(Z2),matmul(M,Z2))
if(nVV > 0) call orthogonalization_matrix(1,nVV,S1,O1) if(nVV > 0) call orthogonalization_matrix(1,nVV,S1,O1)
if(nOO > 0) call orthogonalization_matrix(1,nOO,S2,O2) 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) Z2 = matmul(Z2,O2)
! Define submatrices X1, Y1, X2, & Y2 ! Define submatrices X1, Y1, X2, & Y2