From 992e3dff4b98142130638f597d8cad240e21ba25 Mon Sep 17 00:00:00 2001
From: Abdallah Ammar <abd.ammar.phys@gmail.com>
Date: Wed, 21 Aug 2024 13:51:12 +0200
Subject: [PATCH] optim in AdAt

---
 src/RPA/sort_ppRPA.f90                 |  9 ++++----
 src/utils/orthogonalization_matrix.f90 |  2 +-
 src/utils/utils.f90                    | 30 +++++++++++++++++++-------
 src/utils/wrap_lapack.f90              |  4 ++++
 4 files changed, 32 insertions(+), 13 deletions(-)

diff --git a/src/RPA/sort_ppRPA.f90 b/src/RPA/sort_ppRPA.f90
index e7037ad..8518ecb 100644
--- a/src/RPA/sort_ppRPA.f90
+++ b/src/RPA/sort_ppRPA.f90
@@ -39,9 +39,10 @@ subroutine sort_ppRPA(nOO,nVV,Om,Z,Om1,X1,Y1,Om2,X2,Y2)
   double precision,intent(out)  :: X2(nVV,nOO)
   double precision,intent(out)  :: Y2(nOO,nOO)
 
+
 ! Memory allocation
 
- allocate(M(nOO+nVV,nOO+nVV),Z1(nOO+nVV,nVV),Z2(nOO+nVV,nOO),order1(nVV),order2(nOO))
+  allocate(M(nOO+nVV,nOO+nVV),Z1(nOO+nVV,nVV),Z2(nOO+nVV,nOO),order1(nVV),order2(nOO))
 
 ! Initializatiom
 
@@ -86,7 +87,7 @@ subroutine sort_ppRPA(nOO,nVV,Om,Z,Om1,X1,Y1,Om2,X2,Y2)
 
     end if
 
- end do
+  end do
 
   if(minval(Om1) < 0d0 .or. ab /= nVV) call print_warning('You may have instabilities in pp-RPA!!')
   if(maxval(Om2) > 0d0 .or. ij /= nOO) call print_warning('You may have instabilities in pp-RPA!!')
@@ -111,7 +112,8 @@ subroutine sort_ppRPA(nOO,nVV,Om,Z,Om1,X1,Y1,Om2,X2,Y2)
     call quick_sort(Om2,order2,nOO)
     call set_order(Z2,order2,nOO+nVV,nOO)
 
- end if
+  end if
+
 
 ! Orthogonalize eigenvectors
 
@@ -202,7 +204,6 @@ subroutine sort_ppRPA(nOO,nVV,Om,Z,Om1,X1,Y1,Om2,X2,Y2)
 
   if(nVV > 0) call dgemm ('N', 'N', nOO+nVV, nVV, nOO+nVV, 1d0,  M, nOO+nVV, Z1, nOO+nVV, 0d0, tmp1, nOO+nVV)
   if(nVV > 0) call dgemm ('T', 'N', nVV    , nVV, nOO+nVV, 1d0, Z1, nOO+nVV, tmp1, nOO+nVV, 0d0, S1, nVV)
-
   if(nOO > 0) call dgemm ('N', 'N', nOO+nVV, nOO, nOO+nVV, 1d0,  M, nOO+nVV, -1d0*Z2, nOO+nVV, 0d0, tmp2, nOO+nVV)
   if(nOO > 0) call dgemm ('T', 'N', nOO    , nOO, nOO+nVV, 1d0, Z2, nOO+nVV, tmp2, nOO+nVV, 0d0, S2, nOO)
 
diff --git a/src/utils/orthogonalization_matrix.f90 b/src/utils/orthogonalization_matrix.f90
index 176862a..d7e0089 100644
--- a/src/utils/orthogonalization_matrix.f90
+++ b/src/utils/orthogonalization_matrix.f90
@@ -54,7 +54,7 @@ subroutine orthogonalization_matrix(nBas,S,X)
 
     end do
     
-    call ADAt(nBas,Uvec,Uval,X)
+    call ADAt(nBas, Uvec(1,1), Uval(1), X(1,1))
 
   elseif(ortho_type == 2) then
 
diff --git a/src/utils/utils.f90 b/src/utils/utils.f90
index efd3087..80b8322 100644
--- a/src/utils/utils.f90
+++ b/src/utils/utils.f90
@@ -375,15 +375,29 @@ subroutine ADAt(N,A,D,B)
 
   double precision,intent(out)  :: B(N,N)
 
-  B = 0d0
+  double precision, allocatable :: tmp(:,:)
 
-  do i=1,N
-    do j=1,N
-      do k=1,N
-        B(i,k) = B(i,k) + A(i,j)*D(j)*A(k,j)
-      end do
-    end do
-  end do
+  allocate(tmp(N,N))
+  !$OMP PARALLEL DEFAULT(NONE) PRIVATE(i, j) SHARED(N, A, D, tmp)
+  !$OMP DO
+  do i = 1, N
+    do j = 1, N
+      tmp(i,j) = D(i) * A(j,i)
+    enddo
+  enddo
+  !$OMP END DO
+  !$OMP END PARALLEL
+  call dgemm("N", "N", N, N, N, 1.d0, A(1,1), N, tmp(1,1), N, 0.d0, B(1,1), N)
+  deallocate(tmp)
+
+!  B = 0d0
+!  do i=1,N
+!    do j=1,N
+!      do k=1,N
+!        B(i,k) = B(i,k) + A(i,j)*D(j)*A(k,j)
+!      end do
+!    end do
+!  end do
 
 end subroutine 
 !------------------------------------------------------------------------
diff --git a/src/utils/wrap_lapack.f90 b/src/utils/wrap_lapack.f90
index d63b722..86280fe 100644
--- a/src/utils/wrap_lapack.f90
+++ b/src/utils/wrap_lapack.f90
@@ -31,6 +31,8 @@ subroutine diagonalize_general_matrix(N,A,WR,VR)
 
   call dgeev('V','V',N,A,N,WR,WI,VL,N,VR,N,work,lwork,info)
 
+  deallocate(work, WI, VL)
+
   if(info /= 0) then 
     print*,'Problem in diagonalize_general_matrix (dgeev)!!'
   end if
@@ -67,6 +69,8 @@ subroutine diagonalize_matrix(N,A,e)
   allocate(work(lwork))
 
   call dsyev('V','U',N,A,N,e,work,lwork,info)
+
+  deallocate(work)
  
   if(info /= 0) then 
     print*,'Problem in diagonalize_matrix (dsyev)!!'