From cf361c396bdb4e10c9c4e37ded3d6f77be1abcc5 Mon Sep 17 00:00:00 2001 From: Francois Coppens Date: Thu, 14 Oct 2021 15:12:40 +0200 Subject: [PATCH] Replaced 1-update call to det_update() with call to qmckl_sherman_morrison_splitting() as a first step to make a general call to qmckl_sherman_morrison_smwb32s(). --- src/det.irp.f | 1244 +++++-------------------------------------------- src/qmckl | 2 +- 2 files changed, 127 insertions(+), 1119 deletions(-) diff --git a/src/det.irp.f b/src/det.irp.f index 6218438..a58c8ef 100644 --- a/src/det.irp.f +++ b/src/det.irp.f @@ -20,1109 +20,6 @@ END_PROVIDER END_PROVIDER - subroutine det_update(n,LDS,m,l,S,S_inv,d) - use qmckl - implicit none - - integer, intent(in) :: n,LDS ! Dimension of the vector - real, intent(in) :: m(LDS) ! New vector - integer, intent(in) :: l ! New position in S - - real,intent(inout) :: S(LDS,n) ! Slater matrix - double precision,intent(inout) :: S_inv(LDS,n) ! Inverse Slater matrix - double precision,intent(inout) :: d ! Det(S) - double precision :: zl, lambda, d_inv - - if (d == 0.d0) then - return - endif - -!! -- START QMCKL/SMW ADDITIONS -- !! -!! !! - integer :: i, j - integer (qmckl_exit_code) :: rc - integer (qmckl_context) :: context - integer(kind=8) :: ddim, nupdates, updates_index(1) - real(c_double) :: updates(n,1), breakdown - double precision :: S_inv_sq(n,n) - - S_inv = S_inv / d - - ! open(unit = 2000, file = "Slater_old.dat") - ! open(unit = 3000, file = "Slater_old_inv.dat") - ! do i=1,n - ! do j=1,n - ! write(2000,"(E23.15, 1X)", advance="no") S(j,i) ! write transpose for Octave - ! write(3000,"(E23.15, 1X)", advance="no") S_inv(i,j) - ! end do - ! write(2000,*) - ! write(3000,*) - ! end do - ! flush(2000) - ! flush(3000) - ! close(2000) - ! close(3000) - - do i=1,n - updates(i,1) = m(i) - S(i,l) ! transform repl. upds. into additive upds. - S(i,l) = m(i) ! update S with repl. upds - end do - - zl = 0 - do i=1,n - zl = zl + S_inv(i,l) * updates(i,1) - end do - - d_inv = 1.d0/d - d = d + zl - lambda = d * d_inv - if ( dabs(lambda) < 1.d-3 ) then - d = 0.d0 - return - endif - - context = qmckl_context_create() - ddim = n - nupdates = 1 - updates_index(1) = l - breakdown = 1e-3 - S_inv_sq = S_inv(1:n,1:n) - rc = qmckl_sherman_morrison_splitting(context, ddim, nupdates, updates, updates_index, breakdown, S_inv_sq) - S_inv = S_inv_sq(1:n,1:n) - rc = qmckl_context_destroy(context) - - ! open(unit = 4000, file = "Slater.dat") - ! open(unit = 5000, file = "Slater_inv.dat") - ! do i=1,n - ! do j=1,n - ! write(4000,"(E23.15, 1X)", advance="no") S(j,i) ! write transpose for Octave - ! write(5000,"(E23.15, 1X)", advance="no") S_inv(i,j) - ! end do - ! write(4000,*) - ! write(5000,*) - ! end do - ! flush(4000) - ! flush(5000) - ! close(4000) - ! close(5000) - - S_inv = S_inv * d - - ! select case (n) - ! case default - ! call det_update_general(n,LDS,m,l,S,S_inv,d) - ! BEGIN_TEMPLATE - ! case ($n) - ! call det_update$n(n,LDS,m,l,S,S_inv,d) - ! SUBST [n] - ! 1;; - ! 2;; - ! 3;; - ! 4;; - ! 5;; - ! 6;; - ! 7;; - ! 8;; - ! 9;; - ! 10;; - ! 11;; - ! 12;; - ! 13;; - ! 14;; - ! 15;; - ! 16;; - ! 17;; - ! 18;; - ! 19;; - ! 20;; - ! 21;; - ! 22;; - ! 23;; - ! 24;; - ! 25;; - ! 26;; - ! 27;; - ! 28;; - ! 29;; - ! 30;; - ! 31;; - ! 32;; - ! 33;; - ! 34;; - ! 35;; - ! 36;; - ! 37;; - ! 38;; - ! 39;; - ! 40;; - ! 41;; - ! 42;; - ! 43;; - ! 44;; - ! 45;; - ! 46;; - ! 47;; - ! 48;; - ! 49;; - ! 50;; - ! 51;; - ! 52;; - ! 53;; - ! 54;; - ! 55;; - ! 56;; - ! 57;; - ! 58;; - ! 59;; - ! 60;; - ! 61;; - ! 62;; - ! 63;; - ! 64;; - ! 65;; - ! 66;; - ! 67;; - ! 68;; - ! 69;; - ! 70;; - ! 71;; - ! 72;; - ! 73;; - ! 74;; - ! 75;; - ! 76;; - ! 77;; - ! 78;; - ! 79;; - ! 80;; - ! 81;; - ! 82;; - ! 83;; - ! 84;; - ! 85;; - ! 86;; - ! 87;; - ! 88;; - ! 89;; - ! 90;; - ! 91;; - ! 92;; - ! 93;; - ! 94;; - ! 95;; - ! 96;; - ! 97;; - ! 98;; - ! 99;; - ! 100;; - ! 101;; - ! 102;; - ! 103;; - ! 104;; - ! 105;; - ! 106;; - ! 107;; - ! 108;; - ! 109;; - ! 110;; - ! 111;; - ! 112;; - ! 113;; - ! 114;; - ! 115;; - ! 116;; - ! 117;; - ! 118;; - ! 119;; - ! 120;; - ! 121;; - ! 122;; - ! 123;; - ! 124;; - ! 125;; - ! 126;; - ! 127;; - ! 128;; - ! 129;; - ! 130;; - ! 131;; - ! 132;; - ! 133;; - ! 134;; - ! 135;; - ! 136;; - ! 137;; - ! 138;; - ! 139;; - ! 140;; - ! 141;; - ! 142;; - ! 143;; - ! 144;; - ! 145;; - ! 146;; - ! 147;; - ! 148;; - ! 149;; - ! 150;; - ! END_TEMPLATE - ! end select - -!! !! -!! -- END QMCKL/SMW ADDITIONS -- !! -end - -subroutine det_update2(n,LDS,m,l,S,S_inv,d) - implicit none - - integer, intent(in) :: n,LDS ! Dimension of the vector - real, intent(in) :: m(2) ! New vector - integer, intent(in) :: l ! New position in S - - real,intent(inout) :: S(LDS,2) ! Slater matrix - double precision,intent(inout) :: S_inv(LDS,2) ! Inverse Slater matrix - double precision,intent(inout) :: d ! Det(S) - - S(1,l) = m(1) - S(2,l) = m(2) - S_inv(1,1) = S(1,1) - S_inv(1,2) = S(2,1) - S_inv(2,1) = S(1,2) - S_inv(2,2) = S(2,2) - call invert2(S_inv,LDS,n,d) - -end - -subroutine det_update1(n,LDS,m,l,S,S_inv,d) - implicit none - - integer, intent(in) :: n,LDS ! Dimension of the vector - real, intent(in) :: m(1) ! New vector - integer, intent(in) :: l ! New position in S - - real,intent(inout) :: S(LDS,1) ! Slater matrix - double precision,intent(inout) :: S_inv(LDS,1) ! Inverse Slater matrix - double precision,intent(inout) :: d ! Det(S) - - S(1,l) = m(1) - S_inv(1,1) = S(1,1) - call invert1(S_inv,LDS,n,d) - -end - -subroutine det_update3(n,LDS,m,l,S,S_inv,d) - implicit none - - integer, intent(in) :: n,LDS ! Dimension of the vector - real, intent(in) :: m(3) ! New vector - integer, intent(in) :: l ! New position in S - - real,intent(inout) :: S(LDS,3) ! Slater matrix - double precision,intent(inout) :: S_inv(LDS,3) ! Inverse Slater matrix - double precision,intent(inout) :: d ! Det(S) - - integer :: i - do i=1,3 - S(i,l) = m(i) - enddo - do i=1,3 - S_inv(1,i) = S(i,1) - S_inv(2,i) = S(i,2) - S_inv(3,i) = S(i,3) - enddo - - call invert3(S_inv,LDS,n,d) - -end - -subroutine det_update4(n,LDS,m,l,S,S_inv,d) - implicit none - - integer, intent(in) :: n,LDS ! Dimension of the vector - real, intent(in) :: m(4) ! New vector - integer, intent(in) :: l ! New position in S - - real,intent(inout) :: S(LDS,4) ! Slater matrix - double precision,intent(inout) :: S_inv(LDS,4) ! Inverse Slater matrix - double precision,intent(inout) :: d ! Det(S) - - double precision :: u(4), z(4), w(4), lambda, d_inv - !DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: z, w, u - integer :: i,j - u(1) = m(1) - S(1,l) - u(2) = m(2) - S(2,l) - u(3) = m(3) - S(3,l) - u(4) = m(4) - S(4,l) - z(1) = S_inv(1,1)*u(1) + S_inv(2,1)*u(2) + S_inv(3,1)*u(3) + S_inv(4,1)*u(4) - z(2) = S_inv(1,2)*u(1) + S_inv(2,2)*u(2) + S_inv(3,2)*u(3) + S_inv(4,2)*u(4) - z(3) = S_inv(1,3)*u(1) + S_inv(2,3)*u(2) + S_inv(3,3)*u(3) + S_inv(4,3)*u(4) - z(4) = S_inv(1,4)*u(1) + S_inv(2,4)*u(2) + S_inv(3,4)*u(3) + S_inv(4,4)*u(4) - - d_inv = 1.d0/d - d = d+z(l) - lambda = d_inv*d - if (dabs(lambda) < 1.d-3) then - d = 0.d0 - return - endif - - !DIR$ VECTOR ALIGNED - do i=1,4 - w(i) = S_inv(i,l)*d_inv - S(i,l) = m(i) - enddo - - do i=1,4 - !DIR$ VECTOR ALIGNED - do j=1,4 - S_inv(j,i) = S_inv(j,i)*lambda -z(i)*w(j) - enddo - enddo - -end - -BEGIN_TEMPLATE -! Version for mod(n,4) = 0 -subroutine det_update$n(n,LDS,m,l,S,S_inv,d) - implicit none - - integer, intent(in) :: n,LDS ! Dimension of the vector - real, intent(in) :: m($n) ! New vector - integer, intent(in) :: l ! New position in S - - real,intent(inout) :: S(LDS,$n) ! Slater matrix - double precision,intent(inout) :: S_inv(LDS,$n) ! Inverse Slater matrix - double precision,intent(inout) :: d ! Det(S) - - double precision :: u($n), z($n), w($n), lambda, d_inv - !DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: z, w, u - !DIR$ ASSUME_ALIGNED S : $IRP_ALIGN - !DIR$ ASSUME_ALIGNED S_inv : $IRP_ALIGN - !DIR$ ASSUME (mod(LDS,$IRP_ALIGN/8) == 0) - !DIR$ ASSUME (LDS >= $n) - integer :: i,j - double precision :: zj, zj1, zj2, zj3 - - !DIR$ NOPREFETCH - do i=1,$n - u(i) = m(i) - S(i,l) - enddo - - zj = 0.d0 - !DIR$ VECTOR ALIGNED - !DIR$ NOPREFETCH - do i=1,$n-1,4 - zj = zj + S_inv(i,l)*u(i) + S_inv(i+1,l)*u(i+1) & - + S_inv(i+2,l)*u(i+2) + S_inv(i+3,l)*u(i+3) - enddo - - d_inv = 1.d0/d - d = d+zj - lambda = d*d_inv - if (dabs(lambda) < 1.d-3) then - d = 0.d0 - return - endif - - !DIR$ VECTOR ALIGNED - do j=1,$n,4 - zj = 0.d0 - zj1 = 0.d0 - zj2 = 0.d0 - zj3 = 0.d0 - !DIR$ VECTOR ALIGNED - !DIR$ NOPREFETCH - do i=1,$n - zj = zj + S_inv(i,j )*u(i) - zj1 = zj1 + S_inv(i,j+1)*u(i) - zj2 = zj2 + S_inv(i,j+2)*u(i) - zj3 = zj3 + S_inv(i,j+3)*u(i) - enddo - z(j ) = zj - z(j+1) = zj1 - z(j+2) = zj2 - z(j+3) = zj3 - enddo - - !DIR$ NOPREFETCH - do i=1,$n - w(i) = S_inv(i,l)*d_inv - S(i,l) = m(i) - enddo - - do i=1,$n,4 - zj = z(i ) - zj1 = z(i+1) - zj2 = z(i+2) - zj3 = z(i+3) - !DIR$ VECTOR ALIGNED - !DIR$ NOPREFETCH - do j=1,$n - S_inv(j,i ) = S_inv(j,i )*lambda - w(j)*zj - S_inv(j,i+1) = S_inv(j,i+1)*lambda - w(j)*zj1 - S_inv(j,i+2) = S_inv(j,i+2)*lambda - w(j)*zj2 - S_inv(j,i+3) = S_inv(j,i+3)*lambda - w(j)*zj3 - enddo - enddo - -end - -SUBST [ n ] -8 ;; -12 ;; -16 ;; -20 ;; -24 ;; -28 ;; -32 ;; -36 ;; -40 ;; -44 ;; -48 ;; -52 ;; -56 ;; -60 ;; -64 ;; -68 ;; -72 ;; -76 ;; -80 ;; -84 ;; -88 ;; -92 ;; -96 ;; -100 ;; -104 ;; -108 ;; -112 ;; -116 ;; -120 ;; -124 ;; -128 ;; -132 ;; -136 ;; -140 ;; -144 ;; -148 ;; - -END_TEMPLATE - -BEGIN_TEMPLATE -! Version for mod(n,4) = 1 -subroutine det_update$n(n,LDS,m,l,S,S_inv,d) - implicit none - - integer, intent(in) :: n,LDS ! Dimension of the vector - real, intent(in) :: m($n) ! New vector - integer, intent(in) :: l ! New position in S - - real,intent(inout) :: S(LDS,$n) ! Slater matrix - double precision,intent(inout) :: S_inv(LDS,$n) ! Inverse Slater matrix - double precision,intent(inout) :: d ! Det(S) - - double precision :: u($n), z($n), w($n), lambda, d_inv - !DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: z, w, u - !DIR$ ASSUME_ALIGNED S : $IRP_ALIGN - !DIR$ ASSUME_ALIGNED S_inv : $IRP_ALIGN - !DIR$ ASSUME (mod(LDS,$IRP_ALIGN/8) == 0) - !DIR$ ASSUME (LDS >= $n) - integer :: i,j - double precision :: zj, zj1, zj2, zj3 - - do i=1,$n - u(i) = m(i) - S(i,l) - enddo - - zj = 0.d0 - !DIR$ NOPREFETCH - do i=1,$n-1,4 - zj = zj + S_inv(i,l)*u(i) + S_inv(i+1,l)*u(i+1) & - + S_inv(i+2,l)*u(i+2) + S_inv(i+3,l)*u(i+3) - enddo - zj = zj + S_inv($n,l)*u($n) - - d_inv = 1.d0/d - d = d+zj - lambda = d*d_inv - if (dabs(lambda) < 1.d-3) then - d = 0.d0 - return - endif - - !DIR$ VECTOR ALIGNED - do j=1,$n-1,4 - zj = 0.d0 - zj1 = 0.d0 - zj2 = 0.d0 - zj3 = 0.d0 - !DIR$ VECTOR ALIGNED - !DIR$ NOPREFETCH - do i=1,$n-1 - zj = zj + S_inv(i,j )*u(i) - zj1 = zj1 + S_inv(i,j+1)*u(i) - zj2 = zj2 + S_inv(i,j+2)*u(i) - zj3 = zj3 + S_inv(i,j+3)*u(i) - enddo - z(j ) = zj + S_inv($n,j )*u($n) - z(j+1) = zj1 + S_inv($n,j+1)*u($n) - z(j+2) = zj2 + S_inv($n,j+2)*u($n) - z(j+3) = zj3 + S_inv($n,j+3)*u($n) - enddo - - zj = 0.d0 - !DIR$ VECTOR ALIGNED - !DIR$ NOPREFETCH - do i=1,$n-1 - zj = zj + S_inv(i,$n)*u(i) - enddo - z($n) = zj + S_inv($n,$n)*u($n) - - !DIR$ NOPREFETCH - do i=1,$n - w(i) = S_inv(i,l)*d_inv - S(i,l) = m(i) - enddo - - do i=1,$n-1,4 - zj = z(i ) - zj1 = z(i+1) - zj2 = z(i+2) - zj3 = z(i+3) - !DIR$ VECTOR ALIGNED - !DIR$ NOPREFETCH - do j=1,$n-1 - S_inv(j,i ) = S_inv(j,i )*lambda - w(j)*zj - S_inv(j,i+1) = S_inv(j,i+1)*lambda - w(j)*zj1 - S_inv(j,i+2) = S_inv(j,i+2)*lambda - w(j)*zj2 - S_inv(j,i+3) = S_inv(j,i+3)*lambda - w(j)*zj3 - enddo - S_inv($n,i ) = S_inv($n,i )*lambda - w($n)*zj - S_inv($n,i+1) = S_inv($n,i+1)*lambda - w($n)*zj1 - S_inv($n,i+2) = S_inv($n,i+2)*lambda - w($n)*zj2 - S_inv($n,i+3) = S_inv($n,i+3)*lambda - w($n)*zj3 - enddo - - zj = z($n) - !DIR$ VECTOR ALIGNED - !DIR$ NOPREFETCH - do i=1,$n - S_inv(i,$n) = S_inv(i,$n)*lambda -w(i)*zj - enddo - - -end - -SUBST [ n ] -5 ;; -9 ;; -13 ;; -17 ;; -21 ;; -25 ;; -29 ;; -33 ;; -37 ;; -41 ;; -45 ;; -49 ;; -53 ;; -57 ;; -61 ;; -65 ;; -69 ;; -73 ;; -77 ;; -81 ;; -85 ;; -89 ;; -93 ;; -97 ;; -101 ;; -105 ;; -109 ;; -113 ;; -117 ;; -121 ;; -125 ;; -129 ;; -133 ;; -137 ;; -141 ;; -145 ;; -149 ;; - -END_TEMPLATE - - -BEGIN_TEMPLATE -! Version for mod(n,4) = 2 -subroutine det_update$n(n,LDS,m,l,S,S_inv,d) - implicit none - - integer, intent(in) :: n,LDS ! Dimension of the vector - real, intent(in) :: m($n) ! New vector - integer, intent(in) :: l ! New position in S - - real,intent(inout) :: S(LDS,$n) ! Slater matrix - double precision,intent(inout) :: S_inv(LDS,$n) ! Inverse Slater matrix - double precision,intent(inout) :: d ! Det(S) - - double precision :: u($n), z($n), w($n), lambda, d_inv - !DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: z, w, u - !DIR$ ASSUME_ALIGNED S : $IRP_ALIGN - !DIR$ ASSUME_ALIGNED S_inv : $IRP_ALIGN - !DIR$ ASSUME (mod(LDS,$IRP_ALIGN/8) == 0) - !DIR$ ASSUME (LDS >= $n) - integer :: i,j - - double precision :: zj, zj1, zj2, zj3 - !DIR$ NOPREFETCH - do i=1,$n - u(i) = m(i) - S(i,l) - enddo - - zj = 0.d0 - !DIR$ VECTOR ALIGNED - !DIR$ NOPREFETCH - do i=1,$n-2,4 - zj = zj + S_inv(i,l)*u(i) + S_inv(i+1,l)*u(i+1) & - + S_inv(i+2,l)*u(i+2) + S_inv(i+3,l)*u(i+3) - enddo - i=$n-1 - zj = zj + S_inv(i,l)*u(i) + S_inv(i+1,l)*u(i+1) - - d_inv = 1.d0/d - d = d+zj - lambda = d*d_inv - if (dabs(lambda) < 1.d-3) then - d = 0.d0 - return - endif - - !DIR$ VECTOR ALIGNED - do j=1,$n-2,4 - zj = 0.d0 - zj1 = 0.d0 - zj2 = 0.d0 - zj3 = 0.d0 - !DIR$ VECTOR ALIGNED - do i=1,$n-2 - zj = zj + S_inv(i,j )*u(i) - zj1 = zj1 + S_inv(i,j+1)*u(i) - zj2 = zj2 + S_inv(i,j+2)*u(i) - zj3 = zj3 + S_inv(i,j+3)*u(i) - enddo - z(j ) = zj + S_inv($n-1,j )*u($n-1) - z(j ) = z(j ) + S_inv($n,j )*u($n) - z(j+1) = zj1 + S_inv($n-1,j+1)*u($n-1) - z(j+1) = z(j+1) + S_inv($n,j+1)*u($n) - z(j+2) = zj2 + S_inv($n-1,j+2)*u($n-1) - z(j+2) = z(j+2) + S_inv($n,j+2)*u($n) - z(j+3) = zj3 + S_inv($n-1,j+3)*u($n-1) - z(j+3) = z(j+3) + S_inv($n,j+3)*u($n) - enddo - - j=$n-1 - zj = 0.d0 - zj1 = 0.d0 - !DIR$ VECTOR ALIGNED - !DIR$ NOPREFETCH - do i=1,$n-2 - zj = zj + S_inv(i,j )*u(i) - zj1 = zj1 + S_inv(i,j+1)*u(i) - enddo - z(j ) = zj + S_inv($n-1,j )*u($n-1) - z(j ) = z(j ) + S_inv($n,j )*u($n) - z(j+1) = zj1 + S_inv($n-1,j+1)*u($n-1) - z(j+1) = z(j+1) + S_inv($n,j+1)*u($n) - - !DIR$ NOPREFETCH - do i=1,$n - w(i) = S_inv(i,l)*d_inv - S(i,l) = m(i) - enddo - - do i=1,$n-2,4 - zj = z(i) - zj1 = z(i+1) - zj2 = z(i+2) - zj3 = z(i+3) - !DIR$ VECTOR ALIGNED - do j=1,$n-2 - S_inv(j,i ) = S_inv(j,i )*lambda -zj *w(j) - S_inv(j,i+1) = S_inv(j,i+1)*lambda -zj1*w(j) - S_inv(j,i+2) = S_inv(j,i+2)*lambda -zj2*w(j) - S_inv(j,i+3) = S_inv(j,i+3)*lambda -zj3*w(j) - enddo - S_inv($n-1,i ) = S_inv($n-1,i )*lambda -zj *w($n-1) - S_inv($n ,i ) = S_inv($n ,i )*lambda -zj *w($n ) - S_inv($n-1,i+1) = S_inv($n-1,i+1)*lambda -zj1*w($n-1) - S_inv($n ,i+1) = S_inv($n ,i+1)*lambda -zj1*w($n ) - S_inv($n-1,i+2) = S_inv($n-1,i+2)*lambda -zj2*w($n-1) - S_inv($n ,i+2) = S_inv($n ,i+2)*lambda -zj2*w($n ) - S_inv($n-1,i+3) = S_inv($n-1,i+3)*lambda -zj3*w($n-1) - S_inv($n ,i+3) = S_inv($n ,i+3)*lambda -zj3*w($n ) - enddo - - i=$n-1 - zj = z(i) - zj1= z(i+1) - !DIR$ VECTOR ALIGNED - do j=1,$n-2 - S_inv(j,i ) = S_inv(j,i )*lambda -zj*w(j) - S_inv(j,i+1) = S_inv(j,i+1)*lambda -zj1*w(j) - enddo - S_inv($n-1,i ) = S_inv($n-1,i )*lambda -zj*w($n-1) - S_inv($n-1,i+1) = S_inv($n-1,i+1)*lambda -zj1*w($n-1) - S_inv($n ,i ) = S_inv($n ,i )*lambda -zj*w($n ) - S_inv($n ,i+1) = S_inv($n ,i+1)*lambda -zj1*w($n ) - -end - -SUBST [ n ] -6 ;; -10 ;; -14 ;; -18 ;; -22 ;; -26 ;; -30 ;; -34 ;; -38 ;; -42 ;; -46 ;; -50 ;; -54 ;; -58 ;; -62 ;; -66 ;; -70 ;; -74 ;; -78 ;; -82 ;; -86 ;; -90 ;; -94 ;; -98 ;; -102 ;; -106 ;; -110 ;; -114 ;; -118 ;; -122 ;; -126 ;; -130 ;; -134 ;; -138 ;; -142 ;; -146 ;; -150 ;; - -END_TEMPLATE - -BEGIN_TEMPLATE -! Version for mod(n,4) = 3 -subroutine det_update$n(n,LDS,m,l,S,S_inv,d) - implicit none - - integer, intent(in) :: n,LDS ! Dimension of the vector - real, intent(in) :: m($n) ! New vector - integer, intent(in) :: l ! New position in S - - real,intent(inout) :: S(LDS,$n) ! Slater matrix - double precision,intent(inout) :: S_inv(LDS,$n) ! Inverse Slater matrix - double precision,intent(inout) :: d ! Det(S) - - double precision :: u($n), z($n), w($n), lambda, d_inv - !DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: z, w, u - !DIR$ ASSUME_ALIGNED S : $IRP_ALIGN - !DIR$ ASSUME_ALIGNED S_inv : $IRP_ALIGN - !DIR$ ASSUME (mod(LDS,$IRP_ALIGN/8) == 0) - !DIR$ ASSUME (LDS >= $n) - integer :: i,j - - double precision :: zj, zj1, zj2, zj3 - - do i=1,$n - u(i) = m(i) - S(i,l) - enddo - - zj = 0.d0 - !DIR$ VECTOR ALIGNED - !DIR$ NOPREFETCH - do i=1,$n-3,4 - zj = zj + S_inv(i,l)*u(i) + S_inv(i+1,l)*u(i+1) & - + S_inv(i+2,l)*u(i+2) + S_inv(i+3,l)*u(i+3) - enddo - i=$n-2 - zj = zj + S_inv(i,l)*u(i) + S_inv(i+1,l)*u(i+1) + S_inv(i+2,l)*u(i+2) - - - d_inv = 1.d0/d - d = d+zj - lambda = d*d_inv - if (dabs(lambda) < 1.d-3) then - d = 0.d0 - return - endif - - !DIR$ VECTOR ALIGNED - do j=1,$n-3,4 - zj = 0.d0 - zj1 = 0.d0 - zj2 = 0.d0 - zj3 = 0.d0 - !DIR$ VECTOR ALIGNED - do i=1,$n-3 - zj = zj + S_inv(i,j )*u(i) - zj1 = zj1 + S_inv(i,j+1)*u(i) - zj2 = zj2 + S_inv(i,j+2)*u(i) - zj3 = zj3 + S_inv(i,j+3)*u(i) - enddo - z(j ) = zj + S_inv($n-2,j )*u($n-2) - z(j ) = z(j ) + S_inv($n-1,j )*u($n-1) - z(j ) = z(j ) + S_inv($n,j )*u($n) - z(j+1) = zj1 + S_inv($n-2,j+1)*u($n-2) - z(j+1) = z(j+1) + S_inv($n-1,j+1)*u($n-1) - z(j+1) = z(j+1) + S_inv($n,j+1)*u($n) - z(j+2) = zj2 + S_inv($n-2,j+2)*u($n-2) - z(j+2) = z(j+2) + S_inv($n-1,j+2)*u($n-1) - z(j+2) = z(j+2) + S_inv($n,j+2)*u($n) - z(j+3) = zj3 + S_inv($n-2,j+3)*u($n-2) - z(j+3) = z(j+3) + S_inv($n-1,j+3)*u($n-1) - z(j+3) = z(j+3) + S_inv($n,j+3)*u($n) - enddo - - j=$n-2 - zj = 0.d0 - zj1 = 0.d0 - zj2 = 0.d0 - !DIR$ VECTOR ALIGNED - !DIR$ NOPREFETCH - do i=1,$n-3 - zj = zj + S_inv(i,j )*u(i) - zj1 = zj1 + S_inv(i,j+1)*u(i) - zj2 = zj2 + S_inv(i,j+2)*u(i) - enddo - z(j ) = zj + S_inv($n-2,j )*u($n-2) - z(j ) = z(j ) + S_inv($n-1,j )*u($n-1) - z(j ) = z(j ) + S_inv($n,j )*u($n) - z(j+1) = zj1 + S_inv($n-2,j+1)*u($n-2) - z(j+1) = z(j+1) + S_inv($n-1,j+1)*u($n-1) - z(j+1) = z(j+1) + S_inv($n,j+1)*u($n) - z(j+2) = zj2 + S_inv($n-2,j+2)*u($n-2) - z(j+2) = z(j+2) + S_inv($n-1,j+2)*u($n-1) - z(j+2) = z(j+2) + S_inv($n,j+2)*u($n) - - !DIR$ NOPREFETCH - do i=1,$n - w(i) = S_inv(i,l)*d_inv - S(i,l) = m(i) - enddo - - do i=1,$n-3,4 - zj = z(i) - zj1 = z(i+1) - zj2 = z(i+2) - zj3 = z(i+3) - !DIR$ VECTOR ALIGNED - !DIR$ NOPREFETCH - do j=1,$n-3 - S_inv(j,i ) = S_inv(j,i )*lambda - w(j)*zj - S_inv(j,i+1) = S_inv(j,i+1)*lambda - w(j)*zj1 - S_inv(j,i+2) = S_inv(j,i+2)*lambda - w(j)*zj2 - S_inv(j,i+3) = S_inv(j,i+3)*lambda - w(j)*zj3 - enddo - S_inv($n-2,i ) = S_inv($n-2,i )*lambda -zj *w($n-2) - S_inv($n-1,i ) = S_inv($n-1,i )*lambda -zj *w($n-1) - S_inv($n ,i ) = S_inv($n ,i )*lambda -zj *w($n ) - S_inv($n-2,i+1) = S_inv($n-2,i+1)*lambda -zj1*w($n-2) - S_inv($n-1,i+1) = S_inv($n-1,i+1)*lambda -zj1*w($n-1) - S_inv($n ,i+1) = S_inv($n ,i+1)*lambda -zj1*w($n ) - S_inv($n-2,i+2) = S_inv($n-2,i+2)*lambda -zj2*w($n-2) - S_inv($n-1,i+2) = S_inv($n-1,i+2)*lambda -zj2*w($n-1) - S_inv($n ,i+2) = S_inv($n ,i+2)*lambda -zj2*w($n ) - S_inv($n-2,i+3) = S_inv($n-2,i+3)*lambda -zj3*w($n-2) - S_inv($n-1,i+3) = S_inv($n-1,i+3)*lambda -zj3*w($n-1) - S_inv($n ,i+3) = S_inv($n ,i+3)*lambda -zj3*w($n ) - enddo - - i=$n-2 - zj = z(i) - zj1 = z(i+1) - zj2 = z(i+2) - !DIR$ VECTOR ALIGNED - !DIR$ NOPREFETCH - do j=1,$n - S_inv(j,i ) = S_inv(j,i )*lambda - w(j)*zj - S_inv(j,i+1) = S_inv(j,i+1)*lambda - w(j)*zj1 - S_inv(j,i+2) = S_inv(j,i+2)*lambda - w(j)*zj2 - enddo - - -end - -SUBST [ n ] -7 ;; -11 ;; -15 ;; -19 ;; -23 ;; -27 ;; -31 ;; -35 ;; -39 ;; -43 ;; -47 ;; -51 ;; -55 ;; -59 ;; -63 ;; -67 ;; -71 ;; -75 ;; -79 ;; -83 ;; -87 ;; -91 ;; -95 ;; -99 ;; -103 ;; -107 ;; -111 ;; -115 ;; -119 ;; -123 ;; -127 ;; -131 ;; -135 ;; -139 ;; -143 ;; -147 ;; - -END_TEMPLATE - - - -subroutine det_update_general(n,LDS,m,l,S,S_inv,d) - implicit none - - integer, intent(in) :: n,LDS ! Dimension of the vector - real, intent(in) :: m(LDS) ! New vector - integer, intent(in) :: l ! New position in S - - real,intent(inout) :: S(LDS,n) ! Slater matrix - double precision,intent(inout) :: S_inv(LDS,n) ! Inverse Slater matrix - double precision,intent(inout) :: d ! Det(S) - - double precision :: lambda, d_inv - double precision :: u(3840), z(3840), w(3840) - !DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: z, w, u - !DIR$ ASSUME_ALIGNED S : $IRP_ALIGN - !DIR$ ASSUME_ALIGNED S_inv : $IRP_ALIGN - !DIR$ ASSUME (LDS >= n) - !DIR$ ASSUME (LDS <= 3840) - !DIR$ ASSUME (MOD(LDS,$IRP_ALIGN/8) == 0) - !DIR$ ASSUME (n>150) - - integer :: i,j,n4 - double precision :: zl - - !DIR$ NOPREFETCH - do i=1,n - u(i) = m(i) - S(i,l) - enddo - - zl = 0.d0 - !DIR$ VECTOR ALIGNED - !DIR$ NOPREFETCH - do i=1,n - zl = zl + S_inv(i,l)*u(i) - enddo - - d_inv = 1.d0/d - d = d+zl - lambda = d*d_inv -! - if ( dabs(lambda) < 1.d-3 ) then - d = 0.d0 - return - endif - - double precision :: zj, zj1, zj2, zj3 - - n4 = iand(n,not(3)) - !DIR$ VECTOR ALIGNED - !DIR$ NOPREFETCH - do j=1,n4,4 - zj = 0.d0 - zj1 = 0.d0 - zj2 = 0.d0 - zj3 = 0.d0 - !DIR$ VECTOR ALIGNED - !DIR$ NOPREFETCH - do i=1,n - zj = zj + S_inv(i,j )*u(i) - zj1 = zj1 + S_inv(i,j+1)*u(i) - zj2 = zj2 + S_inv(i,j+2)*u(i) - zj3 = zj3 + S_inv(i,j+3)*u(i) - enddo - z(j ) = zj - z(j+1) = zj1 - z(j+2) = zj2 - z(j+3) = zj3 - enddo - - do j=n4+1,n - zj = 0.d0 - !DIR$ VECTOR ALIGNED - !DIR$ NOPREFETCH - do i=1,n - zj = zj + S_inv(i,j)*u(i) - enddo - z(j ) = zj - enddo - - !DIR$ NOPREFETCH - do i=1,n - w(i) = S_inv(i,l)*d_inv - S(i,l) = m(i) - enddo - - !DIR$ NOPREFETCH - do i=1,n - w(i) = S_inv(i,l)*d_inv - S(i,l) = m(i) - enddo - - do i=1,n4,4 - zj = z(i) - zj1 = z(i+1) - zj2 = z(i+2) - zj3 = z(i+3) - !DIR$ VECTOR ALIGNED - !DIR$ NOPREFETCH - do j=1,n - S_inv(j,i ) = S_inv(j,i )*lambda -zj *w(j) - S_inv(j,i+1) = S_inv(j,i+1)*lambda -zj1*w(j) - S_inv(j,i+2) = S_inv(j,i+2)*lambda -zj2*w(j) - S_inv(j,i+3) = S_inv(j,i+3)*lambda -zj3*w(j) - enddo - enddo - - do i=n4+1,n - zj = z(i) - !DIR$ VECTOR ALIGNED - !DIR$ NOPREFETCH - do j=1,n - S_inv(j,i) = S_inv(j,i)*lambda -zj*w(j) - enddo - enddo - -end - - - subroutine bitstring_to_list( string, list, n_elements, Nint) implicit none BEGIN_DOC @@ -1271,6 +168,8 @@ END_PROVIDER ! slater_matrix_alpha_inv_det: Inverse of the alpha Slater matrix x determinant END_DOC + use qmckl + double precision :: ddet integer :: i,j,k,imo,l integer :: to_do(elec_alpha_num), n_to_do_old, n_to_do @@ -1278,6 +177,12 @@ END_PROVIDER real :: tmp_det(elec_alpha_num_8) integer, save :: ifirst logical :: file_exists + + integer (qmckl_exit_code) :: rc + integer (qmckl_context) :: context + integer(kind=8) :: ddim, nupdates, updates_index(1) + real(c_double) :: updates(elec_alpha_num,1), breakdown + !DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: tmp_inv, tmp_det if (ifirst == 0) then @@ -1373,13 +278,61 @@ END_PROVIDER do l=1,n_to_do_old k = to_do(l) ! select electron to change imo = mo_list_alpha_curr(k) ! select mo to change - ! write(*,*) "k, imo, mo_value(1,imo) = ", k, imo, mo_value(1,imo) - call det_update(elec_alpha_num, elec_alpha_num_8, & - mo_value(1,imo), & - k, & - slater_matrix_alpha, & - slater_matrix_alpha_inv_det, & - ddet) + + slater_matrix_alpha_inv_det = slater_matrix_alpha_inv_det / ddet + + open(unit = 2000, file = "SlaterA_old.dat") + open(unit = 3000, file = "SlaterA_old_inv.dat") + do i=1, elec_alpha_num + do j=1, elec_alpha_num + write(2000,"(E23.15, 1X)", advance="no") slater_matrix_alpha(j,i) ! avoid transpose later + write(3000,"(E23.15, 1X)", advance="no") slater_matrix_alpha_inv_det(i,j) + end do + write(2000,*) + write(3000,*) + end do + flush(2000) + flush(3000) + close(2000) + close(3000) + + do i = 1,elec_alpha_num + updates(i,1) = mo_value(i, imo) - slater_matrix_alpha(i, imo) + slater_matrix_alpha(i,imo) = mo_value(i, imo) + end do + + context = qmckl_context_create() + ddim = elec_alpha_num + nupdates = 1 + updates_index(1) = imo + breakdown = 1e-3 + rc = qmckl_sherman_morrison_splitting(context, & + ddim, & + nupdates, & + updates, & + updates_index, & + breakdown, & + slater_matrix_alpha_inv_det(1:elec_alpha_num, 1:elec_alpha_num), & + ddet) + rc = qmckl_context_destroy(context) + + open(unit = 4000, file = "SlaterA.dat") + open(unit = 5000, file = "SlaterA_inv.dat") + do i=1, elec_alpha_num + do j=1, elec_alpha_num + write(4000,"(E23.15, 1X)", advance="no") slater_matrix_alpha(j,i) ! avoid transpose later + write(5000,"(E23.15, 1X)", advance="no") slater_matrix_alpha_inv_det(i,j) + end do + write(4000,*) + write(5000,*) + end do + flush(4000) + flush(5000) + close(4000) + close(5000) + + slater_matrix_alpha_inv_det = slater_matrix_alpha_inv_det * ddet + if (ddet /= 0.d0) then det_alpha_value_curr = ddet else @@ -1445,6 +398,8 @@ END_PROVIDER ! slater_matrix_beta_inv_det : Inverse of the beta Slater matrix x determinant END_DOC + use qmckl + double precision :: ddet integer :: i,j,k,imo,l integer :: to_do(elec_alpha_num-mo_closed_num), n_to_do_old, n_to_do @@ -1452,6 +407,11 @@ END_PROVIDER real :: tmp_det(elec_alpha_num_8) integer, save :: ifirst + integer (qmckl_exit_code) :: rc + integer (qmckl_context) :: context + integer(kind=8) :: ddim, nupdates, updates_index(1) + real(c_double) :: updates(elec_beta_num,1), breakdown + if (elec_beta_num == 0) then det_beta_value_curr = 0.d0 return @@ -1547,15 +507,64 @@ END_PROVIDER n_to_do_old = n_to_do n_to_do = 0 loopcount = 0 - do l=1,n_to_do_old + do l= 1,n_to_do_old k = to_do(l) imo = mo_list_beta_curr(k) - call det_update(elec_beta_num, elec_beta_num_8, & - mo_value(elec_alpha_num+1,imo), & - k, & - slater_matrix_beta, & - slater_matrix_beta_inv_det, & - ddet) + slater_matrix_beta_inv_det = slater_matrix_beta_inv_det / ddet + + open(unit = 2000, file = "SlaterB_old.dat") + open(unit = 3000, file = "SlaterB_old_inv.dat") + do i=1, elec_beta_num + do j=1, elec_beta_num + write(2000,"(E23.15, 1X)", advance="no") slater_matrix_beta(j,i) ! avoid transpose later + write(3000,"(E23.15, 1X)", advance="no") slater_matrix_beta_inv_det(i,j) + end do + write(2000,*) + write(3000,*) + end do + flush(2000) + flush(3000) + close(2000) + close(3000) + + do i = 1,elec_beta_num + updates(i,1) = mo_value(elec_alpha_num + i, imo) - slater_matrix_beta(i, imo) + slater_matrix_beta(i,imo) = mo_value(elec_alpha_num + i, imo) + end do + + context = qmckl_context_create() + ddim = elec_beta_num + nupdates = 1 + updates_index(1) = imo + breakdown = 1e-3 + rc = qmckl_sherman_morrison_splitting(context, & + ddim, & + nupdates, & + updates, & + updates_index, & + breakdown, & + slater_matrix_beta_inv_det(1:elec_beta_num, 1:elec_beta_num), & + ddet) + rc = qmckl_context_destroy(context) + + open(unit = 4000, file = "SlaterB.dat") + open(unit = 5000, file = "SlaterB_inv.dat") + do i=1, elec_beta_num + do j=1, elec_beta_num + write(4000,"(E23.15, 1X)", advance="no") slater_matrix_beta(j,i) ! avoid transpose later + write(5000,"(E23.15, 1X)", advance="no") slater_matrix_beta_inv_det(i,j) + end do + write(4000,*) + write(5000,*) + end do + flush(4000) + flush(5000) + close(4000) + close(5000) + + slater_matrix_beta_inv_det = slater_matrix_beta_inv_det * ddet + stop + if (ddet /= 0.d0) then det_beta_value_curr = ddet else @@ -2248,4 +1257,3 @@ BEGIN_PROVIDER [ double precision, psidet_lapl ] psidet_lapl = psidet_lapl + psidet_grad_lapl(4,j) enddo END_PROVIDER - diff --git a/src/qmckl b/src/qmckl index 8535e3c..9c2d01b 160000 --- a/src/qmckl +++ b/src/qmckl @@ -1 +1 @@ -Subproject commit 8535e3c0b5a4f4194ca0400b348774c3cc7c74d0 +Subproject commit 9c2d01b33d4c972945ffc198a4aa853f70b4ef1e