From 21dc0f53803d6df767a26d006f65ee570a54a937 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Wed, 3 Jul 2019 08:58:30 +0200 Subject: [PATCH] dgemm --- src/casscf/bielec_natorb.irp.f | 157 +++++++++++++++++++-------------- 1 file changed, 89 insertions(+), 68 deletions(-) diff --git a/src/casscf/bielec_natorb.irp.f b/src/casscf/bielec_natorb.irp.f index 9826d80c..691c7441 100644 --- a/src/casscf/bielec_natorb.irp.f +++ b/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_orb+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 - 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(pp)+=bielec_PxxQ_no(list_act(q),k,l,j)*natorbsCI(q,p) + 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 - bielec_PxxQ_no(list_act(p),k,l,j)=d(p) + pp=n_act_orb-p+1 + bielec_PxxQ_no(j,l,n_core_inact_orb+p,k)=d(pp,1,1) end do end do end do end do ! 2nd quarter - 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 k=1,n_core_inact_act_orb + do l=1,n_core_inact_act_orb + do j=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 - pp=n_act_orb-p+1 do q=1,n_act_orb - d(pp)+=bielec_PxxQ_no(j,k,l,list_act(q))*natorbsCI(q,p) + 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 - 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) + bielec_PxxQ_no(j,k,l,list_act(p))=d(pp,1,1) 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 - pp=n_act_orb-p+1 do q=1,n_act_orb - d(pp)+=bielecCI_no(q,j,k,l)*natorbsCI(q,p) + d(p,1,1)+=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) + pp=n_act_orb-p+1 + bielecCI_no(p,j,k,l)=d(pp,1,1) 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 - pp=n_act_orb-p+1 do q=1,n_act_orb - d(pp)+=bielecCI_no(j,q,k,l)*natorbsCI(q,p) + d(p,1,1)+=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) + pp=n_act_orb-p+1 + bielecCI_no(j,p,k,l)=d(pp,1,1) 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 - pp=n_act_orb-p+1 do q=1,n_act_orb - d(pp)+=bielecCI_no(j,k,q,l)*natorbsCI(q,p) + d(p,1,1)+=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) + pp=n_act_orb-p+1 + bielecCI_no(j,k,p,l)=d(pp,1,1) 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 - pp=n_act_orb-p+1 do q=1,n_act_orb - d(pp)+=bielecCI_no(j,k,l,list_act(q))*natorbsCI(q,p) + 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 - bielecCI_no(j,k,l,list_act(p))=d(p) + pp=n_act_orb-p+1 + bielecCI_no(j,k,l,list_act(p))=d(pp,1,1) end do end do end do