dgemm

src/casscf/bielec_natorb.irp.f

@@ -113,88 +113,106 @@ END_PROVIDER
 
 
 
-BEGIN_PROVIDER [real*8, bielec_PxxQ_no, (mo_num,n_core_inact_orb+n_act_orb,n_core_inact_orb+n_act_orb, mo_num)]
+BEGIN_PROVIDER [real*8, bielec_PxxQ_no, (mo_num,n_core_inact_act_orb,n_core_inact_act_orb, mo_num)]
   BEGIN_DOC
   ! integral (px|xq) in the basis of natural MOs
   ! indices are unshifted orbital numbers
   END_DOC
   implicit none
   integer                        :: i,j,k,l,t,u,p,q,pp
-  real*8                         :: d(n_act_orb)
+  double precision, allocatable  :: f(:,:,:), d(:,:,:)
 
   bielec_PxxQ_no(:,:,:,:) = bielec_PxxQ(:,:,:,:)
 
+  allocate (f(n_act_orb,n_core_inact_act_orb,n_core_inact_act_orb), &
+      d(n_act_orb,n_core_inact_act_orb,n_core_inact_act_orb))
+
   do j=1,mo_num
-    do k=1,n_core_inact_orb+n_act_orb
+    do l=1,n_core_inact_act_orb
-      do l=1,n_core_inact_orb+n_act_orb
+      do k=1,n_core_inact_act_orb
         do p=1,n_act_orb
-          d(p)=0.D0             
+            f(p,k,l) = bielec_PxxQ_no(list_act(p),k,l,j)
+        end do
+      end do
+    end do
+    call dgemm('T','N',n_act_orb,n_core_inact_act_orb**2,n_act_orb,1.d0,  &
+          natorbsCI, size(natorbsCI,1),                              &
+          f, n_act_orb,                                              &
+          0.d0,                                                      &
+          d, n_act_orb)
+    do l=1,n_core_inact_act_orb
+      do k=1,n_core_inact_act_orb
+        do p=1,n_act_orb
+          pp=n_act_orb-p+1
+          bielec_PxxQ_no(list_act(p),k,l,j)=d(pp,k,l)
+        end do
+      end do
+    end do
+  end do
+
+  deallocate (f,d)
+
+  allocate (f(n_act_orb,mo_num,n_core_inact_act_orb), &
+    d(n_act_orb,mo_num,n_core_inact_act_orb))
+
+  ! 3rd quarter
+  do k=1,mo_num
+    do l=1,n_core_inact_act_orb
+      do j=1,mo_num
+        do p=1,n_act_orb
+          f(p,j,l) = bielec_PxxQ_no(j,n_core_inact_orb+p,l,k)
+        end do
+      end do
+    end do
+    call dgemm('T','N',n_act_orb,mo_num*n_core_inact_act_orb,n_act_orb,1.d0,  &
+          natorbsCI, size(natorbsCI,1),                              &
+          f, n_act_orb,                                              &
+          0.d0,                                                      &
+          d, n_act_orb)
+    do l=1,n_core_inact_act_orb
+      do j=1,mo_num
+        do p=1,n_act_orb
+          pp=n_act_orb-p+1
+          bielec_PxxQ_no(j,n_core_inact_orb+p,l,k)=d(pp,j,l)
+        end do
+      end do
+    end do
+  end do
+
+  ! 4th quarter
+  do k=1,mo_num
+    do l=1,n_core_inact_act_orb
+      do j=1,mo_num
+        do p=1,n_act_orb
+          d(p,1,1)=0.D0
+        end do
+        do p=1,n_act_orb
+          do q=1,n_act_orb
+            d(p,1,1)+=bielec_PxxQ_no(j,l,n_core_inact_orb+q,k)*natorbsCI(q,p)
+          end do
         end do
         do p=1,n_act_orb
           pp=n_act_orb-p+1
-          do q=1,n_act_orb
+          bielec_PxxQ_no(j,l,n_core_inact_orb+p,k)=d(pp,1,1)
-            d(pp)+=bielec_PxxQ_no(list_act(q),k,l,j)*natorbsCI(q,p)
-          end do
-        end do
-        do p=1,n_act_orb
-          bielec_PxxQ_no(list_act(p),k,l,j)=d(p)
         end do
       end do
     end do
   end do
   ! 2nd quarter
+  do k=1,n_core_inact_act_orb
+    do l=1,n_core_inact_act_orb
       do j=1,mo_num
-    do k=1,n_core_inact_orb+n_act_orb
-      do l=1,n_core_inact_orb+n_act_orb
         do p=1,n_act_orb
-          d(p)=0.D0
+          d(p,1,1)=0.D0
+        end do
+        do p=1,n_act_orb
+          do q=1,n_act_orb
+            d(p,1,1)+=bielec_PxxQ_no(j,k,l,list_act(q))*natorbsCI(q,p)
+          end do
         end do
         do p=1,n_act_orb
           pp=n_act_orb-p+1
-          do q=1,n_act_orb
+          bielec_PxxQ_no(j,k,l,list_act(p))=d(pp,1,1)
-            d(pp)+=bielec_PxxQ_no(j,k,l,list_act(q))*natorbsCI(q,p)
-          end do
-        end do
-        do p=1,n_act_orb
-          bielec_PxxQ_no(j,k,l,list_act(p))=d(p)
-        end do
-      end do
-    end do
-  end do
-  ! 3rd quarter
-  do j=1,mo_num
-    do k=1,mo_num
-      do l=1,n_core_inact_orb+n_act_orb
-        do p=1,n_act_orb
-          d(p)=0.D0
-        end do
-        do p=1,n_act_orb
-          pp=n_act_orb-p+1
-          do q=1,n_act_orb
-            d(pp)+=bielec_PxxQ_no(j,n_core_inact_orb+q,l,k)*natorbsCI(q,p)
-          end do
-        end do
-        do p=1,n_act_orb
-          bielec_PxxQ_no(j,n_core_inact_orb+p,l,k)=d(p)
-        end do
-      end do
-    end do
-  end do
-  ! 4th quarter
-  do j=1,mo_num
-    do k=1,mo_num
-      do l=1,n_core_inact_orb+n_act_orb
-        do p=1,n_act_orb
-          d(p)=0.D0
-        end do
-        do p=1,n_act_orb
-          pp=n_act_orb-p+1
-          do q=1,n_act_orb
-            d(pp)+=bielec_PxxQ_no(j,l,n_core_inact_orb+q,k)*natorbsCI(q,p)
-          end do
-        end do
-        do p=1,n_act_orb
-          bielec_PxxQ_no(j,l,n_core_inact_orb+p,k)=d(p)
         end do
       end do
     end do
@@ -210,24 +228,27 @@ BEGIN_PROVIDER [real*8, bielecCI_no, (n_act_orb,n_act_orb,n_act_orb, mo_num)]
   END_DOC
   implicit none
   integer                        :: i,j,k,l,t,u,p,q,pp
-  real*8                         :: d(n_act_orb)
+  double precision, allocatable  :: f(:,:,:), d(:,:,:)
   
   bielecCI_no(:,:,:,:) = bielecCI(:,:,:,:)
 
+  allocate (f(n_act_orb,mo_num,n_act_orb), &
+      d(n_act_orb,mo_num,n_act_orb))
+
   do j=1,n_act_orb
     do k=1,n_act_orb
       do l=1,mo_num
         do p=1,n_act_orb
-          d(p)=0.D0
+          d(p,1,1)=0.D0
+        end do
+        do p=1,n_act_orb
+          do q=1,n_act_orb
+            d(p,1,1)+=bielecCI_no(q,j,k,l)*natorbsCI(q,p)
+          end do
         end do
         do p=1,n_act_orb
           pp=n_act_orb-p+1
-          do q=1,n_act_orb
+          bielecCI_no(p,j,k,l)=d(pp,1,1)
-            d(pp)+=bielecCI_no(q,j,k,l)*natorbsCI(q,p)
-          end do
-        end do
-        do p=1,n_act_orb
-          bielecCI_no(p,j,k,l)=d(p)
         end do
       end do
     end do
@@ -237,16 +258,16 @@ BEGIN_PROVIDER [real*8, bielecCI_no, (n_act_orb,n_act_orb,n_act_orb, mo_num)]
     do k=1,n_act_orb
       do l=1,mo_num
         do p=1,n_act_orb
-          d(p)=0.D0
+          d(p,1,1)=0.D0
+        end do
+        do p=1,n_act_orb
+          do q=1,n_act_orb
+            d(p,1,1)+=bielecCI_no(j,q,k,l)*natorbsCI(q,p)
+          end do
         end do
         do p=1,n_act_orb
           pp=n_act_orb-p+1
-          do q=1,n_act_orb
+          bielecCI_no(j,p,k,l)=d(pp,1,1)
-            d(pp)+=bielecCI_no(j,q,k,l)*natorbsCI(q,p)
-          end do
-        end do
-        do p=1,n_act_orb
-          bielecCI_no(j,p,k,l)=d(p)
         end do
       end do
     end do
@@ -256,16 +277,16 @@ BEGIN_PROVIDER [real*8, bielecCI_no, (n_act_orb,n_act_orb,n_act_orb, mo_num)]
     do k=1,n_act_orb
       do l=1,mo_num
         do p=1,n_act_orb
-          d(p)=0.D0
+          d(p,1,1)=0.D0
+        end do
+        do p=1,n_act_orb
+          do q=1,n_act_orb
+            d(p,1,1)+=bielecCI_no(j,k,q,l)*natorbsCI(q,p)
+          end do
         end do
         do p=1,n_act_orb
           pp=n_act_orb-p+1
-          do q=1,n_act_orb
+          bielecCI_no(j,k,p,l)=d(pp,1,1)
-            d(pp)+=bielecCI_no(j,k,q,l)*natorbsCI(q,p)
-          end do
-        end do
-        do p=1,n_act_orb
-          bielecCI_no(j,k,p,l)=d(p)
         end do
       end do 
     end do  
@@ -275,16 +296,16 @@ BEGIN_PROVIDER [real*8, bielecCI_no, (n_act_orb,n_act_orb,n_act_orb, mo_num)]
     do k=1,n_act_orb
       do l=1,n_act_orb
         do p=1,n_act_orb
-          d(p)=0.D0
+          d(p,1,1)=0.D0
+        end do
+        do p=1,n_act_orb
+          do q=1,n_act_orb
+            d(p,1,1)+=bielecCI_no(j,k,l,list_act(q))*natorbsCI(q,p)
+          end do
         end do
         do p=1,n_act_orb
           pp=n_act_orb-p+1
-          do q=1,n_act_orb
+          bielecCI_no(j,k,l,list_act(p))=d(pp,1,1)
-            d(pp)+=bielecCI_no(j,k,l,list_act(q))*natorbsCI(q,p)
-          end do
-        end do
-        do p=1,n_act_orb
-          bielecCI_no(j,k,l,list_act(p))=d(p)
         end do
       end do 
     end do