TREX/qmckl
1
0
Fork 0
mirror of https://github.com/TREX-CoE/qmckl.git synced 2024-10-18 05:51:52 +02:00
Code Issues Projects Releases Wiki Activity

Fix warnings

Browse Source
This commit is contained in:
Anthony Scemama 2022-01-08 15:36:07 +01:00
parent 1539a40dfe
commit 28dc3978f4
3 changed files with 93 additions and 85 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

43
org/qmckl_ao.org
View File

@ -601,7 +601,7 @@ qmckl_set_ao_basis_nucleus_shell_num (qmckl_context context,
if (size_max < nucl_num) { if (size_max < nucl_num) {
return qmckl_failwith( context, return qmckl_failwith( context,
QMCKL_FAILURE, QMCKL_INVALID_ARG_3,
"qmckl_set_ao_basis_nucleus_shell_num", "qmckl_set_ao_basis_nucleus_shell_num",
"input array too small"); "input array too small");
} }
@ -663,7 +663,7 @@ qmckl_set_ao_basis_shell_ang_mom (qmckl_context context,
if (size_max < shell_num) { if (size_max < shell_num) {
return qmckl_failwith( context, return qmckl_failwith( context,
QMCKL_FAILURE, QMCKL_INVALID_ARG_3,
"qmckl_set_ao_basis_shell_ang_mom", "qmckl_set_ao_basis_shell_ang_mom",
"input array too small"); "input array too small");
} }
@ -726,7 +726,7 @@ qmckl_set_ao_basis_shell_prim_num (qmckl_context context,
if (size_max < shell_num) { if (size_max < shell_num) {
return qmckl_failwith( context, return qmckl_failwith( context,
QMCKL_FAILURE, QMCKL_INVALID_ARG_3,
"qmckl_set_ao_basis_shell_prim_num", "qmckl_set_ao_basis_shell_prim_num",
"input array too small"); "input array too small");
} }
@ -789,7 +789,7 @@ qmckl_set_ao_basis_shell_prim_index (qmckl_context context,
if (size_max < shell_num) { if (size_max < shell_num) {
return qmckl_failwith( context, return qmckl_failwith( context,
QMCKL_FAILURE, QMCKL_INVALID_ARG_3,
"qmckl_set_ao_basis_shell_prim_index", "qmckl_set_ao_basis_shell_prim_index",
"input array too small"); "input array too small");
} }
@ -851,7 +851,7 @@ qmckl_set_ao_basis_shell_factor (qmckl_context context,
if (size_max < shell_num) { if (size_max < shell_num) {
return qmckl_failwith( context, return qmckl_failwith( context,
QMCKL_FAILURE, QMCKL_INVALID_ARG_3,
"qmckl_set_ao_basis_shell_factor", "qmckl_set_ao_basis_shell_factor",
"input array too small"); "input array too small");
} }
@ -913,7 +913,7 @@ qmckl_set_ao_basis_exponent (qmckl_context context,
if (size_max < prim_num) { if (size_max < prim_num) {
return qmckl_failwith( context, return qmckl_failwith( context,
QMCKL_FAILURE, QMCKL_INVALID_ARG_3,
"qmckl_set_ao_basis_exponent", "qmckl_set_ao_basis_exponent",
"input array too small"); "input array too small");
} }
@ -975,7 +975,7 @@ qmckl_set_ao_basis_coefficient (qmckl_context context,
if (size_max < prim_num) { if (size_max < prim_num) {
return qmckl_failwith( context, return qmckl_failwith( context,
QMCKL_FAILURE, QMCKL_INVALID_ARG_3,
"qmckl_set_ao_basis_coefficient", "qmckl_set_ao_basis_coefficient",
"input array too small"); "input array too small");
} }
@ -1037,7 +1037,7 @@ qmckl_set_ao_basis_prim_factor (qmckl_context context,
if (size_max < prim_num) { if (size_max < prim_num) {
return qmckl_failwith( context, return qmckl_failwith( context,
QMCKL_FAILURE, QMCKL_INVALID_ARG_3,
"qmckl_set_ao_basis_prim_factor", "qmckl_set_ao_basis_prim_factor",
"input array too small"); "input array too small");
} }
@ -1099,7 +1099,7 @@ qmckl_set_ao_basis_ao_factor (qmckl_context context,
if (size_max < ao_num) { if (size_max < ao_num) {
return qmckl_failwith( context, return qmckl_failwith( context,
QMCKL_FAILURE, QMCKL_INVALID_ARG_3,
"qmckl_set_ao_basis_ao_factor", "qmckl_set_ao_basis_ao_factor",
"input array too small"); "input array too small");
} }
@ -2638,6 +2638,13 @@ qmckl_get_ao_basis_primitive_vgl (qmckl_context context,
NULL); 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; qmckl_exit_code rc;
rc = qmckl_provide_ao_basis_primitive_vgl(context); 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; qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
assert (ctx != NULL); 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) { if (size_max < sze) {
return qmckl_failwith( context, return qmckl_failwith( context,
QMCKL_FAILURE, QMCKL_INVALID_ARG_3,
"qmckl_get_ao_basis_primitive_vgl", "qmckl_get_ao_basis_primitive_vgl",
"input array too small"); "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; 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; qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
assert (ctx != NULL); 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) { if (size_max < sze) {
return qmckl_failwith( context, return qmckl_failwith( context,
QMCKL_FAILURE, QMCKL_INVALID_ARG_3,
"qmckl_get_ao_basis_shell_vgl", "qmckl_get_ao_basis_shell_vgl",
"input array too small"); "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; 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; qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
assert (ctx != NULL); 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) { if (size_max < sze) {
return qmckl_failwith( context, return qmckl_failwith( context,
QMCKL_FAILURE, QMCKL_INVALID_ARG_3,
"qmckl_get_ao_basis_ao_vgl", "qmckl_get_ao_basis_ao_vgl",
"input array too small"); "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; return QMCKL_SUCCESS;
} }

134
org/qmckl_blas.org
View File

@ -22,12 +22,12 @@ int main() {
* Matrix operations * Matrix operations
** ~qmckl_dgemm~ ** ~qmckl_dgemm~
Matrix multiplication: Matrix multiplication:
\[ \[
C_{ij} = \beta C_{ij} + \alpha \sum_{k} A_{ik} \cdot B_{kj} C_{ij} = \beta C_{ij} + \alpha \sum_{k} A_{ik} \cdot B_{kj}
\] \]
# TODO: Add description about the external library dependence. # TODO: Add description about the external library dependence.
@ -48,7 +48,7 @@ int main() {
| ~beta~ | ~double~ | in | Array containing the matrix $B$ | | ~beta~ | ~double~ | in | Array containing the matrix $B$ |
| ~C~ | ~double[][ldc]~ | out | Array containing the matrix $B$ | | ~C~ | ~double[][ldc]~ | out | Array containing the matrix $B$ |
| ~ldc~ | ~int64_t~ | in | Leading dimension of array ~B~ | | ~ldc~ | ~int64_t~ | in | Leading dimension of array ~B~ |
Requirements: Requirements:
- ~context~ is not ~QMCKL_NULL_CONTEXT~ - ~context~ is not ~QMCKL_NULL_CONTEXT~
@ -61,7 +61,7 @@ int main() {
- ~A~ is allocated with at least $m \times k \times 8$ bytes - ~A~ is allocated with at least $m \times k \times 8$ bytes
- ~B~ is allocated with at least $k \times n \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 - ~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") #+CALL: generate_c_header(table=qmckl_dgemm_args,rettyp="qmckl_exit_code",fname="qmckl_dgemm")
#+RESULTS: #+RESULTS:
@ -80,7 +80,7 @@ int main() {
const int64_t ldb, const int64_t ldb,
const double beta, const double beta,
double* const C, double* const C,
const int64_t ldc ); const int64_t ldc );
#+end_src #+end_src
#+begin_src f90 :tangle (eval f) #+begin_src f90 :tangle (eval f)
@ -100,28 +100,28 @@ integer function qmckl_dgemm_f(context, TransA, TransB, &
integer*8 , intent(in) :: ldc integer*8 , intent(in) :: ldc
double precision , intent(out) :: C(ldc,*) double precision , intent(out) :: C(ldc,*)
info = QMCKL_SUCCESS info = QMCKL_SUCCESS
if (context == QMCKL_NULL_CONTEXT) then if (context == QMCKL_NULL_CONTEXT) then
info = QMCKL_INVALID_CONTEXT info = QMCKL_INVALID_CONTEXT
return return
endif endif
if (m <= 0_8) then if (m <= 0_8) then
info = QMCKL_INVALID_ARG_4 info = QMCKL_INVALID_ARG_4
return return
endif endif
if (n <= 0_8) then if (n <= 0_8) then
info = QMCKL_INVALID_ARG_5 info = QMCKL_INVALID_ARG_5
return return
endif endif
if (k <= 0_8) then if (k <= 0_8) then
info = QMCKL_INVALID_ARG_6 info = QMCKL_INVALID_ARG_6
return return
endif endif
if (LDA <= 0) then if (LDA <= 0) then
info = QMCKL_INVALID_ARG_9 info = QMCKL_INVALID_ARG_9
return return
@ -142,7 +142,7 @@ integer function qmckl_dgemm_f(context, TransA, TransB, &
end function qmckl_dgemm_f end function qmckl_dgemm_f
#+end_src #+end_src
*** C interface :noexport: *** C interface :noexport:
#+CALL: generate_c_interface(table=qmckl_dgemm_args,rettyp="qmckl_exit_code",fname="qmckl_dgemm") #+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 function qmckl_dgemm
#+end_src #+end_src
#+CALL: generate_f_interface(table=qmckl_dgemm_args,rettyp="qmckl_exit_code",fname="qmckl_dgemm") #+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 function qmckl_dgemm
end interface end interface
#+end_src #+end_src
*** Test :noexport: *** Test :noexport:
#+begin_src f90 :tangle (eval f_test) #+begin_src f90 :tangle (eval f_test)
@ -217,13 +217,13 @@ integer(qmckl_exit_code) function test_qmckl_dgemm(context) bind(C)
use qmckl use qmckl
implicit none implicit none
integer(qmckl_context), intent(in), value :: context integer(qmckl_context), intent(in), value :: context
double precision, allocatable :: A(:,:), B(:,:), C(:,:), D(:,:) double precision, allocatable :: A(:,:), B(:,:), C(:,:), D(:,:)
integer*8 :: m, n, k, LDA, LDB, LDC integer*8 :: m, n, k, LDA, LDB, LDC
integer*8 :: i,j,l integer*8 :: i,j,l
character :: TransA, TransB character :: TransA, TransB
double precision :: x, alpha, beta double precision :: x, alpha, beta
TransA = 'N' TransA = 'N'
TransB = 'N' TransB = 'N'
m = 1_8 m = 1_8
@ -234,7 +234,7 @@ integer(qmckl_exit_code) function test_qmckl_dgemm(context) bind(C)
LDC = m LDC = m
allocate( A(LDA,k), B(LDB,n) , C(LDC,n), D(LDC,n)) allocate( A(LDA,k), B(LDB,n) , C(LDC,n), D(LDC,n))
A = 0.d0 A = 0.d0
B = 0.d0 B = 0.d0
C = 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) A(i,j) = -10.d0 + dble(i+j)
end do end do
end do end do
do j=1,n do j=1,n
do i=1,k do i=1,k
B(i,j) = -10.d0 + dble(i+j) B(i,j) = -10.d0 + dble(i+j)
end do end do
end do end do
test_qmckl_dgemm = qmckl_dgemm(context, TransA, TransB, m, n, k, alpha, A, LDA, B, LDB, beta, C, LDC) 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 if (test_qmckl_dgemm /= QMCKL_SUCCESS) return
test_qmckl_dgemm = QMCKL_FAILURE test_qmckl_dgemm = QMCKL_FAILURE
x = 0.d0 x = 0.d0
do j=1,n do j=1,n
do i=1,m 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 x = x + (D(i,j) - C(i,j))**2
end do end do
end do end do
if (dabs(x) <= 1.d-12) then if (dabs(x) <= 1.d-12) then
test_qmckl_dgemm = QMCKL_SUCCESS test_qmckl_dgemm = QMCKL_SUCCESS
endif endif
deallocate(A,B,C,D) deallocate(A,B,C,D)
end function test_qmckl_dgemm end function test_qmckl_dgemm
#+end_src #+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); qmckl_exit_code test_qmckl_dgemm(qmckl_context context);
assert(QMCKL_SUCCESS == test_qmckl_dgemm(context)); assert(QMCKL_SUCCESS == test_qmckl_dgemm(context));
#+end_src #+end_src
** ~qmckl_adjugate~ ** ~qmckl_adjugate~
Given a matrix $\mathbf{A}$, the adjugate matrix Given a matrix $\mathbf{A}$, the adjugate matrix
$\text{adj}(\mathbf{A})$ is the transpose of the cofactors matrix $\text{adj}(\mathbf{A})$ is the transpose of the cofactors matrix
of $\mathbf{A}$. of $\mathbf{A}$.
@ -305,7 +305,7 @@ assert(QMCKL_SUCCESS == test_qmckl_dgemm(context));
| ~B~ | ~double[][ldb]~ | out | Adjugate of $A$ | | ~B~ | ~double[][ldb]~ | out | Adjugate of $A$ |
| ~ldb~ | ~int64_t~ | in | Leading dimension of array ~B~ | | ~ldb~ | ~int64_t~ | in | Leading dimension of array ~B~ |
| ~det_l~ | ~double~ | inout | determinant of $A$ | | ~det_l~ | ~double~ | inout | determinant of $A$ |
Requirements: Requirements:
- ~context~ is not ~QMCKL_NULL_CONTEXT~ - ~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 - ~A~ is allocated with at least $m \times m \times 8$ bytes
- ~ldb >= m~ - ~ldb >= m~
- ~B~ is allocated with at least $m \times m \times 8$ bytes - ~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") #+CALL: generate_c_header(table=qmckl_adjugate_args,rettyp="qmckl_exit_code",fname="qmckl_adjugate")
#+RESULTS: #+RESULTS:
#+begin_src c :tangle (eval h_func) :comments org #+begin_src c :tangle (eval h_func) :comments org
qmckl_exit_code qmckl_adjugate ( qmckl_exit_code qmckl_adjugate (
@ -326,13 +326,13 @@ assert(QMCKL_SUCCESS == test_qmckl_dgemm(context));
const int64_t lda, const int64_t lda,
double* const B, double* const B,
const int64_t ldb, const int64_t ldb,
double* det_l ); double* det_l );
#+end_src #+end_src
For small matrices (\le 5\times 5), we use specialized routines For small matrices (\le 5\times 5), we use specialized routines
for performance motivations. For larger sizes, we rely on the for performance motivations. For larger sizes, we rely on the
LAPACK library. LAPACK library.
#+begin_src f90 :tangle (eval f) #+begin_src f90 :tangle (eval f)
integer function qmckl_adjugate_f(context, na, A, LDA, B, ldb, det_l) & integer function qmckl_adjugate_f(context, na, A, LDA, B, ldb, det_l) &
result(info) 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 double precision, intent(inout) :: det_l
integer :: i,j integer :: i,j
info = QMCKL_SUCCESS info = QMCKL_SUCCESS
if (context == QMCKL_NULL_CONTEXT) then if (context == QMCKL_NULL_CONTEXT) then
info = QMCKL_INVALID_CONTEXT info = QMCKL_INVALID_CONTEXT
return return
endif endif
if (na <= 0_8) then if (na <= 0_8) then
info = QMCKL_INVALID_ARG_2 info = QMCKL_INVALID_ARG_2
return return
@ -385,7 +385,7 @@ integer function qmckl_adjugate_f(context, na, A, LDA, B, ldb, det_l) &
case default case default
call adjugate_general(context, na, A, LDA, B, LDB, det_l) call adjugate_general(context, na, A, LDA, B, LDB, det_l)
end select end select
end function qmckl_adjugate_f end function qmckl_adjugate_f
#+end_src #+end_src
@ -399,9 +399,9 @@ subroutine adjugate2(A,LDA,B,LDB,na,det_l)
double precision, intent(inout) :: det_l double precision, intent(inout) :: det_l
double precision :: C(2,2) double precision :: C(2,2)
call cofactor2(A,LDA,C,2_8,na,det_l) call cofactor2(A,LDA,C,2_8,na,det_l)
B(1,1) = C(1,1) B(1,1) = C(1,1)
B(2,1) = C(1,2) B(2,1) = C(1,2)
B(1,2) = C(2,1) 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, intent(inout) :: det_l
double precision :: C(4,3) double precision :: C(4,3)
call cofactor3(A,LDA,C,4_8,na,det_l) call cofactor3(A,LDA,C,4_8,na,det_l)
B(1,1) = C(1,1) B(1,1) = C(1,1)
B(1,2) = C(2,1) B(1,2) = C(2,1)
B(1,3) = C(3,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, intent(inout) :: det_l
double precision :: C(4,4) double precision :: C(4,4)
call cofactor4(A,LDA,B,4_8,na,det_l) call cofactor4(A,LDA,B,4_8,na,det_l)
B(1,1) = C(1,1) B(1,1) = C(1,1)
B(1,2) = C(2,1) B(1,2) = C(2,1)
B(1,3) = C(3,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, intent(inout) :: det_l
double precision :: C(8,5) double precision :: C(8,5)
call cofactor5(A,LDA,C,8_8,na,det_l) call cofactor5(A,LDA,C,8_8,na,det_l)
B(1,1) = C(1,1) B(1,1) = C(1,1)
B(1,2) = C(2,1) B(1,2) = C(2,1)
B(1,3) = C(3,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, intent(in) :: LDA, LDB
integer*8 :: na integer*8 :: na
double precision :: det_l double precision :: det_l
det_l = a(1,1)*a(2,2) - a(1,2)*a(2,1) det_l = a(1,1)*a(2,2) - a(1,2)*a(2,1)
b(1,1) = a(2,2) b(1,1) = a(2,2)
b(2,1) = -a(2,1) 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(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(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(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(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(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(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(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(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) b(3,3) = a(1,1)*a(2,2) - a(1,2)*a(2,1)
end subroutine cofactor3 end subroutine cofactor3
subroutine cofactor4(a,LDA,b,LDB,na,det_l) 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(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(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(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(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(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(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(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(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(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(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(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(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(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(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)) 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 end subroutine cofactor4
subroutine cofactor5(A,LDA,B,LDB,na,det_l) subroutine cofactor5(A,LDA,B,LDB,na,det_l)
@ -818,7 +818,7 @@ end
end function qmckl_adjugate end function qmckl_adjugate
end interface end interface
#+end_src #+end_src
#+begin_src f90 :tangle (eval f) #+begin_src f90 :tangle (eval f)
subroutine adjugate_general(context, na, A, LDA, B, LDB, det_l) 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 :: inf
integer :: ipiv(LDA) integer :: ipiv(LDA)
integer :: lwork integer :: lwork
integer :: i, j integer(8) :: i, j
#+end_src #+end_src
@ -859,7 +859,7 @@ subroutine adjugate_general(context, na, A, LDA, B, LDB, det_l)
#+begin_src f90 :tangle (eval f) #+begin_src f90 :tangle (eval f)
det_l = 1.d0 det_l = 1.d0
j=0 j=0_8
do i=1,na do i=1,na
j = j+min(abs(ipiv(i)-i),1) j = j+min(abs(ipiv(i)-i),1)
det_l = det_l*B(i,i) 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 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 sign of the determinant is computed as $-1^m$ where $m$ is the
number of permutations. number of permutations.
#+begin_src f90 :tangle (eval f) #+begin_src f90 :tangle (eval f)
if (iand(j,1) /= 0) then if (iand(j,1_8) /= 0_8) then
det_l = -det_l det_l = -det_l
endif endif
#+end_src #+end_src
Then, the inverse of $A$ is computed using ~dgetri~: Then, the inverse of $A$ is computed using ~dgetri~:
#+begin_src f90 :tangle (eval f) #+begin_src f90 :tangle (eval f)
lwork = SIZE(work) lwork = SIZE(work)
call dgetri(na, B, LDB, ipiv, work, lwork, inf ) 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 use qmckl
implicit none implicit none
integer(qmckl_context), intent(in), value :: context integer(qmckl_context), intent(in), value :: context
double precision, allocatable :: A(:,:), B(:,:) double precision, allocatable :: A(:,:), B(:,:)
integer*8 :: m, n, k, LDA, LDB integer*8 :: m, n, k, LDA, LDB
integer*8 :: i,j,l integer*8 :: i,j,l
double precision :: x, det_l, det_l_ref double precision :: x, det_l, det_l_ref
LDA = 6 LDA = 6_8
LDB = 6 LDB = 6_8
allocate( A(LDA,6), B(LDB,6)) allocate( A(LDA,6), B(LDB,6))
A = 0.1d0 A = 0.1d0
@ -919,9 +919,9 @@ integer(qmckl_exit_code) function test_qmckl_adjugate(context) bind(C)
A(6,6) = 6.0d0; A(6,6) = 6.0d0;
test_qmckl_adjugate = QMCKL_SUCCESS test_qmckl_adjugate = QMCKL_SUCCESS
#+end_src #+end_src
#+begin_src python :results output :output drawer #+begin_src python :results output :output drawer
import numpy as np import numpy as np
import numpy.linalg as la import numpy.linalg as la
@ -958,7 +958,7 @@ for i in range(N):
print (f" return") print (f" return")
print (f" end if") print (f" end if")
print (f"") print (f"")
print ("#+end_src") print ("#+end_src")
# print(adj(A[0:i+1,0:i+1])) # print(adj(A[0:i+1,0:i+1]))
@ -1158,17 +1158,17 @@ print ("#+end_src")
#+end_example #+end_example
#+begin_src f90 :tangle (eval f_test) #+begin_src f90 :tangle (eval f_test)
deallocate(A,B) deallocate(A,B)
end function test_qmckl_adjugate end function test_qmckl_adjugate
#+end_src #+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_adjugate(qmckl_context context); qmckl_exit_code test_qmckl_adjugate(qmckl_context context);
assert(QMCKL_SUCCESS == test_qmckl_adjugate(context)); assert(QMCKL_SUCCESS == test_qmckl_adjugate(context));
#+end_src #+end_src
* End of files :noexport: * End of files :noexport:
#+begin_src c :comments link :tangle (eval c_test) #+begin_src c :comments link :tangle (eval c_test)

1
org/qmckl_determinant.org
View File

@ -557,6 +557,7 @@ qmckl_exit_code qmckl_finalize_determinant(qmckl_context context) {
"qmckl_finalize_determinant", "qmckl_finalize_determinant",
NULL); NULL);
} }
return rc;
} }
#+end_src #+end_src