TREX/qmckl
1
0
Fork 0
mirror of https://github.com/TREX-CoE/qmckl.git synced 2024-07-01 17:06:11 +02:00
Code Issues Projects Releases Wiki Activity
60356cdbbb
qmckl/org/qmckl_blas.org
Anthony Scemama 60356cdbbb Removed DIR$
2021-12-12 11:28:57 +01:00

1058 lines
39 KiB
Org Mode
Raw Blame History

#+TITLE: BLAS functions
#+SETUPFILE: ../tools/theme.setup
#+INCLUDE: ../tools/lib.org
* Headers :noexport:
#+begin_src elisp :noexport :results none
(org-babel-lob-ingest "../tools/lib.org")
#+end_src
#+begin_src c :comments link :tangle (eval c_test) :noweb yes
#include "qmckl.h"
#include "assert.h"
#ifdef HAVE_CONFIG_H
#include "config.h"
#endif
int main() {
qmckl_context context;
context = qmckl_context_create();
#+end_src
* Matrix operations
** ~qmckl_dgemm~
Matrix multiply: $C_{ij} = \beta C_{ij} + \alpha \sum_{k} A_{ik} \cdot B_{kj}$ using Fortran ~matmul~ function.
TODO: Add description about the external library dependence.
#+NAME: qmckl_dgemm_args
| qmckl_context | context | in | Global state |
| bool | TransA | in | Number of rows of the input matrix |
| bool | TransB | in | Number of rows of the input matrix |
| int64_t | m | in | Number of rows 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 |
| double | alpha | in | Number of columns of the input matrix |
| double | A[][lda] | in | Array containing the $m \times n$ matrix $A$ |
| int64_t | lda | in | Leading dimension of array ~A~ |
| double | B[][ldb] | in | Array containing the $n \times m$ matrix $B$ |
| int64_t | ldb | in | Leading dimension of array ~B~ |
| double | beta | in | Array containing the $n \times m$ matrix $B$ |
| double | C[][ldc] | out | Array containing the $n \times m$ matrix $B$ |
| int64_t | ldc | in | Leading dimension of array ~B~ |
*** Requirements
- ~context~ is not ~QMCKL_NULL_CONTEXT~
- ~m > 0~
- ~n > 0~
- ~k > 0~
- ~lda >= m~
- ~ldb >= n~
- ~ldc >= n~
- ~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
*** C header
#+CALL: generate_c_header(table=qmckl_dgemm_args,rettyp="qmckl_exit_code",fname="qmckl_dgemm")
#+RESULTS:
#+BEGIN_src c :tangle (eval h_func) :comments org
qmckl_exit_code qmckl_dgemm (
const qmckl_context context,
const bool TransA,
const bool TransB,
const int64_t m,
const int64_t n,
const int64_t k,
const double alpha,
const double* A,
const int64_t lda,
const double* B,
const int64_t ldb,
const double beta,
double* const C,
const int64_t ldc );
#+END_src
*** Source
#+begin_src f90 :tangle (eval f)
integer function qmckl_dgemm_f(context, TransA, TransB, m, n, k, alpha, A, LDA, B, LDB, beta, C, LDC) &
result(info)
use qmckl
implicit none
integer(qmckl_context) , intent(in) :: context
logical*8 , intent(in) :: TransA, TransB
integer*8 , intent(in) :: m, n, k
real*8 , intent(in) :: alpha, beta
integer*8 , intent(in) :: lda
real*8 , intent(in) :: A(lda,*)
integer*8 , intent(in) :: ldb
real*8 , intent(in) :: B(ldb,*)
integer*8 , intent(in) :: ldc
real*8 , intent(out) :: C(ldc,*)
real*8, allocatable :: AT(:,:), BT(:,:), CT(:,:)
integer*4 :: qmckl_dgemm_N_N_f
integer*8 :: i,j,l, LDA_2, LDB_2
info = QMCKL_SUCCESS
if (TransA) then
allocate(AT(m,k))
do i = 1, k
do j = 1, m
AT(j,i) = A(i,j)
end do
end do
LDA_2 = M
else
LDA_2 = LDA
endif
if (TransB) then
allocate(BT(k,n))
do i = 1, n
do j = 1, k
BT(j,i) = B(i,j)
end do
end do
LDB_2 = K
else
LDB_2 = LDB
endif
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_2 /= m) then
info = QMCKL_INVALID_ARG_9
return
endif
if (LDB_2 /= k) then
info = QMCKL_INVALID_ARG_10
return
endif
if (LDC /= m) then
info = QMCKL_INVALID_ARG_13
return
endif
if (TransA) then
info = qmckl_dgemm_N_N_f(context, m, n, k, alpha, AT, LDA_2, B, LDB_2, beta, c, LDC)
else if (TransB) then
info = qmckl_dgemm_N_N_f(context, m, n, k, alpha, A, LDA_2, BT, LDB_2, beta, c, LDC)
else if (TransA .and. TransB) then
info = qmckl_dgemm_N_N_f(context, m, n, k, alpha, AT, LDA_2, BT, LDB_2, beta, c, LDC)
else
info = qmckl_dgemm_N_N_f(context, m, n, k, alpha, A, LDA_2, B, LDB_2, beta, c, LDC)
endif
end function qmckl_dgemm_f
integer function qmckl_dgemm_N_N_f(context, m, n, k, alpha, A, LDA, B, LDB, beta, C, LDC) &
result(info)
use qmckl
implicit none
integer(qmckl_context) , intent(in) :: context
integer*8 , intent(in) :: m, n, k
real*8 , intent(in) :: alpha, beta
integer*8 , intent(in) :: lda
real*8 , intent(in) :: A(lda,k)
integer*8 , intent(in) :: ldb
real*8 , intent(in) :: B(ldb,n)
integer*8 , intent(in) :: ldc
real*8 , intent(out) :: C(ldc,n)
integer*8 :: i,j,l, LDA_2, LDB_2
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 /= m) then
info = QMCKL_INVALID_ARG_9
return
endif
if (LDB /= k) then
info = QMCKL_INVALID_ARG_10
return
endif
if (LDC /= m) then
info = QMCKL_INVALID_ARG_13
return
endif
if (alpha == 1.0d0 .and. beta == 0.0d0) then
C = matmul(A,B)
else
C = beta*C + alpha*matmul(A,B)
endif
end function qmckl_dgemm_N_N_f
#+end_src
*** C interface :noexport:
#+CALL: generate_c_interface(table=qmckl_dgemm_args,rettyp="qmckl_exit_code",fname="qmckl_dgemm")
#+RESULTS:
#+BEGIN_src f90 :tangle (eval f) :comments org :exports none
integer(c_int32_t) function qmckl_dgemm &
(context, TransA, TransB, m, n, k, alpha, A, lda, B, ldb, beta, C, ldc) &
bind(C) result(info)
use, intrinsic :: iso_c_binding
implicit none
integer (c_int64_t) , intent(in) , value :: context
logical*8 , intent(in) , value :: TransA
logical*8 , intent(in) , value :: TransB
integer (c_int64_t) , intent(in) , value :: m
integer (c_int64_t) , intent(in) , value :: n
integer (c_int64_t) , intent(in) , value :: k
real (c_double ) , intent(in) , value :: alpha
integer (c_int64_t) , intent(in) , value :: lda
real (c_double ) , intent(in) :: A(lda,*)
integer (c_int64_t) , intent(in) , value :: ldb
real (c_double ) , intent(in) :: B(ldb,*)
real (c_double ) , intent(in) , value :: beta
integer (c_int64_t) , intent(in) , value :: ldc
real (c_double ) , intent(out) :: C(ldc,*)
integer(c_int32_t), external :: qmckl_dgemm_f
info = qmckl_dgemm_f &
(context, TransA, TransB, m, n, k, alpha, A, lda, B, ldb, beta, C, ldc)
end function qmckl_dgemm
#+END_src
#+CALL: generate_f_interface(table=qmckl_dgemm_args,rettyp="qmckl_exit_code",fname="qmckl_dgemm")
#+RESULTS:
#+begin_src f90 :tangle (eval fh_func) :comments org :exports none
interface
integer(c_int32_t) function qmckl_dgemm &
(context, TransA, TransB, m, n, k, alpha, A, lda, B, ldb, beta, C, ldc) &
bind(C)
use, intrinsic :: iso_c_binding
import
implicit none
integer (c_int64_t) , intent(in) , value :: context
logical*8 , intent(in) , value :: TransA
logical*8 , intent(in) , value :: TransB
integer (c_int64_t) , intent(in) , value :: m
integer (c_int64_t) , intent(in) , value :: n
integer (c_int64_t) , intent(in) , value :: k
real (c_double ) , intent(in) , value :: alpha
integer (c_int64_t) , intent(in) , value :: lda
real (c_double ) , intent(in) :: A(lda,*)
integer (c_int64_t) , intent(in) , value :: ldb
real (c_double ) , intent(in) :: B(ldb,*)
real (c_double ) , intent(in) , value :: beta
integer (c_int64_t) , intent(in) , value :: ldc
real (c_double ) , intent(out) :: C(ldc,*)
end function qmckl_dgemm
end interface
#+end_src
*** Test :noexport:
#+begin_src f90 :tangle (eval f_test)
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
logical*8 :: TransA, TransB
double precision :: x, alpha, beta
TransA = .False.
TransB = .False.
m = 1_8
k = 4_8
n = 6_8
LDA = m
LDB = k
LDC = m
allocate( A(LDA,k), B(LDB,n) , C(LDC,n), D(LDC,n))
A = 0.d0
B = 0.d0
C = 0.d0
D = 0.d0
alpha = 1.0d0
beta = 0.0d0
do j=1,k
do i=1,m
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
do l=1,k
D(i,j) = D(i,j) + A(i,l)*B(l,j)
end do
x = x + (D(i,j) - C(i,j))**2
end do
end do
if (dabs(x) <= 1.d-15) 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)
qmckl_exit_code test_qmckl_dgemm(qmckl_context context);
assert(QMCKL_SUCCESS == test_qmckl_dgemm(context));
#+end_src
** ~qmckl_adjoint~
Matrix adjoint. Given a matrix M, returns a matrix M⁻¹ such that:
\[
M · M^{-1} = I
\]
This is a native Fortran implementation hand written (by: A. Scemama)
only for small matrices (<=5x5).
TODO: Add description about the external library dependence.
#+NAME: qmckl_adjoint_args
| qmckl_context | context | in | Global state |
| int64_t | m | in | Number of rows of the input matrix |
| int64_t | n | in | Number of columns of the input matrix |
| int64_t | lda | in | Leading dimension of array ~A~ |
| double | A[][lda] | inout | Array containing the $m \times n$ matrix $A$ |
| double | det_l | inout | determinant of A |
*** Requirements
- ~context~ is not ~QMCKL_NULL_CONTEXT~
- ~m > 0~
- ~n > 0~
- ~lda >= m~
- ~A~ is allocated with at least $m \times n \times 8$ bytes
*** C header
#+CALL: generate_c_header(table=qmckl_adjoint_args,rettyp="qmckl_exit_code",fname="qmckl_adjoint")
#+RESULTS:
#+begin_src c :tangle (eval h_func) :comments org
qmckl_exit_code qmckl_adjoint (
const qmckl_context context,
const int64_t m,
const int64_t n,
const int64_t lda,
double* A,
double det_l );
#+end_src
*** Source
#+begin_src f90 :tangle (eval f)
integer function qmckl_adjoint_f(context, ma, na, LDA, A, det_l) &
result(info)
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) :: ma
integer*8, intent(in) :: na
double precision, intent(inout) :: det_l
integer :: i,j
info = QMCKL_SUCCESS
select case (na)
case default
print *," TODO: Implement general adjoint"
stop 0
case (5)
call adjoint5(a,LDA,na,det_l)
case (4)
call adjoint4(a,LDA,na,det_l)
case (3)
call adjoint3(a,LDA,na,det_l)
case (2)
call adjoint2(a,LDA,na,det_l)
case (1)
call adjoint1(a,LDA,na,det_l)
case (0)
det_l=1.d0
end select
end function qmckl_adjoint_f
subroutine adjoint1(a,LDA,na,det_l)
implicit none
double precision, intent(inout) :: a (LDA,na)
integer*8, intent(in) :: LDA
integer*8, intent(in) :: na
double precision, intent(inout) :: det_l
call cofactor1(a,LDA,na,det_l)
end
subroutine adjoint2(a,LDA,na,det_l)
implicit none
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 :: b(2,2)
call cofactor2(a,LDA,na,det_l)
! Calculate the transpose
b(1,1) = a(1,1)
b(1,2) = a(2,1)
b(2,1) = a(1,2)
b(2,2) = a(2,2)
a(1,1) = b(1,1)
a(1,2) = b(1,2)
a(2,1) = b(2,1)
a(2,2) = b(2,2)
end
subroutine adjoint3(a,LDA,na,det_l)
implicit none
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 :: b(3,3)
call cofactor3(a,LDA,na,det_l)
! Calculate the transpose
b(1,1) = a(1,1)
b(1,2) = a(2,1)
b(1,3) = a(3,1)
b(2,1) = a(1,2)
b(2,2) = a(2,2)
b(2,3) = a(3,2)
b(3,1) = a(1,3)
b(3,2) = a(2,3)
b(3,3) = a(3,3)
! copy
a(1,1) = b(1,1)
a(2,1) = b(2,1)
a(3,1) = b(3,1)
a(1,2) = b(1,2)
a(2,2) = b(2,2)
a(3,2) = b(3,2)
a(1,3) = b(1,3)
a(2,3) = b(2,3)
a(3,3) = b(3,3)
end
subroutine adjoint4(a,LDA,na,det_l)
implicit none
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 :: b(4,4)
call cofactor4(a,LDA,na,det_l)
! Calculate the transpose
b(1,1) = a(1,1)
b(1,2) = a(2,1)
b(1,3) = a(3,1)
b(1,4) = a(4,1)
b(2,1) = a(1,2)
b(2,2) = a(2,2)
b(2,3) = a(3,2)
b(2,4) = a(4,2)
b(3,1) = a(1,3)
b(3,2) = a(2,3)
b(3,3) = a(3,3)
b(3,4) = a(4,3)
b(4,1) = a(1,4)
b(4,2) = a(2,4)
b(4,3) = a(3,4)
b(4,4) = a(4,4)
! copy
a(1,1) = b(1,1)
a(2,1) = b(2,1)
a(3,1) = b(3,1)
a(4,1) = b(4,1)
a(1,2) = b(1,2)
a(2,2) = b(2,2)
a(3,2) = b(3,2)
a(4,2) = b(4,2)
a(1,3) = b(1,3)
a(2,3) = b(2,3)
a(3,3) = b(3,3)
a(4,3) = b(4,3)
a(1,4) = b(1,4)
a(2,4) = b(2,4)
a(3,4) = b(3,4)
a(4,4) = b(4,4)
end
subroutine adjoint5(a,LDA,na,det_l)
implicit none
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 :: b(5,5)
call cofactor5(a,LDA,na,det_l)
! Calculate the transpose
b(1,1) = a(1,1)
b(1,2) = a(2,1)
b(1,3) = a(3,1)
b(1,4) = a(4,1)
b(1,5) = a(5,1)
b(2,1) = a(1,2)
b(2,2) = a(2,2)
b(2,3) = a(3,2)
b(2,4) = a(4,2)
b(2,5) = a(5,2)
b(3,1) = a(1,3)
b(3,2) = a(2,3)
b(3,3) = a(3,3)
b(3,4) = a(4,3)
b(3,5) = a(5,3)
b(4,1) = a(1,4)
b(4,2) = a(2,4)
b(4,3) = a(3,4)
b(4,4) = a(4,4)
b(4,5) = a(5,4)
b(5,1) = a(1,5)
b(5,2) = a(2,5)
b(5,3) = a(3,5)
b(5,4) = a(4,5)
b(5,5) = a(5,5)
! copy
a(1,1) = b(1,1)
a(2,1) = b(2,1)
a(3,1) = b(3,1)
a(4,1) = b(4,1)
a(5,1) = b(5,1)
a(1,2) = b(1,2)
a(2,2) = b(2,2)
a(3,2) = b(3,2)
a(4,2) = b(4,2)
a(5,2) = b(5,2)
a(1,3) = b(1,3)
a(2,3) = b(2,3)
a(3,3) = b(3,3)
a(4,3) = b(4,3)
a(5,3) = b(5,3)
a(1,4) = b(1,4)
a(2,4) = b(2,4)
a(3,4) = b(3,4)
a(4,4) = b(4,4)
a(5,4) = b(5,4)
a(1,5) = b(1,5)
a(2,5) = b(2,5)
a(3,5) = b(3,5)
a(4,5) = b(4,5)
a(5,5) = b(5,5)
end
subroutine cofactor1(a,LDA,na,det_l)
implicit none
double precision, intent(inout) :: a (LDA,na)
integer*8, intent(in) :: LDA
integer*8, intent(in) :: na
double precision, intent(inout) :: det_l
det_l = a(1,1)
a(1,1) = 1.d0
end
subroutine cofactor2(a,LDA,na,det_l)
implicit none
double precision :: a (LDA,na)
integer*8 :: LDA
integer*8 :: na
double precision :: det_l
double precision :: b(2,2)
b(1,1) = a(1,1)
b(2,1) = a(2,1)
b(1,2) = a(1,2)
b(2,2) = a(2,2)
det_l = a(1,1)*a(2,2) - a(1,2)*a(2,1)
a(1,1) = b(2,2)
a(2,1) = -b(2,1)
a(1,2) = -b(1,2)
a(2,2) = b(1,1)
end
subroutine cofactor3(a,LDA,na,det_l)
implicit none
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 :: b(4,3)
integer :: i
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,3)*(a(2,1)*a(3,2)-a(2,2)*a(3,1))
do i=1,4
b(i,1) = a(i,1)
b(i,2) = a(i,2)
b(i,3) = a(i,3)
end do
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(3,1) = b(2,1)*b(3,2) - b(2,2)*b(3,1)
a(1,2) = b(1,3)*b(3,2) - b(1,2)*b(3,3)
a(2,2) = b(1,1)*b(3,3) - b(1,3)*b(3,1)
a(3,2) = b(1,2)*b(3,1) - b(1,1)*b(3,2)
a(1,3) = b(1,2)*b(2,3) - b(1,3)*b(2,2)
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)
end
subroutine cofactor4(a,LDA,na,det_l)
implicit none
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 :: b(4,4)
integer :: i,j
det_l = a(1,1)*(a(2,2)*(a(3,3)*a(4,4)-a(3,4)*a(4,3)) &
-a(2,3)*(a(3,2)*a(4,4)-a(3,4)*a(4,2)) &
+a(2,4)*(a(3,2)*a(4,3)-a(3,3)*a(4,2))) &
-a(1,2)*(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))) &
+a(1,3)*(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))) &
-a(1,4)*(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)))
do i=1,4
b(1,i) = a(1,i)
b(2,i) = a(2,i)
b(3,i) = a(3,i)
b(4,i) = a(4,i)
end do
a(1,1) = b(2,2)*(b(3,3)*b(4,4)-b(3,4)*b(4,3))-b(2,3)*(b(3,2)*b(4,4)-b(3,4)*b(4,2))+b(2,4)*(b(3,2)*b(4,3)-b(3,3)*b(4,2))
a(2,1) = -b(2,1)*(b(3,3)*b(4,4)-b(3,4)*b(4,3))+b(2,3)*(b(3,1)*b(4,4)-b(3,4)*b(4,1))-b(2,4)*(b(3,1)*b(4,3)-b(3,3)*b(4,1))
a(3,1) = b(2,1)*(b(3,2)*b(4,4)-b(3,4)*b(4,2))-b(2,2)*(b(3,1)*b(4,4)-b(3,4)*b(4,1))+b(2,4)*(b(3,1)*b(4,2)-b(3,2)*b(4,1))
a(4,1) = -b(2,1)*(b(3,2)*b(4,3)-b(3,3)*b(4,2))+b(2,2)*(b(3,1)*b(4,3)-b(3,3)*b(4,1))-b(2,3)*(b(3,1)*b(4,2)-b(3,2)*b(4,1))
a(1,2) = -b(1,2)*(b(3,3)*b(4,4)-b(3,4)*b(4,3))+b(1,3)*(b(3,2)*b(4,4)-b(3,4)*b(4,2))-b(1,4)*(b(3,2)*b(4,3)-b(3,3)*b(4,2))
a(2,2) = b(1,1)*(b(3,3)*b(4,4)-b(3,4)*b(4,3))-b(1,3)*(b(3,1)*b(4,4)-b(3,4)*b(4,1))+b(1,4)*(b(3,1)*b(4,3)-b(3,3)*b(4,1))
a(3,2) = -b(1,1)*(b(3,2)*b(4,4)-b(3,4)*b(4,2))+b(1,2)*(b(3,1)*b(4,4)-b(3,4)*b(4,1))-b(1,4)*(b(3,1)*b(4,2)-b(3,2)*b(4,1))
a(4,2) = b(1,1)*(b(3,2)*b(4,3)-b(3,3)*b(4,2))-b(1,2)*(b(3,1)*b(4,3)-b(3,3)*b(4,1))+b(1,3)*(b(3,1)*b(4,2)-b(3,2)*b(4,1))
a(1,3) = b(1,2)*(b(2,3)*b(4,4)-b(2,4)*b(4,3))-b(1,3)*(b(2,2)*b(4,4)-b(2,4)*b(4,2))+b(1,4)*(b(2,2)*b(4,3)-b(2,3)*b(4,2))
a(2,3) = -b(1,1)*(b(2,3)*b(4,4)-b(2,4)*b(4,3))+b(1,3)*(b(2,1)*b(4,4)-b(2,4)*b(4,1))-b(1,4)*(b(2,1)*b(4,3)-b(2,3)*b(4,1))
a(3,3) = b(1,1)*(b(2,2)*b(4,4)-b(2,4)*b(4,2))-b(1,2)*(b(2,1)*b(4,4)-b(2,4)*b(4,1))+b(1,4)*(b(2,1)*b(4,2)-b(2,2)*b(4,1))
a(4,3) = -b(1,1)*(b(2,2)*b(4,3)-b(2,3)*b(4,2))+b(1,2)*(b(2,1)*b(4,3)-b(2,3)*b(4,1))-b(1,3)*(b(2,1)*b(4,2)-b(2,2)*b(4,1))
a(1,4) = -b(1,2)*(b(2,3)*b(3,4)-b(2,4)*b(3,3))+b(1,3)*(b(2,2)*b(3,4)-b(2,4)*b(3,2))-b(1,4)*(b(2,2)*b(3,3)-b(2,3)*b(3,2))
a(2,4) = b(1,1)*(b(2,3)*b(3,4)-b(2,4)*b(3,3))-b(1,3)*(b(2,1)*b(3,4)-b(2,4)*b(3,1))+b(1,4)*(b(2,1)*b(3,3)-b(2,3)*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))
end
subroutine cofactor5(a,LDA,na,det_l)
implicit none
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 :: b(5,5)
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)*( &
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(4,5)*a(5,2))+a(3,5)*(a(4,2)*a(5,4)-a(4,4)*a(5,2)))+a(2,4)*(a(3,2)*( &
a(4,3)*a(5,5)-a(4,5)*a(5,3))-a(3,3)*(a(4,2)*a(5,5)-a(4,5)*a(5,2))+ &
a(3,5)*(a(4,2)*a(5,3)-a(4,3)*a(5,2)))-a(2,5)*(a(3,2)*(a(4,3)*a(5,4)- &
a(4,4)*a(5,3))-a(3,3)*(a(4,2)*a(5,4)-a(4,4)*a(5,2))+a(3,4)*(a(4,2)* &
a(5,3)-a(4,3)*a(5,2))))-a(1,2)*(a(2,1)*(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(2,3)*(a(3,1)*(a(4,4)*a(5,5)-a(4,5)*a(5,4))-a(3,4)*( &
a(4,1)*a(5,5)-a(4,5)*a(5,1))+a(3,5)*(a(4,1)*a(5,4)-a(4,4)*a(5,1)))+ &
a(2,4)*(a(3,1)*(a(4,3)*a(5,5)-a(4,5)*a(5,3))-a(3,3)*(a(4,1)*a(5,5)- &
a(4,5)*a(5,1))+a(3,5)*(a(4,1)*a(5,3)-a(4,3)*a(5,1)))-a(2,5)*(a(3,1)*( &
a(4,3)*a(5,4)-a(4,4)*a(5,3))-a(3,3)*(a(4,1)*a(5,4)-a(4,4)*a(5,1))+ &
a(3,4)*(a(4,1)*a(5,3)-a(4,3)*a(5,1))))+a(1,3)*(a(2,1)*(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(4,5)*a(5,2))+a(3,5)*( &
a(4,2)*a(5,4)-a(4,4)*a(5,2)))-a(2,2)*(a(3,1)*(a(4,4)*a(5,5)-a(4,5)* &
a(5,4))-a(3,4)*(a(4,1)*a(5,5)-a(4,5)*a(5,1))+a(3,5)*(a(4,1)*a(5,4)- &
a(4,4)*a(5,1)))+a(2,4)*(a(3,1)*(a(4,2)*a(5,5)-a(4,5)*a(5,2))-a(3,2)*( &
a(4,1)*a(5,5)-a(4,5)*a(5,1))+a(3,5)*(a(4,1)*a(5,2)-a(4,2)*a(5,1)))- &
a(2,5)*(a(3,1)*(a(4,2)*a(5,4)-a(4,4)*a(5,2))-a(3,2)*(a(4,1)*a(5,4)- &
a(4,4)*a(5,1))+a(3,4)*(a(4,1)*a(5,2)-a(4,2)*a(5,1))))-a(1,4)*(a(2,1)*( &
a(3,2)*(a(4,3)*a(5,5)-a(4,5)*a(5,3))-a(3,3)*(a(4,2)*a(5,5)-a(4,5)* &
a(5,2))+a(3,5)*(a(4,2)*a(5,3)-a(4,3)*a(5,2)))-a(2,2)*(a(3,1)*(a(4,3)* &
a(5,5)-a(4,5)*a(5,3))-a(3,3)*(a(4,1)*a(5,5)-a(4,5)*a(5,1))+a(3,5)*( &
a(4,1)*a(5,3)-a(4,3)*a(5,1)))+a(2,3)*(a(3,1)*(a(4,2)*a(5,5)-a(4,5)* &
a(5,2))-a(3,2)*(a(4,1)*a(5,5)-a(4,5)*a(5,1))+a(3,5)*(a(4,1)*a(5,2)- &
a(4,2)*a(5,1)))-a(2,5)*(a(3,1)*(a(4,2)*a(5,3)-a(4,3)*a(5,2))-a(3,2)*( &
a(4,1)*a(5,3)-a(4,3)*a(5,1))+a(3,3)*(a(4,1)*a(5,2)-a(4,2)*a(5,1))))+ &
a(1,5)*(a(2,1)*(a(3,2)*(a(4,3)*a(5,4)-a(4,4)*a(5,3))-a(3,3)*(a(4,2)* &
a(5,4)-a(4,4)*a(5,2))+a(3,4)*(a(4,2)*a(5,3)-a(4,3)*a(5,2)))-a(2,2)*( &
a(3,1)*(a(4,3)*a(5,4)-a(4,4)*a(5,3))-a(3,3)*(a(4,1)*a(5,4)-a(4,4)* &
a(5,1))+a(3,4)*(a(4,1)*a(5,3)-a(4,3)*a(5,1)))+a(2,3)*(a(3,1)*(a(4,2)* &
a(5,4)-a(4,4)*a(5,2))-a(3,2)*(a(4,1)*a(5,4)-a(4,4)*a(5,1))+a(3,4)*( &
a(4,1)*a(5,2)-a(4,2)*a(5,1)))-a(2,4)*(a(3,1)*(a(4,2)*a(5,3)-a(4,3)* &
a(5,2))-a(3,2)*(a(4,1)*a(5,3)-a(4,3)*a(5,1))+a(3,3)*(a(4,1)*a(5,2)- &
a(4,2)*a(5,1))))
do i=1,5
b(1,i) = a(1,i)
b(2,i) = a(2,i)
b(3,i) = a(3,i)
b(4,i) = a(4,i)
b(5,i) = a(5,i)
end do
a(1,1) = &
(b(2,2)*(b(3,3)*(b(4,4)*b(5,5)-b(4,5)*b(5,4))-b(3,4)*(b(4,3)*b(5,5)-b(4,5)*b(5,3))+b(3,5)*(b(4,3)*b(5,4)-b(4,4)*b(5,3)))-b(2,3)* &
(b(3,2)*(b(4,4)*b(5,5)-b(4,5)*b(5,4))-b(3,4)*(b(4,2)*b(5,5)-b(4,5)*b(5,2))+b(3,5)*(b(4,2)*b(5,4)-b(4,4)*b(5,2)))+b(2,4)* &
(b(3,2)*(b(4,3)*b(5,5)-b(4,5)*b(5,3))-b(3,3)*(b(4,2)*b(5,5)-b(4,5)*b(5,2))+b(3,5)*(b(4,2)*b(5,3)-b(4,3)*b(5,2)))-b(2,5)* &
(b(3,2)*(b(4,3)*b(5,4)-b(4,4)*b(5,3))-b(3,3)*(b(4,2)*b(5,4)-b(4,4)*b(5,2))+b(3,4)*(b(4,2)*b(5,3)-b(4,3)*b(5,2))))
a(2,1) = &
(-b(2,1)*(b(3,3)*(b(4,4)*b(5,5)-b(4,5)*b(5,4))-b(3,4)*(b(4,3)*b(5,5)-b(4,5)*b(5,3))+b(3,5)*(b(4,3)*b(5,4)-b(4,4)*b(5,3)))+b(2,3)* &
(b(3,1)*(b(4,4)*b(5,5)-b(4,5)*b(5,4))-b(3,4)*(b(4,1)*b(5,5)-b(4,5)*b(5,1))+b(3,5)*(b(4,1)*b(5,4)-b(4,4)*b(5,1)))-b(2,4)* &
(b(3,1)*(b(4,3)*b(5,5)-b(4,5)*b(5,3))-b(3,3)*(b(4,1)*b(5,5)-b(4,5)*b(5,1))+b(3,5)*(b(4,1)*b(5,3)-b(4,3)*b(5,1)))+b(2,5)* &
(b(3,1)*(b(4,3)*b(5,4)-b(4,4)*b(5,3))-b(3,3)*(b(4,1)*b(5,4)-b(4,4)*b(5,1))+b(3,4)*(b(4,1)*b(5,3)-b(4,3)*b(5,1))))
a(3,1) = &
(b(2,1)*(b(3,2)*(b(4,4)*b(5,5)-b(4,5)*b(5,4))-b(3,4)*(b(4,2)*b(5,5)-b(4,5)*b(5,2))+b(3,5)*(b(4,2)*b(5,4)-b(4,4)*b(5,2)))-b(2,2)* &
(b(3,1)*(b(4,4)*b(5,5)-b(4,5)*b(5,4))-b(3,4)*(b(4,1)*b(5,5)-b(4,5)*b(5,1))+b(3,5)*(b(4,1)*b(5,4)-b(4,4)*b(5,1)))+b(2,4)* &
(b(3,1)*(b(4,2)*b(5,5)-b(4,5)*b(5,2))-b(3,2)*(b(4,1)*b(5,5)-b(4,5)*b(5,1))+b(3,5)*(b(4,1)*b(5,2)-b(4,2)*b(5,1)))-b(2,5)* &
(b(3,1)*(b(4,2)*b(5,4)-b(4,4)*b(5,2))-b(3,2)*(b(4,1)*b(5,4)-b(4,4)*b(5,1))+b(3,4)*(b(4,1)*b(5,2)-b(4,2)*b(5,1))))
a(4,1) = &
(-b(2,1)*(b(3,2)*(b(4,3)*b(5,5)-b(4,5)*b(5,3))-b(3,3)*(b(4,2)*b(5,5)-b(4,5)*b(5,2))+b(3,5)*(b(4,2)*b(5,3)-b(4,3)*b(5,2)))+b(2,2)* &
(b(3,1)*(b(4,3)*b(5,5)-b(4,5)*b(5,3))-b(3,3)*(b(4,1)*b(5,5)-b(4,5)*b(5,1))+b(3,5)*(b(4,1)*b(5,3)-b(4,3)*b(5,1)))-b(2,3)* &
(b(3,1)*(b(4,2)*b(5,5)-b(4,5)*b(5,2))-b(3,2)*(b(4,1)*b(5,5)-b(4,5)*b(5,1))+b(3,5)*(b(4,1)*b(5,2)-b(4,2)*b(5,1)))+b(2,5)* &
(b(3,1)*(b(4,2)*b(5,3)-b(4,3)*b(5,2))-b(3,2)*(b(4,1)*b(5,3)-b(4,3)*b(5,1))+b(3,3)*(b(4,1)*b(5,2)-b(4,2)*b(5,1))))
a(5,1) = &
(b(2,1)*(b(3,2)*(b(4,3)*b(5,4)-b(4,4)*b(5,3))-b(3,3)*(b(4,2)*b(5,4)-b(4,4)*b(5,2))+b(3,4)*(b(4,2)*b(5,3)-b(4,3)*b(5,2)))-b(2,2)* &
(b(3,1)*(b(4,3)*b(5,4)-b(4,4)*b(5,3))-b(3,3)*(b(4,1)*b(5,4)-b(4,4)*b(5,1))+b(3,4)*(b(4,1)*b(5,3)-b(4,3)*b(5,1)))+b(2,3)* &
(b(3,1)*(b(4,2)*b(5,4)-b(4,4)*b(5,2))-b(3,2)*(b(4,1)*b(5,4)-b(4,4)*b(5,1))+b(3,4)*(b(4,1)*b(5,2)-b(4,2)*b(5,1)))-b(2,4)* &
(b(3,1)*(b(4,2)*b(5,3)-b(4,3)*b(5,2))-b(3,2)*(b(4,1)*b(5,3)-b(4,3)*b(5,1))+b(3,3)*(b(4,1)*b(5,2)-b(4,2)*b(5,1))))
a(1,2) = &
(-b(1,2)*(b(3,3)*(b(4,4)*b(5,5)-b(4,5)*b(5,4))-b(3,4)*(b(4,3)*b(5,5)-b(4,5)*b(5,3))+b(3,5)*(b(4,3)*b(5,4)-b(4,4)*b(5,3)))+b(1,3)* &
(b(3,2)*(b(4,4)*b(5,5)-b(4,5)*b(5,4))-b(3,4)*(b(4,2)*b(5,5)-b(4,5)*b(5,2))+b(3,5)*(b(4,2)*b(5,4)-b(4,4)*b(5,2)))-b(1,4)* &
(b(3,2)*(b(4,3)*b(5,5)-b(4,5)*b(5,3))-b(3,3)*(b(4,2)*b(5,5)-b(4,5)*b(5,2))+b(3,5)*(b(4,2)*b(5,3)-b(4,3)*b(5,2)))+b(1,5)* &
(b(3,2)*(b(4,3)*b(5,4)-b(4,4)*b(5,3))-b(3,3)*(b(4,2)*b(5,4)-b(4,4)*b(5,2))+b(3,4)*(b(4,2)*b(5,3)-b(4,3)*b(5,2))))
a(2,2) = &
(b(1,1)*(b(3,3)*(b(4,4)*b(5,5)-b(4,5)*b(5,4))-b(3,4)*(b(4,3)*b(5,5)-b(4,5)*b(5,3))+b(3,5)*(b(4,3)*b(5,4)-b(4,4)*b(5,3)))-b(1,3)* &
(b(3,1)*(b(4,4)*b(5,5)-b(4,5)*b(5,4))-b(3,4)*(b(4,1)*b(5,5)-b(4,5)*b(5,1))+b(3,5)*(b(4,1)*b(5,4)-b(4,4)*b(5,1)))+b(1,4)* &
(b(3,1)*(b(4,3)*b(5,5)-b(4,5)*b(5,3))-b(3,3)*(b(4,1)*b(5,5)-b(4,5)*b(5,1))+b(3,5)*(b(4,1)*b(5,3)-b(4,3)*b(5,1)))-b(1,5)* &
(b(3,1)*(b(4,3)*b(5,4)-b(4,4)*b(5,3))-b(3,3)*(b(4,1)*b(5,4)-b(4,4)*b(5,1))+b(3,4)*(b(4,1)*b(5,3)-b(4,3)*b(5,1))))
a(3,2) = &
(-b(1,1)*(b(3,2)*(b(4,4)*b(5,5)-b(4,5)*b(5,4))-b(3,4)*(b(4,2)*b(5,5)-b(4,5)*b(5,2))+b(3,5)*(b(4,2)*b(5,4)-b(4,4)*b(5,2)))+b(1,2)* &
(b(3,1)*(b(4,4)*b(5,5)-b(4,5)*b(5,4))-b(3,4)*(b(4,1)*b(5,5)-b(4,5)*b(5,1))+b(3,5)*(b(4,1)*b(5,4)-b(4,4)*b(5,1)))-b(1,4)* &
(b(3,1)*(b(4,2)*b(5,5)-b(4,5)*b(5,2))-b(3,2)*(b(4,1)*b(5,5)-b(4,5)*b(5,1))+b(3,5)*(b(4,1)*b(5,2)-b(4,2)*b(5,1)))+b(1,5)* &
(b(3,1)*(b(4,2)*b(5,4)-b(4,4)*b(5,2))-b(3,2)*(b(4,1)*b(5,4)-b(4,4)*b(5,1))+b(3,4)*(b(4,1)*b(5,2)-b(4,2)*b(5,1))))
a(4,2) = &
(b(1,1)*(b(3,2)*(b(4,3)*b(5,5)-b(4,5)*b(5,3))-b(3,3)*(b(4,2)*b(5,5)-b(4,5)*b(5,2))+b(3,5)*(b(4,2)*b(5,3)-b(4,3)*b(5,2)))-b(1,2)* &
(b(3,1)*(b(4,3)*b(5,5)-b(4,5)*b(5,3))-b(3,3)*(b(4,1)*b(5,5)-b(4,5)*b(5,1))+b(3,5)*(b(4,1)*b(5,3)-b(4,3)*b(5,1)))+b(1,3)* &
(b(3,1)*(b(4,2)*b(5,5)-b(4,5)*b(5,2))-b(3,2)*(b(4,1)*b(5,5)-b(4,5)*b(5,1))+b(3,5)*(b(4,1)*b(5,2)-b(4,2)*b(5,1)))-b(1,5)* &
(b(3,1)*(b(4,2)*b(5,3)-b(4,3)*b(5,2))-b(3,2)*(b(4,1)*b(5,3)-b(4,3)*b(5,1))+b(3,3)*(b(4,1)*b(5,2)-b(4,2)*b(5,1))))
a(5,2) = &
(-b(1,1)*(b(3,2)*(b(4,3)*b(5,4)-b(4,4)*b(5,3))-b(3,3)*(b(4,2)*b(5,4)-b(4,4)*b(5,2))+b(3,4)*(b(4,2)*b(5,3)-b(4,3)*b(5,2)))+b(1,2)* &
(b(3,1)*(b(4,3)*b(5,4)-b(4,4)*b(5,3))-b(3,3)*(b(4,1)*b(5,4)-b(4,4)*b(5,1))+b(3,4)*(b(4,1)*b(5,3)-b(4,3)*b(5,1)))-b(1,3)* &
(b(3,1)*(b(4,2)*b(5,4)-b(4,4)*b(5,2))-b(3,2)*(b(4,1)*b(5,4)-b(4,4)*b(5,1))+b(3,4)*(b(4,1)*b(5,2)-b(4,2)*b(5,1)))+b(1,4)* &
(b(3,1)*(b(4,2)*b(5,3)-b(4,3)*b(5,2))-b(3,2)*(b(4,1)*b(5,3)-b(4,3)*b(5,1))+b(3,3)*(b(4,1)*b(5,2)-b(4,2)*b(5,1))))
a(1,3) = &
(b(1,2)*(b(2,3)*(b(4,4)*b(5,5)-b(4,5)*b(5,4))-b(2,4)*(b(4,3)*b(5,5)-b(4,5)*b(5,3))+b(2,5)*(b(4,3)*b(5,4)-b(4,4)*b(5,3)))-b(1,3)* &
(b(2,2)*(b(4,4)*b(5,5)-b(4,5)*b(5,4))-b(2,4)*(b(4,2)*b(5,5)-b(4,5)*b(5,2))+b(2,5)*(b(4,2)*b(5,4)-b(4,4)*b(5,2)))+b(1,4)* &
(b(2,2)*(b(4,3)*b(5,5)-b(4,5)*b(5,3))-b(2,3)*(b(4,2)*b(5,5)-b(4,5)*b(5,2))+b(2,5)*(b(4,2)*b(5,3)-b(4,3)*b(5,2)))-b(1,5)* &
(b(2,2)*(b(4,3)*b(5,4)-b(4,4)*b(5,3))-b(2,3)*(b(4,2)*b(5,4)-b(4,4)*b(5,2))+b(2,4)*(b(4,2)*b(5,3)-b(4,3)*b(5,2))))
a(2,3) = &
(-b(1,1)*(b(2,3)*(b(4,4)*b(5,5)-b(4,5)*b(5,4))-b(2,4)*(b(4,3)*b(5,5)-b(4,5)*b(5,3))+b(2,5)*(b(4,3)*b(5,4)-b(4,4)*b(5,3)))+b(1,3)* &
(b(2,1)*(b(4,4)*b(5,5)-b(4,5)*b(5,4))-b(2,4)*(b(4,1)*b(5,5)-b(4,5)*b(5,1))+b(2,5)*(b(4,1)*b(5,4)-b(4,4)*b(5,1)))-b(1,4)* &
(b(2,1)*(b(4,3)*b(5,5)-b(4,5)*b(5,3))-b(2,3)*(b(4,1)*b(5,5)-b(4,5)*b(5,1))+b(2,5)*(b(4,1)*b(5,3)-b(4,3)*b(5,1)))+b(1,5)* &
(b(2,1)*(b(4,3)*b(5,4)-b(4,4)*b(5,3))-b(2,3)*(b(4,1)*b(5,4)-b(4,4)*b(5,1))+b(2,4)*(b(4,1)*b(5,3)-b(4,3)*b(5,1))))
a(3,3) = &
(b(1,1)*(b(2,2)*(b(4,4)*b(5,5)-b(4,5)*b(5,4))-b(2,4)*(b(4,2)*b(5,5)-b(4,5)*b(5,2))+b(2,5)*(b(4,2)*b(5,4)-b(4,4)*b(5,2)))-b(1,2)* &
(b(2,1)*(b(4,4)*b(5,5)-b(4,5)*b(5,4))-b(2,4)*(b(4,1)*b(5,5)-b(4,5)*b(5,1))+b(2,5)*(b(4,1)*b(5,4)-b(4,4)*b(5,1)))+b(1,4)* &
(b(2,1)*(b(4,2)*b(5,5)-b(4,5)*b(5,2))-b(2,2)*(b(4,1)*b(5,5)-b(4,5)*b(5,1))+b(2,5)*(b(4,1)*b(5,2)-b(4,2)*b(5,1)))-b(1,5)* &
(b(2,1)*(b(4,2)*b(5,4)-b(4,4)*b(5,2))-b(2,2)*(b(4,1)*b(5,4)-b(4,4)*b(5,1))+b(2,4)*(b(4,1)*b(5,2)-b(4,2)*b(5,1))))
a(4,3) = &
(-b(1,1)*(b(2,2)*(b(4,3)*b(5,5)-b(4,5)*b(5,3))-b(2,3)*(b(4,2)*b(5,5)-b(4,5)*b(5,2))+b(2,5)*(b(4,2)*b(5,3)-b(4,3)*b(5,2)))+b(1,2)* &
(b(2,1)*(b(4,3)*b(5,5)-b(4,5)*b(5,3))-b(2,3)*(b(4,1)*b(5,5)-b(4,5)*b(5,1))+b(2,5)*(b(4,1)*b(5,3)-b(4,3)*b(5,1)))-b(1,3)* &
(b(2,1)*(b(4,2)*b(5,5)-b(4,5)*b(5,2))-b(2,2)*(b(4,1)*b(5,5)-b(4,5)*b(5,1))+b(2,5)*(b(4,1)*b(5,2)-b(4,2)*b(5,1)))+b(1,5)* &
(b(2,1)*(b(4,2)*b(5,3)-b(4,3)*b(5,2))-b(2,2)*(b(4,1)*b(5,3)-b(4,3)*b(5,1))+b(2,3)*(b(4,1)*b(5,2)-b(4,2)*b(5,1))))
a(5,3) = &
(b(1,1)*(b(2,2)*(b(4,3)*b(5,4)-b(4,4)*b(5,3))-b(2,3)*(b(4,2)*b(5,4)-b(4,4)*b(5,2))+b(2,4)*(b(4,2)*b(5,3)-b(4,3)*b(5,2)))-b(1,2)* &
(b(2,1)*(b(4,3)*b(5,4)-b(4,4)*b(5,3))-b(2,3)*(b(4,1)*b(5,4)-b(4,4)*b(5,1))+b(2,4)*(b(4,1)*b(5,3)-b(4,3)*b(5,1)))+b(1,3)* &
(b(2,1)*(b(4,2)*b(5,4)-b(4,4)*b(5,2))-b(2,2)*(b(4,1)*b(5,4)-b(4,4)*b(5,1))+b(2,4)*(b(4,1)*b(5,2)-b(4,2)*b(5,1)))-b(1,4)* &
(b(2,1)*(b(4,2)*b(5,3)-b(4,3)*b(5,2))-b(2,2)*(b(4,1)*b(5,3)-b(4,3)*b(5,1))+b(2,3)*(b(4,1)*b(5,2)-b(4,2)*b(5,1))))
a(1,4) = &
(-b(1,2)*(b(2,3)*(b(3,4)*b(5,5)-b(3,5)*b(5,4))-b(2,4)*(b(3,3)*b(5,5)-b(3,5)*b(5,3))+b(2,5)*(b(3,3)*b(5,4)-b(3,4)*b(5,3)))+b(1,3)* &
(b(2,2)*(b(3,4)*b(5,5)-b(3,5)*b(5,4))-b(2,4)*(b(3,2)*b(5,5)-b(3,5)*b(5,2))+b(2,5)*(b(3,2)*b(5,4)-b(3,4)*b(5,2)))-b(1,4)* &
(b(2,2)*(b(3,3)*b(5,5)-b(3,5)*b(5,3))-b(2,3)*(b(3,2)*b(5,5)-b(3,5)*b(5,2))+b(2,5)*(b(3,2)*b(5,3)-b(3,3)*b(5,2)))+b(1,5)* &
(b(2,2)*(b(3,3)*b(5,4)-b(3,4)*b(5,3))-b(2,3)*(b(3,2)*b(5,4)-b(3,4)*b(5,2))+b(2,4)*(b(3,2)*b(5,3)-b(3,3)*b(5,2))))
a(2,4) = &
(b(1,1)*(b(2,3)*(b(3,4)*b(5,5)-b(3,5)*b(5,4))-b(2,4)*(b(3,3)*b(5,5)-b(3,5)*b(5,3))+b(2,5)*(b(3,3)*b(5,4)-b(3,4)*b(5,3)))-b(1,3)* &
(b(2,1)*(b(3,4)*b(5,5)-b(3,5)*b(5,4))-b(2,4)*(b(3,1)*b(5,5)-b(3,5)*b(5,1))+b(2,5)*(b(3,1)*b(5,4)-b(3,4)*b(5,1)))+b(1,4)* &
(b(2,1)*(b(3,3)*b(5,5)-b(3,5)*b(5,3))-b(2,3)*(b(3,1)*b(5,5)-b(3,5)*b(5,1))+b(2,5)*(b(3,1)*b(5,3)-b(3,3)*b(5,1)))-b(1,5)* &
(b(2,1)*(b(3,3)*b(5,4)-b(3,4)*b(5,3))-b(2,3)*(b(3,1)*b(5,4)-b(3,4)*b(5,1))+b(2,4)*(b(3,1)*b(5,3)-b(3,3)*b(5,1))))
a(3,4) = &
(-b(1,1)*(b(2,2)*(b(3,4)*b(5,5)-b(3,5)*b(5,4))-b(2,4)*(b(3,2)*b(5,5)-b(3,5)*b(5,2))+b(2,5)*(b(3,2)*b(5,4)-b(3,4)*b(5,2)))+b(1,2)* &
(b(2,1)*(b(3,4)*b(5,5)-b(3,5)*b(5,4))-b(2,4)*(b(3,1)*b(5,5)-b(3,5)*b(5,1))+b(2,5)*(b(3,1)*b(5,4)-b(3,4)*b(5,1)))-b(1,4)* &
(b(2,1)*(b(3,2)*b(5,5)-b(3,5)*b(5,2))-b(2,2)*(b(3,1)*b(5,5)-b(3,5)*b(5,1))+b(2,5)*(b(3,1)*b(5,2)-b(3,2)*b(5,1)))+b(1,5)* &
(b(2,1)*(b(3,2)*b(5,4)-b(3,4)*b(5,2))-b(2,2)*(b(3,1)*b(5,4)-b(3,4)*b(5,1))+b(2,4)*(b(3,1)*b(5,2)-b(3,2)*b(5,1))))
a(4,4) = &
(b(1,1)*(b(2,2)*(b(3,3)*b(5,5)-b(3,5)*b(5,3))-b(2,3)*(b(3,2)*b(5,5)-b(3,5)*b(5,2))+b(2,5)*(b(3,2)*b(5,3)-b(3,3)*b(5,2)))-b(1,2)* &
(b(2,1)*(b(3,3)*b(5,5)-b(3,5)*b(5,3))-b(2,3)*(b(3,1)*b(5,5)-b(3,5)*b(5,1))+b(2,5)*(b(3,1)*b(5,3)-b(3,3)*b(5,1)))+b(1,3)* &
(b(2,1)*(b(3,2)*b(5,5)-b(3,5)*b(5,2))-b(2,2)*(b(3,1)*b(5,5)-b(3,5)*b(5,1))+b(2,5)*(b(3,1)*b(5,2)-b(3,2)*b(5,1)))-b(1,5)* &
(b(2,1)*(b(3,2)*b(5,3)-b(3,3)*b(5,2))-b(2,2)*(b(3,1)*b(5,3)-b(3,3)*b(5,1))+b(2,3)*(b(3,1)*b(5,2)-b(3,2)*b(5,1))))
a(5,4) = &
(-b(1,1)*(b(2,2)*(b(3,3)*b(5,4)-b(3,4)*b(5,3))-b(2,3)*(b(3,2)*b(5,4)-b(3,4)*b(5,2))+b(2,4)*(b(3,2)*b(5,3)-b(3,3)*b(5,2)))+b(1,2)* &
(b(2,1)*(b(3,3)*b(5,4)-b(3,4)*b(5,3))-b(2,3)*(b(3,1)*b(5,4)-b(3,4)*b(5,1))+b(2,4)*(b(3,1)*b(5,3)-b(3,3)*b(5,1)))-b(1,3)* &
(b(2,1)*(b(3,2)*b(5,4)-b(3,4)*b(5,2))-b(2,2)*(b(3,1)*b(5,4)-b(3,4)*b(5,1))+b(2,4)*(b(3,1)*b(5,2)-b(3,2)*b(5,1)))+b(1,4)* &
(b(2,1)*(b(3,2)*b(5,3)-b(3,3)*b(5,2))-b(2,2)*(b(3,1)*b(5,3)-b(3,3)*b(5,1))+b(2,3)*(b(3,1)*b(5,2)-b(3,2)*b(5,1))))
a(1,5) = &
(b(1,2)*(b(2,3)*(b(3,4)*b(4,5)-b(3,5)*b(4,4))-b(2,4)*(b(3,3)*b(4,5)-b(3,5)*b(4,3))+b(2,5)*(b(3,3)*b(4,4)-b(3,4)*b(4,3)))-b(1,3)* &
(b(2,2)*(b(3,4)*b(4,5)-b(3,5)*b(4,4))-b(2,4)*(b(3,2)*b(4,5)-b(3,5)*b(4,2))+b(2,5)*(b(3,2)*b(4,4)-b(3,4)*b(4,2)))+b(1,4)* &
(b(2,2)*(b(3,3)*b(4,5)-b(3,5)*b(4,3))-b(2,3)*(b(3,2)*b(4,5)-b(3,5)*b(4,2))+b(2,5)*(b(3,2)*b(4,3)-b(3,3)*b(4,2)))-b(1,5)* &
(b(2,2)*(b(3,3)*b(4,4)-b(3,4)*b(4,3))-b(2,3)*(b(3,2)*b(4,4)-b(3,4)*b(4,2))+b(2,4)*(b(3,2)*b(4,3)-b(3,3)*b(4,2))))
a(2,5) = &
(-b(1,1)*(b(2,3)*(b(3,4)*b(4,5)-b(3,5)*b(4,4))-b(2,4)*(b(3,3)*b(4,5)-b(3,5)*b(4,3))+b(2,5)*(b(3,3)*b(4,4)-b(3,4)*b(4,3)))+b(1,3)* &
(b(2,1)*(b(3,4)*b(4,5)-b(3,5)*b(4,4))-b(2,4)*(b(3,1)*b(4,5)-b(3,5)*b(4,1))+b(2,5)*(b(3,1)*b(4,4)-b(3,4)*b(4,1)))-b(1,4)* &
(b(2,1)*(b(3,3)*b(4,5)-b(3,5)*b(4,3))-b(2,3)*(b(3,1)*b(4,5)-b(3,5)*b(4,1))+b(2,5)*(b(3,1)*b(4,3)-b(3,3)*b(4,1)))+b(1,5)* &
(b(2,1)*(b(3,3)*b(4,4)-b(3,4)*b(4,3))-b(2,3)*(b(3,1)*b(4,4)-b(3,4)*b(4,1))+b(2,4)*(b(3,1)*b(4,3)-b(3,3)*b(4,1))))
a(3,5) = &
(b(1,1)*(b(2,2)*(b(3,4)*b(4,5)-b(3,5)*b(4,4))-b(2,4)*(b(3,2)*b(4,5)-b(3,5)*b(4,2))+b(2,5)*(b(3,2)*b(4,4)-b(3,4)*b(4,2)))-b(1,2)* &
(b(2,1)*(b(3,4)*b(4,5)-b(3,5)*b(4,4))-b(2,4)*(b(3,1)*b(4,5)-b(3,5)*b(4,1))+b(2,5)*(b(3,1)*b(4,4)-b(3,4)*b(4,1)))+b(1,4)* &
(b(2,1)*(b(3,2)*b(4,5)-b(3,5)*b(4,2))-b(2,2)*(b(3,1)*b(4,5)-b(3,5)*b(4,1))+b(2,5)*(b(3,1)*b(4,2)-b(3,2)*b(4,1)))-b(1,5)* &
(b(2,1)*(b(3,2)*b(4,4)-b(3,4)*b(4,2))-b(2,2)*(b(3,1)*b(4,4)-b(3,4)*b(4,1))+b(2,4)*(b(3,1)*b(4,2)-b(3,2)*b(4,1))))
a(4,5) = &
(-b(1,1)*(b(2,2)*(b(3,3)*b(4,5)-b(3,5)*b(4,3))-b(2,3)*(b(3,2)*b(4,5)-b(3,5)*b(4,2))+b(2,5)*(b(3,2)*b(4,3)-b(3,3)*b(4,2)))+b(1,2)* &
(b(2,1)*(b(3,3)*b(4,5)-b(3,5)*b(4,3))-b(2,3)*(b(3,1)*b(4,5)-b(3,5)*b(4,1))+b(2,5)*(b(3,1)*b(4,3)-b(3,3)*b(4,1)))-b(1,3)* &
(b(2,1)*(b(3,2)*b(4,5)-b(3,5)*b(4,2))-b(2,2)*(b(3,1)*b(4,5)-b(3,5)*b(4,1))+b(2,5)*(b(3,1)*b(4,2)-b(3,2)*b(4,1)))+b(1,5)* &
(b(2,1)*(b(3,2)*b(4,3)-b(3,3)*b(4,2))-b(2,2)*(b(3,1)*b(4,3)-b(3,3)*b(4,1))+b(2,3)*(b(3,1)*b(4,2)-b(3,2)*b(4,1))))
a(5,5) = &
(b(1,1)*(b(2,2)*(b(3,3)*b(4,4)-b(3,4)*b(4,3))-b(2,3)*(b(3,2)*b(4,4)-b(3,4)*b(4,2))+b(2,4)*(b(3,2)*b(4,3)-b(3,3)*b(4,2)))-b(1,2)* &
(b(2,1)*(b(3,3)*b(4,4)-b(3,4)*b(4,3))-b(2,3)*(b(3,1)*b(4,4)-b(3,4)*b(4,1))+b(2,4)*(b(3,1)*b(4,3)-b(3,3)*b(4,1)))+b(1,3)* &
(b(2,1)*(b(3,2)*b(4,4)-b(3,4)*b(4,2))-b(2,2)*(b(3,1)*b(4,4)-b(3,4)*b(4,1))+b(2,4)*(b(3,1)*b(4,2)-b(3,2)*b(4,1)))-b(1,4)* &
(b(2,1)*(b(3,2)*b(4,3)-b(3,3)*b(4,2))-b(2,2)*(b(3,1)*b(4,3)-b(3,3)*b(4,1))+b(2,3)*(b(3,1)*b(4,2)-b(3,2)*b(4,1))))
end
#+end_src
*** C interface :noexport:
#+CALL: generate_c_interface(table=qmckl_adjoint_args,rettyp="qmckl_exit_code",fname="qmckl_adjoint")
#+RESULTS:
#+begin_src f90 :tangle (eval f) :comments org :exports none
integer(c_int32_t) function qmckl_adjoint &
(context, m, n, lda, A, det_l) &
bind(C) result(info)
use, intrinsic :: iso_c_binding
implicit none
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 :: lda
real (c_double ) , intent(inout) :: A(lda,*)
real (c_double ) , intent(inout) :: det_l
integer(c_int32_t), external :: qmckl_adjoint_f
info = qmckl_adjoint_f &
(context, m, n, lda, A, det_l)
end function qmckl_adjoint
#+end_src
#+CALL: generate_f_interface(table=qmckl_adjoint_args,rettyp="qmckl_exit_code",fname="qmckl_adjoint")
#+RESULTS:
#+begin_src f90 :tangle (eval fh_func) :comments org :exports none
interface
integer(c_int32_t) function qmckl_adjoint &
(context, m, n, lda, A, det_l) &
bind(C)
use, intrinsic :: iso_c_binding
import
implicit none
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 :: lda
real (c_double ) , intent(inout) :: A(lda,*)
real (c_double ) , intent(inout) :: det_l
end function qmckl_adjoint
end interface
#+end_src
*** Test :noexport:
#+begin_src f90 :tangle (eval f_test)
integer(qmckl_exit_code) function test_qmckl_adjoint(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_adjoint = qmckl_adjoint(context, m, k, LDA, A, det_l)
if (test_qmckl_adjoint /= QMCKL_SUCCESS) return
test_qmckl_adjoint = 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_adjoint = QMCKL_SUCCESS
endif
deallocate(A,C)
end function test_qmckl_adjoint
#+end_src
#+begin_src c :comments link :tangle (eval c_test)
qmckl_exit_code test_qmckl_adjoint(qmckl_context context);
assert(QMCKL_SUCCESS == test_qmckl_adjoint(context));
#+end_src
* End of files :noexport:
#+begin_src c :comments link :tangle (eval c_test)
assert (qmckl_context_destroy(context) == QMCKL_SUCCESS);
return 0;
}
#+end_src
# -*- mode: org -*-
# vim: syntax=c
Reference in New Issue View Git Blame Copy Permalink