diff --git a/org/qmckl_ao.org b/org/qmckl_ao.org index 26f87f6..6c04a00 100644 --- a/org/qmckl_ao.org +++ b/org/qmckl_ao.org @@ -601,7 +601,7 @@ qmckl_set_ao_basis_nucleus_shell_num (qmckl_context context, if (size_max < nucl_num) { return qmckl_failwith( context, - QMCKL_FAILURE, + QMCKL_INVALID_ARG_3, "qmckl_set_ao_basis_nucleus_shell_num", "input array too small"); } @@ -663,7 +663,7 @@ qmckl_set_ao_basis_shell_ang_mom (qmckl_context context, if (size_max < shell_num) { return qmckl_failwith( context, - QMCKL_FAILURE, + QMCKL_INVALID_ARG_3, "qmckl_set_ao_basis_shell_ang_mom", "input array too small"); } @@ -726,7 +726,7 @@ qmckl_set_ao_basis_shell_prim_num (qmckl_context context, if (size_max < shell_num) { return qmckl_failwith( context, - QMCKL_FAILURE, + QMCKL_INVALID_ARG_3, "qmckl_set_ao_basis_shell_prim_num", "input array too small"); } @@ -789,7 +789,7 @@ qmckl_set_ao_basis_shell_prim_index (qmckl_context context, if (size_max < shell_num) { return qmckl_failwith( context, - QMCKL_FAILURE, + QMCKL_INVALID_ARG_3, "qmckl_set_ao_basis_shell_prim_index", "input array too small"); } @@ -851,7 +851,7 @@ qmckl_set_ao_basis_shell_factor (qmckl_context context, if (size_max < shell_num) { return qmckl_failwith( context, - QMCKL_FAILURE, + QMCKL_INVALID_ARG_3, "qmckl_set_ao_basis_shell_factor", "input array too small"); } @@ -913,7 +913,7 @@ qmckl_set_ao_basis_exponent (qmckl_context context, if (size_max < prim_num) { return qmckl_failwith( context, - QMCKL_FAILURE, + QMCKL_INVALID_ARG_3, "qmckl_set_ao_basis_exponent", "input array too small"); } @@ -975,7 +975,7 @@ qmckl_set_ao_basis_coefficient (qmckl_context context, if (size_max < prim_num) { return qmckl_failwith( context, - QMCKL_FAILURE, + QMCKL_INVALID_ARG_3, "qmckl_set_ao_basis_coefficient", "input array too small"); } @@ -1037,7 +1037,7 @@ qmckl_set_ao_basis_prim_factor (qmckl_context context, if (size_max < prim_num) { return qmckl_failwith( context, - QMCKL_FAILURE, + QMCKL_INVALID_ARG_3, "qmckl_set_ao_basis_prim_factor", "input array too small"); } @@ -1099,7 +1099,7 @@ qmckl_set_ao_basis_ao_factor (qmckl_context context, if (size_max < ao_num) { return qmckl_failwith( context, - QMCKL_FAILURE, + QMCKL_INVALID_ARG_3, "qmckl_set_ao_basis_ao_factor", "input array too small"); } @@ -2638,6 +2638,13 @@ qmckl_get_ao_basis_primitive_vgl (qmckl_context context, NULL); } + if (size_max <= 0) { + return qmckl_failwith( context, + QMCKL_INVALID_ARG_3, + "qmckl_get_ao_basis_primitive_vgl", + "size_max <= 0"); + } + qmckl_exit_code rc; rc = qmckl_provide_ao_basis_primitive_vgl(context); @@ -2646,14 +2653,14 @@ qmckl_get_ao_basis_primitive_vgl (qmckl_context context, qmckl_context_struct* const ctx = (qmckl_context_struct* const) context; assert (ctx != NULL); - size_t sze = ctx->ao_basis.prim_num * 5 * ctx->electron.num; + int64_t sze = ctx->ao_basis.prim_num * 5 * ctx->electron.num; if (size_max < sze) { return qmckl_failwith( context, - QMCKL_FAILURE, + QMCKL_INVALID_ARG_3, "qmckl_get_ao_basis_primitive_vgl", "input array too small"); } - memcpy(primitive_vgl, ctx->ao_basis.primitive_vgl, sze * sizeof(double)); + memcpy(primitive_vgl, ctx->ao_basis.primitive_vgl, (size_t) sze * sizeof(double)); return QMCKL_SUCCESS; } @@ -2707,14 +2714,14 @@ qmckl_get_ao_basis_shell_vgl (qmckl_context context, qmckl_context_struct* const ctx = (qmckl_context_struct* const) context; assert (ctx != NULL); - size_t sze = ctx->ao_basis.shell_num * 5 * ctx->electron.num; + int64_t sze = ctx->ao_basis.shell_num * 5 * ctx->electron.num; if (size_max < sze) { return qmckl_failwith( context, - QMCKL_FAILURE, + QMCKL_INVALID_ARG_3, "qmckl_get_ao_basis_shell_vgl", "input array too small"); } - memcpy(shell_vgl, ctx->ao_basis.shell_vgl, sze * sizeof(double)); + memcpy(shell_vgl, ctx->ao_basis.shell_vgl, (size_t)sze * sizeof(double)); return QMCKL_SUCCESS; } @@ -2770,14 +2777,14 @@ qmckl_get_ao_basis_ao_vgl (qmckl_context context, qmckl_context_struct* const ctx = (qmckl_context_struct* const) context; assert (ctx != NULL); - size_t sze = ctx->ao_basis.ao_num * 5 * ctx->electron.num; + int64_t sze = ctx->ao_basis.ao_num * 5 * ctx->electron.num; if (size_max < sze) { return qmckl_failwith( context, - QMCKL_FAILURE, + QMCKL_INVALID_ARG_3, "qmckl_get_ao_basis_ao_vgl", "input array too small"); } - memcpy(ao_vgl, ctx->ao_basis.ao_vgl, sze * sizeof(double)); + memcpy(ao_vgl, ctx->ao_basis.ao_vgl, (size_t) sze * sizeof(double)); return QMCKL_SUCCESS; } diff --git a/org/qmckl_blas.org b/org/qmckl_blas.org index a07c816..a424082 100644 --- a/org/qmckl_blas.org +++ b/org/qmckl_blas.org @@ -22,12 +22,12 @@ int main() { * Matrix operations ** ~qmckl_dgemm~ - + Matrix multiplication: \[ C_{ij} = \beta C_{ij} + \alpha \sum_{k} A_{ik} \cdot B_{kj} - \] + \] # TODO: Add description about the external library dependence. @@ -48,7 +48,7 @@ int main() { | ~beta~ | ~double~ | in | Array containing the matrix $B$ | | ~C~ | ~double[][ldc]~ | out | Array containing the matrix $B$ | | ~ldc~ | ~int64_t~ | in | Leading dimension of array ~B~ | - + Requirements: - ~context~ is not ~QMCKL_NULL_CONTEXT~ @@ -61,7 +61,7 @@ int main() { - ~A~ is allocated with at least $m \times k \times 8$ bytes - ~B~ is allocated with at least $k \times n \times 8$ bytes - ~C~ is allocated with at least $m \times n \times 8$ bytes - + #+CALL: generate_c_header(table=qmckl_dgemm_args,rettyp="qmckl_exit_code",fname="qmckl_dgemm") #+RESULTS: @@ -80,7 +80,7 @@ int main() { const int64_t ldb, const double beta, double* const C, - const int64_t ldc ); + const int64_t ldc ); #+end_src #+begin_src f90 :tangle (eval f) @@ -100,28 +100,28 @@ integer function qmckl_dgemm_f(context, TransA, TransB, & integer*8 , intent(in) :: ldc double precision , intent(out) :: C(ldc,*) - info = QMCKL_SUCCESS + info = QMCKL_SUCCESS if (context == QMCKL_NULL_CONTEXT) then info = QMCKL_INVALID_CONTEXT return endif - + if (m <= 0_8) then info = QMCKL_INVALID_ARG_4 return endif - + if (n <= 0_8) then info = QMCKL_INVALID_ARG_5 return endif - + if (k <= 0_8) then info = QMCKL_INVALID_ARG_6 return endif - + if (LDA <= 0) then info = QMCKL_INVALID_ARG_9 return @@ -142,7 +142,7 @@ integer function qmckl_dgemm_f(context, TransA, TransB, & end function qmckl_dgemm_f #+end_src - + *** C interface :noexport: #+CALL: generate_c_interface(table=qmckl_dgemm_args,rettyp="qmckl_exit_code",fname="qmckl_dgemm") @@ -177,7 +177,7 @@ end function qmckl_dgemm_f end function qmckl_dgemm #+end_src - + #+CALL: generate_f_interface(table=qmckl_dgemm_args,rettyp="qmckl_exit_code",fname="qmckl_dgemm") @@ -209,7 +209,7 @@ end function qmckl_dgemm_f end function qmckl_dgemm end interface #+end_src - + *** Test :noexport: #+begin_src f90 :tangle (eval f_test) @@ -217,13 +217,13 @@ integer(qmckl_exit_code) function test_qmckl_dgemm(context) bind(C) use qmckl implicit none integer(qmckl_context), intent(in), value :: context - + double precision, allocatable :: A(:,:), B(:,:), C(:,:), D(:,:) integer*8 :: m, n, k, LDA, LDB, LDC integer*8 :: i,j,l character :: TransA, TransB double precision :: x, alpha, beta - + TransA = 'N' TransB = 'N' m = 1_8 @@ -234,7 +234,7 @@ integer(qmckl_exit_code) function test_qmckl_dgemm(context) bind(C) LDC = m allocate( A(LDA,k), B(LDB,n) , C(LDC,n), D(LDC,n)) - + A = 0.d0 B = 0.d0 C = 0.d0 @@ -246,19 +246,19 @@ integer(qmckl_exit_code) function test_qmckl_dgemm(context) bind(C) A(i,j) = -10.d0 + dble(i+j) end do end do - + do j=1,n do i=1,k B(i,j) = -10.d0 + dble(i+j) end do end do - + test_qmckl_dgemm = qmckl_dgemm(context, TransA, TransB, m, n, k, alpha, A, LDA, B, LDB, beta, C, LDC) - + if (test_qmckl_dgemm /= QMCKL_SUCCESS) return - + test_qmckl_dgemm = QMCKL_FAILURE - + x = 0.d0 do j=1,n do i=1,m @@ -268,23 +268,23 @@ integer(qmckl_exit_code) function test_qmckl_dgemm(context) bind(C) x = x + (D(i,j) - C(i,j))**2 end do end do - + if (dabs(x) <= 1.d-12) then test_qmckl_dgemm = QMCKL_SUCCESS endif - + deallocate(A,B,C,D) end function test_qmckl_dgemm #+end_src - #+begin_src c :comments link :tangle (eval c_test) + #+begin_src c :comments link :tangle (eval c_test) qmckl_exit_code test_qmckl_dgemm(qmckl_context context); assert(QMCKL_SUCCESS == test_qmckl_dgemm(context)); #+end_src ** ~qmckl_adjugate~ - + Given a matrix $\mathbf{A}$, the adjugate matrix $\text{adj}(\mathbf{A})$ is the transpose of the cofactors matrix of $\mathbf{A}$. @@ -305,7 +305,7 @@ assert(QMCKL_SUCCESS == test_qmckl_dgemm(context)); | ~B~ | ~double[][ldb]~ | out | Adjugate of $A$ | | ~ldb~ | ~int64_t~ | in | Leading dimension of array ~B~ | | ~det_l~ | ~double~ | inout | determinant of $A$ | - + Requirements: - ~context~ is not ~QMCKL_NULL_CONTEXT~ @@ -314,9 +314,9 @@ assert(QMCKL_SUCCESS == test_qmckl_dgemm(context)); - ~A~ is allocated with at least $m \times m \times 8$ bytes - ~ldb >= m~ - ~B~ is allocated with at least $m \times m \times 8$ bytes - + #+CALL: generate_c_header(table=qmckl_adjugate_args,rettyp="qmckl_exit_code",fname="qmckl_adjugate") - + #+RESULTS: #+begin_src c :tangle (eval h_func) :comments org qmckl_exit_code qmckl_adjugate ( @@ -326,13 +326,13 @@ assert(QMCKL_SUCCESS == test_qmckl_dgemm(context)); const int64_t lda, double* const B, const int64_t ldb, - double* det_l ); + double* det_l ); #+end_src - + For small matrices (\le 5\times 5), we use specialized routines for performance motivations. For larger sizes, we rely on the LAPACK library. - + #+begin_src f90 :tangle (eval f) integer function qmckl_adjugate_f(context, na, A, LDA, B, ldb, det_l) & result(info) @@ -347,14 +347,14 @@ integer function qmckl_adjugate_f(context, na, A, LDA, B, ldb, det_l) & double precision, intent(inout) :: det_l integer :: i,j - + info = QMCKL_SUCCESS - + if (context == QMCKL_NULL_CONTEXT) then info = QMCKL_INVALID_CONTEXT return endif - + if (na <= 0_8) then info = QMCKL_INVALID_ARG_2 return @@ -385,7 +385,7 @@ integer function qmckl_adjugate_f(context, na, A, LDA, B, ldb, det_l) & case default call adjugate_general(context, na, A, LDA, B, LDB, det_l) end select - + end function qmckl_adjugate_f #+end_src @@ -399,9 +399,9 @@ subroutine adjugate2(A,LDA,B,LDB,na,det_l) double precision, intent(inout) :: det_l double precision :: C(2,2) - + call cofactor2(A,LDA,C,2_8,na,det_l) - + B(1,1) = C(1,1) B(2,1) = C(1,2) B(1,2) = C(2,1) @@ -418,9 +418,9 @@ subroutine adjugate3(a,LDA,B,LDB,na,det_l) double precision, intent(inout) :: det_l double precision :: C(4,3) - + call cofactor3(A,LDA,C,4_8,na,det_l) - + B(1,1) = C(1,1) B(1,2) = C(2,1) B(1,3) = C(3,1) @@ -442,9 +442,9 @@ subroutine adjugate4(a,LDA,B,LDB,na,det_l) double precision, intent(inout) :: det_l double precision :: C(4,4) - + call cofactor4(A,LDA,B,4_8,na,det_l) - + B(1,1) = C(1,1) B(1,2) = C(2,1) B(1,3) = C(3,1) @@ -473,9 +473,9 @@ subroutine adjugate5(A,LDA,B,LDB,na,det_l) double precision, intent(inout) :: det_l double precision :: C(8,5) - + call cofactor5(A,LDA,C,8_8,na,det_l) - + B(1,1) = C(1,1) B(1,2) = C(2,1) B(1,3) = C(3,1) @@ -511,7 +511,7 @@ subroutine cofactor2(a,LDA,b,LDB,na,det_l) integer*8, intent(in) :: LDA, LDB integer*8 :: na double precision :: det_l - + det_l = a(1,1)*a(2,2) - a(1,2)*a(2,1) b(1,1) = a(2,2) b(2,1) = -a(2,1) @@ -535,15 +535,15 @@ subroutine cofactor3(a,LDA,b,LDB,na,det_l) b(1,1) = a(2,2)*a(3,3) - a(2,3)*a(3,2) b(2,1) = a(2,3)*a(3,1) - a(2,1)*a(3,3) b(3,1) = a(2,1)*a(3,2) - a(2,2)*a(3,1) - + b(1,2) = a(1,3)*a(3,2) - a(1,2)*a(3,3) b(2,2) = a(1,1)*a(3,3) - a(1,3)*a(3,1) b(3,2) = a(1,2)*a(3,1) - a(1,1)*a(3,2) - + b(1,3) = a(1,2)*a(2,3) - a(1,3)*a(2,2) b(2,3) = a(1,3)*a(2,1) - a(1,1)*a(2,3) b(3,3) = a(1,1)*a(2,2) - a(1,2)*a(2,1) - + end subroutine cofactor3 subroutine cofactor4(a,LDA,b,LDB,na,det_l) @@ -571,22 +571,22 @@ subroutine cofactor4(a,LDA,b,LDB,na,det_l) b(2,1) = -a(2,1)*(a(3,3)*a(4,4)-a(3,4)*a(4,3))+a(2,3)*(a(3,1)*a(4,4)-a(3,4)*a(4,1))-a(2,4)*(a(3,1)*a(4,3)-a(3,3)*a(4,1)) b(3,1) = a(2,1)*(a(3,2)*a(4,4)-a(3,4)*a(4,2))-a(2,2)*(a(3,1)*a(4,4)-a(3,4)*a(4,1))+a(2,4)*(a(3,1)*a(4,2)-a(3,2)*a(4,1)) b(4,1) = -a(2,1)*(a(3,2)*a(4,3)-a(3,3)*a(4,2))+a(2,2)*(a(3,1)*a(4,3)-a(3,3)*a(4,1))-a(2,3)*(a(3,1)*a(4,2)-a(3,2)*a(4,1)) - + b(1,2) = -a(1,2)*(a(3,3)*a(4,4)-a(3,4)*a(4,3))+a(1,3)*(a(3,2)*a(4,4)-a(3,4)*a(4,2))-a(1,4)*(a(3,2)*a(4,3)-a(3,3)*a(4,2)) b(2,2) = a(1,1)*(a(3,3)*a(4,4)-a(3,4)*a(4,3))-a(1,3)*(a(3,1)*a(4,4)-a(3,4)*a(4,1))+a(1,4)*(a(3,1)*a(4,3)-a(3,3)*a(4,1)) b(3,2) = -a(1,1)*(a(3,2)*a(4,4)-a(3,4)*a(4,2))+a(1,2)*(a(3,1)*a(4,4)-a(3,4)*a(4,1))-a(1,4)*(a(3,1)*a(4,2)-a(3,2)*a(4,1)) b(4,2) = a(1,1)*(a(3,2)*a(4,3)-a(3,3)*a(4,2))-a(1,2)*(a(3,1)*a(4,3)-a(3,3)*a(4,1))+a(1,3)*(a(3,1)*a(4,2)-a(3,2)*a(4,1)) - + b(1,3) = a(1,2)*(a(2,3)*a(4,4)-a(2,4)*a(4,3))-a(1,3)*(a(2,2)*a(4,4)-a(2,4)*a(4,2))+a(1,4)*(a(2,2)*a(4,3)-a(2,3)*a(4,2)) b(2,3) = -a(1,1)*(a(2,3)*a(4,4)-a(2,4)*a(4,3))+a(1,3)*(a(2,1)*a(4,4)-a(2,4)*a(4,1))-a(1,4)*(a(2,1)*a(4,3)-a(2,3)*a(4,1)) b(3,3) = a(1,1)*(a(2,2)*a(4,4)-a(2,4)*a(4,2))-a(1,2)*(a(2,1)*a(4,4)-a(2,4)*a(4,1))+a(1,4)*(a(2,1)*a(4,2)-a(2,2)*a(4,1)) b(4,3) = -a(1,1)*(a(2,2)*a(4,3)-a(2,3)*a(4,2))+a(1,2)*(a(2,1)*a(4,3)-a(2,3)*a(4,1))-a(1,3)*(a(2,1)*a(4,2)-a(2,2)*a(4,1)) - + b(1,4) = -a(1,2)*(a(2,3)*a(3,4)-a(2,4)*a(3,3))+a(1,3)*(a(2,2)*a(3,4)-a(2,4)*a(3,2))-a(1,4)*(a(2,2)*a(3,3)-a(2,3)*a(3,2)) b(2,4) = a(1,1)*(a(2,3)*a(3,4)-a(2,4)*a(3,3))-a(1,3)*(a(2,1)*a(3,4)-a(2,4)*a(3,1))+a(1,4)*(a(2,1)*a(3,3)-a(2,3)*a(3,1)) b(3,4) = -a(1,1)*(a(2,2)*a(3,4)-a(2,4)*a(3,2))+a(1,2)*(a(2,1)*a(3,4)-a(2,4)*a(3,1))-a(1,4)*(a(2,1)*a(3,2)-a(2,2)*a(3,1)) b(4,4) = a(1,1)*(a(2,2)*a(3,3)-a(2,3)*a(3,2))-a(1,2)*(a(2,1)*a(3,3)-a(2,3)*a(3,1))+a(1,3)*(a(2,1)*a(3,2)-a(2,2)*a(3,1)) - + end subroutine cofactor4 subroutine cofactor5(A,LDA,B,LDB,na,det_l) @@ -818,7 +818,7 @@ end end function qmckl_adjugate end interface #+end_src - + #+begin_src f90 :tangle (eval f) subroutine adjugate_general(context, na, A, LDA, B, LDB, det_l) @@ -836,7 +836,7 @@ subroutine adjugate_general(context, na, A, LDA, B, LDB, det_l) integer :: inf integer :: ipiv(LDA) integer :: lwork - integer :: i, j + integer(8) :: i, j #+end_src @@ -859,7 +859,7 @@ subroutine adjugate_general(context, na, A, LDA, B, LDB, det_l) #+begin_src f90 :tangle (eval f) det_l = 1.d0 - j=0 + j=0_8 do i=1,na j = j+min(abs(ipiv(i)-i),1) det_l = det_l*B(i,i) @@ -869,15 +869,15 @@ subroutine adjugate_general(context, na, A, LDA, B, LDB, det_l) As ~dgetrf~ returns $PLU=A$ where $P$ is a permutation matrix, the sign of the determinant is computed as $-1^m$ where $m$ is the number of permutations. - + #+begin_src f90 :tangle (eval f) - if (iand(j,1) /= 0) then + if (iand(j,1_8) /= 0_8) then det_l = -det_l endif #+end_src Then, the inverse of $A$ is computed using ~dgetri~: - + #+begin_src f90 :tangle (eval f) lwork = SIZE(work) call dgetri(na, B, LDB, ipiv, work, lwork, inf ) @@ -899,15 +899,15 @@ integer(qmckl_exit_code) function test_qmckl_adjugate(context) bind(C) use qmckl implicit none integer(qmckl_context), intent(in), value :: context - + double precision, allocatable :: A(:,:), B(:,:) integer*8 :: m, n, k, LDA, LDB integer*8 :: i,j,l double precision :: x, det_l, det_l_ref - LDA = 6 - LDB = 6 - + LDA = 6_8 + LDB = 6_8 + allocate( A(LDA,6), B(LDB,6)) A = 0.1d0 @@ -919,9 +919,9 @@ integer(qmckl_exit_code) function test_qmckl_adjugate(context) bind(C) A(6,6) = 6.0d0; test_qmckl_adjugate = QMCKL_SUCCESS - + #+end_src - + #+begin_src python :results output :output drawer import numpy as np import numpy.linalg as la @@ -958,7 +958,7 @@ for i in range(N): print (f" return") print (f" end if") print (f"") - + print ("#+end_src") # print(adj(A[0:i+1,0:i+1])) @@ -1158,17 +1158,17 @@ print ("#+end_src") #+end_example #+begin_src f90 :tangle (eval f_test) - + deallocate(A,B) end function test_qmckl_adjugate #+end_src - + #+begin_src c :comments link :tangle (eval c_test) qmckl_exit_code test_qmckl_adjugate(qmckl_context context); assert(QMCKL_SUCCESS == test_qmckl_adjugate(context)); #+end_src - + * End of files :noexport: #+begin_src c :comments link :tangle (eval c_test) diff --git a/org/qmckl_determinant.org b/org/qmckl_determinant.org index 13177a8..12cebf0 100644 --- a/org/qmckl_determinant.org +++ b/org/qmckl_determinant.org @@ -557,6 +557,7 @@ qmckl_exit_code qmckl_finalize_determinant(qmckl_context context) { "qmckl_finalize_determinant", NULL); } + return rc; } #+end_src