diff --git a/org/qmckl_blas.org b/org/qmckl_blas.org index 293af48..46ae91a 100644 --- a/org/qmckl_blas.org +++ b/org/qmckl_blas.org @@ -365,8 +365,10 @@ integer function qmckl_invert_f(context, ma, na, LDA, A, det_l) & integer*8, intent(in) :: ma integer*8, intent(in) :: na double precision, intent(inout) :: det_l - integer :: i,j + + info = QMCKL_SUCCESS + select case (na) case default !DIR$ forceinline @@ -742,6 +744,76 @@ end #+end_src *** Test :noexport: + #+begin_src f90 :tangle (eval f_test) +integer(qmckl_exit_code) function test_qmckl_invert(context) bind(C) + use qmckl + implicit none + integer(qmckl_context), intent(in), value :: context + + double precision, allocatable :: A(:,:), C(:,:) + integer*8 :: m, n, k, LDA, LDB, LDC + integer*8 :: i,j,l + double precision :: x, det_l, det_l_ref + + m = 4_8 + k = 4_8 + LDA = m + LDB = m + LDC = m + + allocate( A(LDA,k), C(LDC,k)) + + A = 0.10d0 + C = 0.d0 + A(1,1) = 1.0d0; + A(2,2) = 2.0d0; + A(3,3) = 3.0d0; + A(4,4) = 4.0d0; + + ! Exact inverse (Mathematica) + C(1,1) = 1.0102367161391992d0 + C(2,2) = 0.5036819224578257d0 + C(3,3) = 0.33511197860555897d0 + C(4,4) = 0.2510382472105688d0 + C(1,2) = -0.047782608144589914d0 + C(1,3) = -0.031305846715420985d0 + C(1,4) = -0.023278706531979707d0 + C(2,3) = -0.014829085286252043d0 + C(2,4) = -0.011026755725674596d0 + C(3,4) = -0.007224426165097149d0 + C(2,1) = -0.047782608144589914d0 + C(3,1) = -0.031305846715420985d0 + C(4,1) = -0.023278706531979707d0 + C(3,2) = -0.014829085286252043d0 + C(4,2) = -0.011026755725674596d0 + C(4,3) = -0.007224426165097149d0 + det_l_ref = 23.6697d0 + + test_qmckl_invert = qmckl_invert(context, m, k, LDA, A, det_l) + + if (test_qmckl_invert /= QMCKL_SUCCESS) return + + test_qmckl_invert = QMCKL_FAILURE + + x = 0.d0 + do j=1,m + do i=1,k + x = x + (A(i,j) - (C(i,j) * det_l_ref))**2 + end do + end do + + if (dabs(x) <= 1.d-15 .and. (dabs(det_l_ref - det_l)) <= 1.d-15) then + test_qmckl_invert = QMCKL_SUCCESS + endif + + deallocate(A,C) +end function test_qmckl_invert + #+end_src + + #+begin_src c :comments link :tangle (eval c_test) +qmckl_exit_code test_qmckl_invert(qmckl_context context); +assert(QMCKL_SUCCESS == test_qmckl_invert(context)); + #+end_src * End of files :noexport: diff --git a/org/qmckl_mo.org b/org/qmckl_mo.org index 8fa209f..012caae 100644 --- a/org/qmckl_mo.org +++ b/org/qmckl_mo.org @@ -482,7 +482,7 @@ integer function qmckl_compute_mo_basis_gaussian_vgl_f(context, & ! Don't compute exponentials when the result will be almost zero. ! TODO : Use numerical precision here cutoff = -dlog(1.d-15) - M = 1_8 + M = elec_num N = mo_num * 1_8 K = ao_num * 1_8 LDA = M @@ -490,38 +490,36 @@ integer function qmckl_compute_mo_basis_gaussian_vgl_f(context, & LDC = M do iwalk = 1, walk_num - do ielec = 1, elec_num - ! Value - info_qmckl_dgemm_value = qmckl_dgemm(context,TransA, TransB, M, N, K, alpha, & - ao_vgl(:, ielec, iwalk, 1), LDA, & - coef_normalized(1:ao_num,1:mo_num),size(coef_normalized,1) * 1_8, & - beta, & - mo_vgl(:,ielec,iwalk,1),LDC) - ! Grad_x - info_qmckl_dgemm_Gx = qmckl_dgemm(context,TransA, TransB, M, N, K, alpha, & - ao_vgl(:, ielec, iwalk, 2), LDA, & - coef_normalized(1:ao_num,1:mo_num),size(coef_normalized,1) * 1_8, & - beta, & - mo_vgl(:,ielec,iwalk,2),LDC) - ! Grad_y - info_qmckl_dgemm_Gy = qmckl_dgemm(context,TransA, TransB, M, N, K, alpha, & - ao_vgl(:, ielec, iwalk, 3), LDA, & - coef_normalized(1:ao_num,1:mo_num),size(coef_normalized,1) * 1_8, & - beta, & - mo_vgl(:,ielec,iwalk,3),LDC) - ! Grad_z - info_qmckl_dgemm_Gz = qmckl_dgemm(context,TransA, TransB, M, N, K, alpha, & - ao_vgl(:, ielec, iwalk, 4), LDA, & - coef_normalized(1:ao_num,1:mo_num),size(coef_normalized,1) * 1_8, & - beta, & - mo_vgl(:,ielec,iwalk,4),LDC) - ! Lapl_z - info_qmckl_dgemm_lap = qmckl_dgemm(context, TransA, TransB, M, N, K, alpha, & - ao_vgl(:, ielec, iwalk, 5), LDA, & - coef_normalized(1:ao_num,1:mo_num),size(coef_normalized,1) * 1_8, & - beta, & - mo_vgl(:,ielec,iwalk,5),LDC) - end do + ! Value + info_qmckl_dgemm_value = qmckl_dgemm(context,TransA, TransB, M, N, K, alpha, & + ao_vgl(:, :, iwalk, 1), LDA, & + coef_normalized(1:ao_num,1:mo_num),size(coef_normalized,1) * 1_8, & + beta, & + mo_vgl(:,:,iwalk,1),LDC) + ! Grad_x + info_qmckl_dgemm_Gx = qmckl_dgemm(context,TransA, TransB, M, N, K, alpha, & + ao_vgl(:, ielec, iwalk, 2), LDA, & + coef_normalized(1:ao_num,1:mo_num),size(coef_normalized,1) * 1_8, & + beta, & + mo_vgl(:,ielec,iwalk,2),LDC) + ! Grad_y + info_qmckl_dgemm_Gy = qmckl_dgemm(context,TransA, TransB, M, N, K, alpha, & + ao_vgl(:, ielec, iwalk, 3), LDA, & + coef_normalized(1:ao_num,1:mo_num),size(coef_normalized,1) * 1_8, & + beta, & + mo_vgl(:,ielec,iwalk,3),LDC) + ! Grad_z + info_qmckl_dgemm_Gz = qmckl_dgemm(context,TransA, TransB, M, N, K, alpha, & + ao_vgl(:, ielec, iwalk, 4), LDA, & + coef_normalized(1:ao_num,1:mo_num),size(coef_normalized,1) * 1_8, & + beta, & + mo_vgl(:,ielec,iwalk,4),LDC) + ! Lapl_z + info_qmckl_dgemm_lap = qmckl_dgemm(context, TransA, TransB, M, N, K, alpha, & + ao_vgl(:, ielec, iwalk, 5), LDA, & + coef_normalized(1:ao_num,1:mo_num),size(coef_normalized,1) * 1_8, & + beta, & + mo_vgl(:,ielec,iwalk,5),LDC) end do if(info_qmckl_dgemm_value .eq. QMCKL_SUCCESS .and. & @@ -778,10 +776,10 @@ rc = qmckl_get_mo_basis_vgl(context, &(mo_vgl[0][0][0][0])); assert (rc == QMCKL_SUCCESS); // Test overlap of MO -//double point_x[100]; -//double point_y[100]; -//double point_z[100]; -//int32_t npoints=100; +//double point_x[10]; +//double point_y[10]; +//double point_z[10]; +//int32_t npoints=10; //// obtain points //double dr = 20./(npoints-1); //double dr3 = dr*dr*dr; @@ -795,10 +793,11 @@ assert (rc == QMCKL_SUCCESS); //double ovlmo1 = 0.0; //// Calculate overlap //for (int i=0;i