TREX/qmckl
1
0
Fork 0
mirror of https://github.com/TREX-CoE/qmckl.git synced 2024-12-23 04:44:03 +01:00
Code Issues Projects Releases Wiki Activity

Merge branch 'master' into work

Browse Source
This commit is contained in:
Anthony Scemama 2022-01-05 08:58:30 +01:00
commit 74f6e85176
6 changed files with 793 additions and 511 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

9
.github/workflows/vfc_test_workflow.yml vendored
View File

@ -25,7 +25,7 @@ jobs:
run: | run: |
ln -s /usr/bin/python3 /usr/bin/python ln -s /usr/bin/python3 /usr/bin/python
apt update apt update
apt -y install emacs pkg-config wget libhdf5-dev apt -y install emacs pkg-config wget libhdf5-dev libblas-dev liblapack-dev
- name: Install trexio - name: Install trexio
run: | run: |
@ -33,9 +33,10 @@ jobs:
wget https://github.com/TREX-CoE/trexio/releases/download/v${VERSION}/trexio-${VERSION}.tar.gz wget https://github.com/TREX-CoE/trexio/releases/download/v${VERSION}/trexio-${VERSION}.tar.gz
tar -zxf trexio-${VERSION}.tar.gz tar -zxf trexio-${VERSION}.tar.gz
cd trexio-${VERSION} cd trexio-${VERSION}
./configure --prefix=/usr ./configure --prefix=/usr CC="gcc-7" FC="gfortran-7"
make -j 8 # modify LDFLAGS to include -lhdf5_hl because autoconf sometime fails to detect the HL component
sudo make install make LDFLAGS="-L/usr/lib/x86_64-linux-gnu/hdf5/serial -lhdf5_hl"
make install
- name: Run tests - name: Run tests
run: vfc_ci test -g -r run: vfc_ci test -g -r

2
autogen.sh
View File

@ -1,5 +1,5 @@
#!/bin/bash #!/bin/bash
export srcdir="." export srcdir="."
${PYTHON} ${srcdir}/tools/build_makefile.py python ${srcdir}/tools/build_makefile.py
autoreconf -i -Wall --no-recursive autoreconf -i -Wall --no-recursive

6
configure.ac
View File

@ -124,10 +124,10 @@ PKG_CFLAGS="$PKG_CFLAGS $TREXIO_CFLAGS"
PKG_LIBS="$PKG_LIBS $TREXIO_LIBS" PKG_LIBS="$PKG_LIBS $TREXIO_LIBS"
## BLAS ## BLAS
#AX_BLAS([], [AC_MSG_ERROR([BLAS was not found.])]) AX_BLAS([], [AC_MSG_ERROR([BLAS was not found.])])
## LAPACK ## LAPACK
#AX_LAPACK([], [AC_MSG_ERROR([LAPACK was not found.])]) AX_LAPACK([], [AC_MSG_ERROR([LAPACK was not found.])])
# Options. # Options.
@ -252,7 +252,7 @@ fi
#mkl-dynamic-lp64-seq #mkl-dynamic-lp64-seq
PKG_LIBS="$PKG_LIBS $LIBS" PKG_LIBS="$PKG_LIBS $LIBS"
LIBS="$LAPACK_LIBS $BLAS_LIBS $PKG_LIBS" LIBS="$BLAS_LIBS $LAPACK_LIBS $BLAS_LIBS $PKG_LIBS"
CFLAGS="$CFLAGS $PKG_CFLAGS" CFLAGS="$CFLAGS $PKG_CFLAGS"
AC_SUBST([PKG_LIBS]) AC_SUBST([PKG_LIBS])
AC_SUBST([PKG_CFLAGS]) AC_SUBST([PKG_CFLAGS])

509
org/qmckl_blas.org
View File

@ -35,12 +35,12 @@ int main() {
| int64_t | n | in | Number of columns of the input matrix | | int64_t | n | in | Number of columns of the input matrix |
| int64_t | k | in | Number of columns of the input matrix | | int64_t | k | in | Number of columns of the input matrix |
| double | alpha | in | Number of columns of the input matrix | | double | alpha | in | Number of columns of the input matrix |
| double | A[][lda] | in | Array containing the $m \times n$ matrix $A$ | | double | A[][lda] | in | Array containing the matrix $A$ |
| int64_t | lda | in | Leading dimension of array ~A~ | | int64_t | lda | in | Leading dimension of array ~A~ |
| double | B[][ldb] | in | Array containing the $n \times m$ matrix $B$ | | double | B[][ldb] | in | Array containing the matrix $B$ |
| int64_t | ldb | in | Leading dimension of array ~B~ | | int64_t | ldb | in | Leading dimension of array ~B~ |
| double | beta | in | Array containing the $n \times m$ matrix $B$ | | double | beta | in | Array containing the matrix $B$ |
| double | C[][ldc] | out | Array containing the $n \times m$ matrix $B$ | | double | C[][ldc] | out | Array containing the matrix $B$ |
| int64_t | ldc | in | Leading dimension of array ~B~ | | int64_t | ldc | in | Leading dimension of array ~B~ |
*** Requirements *** Requirements
@ -85,7 +85,7 @@ integer function qmckl_dgemm_f(context, TransA, TransB, m, n, k, alpha, A, LDA,
result(info) result(info)
use qmckl use qmckl
implicit none implicit none
integer(qmckl_context) , intent(in) :: context integer(qmckl_context), intent(in) :: context
logical*8 , intent(in) :: TransA, TransB logical*8 , intent(in) :: TransA, TransB
integer*8 , intent(in) :: m, n, k integer*8 , intent(in) :: m, n, k
real*8 , intent(in) :: alpha, beta real*8 , intent(in) :: alpha, beta
@ -95,9 +95,9 @@ integer function qmckl_dgemm_f(context, TransA, TransB, m, n, k, alpha, A, LDA,
real*8 , intent(in) :: B(ldb,*) real*8 , intent(in) :: B(ldb,*)
integer*8 , intent(in) :: ldc integer*8 , intent(in) :: ldc
real*8 , intent(out) :: C(ldc,*) real*8 , intent(out) :: C(ldc,*)
real*8, allocatable :: AT(:,:), BT(:,:), CT(:,:) real*8, allocatable :: AT(:,:), BT(:,:), CT(:,:)
integer*4 :: qmckl_dgemm_N_N_f integer*4 :: qmckl_dgemm_N_N_f
integer*8 :: i,j,l, LDA_2, LDB_2 integer*8 :: i,j,l, LDA_2, LDB_2
info = QMCKL_SUCCESS info = QMCKL_SUCCESS
@ -170,6 +170,7 @@ integer function qmckl_dgemm_f(context, TransA, TransB, m, n, k, alpha, A, LDA,
else else
info = qmckl_dgemm_N_N_f(context, m, n, k, alpha, A, LDA_2, B, LDB_2, beta, c, LDC) info = qmckl_dgemm_N_N_f(context, m, n, k, alpha, A, LDA_2, B, LDB_2, beta, c, LDC)
endif endif
end function qmckl_dgemm_f end function qmckl_dgemm_f
integer function qmckl_dgemm_N_N_f(context, m, n, k, alpha, A, LDA, B, LDB, beta, C, LDC) & integer function qmckl_dgemm_N_N_f(context, m, n, k, alpha, A, LDA, B, LDB, beta, C, LDC) &
@ -196,32 +197,32 @@ integer function qmckl_dgemm_N_N_f(context, m, n, k, alpha, A, LDA, B, LDB, beta
endif endif
if (m <= 0_8) then if (m <= 0_8) then
info = QMCKL_INVALID_ARG_4 info = QMCKL_INVALID_ARG_2
return return
endif endif
if (n <= 0_8) then if (n <= 0_8) then
info = QMCKL_INVALID_ARG_5 info = QMCKL_INVALID_ARG_3
return return
endif endif
if (k <= 0_8) then if (k <= 0_8) then
info = QMCKL_INVALID_ARG_6 info = QMCKL_INVALID_ARG_4
return return
endif endif
if (LDA /= m) then if (LDA /= m) then
info = QMCKL_INVALID_ARG_9 info = QMCKL_INVALID_ARG_7
return return
endif endif
if (LDB /= k) then if (LDB /= k) then
info = QMCKL_INVALID_ARG_10 info = QMCKL_INVALID_ARG_8
return return
endif endif
if (LDC /= m) then if (LDC /= m) then
info = QMCKL_INVALID_ARG_13 info = QMCKL_INVALID_ARG_11
return return
endif endif
@ -230,6 +231,7 @@ integer function qmckl_dgemm_N_N_f(context, m, n, k, alpha, A, LDA, B, LDB, beta
else else
C = beta*C + alpha*matmul(A,B) C = beta*C + alpha*matmul(A,B)
endif endif
end function qmckl_dgemm_N_N_f end function qmckl_dgemm_N_N_f
#+end_src #+end_src
@ -358,11 +360,12 @@ integer(qmckl_exit_code) function test_qmckl_dgemm(context) bind(C)
end do end do
end do end do
if (dabs(x) <= 1.d-15) 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
@ -371,45 +374,40 @@ 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_adjoint~ ** ~qmckl_adjugate~
Matrix adjoint. Given a matrix M, returns a matrix M⁻¹ such that:
Given a matrix $\mathbf{A}$, the adjugate matrix
$\text{adj}(\mathbf{A})$ is the transpose of the cofactors matrix
of $\mathbf{A}$.
\[ \[
M · M^{-1} = I \text{adj}(\mathbf{A}) = \text{det}(\mathbf{A}) \, \mathbf{A}^{-1}
\] \]
This is a native Fortran implementation hand written (by: A. Scemama) See also: https://en.wikipedia.org/wiki/Adjugate_matrix
only for small matrices (<=5x5).
TODO: Add description about the external library dependence. #+NAME: qmckl_adjugate_args
#+NAME: qmckl_adjoint_args
| qmckl_context | context | in | Global state | | qmckl_context | context | in | Global state |
| int64_t | m | in | Number of rows of the input matrix | | int64_t | n | in | Number of rows and columns of the input matrix |
| int64_t | n | in | Number of columns of the input matrix |
| int64_t | lda | in | Leading dimension of array ~A~ | | int64_t | lda | in | Leading dimension of array ~A~ |
| double | A[][lda] | inout | Array containing the $m \times n$ matrix $A$ | | double | A[][lda] | inout | Array containing the $n \times n$ matrix $A$ |
| double | det_l | inout | determinant of A | | double | det_l | inout | determinant of $A$ |
*** Requirements *** Requirements
- ~context~ is not ~QMCKL_NULL_CONTEXT~ - ~context~ is not ~QMCKL_NULL_CONTEXT~
- ~m > 0~
- ~n > 0~ - ~n > 0~
- ~lda >= m~ - ~lda >= m~
- ~A~ is allocated with at least $m \times n \times 8$ bytes - ~A~ is allocated with at least $m \times m \times 8$ bytes
*** C header *** C header
#+CALL: generate_c_header(table=qmckl_adjoint_args,rettyp="qmckl_exit_code",fname="qmckl_adjoint") #+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_adjoint ( qmckl_exit_code qmckl_adjugate (
const qmckl_context context, const qmckl_context context,
const int64_t m,
const int64_t n, const int64_t n,
const int64_t lda, const int64_t lda,
double* A, double* A,
@ -417,42 +415,49 @@ M · M^{-1} = I
#+end_src #+end_src
*** Source *** Source
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) #+begin_src f90 :tangle (eval f)
integer function qmckl_adjoint_f(context, ma, na, LDA, A, det_l) & integer function qmckl_adjugate_f(context, na, LDA, A, det_l) &
result(info) result(info)
use qmckl use qmckl
implicit none implicit none
integer(qmckl_context) , intent(in) :: context integer(qmckl_context) , intent(in) :: context
double precision, intent(inout) :: A (LDA,na) double precision, intent(inout) :: A (LDA,na)
integer*8, intent(in) :: LDA integer*8, intent(in) :: LDA
integer*8, intent(in) :: ma
integer*8, intent(in) :: na integer*8, intent(in) :: na
double precision, intent(inout) :: det_l double precision, intent(inout) :: det_l
integer :: i,j integer :: i,j
!TODO CHECK ARGUMENTS
info = QMCKL_SUCCESS info = QMCKL_SUCCESS
select case (na) select case (na)
case default case default
print *," TODO: Implement general adjoint" call adjugate_general(context, na, LDA, A, det_l)
stop 0
case (5) case (5)
call adjoint5(a,LDA,na,det_l) call adjugate5(a,LDA,na,det_l)
case (4) case (4)
call adjoint4(a,LDA,na,det_l) call adjugate4(a,LDA,na,det_l)
case (3) case (3)
call adjoint3(a,LDA,na,det_l) call adjugate3(a,LDA,na,det_l)
case (2) case (2)
call adjoint2(a,LDA,na,det_l) call adjugate2(a,LDA,na,det_l)
case (1) case (1)
call adjoint1(a,LDA,na,det_l) call adjugate1(a,LDA,na,det_l)
case (0) case (0)
det_l=1.d0 det_l=1.d0
end select end select
end function qmckl_adjoint_f
subroutine adjoint1(a,LDA,na,det_l) end function qmckl_adjugate_f
#+end_src
#+begin_src f90 :tangle (eval f)
subroutine adjugate1(a,LDA,na,det_l)
implicit none implicit none
double precision, intent(inout) :: a (LDA,na) double precision, intent(inout) :: a (LDA,na)
integer*8, intent(in) :: LDA integer*8, intent(in) :: LDA
@ -460,9 +465,9 @@ subroutine adjoint1(a,LDA,na,det_l)
double precision, intent(inout) :: det_l double precision, intent(inout) :: det_l
call cofactor1(a,LDA,na,det_l) call cofactor1(a,LDA,na,det_l)
end end subroutine adjugate1
subroutine adjoint2(a,LDA,na,det_l) subroutine adjugate2(a,LDA,na,det_l)
implicit none implicit none
double precision, intent(inout) :: a (LDA,na) double precision, intent(inout) :: a (LDA,na)
integer*8, intent(in) :: LDA integer*8, intent(in) :: LDA
@ -481,9 +486,9 @@ subroutine adjoint2(a,LDA,na,det_l)
a(1,2) = b(1,2) a(1,2) = b(1,2)
a(2,1) = b(2,1) a(2,1) = b(2,1)
a(2,2) = b(2,2) a(2,2) = b(2,2)
end end subroutine adjugate2
subroutine adjoint3(a,LDA,na,det_l) subroutine adjugate3(a,LDA,na,det_l)
implicit none implicit none
double precision, intent(inout) :: a (LDA,na) double precision, intent(inout) :: a (LDA,na)
integer*8, intent(in) :: LDA integer*8, intent(in) :: LDA
@ -513,9 +518,9 @@ subroutine adjoint3(a,LDA,na,det_l)
a(1,3) = b(1,3) a(1,3) = b(1,3)
a(2,3) = b(2,3) a(2,3) = b(2,3)
a(3,3) = b(3,3) a(3,3) = b(3,3)
end end subroutine adjugate3
subroutine adjoint4(a,LDA,na,det_l) subroutine adjugate4(a,LDA,na,det_l)
implicit none implicit none
double precision, intent(inout) :: a (LDA,na) double precision, intent(inout) :: a (LDA,na)
integer*8, intent(in) :: LDA integer*8, intent(in) :: LDA
@ -559,9 +564,9 @@ subroutine adjoint4(a,LDA,na,det_l)
a(2,4) = b(2,4) a(2,4) = b(2,4)
a(3,4) = b(3,4) a(3,4) = b(3,4)
a(4,4) = b(4,4) a(4,4) = b(4,4)
end end subroutine adjugate4
subroutine adjoint5(a,LDA,na,det_l) subroutine adjugate5(a,LDA,na,det_l)
implicit none implicit none
double precision, intent(inout) :: a (LDA,na) double precision, intent(inout) :: a (LDA,na)
integer*8, intent(in) :: LDA integer*8, intent(in) :: LDA
@ -623,7 +628,7 @@ subroutine adjoint5(a,LDA,na,det_l)
a(3,5) = b(3,5) a(3,5) = b(3,5)
a(4,5) = b(4,5) a(4,5) = b(4,5)
a(5,5) = b(5,5) a(5,5) = b(5,5)
end end subroutine adjugate5
subroutine cofactor1(a,LDA,na,det_l) subroutine cofactor1(a,LDA,na,det_l)
implicit none implicit none
@ -634,7 +639,7 @@ subroutine cofactor1(a,LDA,na,det_l)
det_l = a(1,1) det_l = a(1,1)
a(1,1) = 1.d0 a(1,1) = 1.d0
end end subroutine cofactor1
subroutine cofactor2(a,LDA,na,det_l) subroutine cofactor2(a,LDA,na,det_l)
implicit none implicit none
@ -653,7 +658,7 @@ subroutine cofactor2(a,LDA,na,det_l)
a(2,1) = -b(2,1) a(2,1) = -b(2,1)
a(1,2) = -b(1,2) a(1,2) = -b(1,2)
a(2,2) = b(1,1) a(2,2) = b(1,1)
end end subroutine cofactor2
subroutine cofactor3(a,LDA,na,det_l) subroutine cofactor3(a,LDA,na,det_l)
implicit none implicit none
@ -666,11 +671,13 @@ subroutine cofactor3(a,LDA,na,det_l)
det_l = a(1,1)*(a(2,2)*a(3,3)-a(2,3)*a(3,2)) & det_l = 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,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)) +a(1,3)*(a(2,1)*a(3,2)-a(2,2)*a(3,1))
do i=1,4 do i=1,4
b(i,1) = a(i,1) b(i,1) = a(i,1)
b(i,2) = a(i,2) b(i,2) = a(i,2)
b(i,3) = a(i,3) b(i,3) = a(i,3)
end do end do
a(1,1) = b(2,2)*b(3,3) - b(2,3)*b(3,2) a(1,1) = b(2,2)*b(3,3) - b(2,3)*b(3,2)
a(2,1) = b(2,3)*b(3,1) - b(2,1)*b(3,3) a(2,1) = b(2,3)*b(3,1) - b(2,1)*b(3,3)
a(3,1) = b(2,1)*b(3,2) - b(2,2)*b(3,1) a(3,1) = b(2,1)*b(3,2) - b(2,2)*b(3,1)
@ -683,7 +690,7 @@ subroutine cofactor3(a,LDA,na,det_l)
a(2,3) = b(1,3)*b(2,1) - b(1,1)*b(2,3) a(2,3) = b(1,3)*b(2,1) - b(1,1)*b(2,3)
a(3,3) = b(1,1)*b(2,2) - b(1,2)*b(2,1) a(3,3) = b(1,1)*b(2,2) - b(1,2)*b(2,1)
end end subroutine cofactor3
subroutine cofactor4(a,LDA,na,det_l) subroutine cofactor4(a,LDA,na,det_l)
implicit none implicit none
@ -732,7 +739,7 @@ subroutine cofactor4(a,LDA,na,det_l)
a(3,4) = -b(1,1)*(b(2,2)*b(3,4)-b(2,4)*b(3,2))+b(1,2)*(b(2,1)*b(3,4)-b(2,4)*b(3,1))-b(1,4)*(b(2,1)*b(3,2)-b(2,2)*b(3,1)) a(3,4) = -b(1,1)*(b(2,2)*b(3,4)-b(2,4)*b(3,2))+b(1,2)*(b(2,1)*b(3,4)-b(2,4)*b(3,1))-b(1,4)*(b(2,1)*b(3,2)-b(2,2)*b(3,1))
a(4,4) = b(1,1)*(b(2,2)*b(3,3)-b(2,3)*b(3,2))-b(1,2)*(b(2,1)*b(3,3)-b(2,3)*b(3,1))+b(1,3)*(b(2,1)*b(3,2)-b(2,2)*b(3,1)) a(4,4) = b(1,1)*(b(2,2)*b(3,3)-b(2,3)*b(3,2))-b(1,2)*(b(2,1)*b(3,3)-b(2,3)*b(3,1))+b(1,3)*(b(2,1)*b(3,2)-b(2,2)*b(3,1))
end end subroutine cofactor4
subroutine cofactor5(a,LDA,na,det_l) subroutine cofactor5(a,LDA,na,det_l)
implicit none implicit none
@ -742,6 +749,7 @@ subroutine cofactor5(a,LDA,na,det_l)
double precision, intent(inout) :: det_l double precision, intent(inout) :: det_l
double precision :: b(5,5) double precision :: b(5,5)
integer :: i,j integer :: i,j
det_l = a(1,1)*(a(2,2)*(a(3,3)*(a(4,4)*a(5,5)-a(4,5)*a(5,4))-a(3,4)*( & det_l = a(1,1)*(a(2,2)*(a(3,3)*(a(4,4)*a(5,5)-a(4,5)*a(5,4))-a(3,4)*( &
a(4,3)*a(5,5)-a(4,5)*a(5,3))+a(3,5)*(a(4,3)*a(5,4)-a(4,4)*a(5,3)))- & a(4,3)*a(5,5)-a(4,5)*a(5,3))+a(3,5)*(a(4,3)*a(5,4)-a(4,4)*a(5,3)))- &
a(2,3)*(a(3,2)*(a(4,4)*a(5,5)-a(4,5)*a(5,4))-a(3,4)*(a(4,2)*a(5,5)- & a(2,3)*(a(3,2)*(a(4,4)*a(5,5)-a(4,5)*a(5,4))-a(3,4)*(a(4,2)*a(5,5)- &
@ -921,126 +929,399 @@ subroutine cofactor5(a,LDA,na,det_l)
end end
#+end_src #+end_src
#+begin_src f90 :tangle (eval f)
subroutine adjugate_general(context, na, LDA, A, det_l)
use qmckl
implicit none
integer(qmckl_context) , intent(in) :: context
double precision, intent(inout) :: A (LDA,na)
integer*8, intent(in) :: LDA
integer*8, intent(in) :: na
double precision, intent(inout) :: det_l
double precision :: work(LDA*max(na,64))
integer :: inf
integer :: ipiv(LDA)
integer :: lwork
integer :: i, j
#+end_src
For larger matrices, we first compute the LU factorization $LU=A$
using the ~dgetrf~ routine.
#+begin_src f90 :tangle (eval f)
call dgetrf(na, na, a, LDA, ipiv, inf )
#+end_src
By convention, the determinant of $L$ is equal to one, so the
determinant of $A$ is equal to the determinant of $U$, which is
simply computed as the product of its diagonal elements.
#+begin_src f90 :tangle (eval f)
det_l = 1.d0
j=0
do i=1,na
j = j+min(abs(ipiv(i)-i),1)
det_l = det_l*a(i,i)
enddo
#+end_src
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
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, a, LDA, ipiv, work, lwork, inf )
#+end_src
and the adjugate matrix is computed as the product of the
determinant with the inverse:
#+begin_src f90 :tangle (eval f)
a(:,:) = a(:,:)*det_l
end subroutine adjugate_general
#+end_src
*** C interface :noexport: *** C interface :noexport:
#+CALL: generate_c_interface(table=qmckl_adjoint_args,rettyp="qmckl_exit_code",fname="qmckl_adjoint") #+CALL: generate_c_interface(table=qmckl_adjugate_args,rettyp="qmckl_exit_code",fname="qmckl_adjugate")
#+RESULTS: #+RESULTS:
#+begin_src f90 :tangle (eval f) :comments org :exports none #+begin_src f90 :tangle (eval f) :comments org :exports none
integer(c_int32_t) function qmckl_adjoint & integer(c_int32_t) function qmckl_adjugate &
(context, m, n, lda, A, det_l) & (context, n, lda, A, det_l) &
bind(C) result(info) bind(C) result(info)
use, intrinsic :: iso_c_binding use, intrinsic :: iso_c_binding
implicit none implicit none
integer (c_int64_t) , intent(in) , value :: context integer (c_int64_t) , intent(in) , value :: context
integer (c_int64_t) , intent(in) , value :: m
integer (c_int64_t) , intent(in) , value :: n integer (c_int64_t) , intent(in) , value :: n
integer (c_int64_t) , intent(in) , value :: lda integer (c_int64_t) , intent(in) , value :: lda
real (c_double ) , intent(inout) :: A(lda,*) real (c_double ) , intent(inout) :: A(lda,*)
real (c_double ) , intent(inout) :: det_l real (c_double ) , intent(inout) :: det_l
integer(c_int32_t), external :: qmckl_adjoint_f integer(c_int32_t), external :: qmckl_adjugate_f
info = qmckl_adjoint_f & info = qmckl_adjugate_f &
(context, m, n, lda, A, det_l) (context, n, lda, A, det_l)
end function qmckl_adjoint end function qmckl_adjugate
#+end_src #+end_src
#+CALL: generate_f_interface(table=qmckl_adjoint_args,rettyp="qmckl_exit_code",fname="qmckl_adjoint") #+CALL: generate_f_interface(table=qmckl_adjugate_args,rettyp="qmckl_exit_code",fname="qmckl_adjugate")
#+RESULTS: #+RESULTS:
#+begin_src f90 :tangle (eval fh_func) :comments org :exports none #+begin_src f90 :tangle (eval fh_func) :comments org :exports none
interface interface
integer(c_int32_t) function qmckl_adjoint & integer(c_int32_t) function qmckl_adjugate &
(context, m, n, lda, A, det_l) & (context, n, lda, A, det_l) &
bind(C) bind(C)
use, intrinsic :: iso_c_binding use, intrinsic :: iso_c_binding
import import
implicit none implicit none
integer (c_int64_t) , intent(in) , value :: context integer (c_int64_t) , intent(in) , value :: context
integer (c_int64_t) , intent(in) , value :: m
integer (c_int64_t) , intent(in) , value :: n integer (c_int64_t) , intent(in) , value :: n
integer (c_int64_t) , intent(in) , value :: lda integer (c_int64_t) , intent(in) , value :: lda
real (c_double ) , intent(inout) :: A(lda,*) real (c_double ) , intent(inout) :: A(lda,*)
real (c_double ) , intent(inout) :: det_l real (c_double ) , intent(inout) :: det_l
end function qmckl_adjoint end function qmckl_adjugate
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)
integer(qmckl_exit_code) function test_qmckl_adjoint(context) bind(C) 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(:,:), C(:,:) double precision, allocatable :: A(:,:), B(:,:)
integer*8 :: m, n, k, LDA, LDB, LDC 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
m = 4_8 LDA = 6
k = 4_8 LDB = 6
LDA = m
LDB = m
LDC = m
allocate( A(LDA,k), C(LDC,k)) allocate( A(LDA,6), B(LDB,6))
A = 0.10d0 A = 0.1d0
C = 0.d0
A(1,1) = 1.0d0; A(1,1) = 1.0d0;
A(2,2) = 2.0d0; A(2,2) = 2.0d0;
A(3,3) = 3.0d0; A(3,3) = 3.0d0;
A(4,4) = 4.0d0; A(4,4) = 4.0d0;
A(5,5) = 5.0d0;
A(6,6) = 6.0d0;
! Exact inverse (Mathematica) test_qmckl_adjugate = QMCKL_SUCCESS
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_adjoint = qmckl_adjoint(context, m, k, LDA, A, det_l) #+end_src
if (test_qmckl_adjoint /= QMCKL_SUCCESS) return #+begin_src python :results output :output drawer
import numpy as np
import numpy.linalg as la
N=6
test_qmckl_adjoint = QMCKL_FAILURE A = np.zeros( (N,N) )
A += 0.1
for i in range(N):
A[i][i] = i+1
def adj(A):
return la.det(A) * la.inv(A)
print ("#+begin_src f90 :tangle (eval f_test)")
for i in range(N):
print ("! N = ", i+1)
print (f" B(:,:) = A(:,:)")
print (f" test_qmckl_adjugate = qmckl_adjugate(context, {i+1}_8, LDB, B, det_l)")
print (f" if (test_qmckl_adjugate /= QMCKL_SUCCESS) return")
print (f" if (dabs((det_l - ({la.det(A[0:i+1,0:i+1])}d0))/det_l) > 1.d-13) then")
print (f" print *, 'N = {i+1}: det = ', det_l, {la.det(A[0:i+1,0:i+1])}d0")
print (f" test_qmckl_adjugate = {i+1}")
print (f" return")
print (f" end if")
print (f"")
print (f" x = 0.d0")
for j in range(i+1):
C = adj(A[0:i+1,0:i+1])
for k in range(i+1):
print (f" x = x + (B({j+1},{k+1}) - ({C[j,k]}) )**2")
print (f" if (dabs(x / det_l) > 1.d-12) then")
print (f" print *, 'N = {i+1}: x = ', x")
print (f" test_qmckl_adjugate = {i+1}")
print (f" return")
print (f" end if")
print (f"")
print ("#+end_src")
# print(adj(A[0:i+1,0:i+1]))
#+end_src
#+RESULTS:
#+begin_src f90 :tangle (eval f_test)
! N = 1
B(:,:) = A(:,:)
test_qmckl_adjugate = qmckl_adjugate(context, 1_8, LDB, B, det_l)
if (test_qmckl_adjugate /= QMCKL_SUCCESS) return
if (dabs((det_l - (1.0d0))/det_l) > 1.d-13) then
print *, 'N = 1: det = ', det_l, 1.0d0
test_qmckl_adjugate = 1
return
end if
x = 0.d0 x = 0.d0
do j=1,m x = x + (B(1,1) - (1.0) )**2
do i=1,k if (dabs(x / det_l) > 1.d-12) then
x = x + (A(i,j) - (C(i,j) * det_l_ref))**2 print *, 'N = 1: x = ', x
end do test_qmckl_adjugate = 1
end do return
end if
if (dabs(x) <= 1.d-15 .and. (dabs(det_l_ref - det_l)) <= 1.d-15) then ! N = 2
test_qmckl_adjoint = QMCKL_SUCCESS B(:,:) = A(:,:)
endif test_qmckl_adjugate = qmckl_adjugate(context, 2_8, LDB, B, det_l)
if (test_qmckl_adjugate /= QMCKL_SUCCESS) return
if (dabs((det_l - (1.99d0))/det_l) > 1.d-13) then
print *, 'N = 2: det = ', det_l, 1.99d0
test_qmckl_adjugate = 2
return
end if
deallocate(A,C) x = 0.d0
end function test_qmckl_adjoint x = x + (B(1,1) - (1.9999999999999998) )**2
x = x + (B(1,2) - (-0.09999999999999999) )**2
x = x + (B(2,1) - (-0.09999999999999999) )**2
x = x + (B(2,2) - (0.9999999999999999) )**2
if (dabs(x / det_l) > 1.d-12) then
print *, 'N = 2: x = ', x
test_qmckl_adjugate = 2
return
end if
! N = 3
B(:,:) = A(:,:)
test_qmckl_adjugate = qmckl_adjugate(context, 3_8, LDB, B, det_l)
if (test_qmckl_adjugate /= QMCKL_SUCCESS) return
if (dabs((det_l - (5.942000000000001d0))/det_l) > 1.d-13) then
print *, 'N = 3: det = ', det_l, 5.942000000000001d0
test_qmckl_adjugate = 3
return
end if
x = 0.d0
x = x + (B(1,1) - (5.990000000000001) )**2
x = x + (B(1,2) - (-0.29000000000000004) )**2
x = x + (B(1,3) - (-0.19000000000000003) )**2
x = x + (B(2,1) - (-0.29000000000000004) )**2
x = x + (B(2,2) - (2.9900000000000007) )**2
x = x + (B(2,3) - (-0.09000000000000001) )**2
x = x + (B(3,1) - (-0.19000000000000003) )**2
x = x + (B(3,2) - (-0.09) )**2
x = x + (B(3,3) - (1.9900000000000002) )**2
if (dabs(x / det_l) > 1.d-12) then
print *, 'N = 3: x = ', x
test_qmckl_adjugate = 3
return
end if
! N = 4
B(:,:) = A(:,:)
test_qmckl_adjugate = qmckl_adjugate(context, 4_8, LDB, B, det_l)
if (test_qmckl_adjugate /= QMCKL_SUCCESS) return
if (dabs((det_l - (23.669700000000006d0))/det_l) > 1.d-13) then
print *, 'N = 4: det = ', det_l, 23.669700000000006d0
test_qmckl_adjugate = 4
return
end if
x = 0.d0
x = x + (B(1,1) - (23.91200000000001) )**2
x = x + (B(1,2) - (-1.1310000000000004) )**2
x = x + (B(1,3) - (-0.7410000000000001) )**2
x = x + (B(1,4) - (-0.5510000000000002) )**2
x = x + (B(2,1) - (-1.1310000000000002) )**2
x = x + (B(2,2) - (11.922000000000002) )**2
x = x + (B(2,3) - (-0.351) )**2
x = x + (B(2,4) - (-0.261) )**2
x = x + (B(3,1) - (-0.7410000000000002) )**2
x = x + (B(3,2) - (-0.351) )**2
x = x + (B(3,3) - (7.932000000000001) )**2
x = x + (B(3,4) - (-0.17100000000000004) )**2
x = x + (B(4,1) - (-0.5510000000000002) )**2
x = x + (B(4,2) - (-0.261) )**2
x = x + (B(4,3) - (-0.17100000000000004) )**2
x = x + (B(4,4) - (5.942000000000001) )**2
if (dabs(x / det_l) > 1.d-12) then
print *, 'N = 4: x = ', x
test_qmckl_adjugate = 4
return
end if
! N = 5
B(:,:) = A(:,:)
test_qmckl_adjugate = qmckl_adjugate(context, 5_8, LDB, B, det_l)
if (test_qmckl_adjugate /= QMCKL_SUCCESS) return
if (dabs((det_l - (117.91554000000008d0))/det_l) > 1.d-13) then
print *, 'N = 5: det = ', det_l, 117.91554000000008d0
test_qmckl_adjugate = 5
return
end if
x = 0.d0
x = x + (B(1,1) - (119.31770000000006) )**2
x = x + (B(1,2) - (-5.541900000000004) )**2
x = x + (B(1,3) - (-3.6309000000000022) )**2
x = x + (B(1,4) - (-2.6999000000000017) )**2
x = x + (B(1,5) - (-2.1489000000000016) )**2
x = x + (B(2,1) - (-5.541900000000004) )**2
x = x + (B(2,2) - (59.435700000000026) )**2
x = x + (B(2,3) - (-1.7199000000000007) )**2
x = x + (B(2,4) - (-1.2789000000000006) )**2
x = x + (B(2,5) - (-1.0179000000000005) )**2
x = x + (B(3,1) - (-3.6309000000000027) )**2
x = x + (B(3,2) - (-1.7199000000000007) )**2
x = x + (B(3,3) - (39.53370000000002) )**2
x = x + (B(3,4) - (-0.8379000000000005) )**2
x = x + (B(3,5) - (-0.6669000000000004) )**2
x = x + (B(4,1) - (-2.699900000000002) )**2
x = x + (B(4,2) - (-1.2789000000000006) )**2
x = x + (B(4,3) - (-0.8379000000000004) )**2
x = x + (B(4,4) - (29.611700000000017) )**2
x = x + (B(4,5) - (-0.4959000000000003) )**2
x = x + (B(5,1) - (-2.1489000000000016) )**2
x = x + (B(5,2) - (-1.0179000000000005) )**2
x = x + (B(5,3) - (-0.6669000000000004) )**2
x = x + (B(5,4) - (-0.4959000000000003) )**2
x = x + (B(5,5) - (23.669700000000013) )**2
if (dabs(x / det_l) > 1.d-12) then
print *, 'N = 5: x = ', x
test_qmckl_adjugate = 5
return
end if
! N = 6
B(:,:) = A(:,:)
test_qmckl_adjugate = qmckl_adjugate(context, 6_8, LDB, B, det_l)
if (test_qmckl_adjugate /= QMCKL_SUCCESS) return
if (dabs((det_l - (705.1783350000001d0))/det_l) > 1.d-13) then
print *, 'N = 6: det = ', det_l, 705.1783350000001d0
test_qmckl_adjugate = 6
return
end if
x = 0.d0
x = x + (B(1,1) - (714.5040400000001) )**2
x = x + (B(1,2) - (-32.697210000000005) )**2
x = x + (B(1,3) - (-21.422310000000003) )**2
x = x + (B(1,4) - (-15.929410000000006) )**2
x = x + (B(1,5) - (-12.678510000000003) )**2
x = x + (B(1,6) - (-10.529610000000003) )**2
x = x + (B(2,1) - (-32.69721) )**2
x = x + (B(2,2) - (355.65834) )**2
x = x + (B(2,3) - (-10.147409999999997) )**2
x = x + (B(2,4) - (-7.54551) )**2
x = x + (B(2,5) - (-6.005610000000001) )**2
x = x + (B(2,6) - (-4.987709999999999) )**2
x = x + (B(3,1) - (-21.422310000000003) )**2
x = x + (B(3,2) - (-10.147409999999999) )**2
x = x + (B(3,3) - (236.51663999999997) )**2
x = x + (B(3,4) - (-4.943610000000001) )**2
x = x + (B(3,5) - (-3.93471) )**2
x = x + (B(3,6) - (-3.267810000000001) )**2
x = x + (B(4,1) - (-15.929410000000003) )**2
x = x + (B(4,2) - (-7.54551) )**2
x = x + (B(4,3) - (-4.9436100000000005) )**2
x = x + (B(4,4) - (177.13894000000002) )**2
x = x + (B(4,5) - (-2.92581) )**2
x = x + (B(4,6) - (-2.42991) )**2
x = x + (B(5,1) - (-12.678510000000001) )**2
x = x + (B(5,2) - (-6.005609999999999) )**2
x = x + (B(5,3) - (-3.9347100000000004) )**2
x = x + (B(5,4) - (-2.92581) )**2
x = x + (B(5,5) - (141.58524) )**2
x = x + (B(5,6) - (-1.93401) )**2
x = x + (B(6,1) - (-10.529610000000003) )**2
x = x + (B(6,2) - (-4.98771) )**2
x = x + (B(6,3) - (-3.2678100000000003) )**2
x = x + (B(6,4) - (-2.42991) )**2
x = x + (B(6,5) - (-1.9340100000000005) )**2
x = x + (B(6,6) - (117.91554000000001) )**2
if (dabs(x / det_l) > 1.d-12) then
print *, 'N = 6: x = ', x
test_qmckl_adjugate = 6
return
end if
#+end_src
#+begin_src f90 :tangle (eval f_test)
deallocate(A,B)
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_adjoint(qmckl_context context); qmckl_exit_code test_qmckl_adjugate(qmckl_context context);
assert(QMCKL_SUCCESS == test_qmckl_adjoint(context)); assert(QMCKL_SUCCESS == test_qmckl_adjugate(context));
#+end_src #+end_src
* End of files :noexport: * End of files :noexport:

4
org/qmckl_determinant.org
View File

@ -1817,7 +1817,7 @@ integer function qmckl_compute_det_inv_matrix_alpha_f(context, &
do iwalk = 1, walk_num do iwalk = 1, walk_num
! Value ! Value
matA(1:alpha_num,1:alpha_num) = det_vgl_alpha(1:alpha_num, 1:alpha_num, 1, iwalk, idet) matA(1:alpha_num,1:alpha_num) = det_vgl_alpha(1:alpha_num, 1:alpha_num, 1, iwalk, idet)
res = qmckl_adjoint(context, alpha_num, alpha_num, LDA, matA, det_l) res = qmckl_adjugate(context, alpha_num, LDA, matA, det_l)
det_adj_matrix_alpha(1:alpha_num, 1:alpha_num, iwalk, idet) = matA det_adj_matrix_alpha(1:alpha_num, 1:alpha_num, iwalk, idet) = matA
det_inv_matrix_alpha(1:alpha_num, 1:alpha_num, iwalk, idet) = matA/det_l det_inv_matrix_alpha(1:alpha_num, 1:alpha_num, iwalk, idet) = matA/det_l
det_value_alpha(iwalk, idet) = det_l det_value_alpha(iwalk, idet) = det_l
@ -1948,7 +1948,7 @@ integer function qmckl_compute_det_inv_matrix_beta_f(context, &
do iwalk = 1, walk_num do iwalk = 1, walk_num
! Value ! Value
matA(1:beta_num,1:beta_num) = det_vgl_beta(1:beta_num, 1:beta_num, 1, iwalk, idet) matA(1:beta_num,1:beta_num) = det_vgl_beta(1:beta_num, 1:beta_num, 1, iwalk, idet)
res = qmckl_adjoint(context, beta_num, beta_num, LDA, matA, det_l) res = qmckl_adjugate(context, beta_num, LDA, matA, det_l)
det_adj_matrix_beta(1:beta_num, 1:beta_num, iwalk, idet) = matA det_adj_matrix_beta(1:beta_num, 1:beta_num, iwalk, idet) = matA
det_inv_matrix_beta(1:beta_num, 1:beta_num, iwalk, idet) = matA/det_l det_inv_matrix_beta(1:beta_num, 1:beta_num, iwalk, idet) = matA/det_l
det_value_beta(iwalk, idet) = det_l det_value_beta(iwalk, idet) = det_l

4
tools/build_makefile.py
View File

@ -250,8 +250,8 @@ def main():
tmp = "EXTRA_DIST += " tmp = "EXTRA_DIST += "
r = subprocess.check_output("git ls-tree --name-only -r HEAD".split()) r = subprocess.check_output("git ls-tree --name-only -r HEAD".split())
for line in r.splitlines(): for line in r.splitlines():
if "share/qmckl/test_data/" in line: if b"share/qmckl/test_data/" in line:
tmp += " \\\n " + line tmp += " \\\n " + line.decode('utf8')
tmp += "\n" tmp += "\n"
output += tmp.split("\n") output += tmp.split("\n")