TREX/qmckl
1
0
Fork 0
mirror of https://github.com/TREX-CoE/qmckl.git synced 2025-01-03 18:16:28 +01:00
Code Issues Projects Releases Wiki Activity

Added test for invert. #41

Browse Source
This commit is contained in:
v1j4y 2021-10-07 14:01:40 +02:00
parent 2735b31c12
commit d24f268369
2 changed files with 110 additions and 39 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

74
org/qmckl_blas.org
View File

@ -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:

19
org/qmckl_mo.org
View File

@ -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,13 +490,12 @@ 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, &
ao_vgl(:, :, iwalk, 1), LDA, &
coef_normalized(1:ao_num,1:mo_num),size(coef_normalized,1) * 1_8, &
beta, &
mo_vgl(:,ielec,iwalk,1),LDC)
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, &
@ -522,7 +521,6 @@ integer function qmckl_compute_mo_basis_gaussian_vgl_f(context, &
beta, &
mo_vgl(:,ielec,iwalk,5),LDC)
end do
end do
if(info_qmckl_dgemm_value .eq. QMCKL_SUCCESS .and. &
info_qmckl_dgemm_Gx .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<npoints;++i) {
// printf(".");
// fflush(stdout);
// for (int j=0;j<npoints;++j) {
// printf(" .. ");
// for (int k=0;k<npoints;++k) {
// printf(" . ");
// // Set point
// elec_coord[0] = point_x[i];
// elec_coord[1] = point_y[j];