mirror of
https://github.com/TREX-CoE/qmckl.git
synced 2024-08-18 03:11:44 +02:00
1187 lines
48 KiB
Org Mode
1187 lines
48 KiB
Org Mode
#+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 multiplication:
|
|
|
|
\[
|
|
C_{ij} = \beta C_{ij} + \alpha \sum_{k} A_{ik} \cdot B_{kj}
|
|
\]
|
|
|
|
# TODO: Add description about the external library dependence.
|
|
|
|
#+NAME: qmckl_dgemm_args
|
|
| Variable | Type | In/Out | Description |
|
|
|-----------+-----------------+--------+---------------------------------------|
|
|
| ~context~ | ~qmckl_context~ | in | Global state |
|
|
| ~TransA~ | ~bool~ | in | Number of rows of the input matrix |
|
|
| ~TransB~ | ~bool~ | in | Number of rows of the input matrix |
|
|
| ~m~ | ~int64_t~ | in | Number of rows of the input matrix |
|
|
| ~n~ | ~int64_t~ | in | Number of columns of the input matrix |
|
|
| ~k~ | ~int64_t~ | in | Number of columns of the input matrix |
|
|
| ~alpha~ | ~double~ | in | Number of columns of the input matrix |
|
|
| ~A~ | ~double[][lda]~ | in | Array containing the matrix $A$ |
|
|
| ~lda~ | ~int64_t~ | in | Leading dimension of array ~A~ |
|
|
| ~B~ | ~double[][ldb]~ | in | Array containing the matrix $B$ |
|
|
| ~ldb~ | ~int64_t~ | in | Leading dimension of array ~B~ |
|
|
| ~beta~ | ~double~ | in | Array containing the matrix $B$ |
|
|
| ~C~ | ~double[][ldc]~ | out | Array containing the matrix $B$ |
|
|
| ~ldc~ | ~int64_t~ | in | Leading dimension of array ~B~ |
|
|
|
|
Requirements:
|
|
|
|
- ~context~ is not ~QMCKL_NULL_CONTEXT~
|
|
- ~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
|
|
|
|
#+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
|
|
|
|
#+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
|
|
double precision , intent(in) :: alpha, beta
|
|
integer*8 , intent(in) :: lda
|
|
double precision , intent(in) :: A(lda,*)
|
|
integer*8 , intent(in) :: ldb
|
|
double precision , intent(in) :: B(ldb,*)
|
|
integer*8 , intent(in) :: ldc
|
|
double precision , intent(out) :: C(ldc,*)
|
|
|
|
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_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
|
|
|
|
call dgemm(transA, transB, int(m,4), int(n,4), int(k,4), &
|
|
alpha, A, int(LDA,4), B, int(LDB,4), beta, C, int(LDC,4))
|
|
|
|
end function qmckl_dgemm_f
|
|
#+end_src
|
|
|
|
*** C interface :noexport:
|
|
|
|
#+CALL: generate_c_interface(table=qmckl_dgemm_args,rettyp="qmckl_exit_code",fname="qmckl_dgemm")
|
|
|
|
#+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
|
|
real (c_double ) , intent(in) :: A(lda,*)
|
|
integer (c_int64_t) , intent(in) , value :: lda
|
|
real (c_double ) , intent(in) :: B(ldb,*)
|
|
integer (c_int64_t) , intent(in) , value :: ldb
|
|
real (c_double ) , intent(in) , value :: beta
|
|
real (c_double ) , intent(out) :: C(ldc,*)
|
|
integer (c_int64_t) , intent(in) , value :: 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
|
|
real (c_double ) , intent(in) :: A(lda,*)
|
|
integer (c_int64_t) , intent(in) , value :: lda
|
|
real (c_double ) , intent(in) :: B(ldb,*)
|
|
integer (c_int64_t) , intent(in) , value :: ldb
|
|
real (c_double ) , intent(in) , value :: beta
|
|
real (c_double ) , intent(out) :: C(ldc,*)
|
|
integer (c_int64_t) , intent(in) , value :: 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-12) then
|
|
test_qmckl_dgemm = QMCKL_SUCCESS
|
|
endif
|
|
|
|
deallocate(A,B,C,D)
|
|
|
|
end function test_qmckl_dgemm
|
|
#+end_src
|
|
|
|
#+begin_src c :comments link :tangle (eval c_test)
|
|
qmckl_exit_code test_qmckl_dgemm(qmckl_context context);
|
|
assert(QMCKL_SUCCESS == test_qmckl_dgemm(context));
|
|
#+end_src
|
|
|
|
** ~qmckl_adjugate~
|
|
|
|
Given a matrix $\mathbf{A}$, the adjugate matrix
|
|
$\text{adj}(\mathbf{A})$ is the transpose of the cofactors matrix
|
|
of $\mathbf{A}$.
|
|
|
|
\[
|
|
\mathbf{B} = \text{adj}(\mathbf{A}) = \text{det}(\mathbf{A}) \, \mathbf{A}^{-1}
|
|
\]
|
|
|
|
See also: https://en.wikipedia.org/wiki/Adjugate_matrix
|
|
|
|
#+NAME: qmckl_adjugate_args
|
|
| Variable | Type | In/Out | Description |
|
|
|-----------+-----------------+--------+------------------------------------------------|
|
|
| ~context~ | ~qmckl_context~ | in | Global state |
|
|
| ~n~ | ~int64_t~ | in | Number of rows and columns of the input matrix |
|
|
| ~A~ | ~double[][lda]~ | in | Array containing the $n \times n$ matrix $A$ |
|
|
| ~lda~ | ~int64_t~ | in | Leading dimension of array ~A~ |
|
|
| ~B~ | ~double[][ldb]~ | out | Adjugate of $A$ |
|
|
| ~ldb~ | ~int64_t~ | in | Leading dimension of array ~B~ |
|
|
| ~det_l~ | ~double~ | inout | determinant of $A$ |
|
|
|
|
Requirements:
|
|
|
|
- ~context~ is not ~QMCKL_NULL_CONTEXT~
|
|
- ~n > 0~
|
|
- ~lda >= m~
|
|
- ~A~ is allocated with at least $m \times m \times 8$ bytes
|
|
- ~ldb >= m~
|
|
- ~B~ is allocated with at least $m \times m \times 8$ bytes
|
|
|
|
#+CALL: generate_c_header(table=qmckl_adjugate_args,rettyp="qmckl_exit_code",fname="qmckl_adjugate")
|
|
|
|
#+RESULTS:
|
|
#+begin_src c :tangle (eval h_func) :comments org
|
|
qmckl_exit_code qmckl_adjugate (
|
|
const qmckl_context context,
|
|
const int64_t n,
|
|
const double* A,
|
|
const int64_t lda,
|
|
double* const B,
|
|
const int64_t ldb,
|
|
double* det_l );
|
|
#+end_src
|
|
|
|
For small matrices (\le 5\times 5), we use specialized routines
|
|
for performance motivations. For larger sizes, we rely on the
|
|
LAPACK library.
|
|
|
|
#+begin_src f90 :tangle (eval f)
|
|
integer function qmckl_adjugate_f(context, na, A, LDA, B, ldb, det_l) &
|
|
result(info)
|
|
use qmckl
|
|
implicit none
|
|
integer(qmckl_context) , intent(in) :: context
|
|
double precision, intent(in) :: A (LDA,*)
|
|
integer*8, intent(in) :: LDA
|
|
double precision, intent(out) :: B (LDB,*)
|
|
integer*8, intent(in) :: LDB
|
|
integer*8, intent(in) :: na
|
|
double precision, intent(inout) :: det_l
|
|
|
|
integer :: i,j
|
|
|
|
info = QMCKL_SUCCESS
|
|
|
|
if (context == QMCKL_NULL_CONTEXT) then
|
|
info = QMCKL_INVALID_CONTEXT
|
|
return
|
|
endif
|
|
|
|
if (na <= 0_8) then
|
|
info = QMCKL_INVALID_ARG_2
|
|
return
|
|
endif
|
|
|
|
if (LDA <= 0_8) then
|
|
info = QMCKL_INVALID_ARG_4
|
|
return
|
|
endif
|
|
|
|
if (LDA < na) then
|
|
info = QMCKL_INVALID_ARG_4
|
|
return
|
|
endif
|
|
|
|
select case (na)
|
|
case (5)
|
|
call adjugate5(A,LDA,B,LDB,na,det_l)
|
|
case (4)
|
|
call adjugate4(A,LDA,B,LDB,na,det_l)
|
|
case (3)
|
|
call adjugate3(A,LDA,B,LDB,na,det_l)
|
|
case (2)
|
|
call adjugate2(A,LDA,B,LDB,na,det_l)
|
|
case (1)
|
|
det_l = a(1,1)
|
|
b(1,1) = 1.d0
|
|
case default
|
|
call adjugate_general(context, na, A, LDA, B, LDB, det_l)
|
|
end select
|
|
|
|
end function qmckl_adjugate_f
|
|
#+end_src
|
|
|
|
#+begin_src f90 :tangle (eval f) :exports none
|
|
subroutine adjugate2(A,LDA,B,LDB,na,det_l)
|
|
implicit none
|
|
double precision, intent(in) :: A (LDA,na)
|
|
double precision, intent(out) :: B (LDA,na)
|
|
integer*8, intent(in) :: LDA, LDB
|
|
integer*8, intent(in) :: na
|
|
double precision, intent(inout) :: det_l
|
|
|
|
double precision :: C(2,2)
|
|
|
|
call cofactor2(A,LDA,C,2,na,det_l)
|
|
|
|
B(1,1) = C(1,1)
|
|
B(2,1) = C(1,2)
|
|
B(1,2) = C(2,1)
|
|
B(2,2) = C(2,2)
|
|
|
|
end subroutine adjugate2
|
|
|
|
subroutine adjugate3(a,LDA,B,LDB,na,det_l)
|
|
implicit none
|
|
double precision, intent(in) :: A (LDA,na)
|
|
double precision, intent(out) :: B (LDA,na)
|
|
integer*8, intent(in) :: LDA, LDB
|
|
integer*8, intent(in) :: na
|
|
double precision, intent(inout) :: det_l
|
|
|
|
double precision :: C(4,3)
|
|
|
|
call cofactor3(A,LDA,C,4,na,det_l)
|
|
|
|
B(1,1) = C(1,1)
|
|
B(1,2) = C(2,1)
|
|
B(1,3) = C(3,1)
|
|
B(2,1) = C(1,2)
|
|
B(2,2) = C(2,2)
|
|
B(2,3) = C(3,2)
|
|
B(3,1) = C(1,3)
|
|
B(3,2) = C(2,3)
|
|
B(3,3) = C(3,3)
|
|
|
|
end subroutine adjugate3
|
|
|
|
subroutine adjugate4(a,LDA,B,LDB,na,det_l)
|
|
implicit none
|
|
double precision, intent(in) :: A (LDA,na)
|
|
double precision, intent(out) :: B (LDA,na)
|
|
integer*8, intent(in) :: LDA, LDB
|
|
integer*8, intent(in) :: na
|
|
double precision, intent(inout) :: det_l
|
|
|
|
double precision :: C(4,4)
|
|
|
|
call cofactor4(A,LDA,B,4,na,det_l)
|
|
|
|
B(1,1) = C(1,1)
|
|
B(1,2) = C(2,1)
|
|
B(1,3) = C(3,1)
|
|
B(1,4) = C(4,1)
|
|
B(2,1) = C(1,2)
|
|
B(2,2) = C(2,2)
|
|
B(2,3) = C(3,2)
|
|
B(2,4) = C(4,2)
|
|
B(3,1) = C(1,3)
|
|
B(3,2) = C(2,3)
|
|
B(3,3) = C(3,3)
|
|
B(3,4) = C(4,3)
|
|
B(4,1) = C(1,4)
|
|
B(4,2) = C(2,4)
|
|
B(4,3) = C(3,4)
|
|
B(4,4) = C(4,4)
|
|
|
|
end subroutine adjugate4
|
|
|
|
subroutine adjugate5(A,LDA,B,LDB,na,det_l)
|
|
implicit none
|
|
double precision, intent(in) :: A (LDA,na)
|
|
double precision, intent(out) :: B (LDA,na)
|
|
integer*8, intent(in) :: LDA, LDB
|
|
integer*8, intent(in) :: na
|
|
double precision, intent(inout) :: det_l
|
|
|
|
double precision :: C(8,5)
|
|
|
|
call cofactor5(A,LDA,C,8,na,det_l)
|
|
|
|
B(1,1) = C(1,1)
|
|
B(1,2) = C(2,1)
|
|
B(1,3) = C(3,1)
|
|
B(1,4) = C(4,1)
|
|
B(1,5) = C(5,1)
|
|
B(2,1) = C(1,2)
|
|
B(2,2) = C(2,2)
|
|
B(2,3) = C(3,2)
|
|
B(2,4) = C(4,2)
|
|
B(2,5) = C(5,2)
|
|
B(3,1) = C(1,3)
|
|
B(3,2) = C(2,3)
|
|
B(3,3) = C(3,3)
|
|
B(3,4) = C(4,3)
|
|
B(3,5) = C(5,3)
|
|
B(4,1) = C(1,4)
|
|
B(4,2) = C(2,4)
|
|
B(4,3) = C(3,4)
|
|
B(4,4) = C(4,4)
|
|
B(4,5) = C(5,4)
|
|
B(5,1) = C(1,5)
|
|
B(5,2) = C(2,5)
|
|
B(5,3) = C(3,5)
|
|
B(5,4) = C(4,5)
|
|
B(5,5) = C(5,5)
|
|
|
|
end subroutine adjugate5
|
|
|
|
subroutine cofactor2(a,LDA,b,LDB,na,det_l)
|
|
implicit none
|
|
double precision, intent(in) :: A (LDA,na)
|
|
double precision, intent(out) :: B (LDA,na)
|
|
integer*8, intent(in) :: LDA, LDB
|
|
integer*8 :: na
|
|
double precision :: det_l
|
|
|
|
det_l = a(1,1)*a(2,2) - a(1,2)*a(2,1)
|
|
b(1,1) = a(2,2)
|
|
b(2,1) = -a(2,1)
|
|
b(1,2) = -a(1,2)
|
|
b(2,2) = a(1,1)
|
|
end subroutine cofactor2
|
|
|
|
subroutine cofactor3(a,LDA,b,LDB,na,det_l)
|
|
implicit none
|
|
double precision, intent(in) :: A (LDA,na)
|
|
double precision, intent(out) :: B (LDA,na)
|
|
integer*8, intent(in) :: LDA, LDB
|
|
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))
|
|
|
|
b(1,1) = a(2,2)*a(3,3) - a(2,3)*a(3,2)
|
|
b(2,1) = a(2,3)*a(3,1) - a(2,1)*a(3,3)
|
|
b(3,1) = a(2,1)*a(3,2) - a(2,2)*a(3,1)
|
|
|
|
b(1,2) = a(1,3)*a(3,2) - a(1,2)*a(3,3)
|
|
b(2,2) = a(1,1)*a(3,3) - a(1,3)*a(3,1)
|
|
b(3,2) = a(1,2)*a(3,1) - a(1,1)*a(3,2)
|
|
|
|
b(1,3) = a(1,2)*a(2,3) - a(1,3)*a(2,2)
|
|
b(2,3) = a(1,3)*a(2,1) - a(1,1)*a(2,3)
|
|
b(3,3) = a(1,1)*a(2,2) - a(1,2)*a(2,1)
|
|
|
|
end subroutine cofactor3
|
|
|
|
subroutine cofactor4(a,LDA,b,LDB,na,det_l)
|
|
implicit none
|
|
double precision, intent(in) :: A (LDA,na)
|
|
double precision, intent(out) :: B (LDA,na)
|
|
integer*8, intent(in) :: LDA, LDB
|
|
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)))
|
|
|
|
b(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))
|
|
b(2,1) = -a(2,1)*(a(3,3)*a(4,4)-a(3,4)*a(4,3))+a(2,3)*(a(3,1)*a(4,4)-a(3,4)*a(4,1))-a(2,4)*(a(3,1)*a(4,3)-a(3,3)*a(4,1))
|
|
b(3,1) = a(2,1)*(a(3,2)*a(4,4)-a(3,4)*a(4,2))-a(2,2)*(a(3,1)*a(4,4)-a(3,4)*a(4,1))+a(2,4)*(a(3,1)*a(4,2)-a(3,2)*a(4,1))
|
|
b(4,1) = -a(2,1)*(a(3,2)*a(4,3)-a(3,3)*a(4,2))+a(2,2)*(a(3,1)*a(4,3)-a(3,3)*a(4,1))-a(2,3)*(a(3,1)*a(4,2)-a(3,2)*a(4,1))
|
|
|
|
b(1,2) = -a(1,2)*(a(3,3)*a(4,4)-a(3,4)*a(4,3))+a(1,3)*(a(3,2)*a(4,4)-a(3,4)*a(4,2))-a(1,4)*(a(3,2)*a(4,3)-a(3,3)*a(4,2))
|
|
b(2,2) = a(1,1)*(a(3,3)*a(4,4)-a(3,4)*a(4,3))-a(1,3)*(a(3,1)*a(4,4)-a(3,4)*a(4,1))+a(1,4)*(a(3,1)*a(4,3)-a(3,3)*a(4,1))
|
|
b(3,2) = -a(1,1)*(a(3,2)*a(4,4)-a(3,4)*a(4,2))+a(1,2)*(a(3,1)*a(4,4)-a(3,4)*a(4,1))-a(1,4)*(a(3,1)*a(4,2)-a(3,2)*a(4,1))
|
|
b(4,2) = a(1,1)*(a(3,2)*a(4,3)-a(3,3)*a(4,2))-a(1,2)*(a(3,1)*a(4,3)-a(3,3)*a(4,1))+a(1,3)*(a(3,1)*a(4,2)-a(3,2)*a(4,1))
|
|
|
|
b(1,3) = a(1,2)*(a(2,3)*a(4,4)-a(2,4)*a(4,3))-a(1,3)*(a(2,2)*a(4,4)-a(2,4)*a(4,2))+a(1,4)*(a(2,2)*a(4,3)-a(2,3)*a(4,2))
|
|
b(2,3) = -a(1,1)*(a(2,3)*a(4,4)-a(2,4)*a(4,3))+a(1,3)*(a(2,1)*a(4,4)-a(2,4)*a(4,1))-a(1,4)*(a(2,1)*a(4,3)-a(2,3)*a(4,1))
|
|
b(3,3) = a(1,1)*(a(2,2)*a(4,4)-a(2,4)*a(4,2))-a(1,2)*(a(2,1)*a(4,4)-a(2,4)*a(4,1))+a(1,4)*(a(2,1)*a(4,2)-a(2,2)*a(4,1))
|
|
b(4,3) = -a(1,1)*(a(2,2)*a(4,3)-a(2,3)*a(4,2))+a(1,2)*(a(2,1)*a(4,3)-a(2,3)*a(4,1))-a(1,3)*(a(2,1)*a(4,2)-a(2,2)*a(4,1))
|
|
|
|
b(1,4) = -a(1,2)*(a(2,3)*a(3,4)-a(2,4)*a(3,3))+a(1,3)*(a(2,2)*a(3,4)-a(2,4)*a(3,2))-a(1,4)*(a(2,2)*a(3,3)-a(2,3)*a(3,2))
|
|
b(2,4) = a(1,1)*(a(2,3)*a(3,4)-a(2,4)*a(3,3))-a(1,3)*(a(2,1)*a(3,4)-a(2,4)*a(3,1))+a(1,4)*(a(2,1)*a(3,3)-a(2,3)*a(3,1))
|
|
b(3,4) = -a(1,1)*(a(2,2)*a(3,4)-a(2,4)*a(3,2))+a(1,2)*(a(2,1)*a(3,4)-a(2,4)*a(3,1))-a(1,4)*(a(2,1)*a(3,2)-a(2,2)*a(3,1))
|
|
b(4,4) = a(1,1)*(a(2,2)*a(3,3)-a(2,3)*a(3,2))-a(1,2)*(a(2,1)*a(3,3)-a(2,3)*a(3,1))+a(1,3)*(a(2,1)*a(3,2)-a(2,2)*a(3,1))
|
|
|
|
end subroutine cofactor4
|
|
|
|
subroutine cofactor5(A,LDA,B,LDB,na,det_l)
|
|
implicit none
|
|
double precision, intent(in) :: A (LDA,na)
|
|
double precision, intent(out) :: B (LDA,na)
|
|
integer*8, intent(in) :: LDA, LDB
|
|
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))))
|
|
|
|
b(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))))
|
|
b(2,1) = &
|
|
(-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))))
|
|
b(3,1) = &
|
|
(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))))
|
|
b(4,1) = &
|
|
(-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))))
|
|
b(5,1) = &
|
|
(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))))
|
|
|
|
b(1,2) = &
|
|
(-a(1,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(1,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(1,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(1,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))))
|
|
b(2,2) = &
|
|
(a(1,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(1,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(1,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(1,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))))
|
|
b(3,2) = &
|
|
(-a(1,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(1,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(1,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(1,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))))
|
|
b(4,2) = &
|
|
(a(1,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(1,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(1,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(1,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))))
|
|
b(5,2) = &
|
|
(-a(1,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(1,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(1,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(1,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))))
|
|
|
|
b(1,3) = &
|
|
(a(1,2)*(a(2,3)*(a(4,4)*a(5,5)-a(4,5)*a(5,4))-a(2,4)*(a(4,3)*a(5,5)-a(4,5)*a(5,3))+a(2,5)*(a(4,3)*a(5,4)-a(4,4)*a(5,3)))-a(1,3)* &
|
|
(a(2,2)*(a(4,4)*a(5,5)-a(4,5)*a(5,4))-a(2,4)*(a(4,2)*a(5,5)-a(4,5)*a(5,2))+a(2,5)*(a(4,2)*a(5,4)-a(4,4)*a(5,2)))+a(1,4)* &
|
|
(a(2,2)*(a(4,3)*a(5,5)-a(4,5)*a(5,3))-a(2,3)*(a(4,2)*a(5,5)-a(4,5)*a(5,2))+a(2,5)*(a(4,2)*a(5,3)-a(4,3)*a(5,2)))-a(1,5)* &
|
|
(a(2,2)*(a(4,3)*a(5,4)-a(4,4)*a(5,3))-a(2,3)*(a(4,2)*a(5,4)-a(4,4)*a(5,2))+a(2,4)*(a(4,2)*a(5,3)-a(4,3)*a(5,2))))
|
|
b(2,3) = &
|
|
(-a(1,1)*(a(2,3)*(a(4,4)*a(5,5)-a(4,5)*a(5,4))-a(2,4)*(a(4,3)*a(5,5)-a(4,5)*a(5,3))+a(2,5)*(a(4,3)*a(5,4)-a(4,4)*a(5,3)))+a(1,3)* &
|
|
(a(2,1)*(a(4,4)*a(5,5)-a(4,5)*a(5,4))-a(2,4)*(a(4,1)*a(5,5)-a(4,5)*a(5,1))+a(2,5)*(a(4,1)*a(5,4)-a(4,4)*a(5,1)))-a(1,4)* &
|
|
(a(2,1)*(a(4,3)*a(5,5)-a(4,5)*a(5,3))-a(2,3)*(a(4,1)*a(5,5)-a(4,5)*a(5,1))+a(2,5)*(a(4,1)*a(5,3)-a(4,3)*a(5,1)))+a(1,5)* &
|
|
(a(2,1)*(a(4,3)*a(5,4)-a(4,4)*a(5,3))-a(2,3)*(a(4,1)*a(5,4)-a(4,4)*a(5,1))+a(2,4)*(a(4,1)*a(5,3)-a(4,3)*a(5,1))))
|
|
b(3,3) = &
|
|
(a(1,1)*(a(2,2)*(a(4,4)*a(5,5)-a(4,5)*a(5,4))-a(2,4)*(a(4,2)*a(5,5)-a(4,5)*a(5,2))+a(2,5)*(a(4,2)*a(5,4)-a(4,4)*a(5,2)))-a(1,2)* &
|
|
(a(2,1)*(a(4,4)*a(5,5)-a(4,5)*a(5,4))-a(2,4)*(a(4,1)*a(5,5)-a(4,5)*a(5,1))+a(2,5)*(a(4,1)*a(5,4)-a(4,4)*a(5,1)))+a(1,4)* &
|
|
(a(2,1)*(a(4,2)*a(5,5)-a(4,5)*a(5,2))-a(2,2)*(a(4,1)*a(5,5)-a(4,5)*a(5,1))+a(2,5)*(a(4,1)*a(5,2)-a(4,2)*a(5,1)))-a(1,5)* &
|
|
(a(2,1)*(a(4,2)*a(5,4)-a(4,4)*a(5,2))-a(2,2)*(a(4,1)*a(5,4)-a(4,4)*a(5,1))+a(2,4)*(a(4,1)*a(5,2)-a(4,2)*a(5,1))))
|
|
b(4,3) = &
|
|
(-a(1,1)*(a(2,2)*(a(4,3)*a(5,5)-a(4,5)*a(5,3))-a(2,3)*(a(4,2)*a(5,5)-a(4,5)*a(5,2))+a(2,5)*(a(4,2)*a(5,3)-a(4,3)*a(5,2)))+a(1,2)* &
|
|
(a(2,1)*(a(4,3)*a(5,5)-a(4,5)*a(5,3))-a(2,3)*(a(4,1)*a(5,5)-a(4,5)*a(5,1))+a(2,5)*(a(4,1)*a(5,3)-a(4,3)*a(5,1)))-a(1,3)* &
|
|
(a(2,1)*(a(4,2)*a(5,5)-a(4,5)*a(5,2))-a(2,2)*(a(4,1)*a(5,5)-a(4,5)*a(5,1))+a(2,5)*(a(4,1)*a(5,2)-a(4,2)*a(5,1)))+a(1,5)* &
|
|
(a(2,1)*(a(4,2)*a(5,3)-a(4,3)*a(5,2))-a(2,2)*(a(4,1)*a(5,3)-a(4,3)*a(5,1))+a(2,3)*(a(4,1)*a(5,2)-a(4,2)*a(5,1))))
|
|
b(5,3) = &
|
|
(a(1,1)*(a(2,2)*(a(4,3)*a(5,4)-a(4,4)*a(5,3))-a(2,3)*(a(4,2)*a(5,4)-a(4,4)*a(5,2))+a(2,4)*(a(4,2)*a(5,3)-a(4,3)*a(5,2)))-a(1,2)* &
|
|
(a(2,1)*(a(4,3)*a(5,4)-a(4,4)*a(5,3))-a(2,3)*(a(4,1)*a(5,4)-a(4,4)*a(5,1))+a(2,4)*(a(4,1)*a(5,3)-a(4,3)*a(5,1)))+a(1,3)* &
|
|
(a(2,1)*(a(4,2)*a(5,4)-a(4,4)*a(5,2))-a(2,2)*(a(4,1)*a(5,4)-a(4,4)*a(5,1))+a(2,4)*(a(4,1)*a(5,2)-a(4,2)*a(5,1)))-a(1,4)* &
|
|
(a(2,1)*(a(4,2)*a(5,3)-a(4,3)*a(5,2))-a(2,2)*(a(4,1)*a(5,3)-a(4,3)*a(5,1))+a(2,3)*(a(4,1)*a(5,2)-a(4,2)*a(5,1))))
|
|
|
|
b(1,4) = &
|
|
(-a(1,2)*(a(2,3)*(a(3,4)*a(5,5)-a(3,5)*a(5,4))-a(2,4)*(a(3,3)*a(5,5)-a(3,5)*a(5,3))+a(2,5)*(a(3,3)*a(5,4)-a(3,4)*a(5,3)))+a(1,3)* &
|
|
(a(2,2)*(a(3,4)*a(5,5)-a(3,5)*a(5,4))-a(2,4)*(a(3,2)*a(5,5)-a(3,5)*a(5,2))+a(2,5)*(a(3,2)*a(5,4)-a(3,4)*a(5,2)))-a(1,4)* &
|
|
(a(2,2)*(a(3,3)*a(5,5)-a(3,5)*a(5,3))-a(2,3)*(a(3,2)*a(5,5)-a(3,5)*a(5,2))+a(2,5)*(a(3,2)*a(5,3)-a(3,3)*a(5,2)))+a(1,5)* &
|
|
(a(2,2)*(a(3,3)*a(5,4)-a(3,4)*a(5,3))-a(2,3)*(a(3,2)*a(5,4)-a(3,4)*a(5,2))+a(2,4)*(a(3,2)*a(5,3)-a(3,3)*a(5,2))))
|
|
b(2,4) = &
|
|
(a(1,1)*(a(2,3)*(a(3,4)*a(5,5)-a(3,5)*a(5,4))-a(2,4)*(a(3,3)*a(5,5)-a(3,5)*a(5,3))+a(2,5)*(a(3,3)*a(5,4)-a(3,4)*a(5,3)))-a(1,3)* &
|
|
(a(2,1)*(a(3,4)*a(5,5)-a(3,5)*a(5,4))-a(2,4)*(a(3,1)*a(5,5)-a(3,5)*a(5,1))+a(2,5)*(a(3,1)*a(5,4)-a(3,4)*a(5,1)))+a(1,4)* &
|
|
(a(2,1)*(a(3,3)*a(5,5)-a(3,5)*a(5,3))-a(2,3)*(a(3,1)*a(5,5)-a(3,5)*a(5,1))+a(2,5)*(a(3,1)*a(5,3)-a(3,3)*a(5,1)))-a(1,5)* &
|
|
(a(2,1)*(a(3,3)*a(5,4)-a(3,4)*a(5,3))-a(2,3)*(a(3,1)*a(5,4)-a(3,4)*a(5,1))+a(2,4)*(a(3,1)*a(5,3)-a(3,3)*a(5,1))))
|
|
b(3,4) = &
|
|
(-a(1,1)*(a(2,2)*(a(3,4)*a(5,5)-a(3,5)*a(5,4))-a(2,4)*(a(3,2)*a(5,5)-a(3,5)*a(5,2))+a(2,5)*(a(3,2)*a(5,4)-a(3,4)*a(5,2)))+a(1,2)* &
|
|
(a(2,1)*(a(3,4)*a(5,5)-a(3,5)*a(5,4))-a(2,4)*(a(3,1)*a(5,5)-a(3,5)*a(5,1))+a(2,5)*(a(3,1)*a(5,4)-a(3,4)*a(5,1)))-a(1,4)* &
|
|
(a(2,1)*(a(3,2)*a(5,5)-a(3,5)*a(5,2))-a(2,2)*(a(3,1)*a(5,5)-a(3,5)*a(5,1))+a(2,5)*(a(3,1)*a(5,2)-a(3,2)*a(5,1)))+a(1,5)* &
|
|
(a(2,1)*(a(3,2)*a(5,4)-a(3,4)*a(5,2))-a(2,2)*(a(3,1)*a(5,4)-a(3,4)*a(5,1))+a(2,4)*(a(3,1)*a(5,2)-a(3,2)*a(5,1))))
|
|
b(4,4) = &
|
|
(a(1,1)*(a(2,2)*(a(3,3)*a(5,5)-a(3,5)*a(5,3))-a(2,3)*(a(3,2)*a(5,5)-a(3,5)*a(5,2))+a(2,5)*(a(3,2)*a(5,3)-a(3,3)*a(5,2)))-a(1,2)* &
|
|
(a(2,1)*(a(3,3)*a(5,5)-a(3,5)*a(5,3))-a(2,3)*(a(3,1)*a(5,5)-a(3,5)*a(5,1))+a(2,5)*(a(3,1)*a(5,3)-a(3,3)*a(5,1)))+a(1,3)* &
|
|
(a(2,1)*(a(3,2)*a(5,5)-a(3,5)*a(5,2))-a(2,2)*(a(3,1)*a(5,5)-a(3,5)*a(5,1))+a(2,5)*(a(3,1)*a(5,2)-a(3,2)*a(5,1)))-a(1,5)* &
|
|
(a(2,1)*(a(3,2)*a(5,3)-a(3,3)*a(5,2))-a(2,2)*(a(3,1)*a(5,3)-a(3,3)*a(5,1))+a(2,3)*(a(3,1)*a(5,2)-a(3,2)*a(5,1))))
|
|
b(5,4) = &
|
|
(-a(1,1)*(a(2,2)*(a(3,3)*a(5,4)-a(3,4)*a(5,3))-a(2,3)*(a(3,2)*a(5,4)-a(3,4)*a(5,2))+a(2,4)*(a(3,2)*a(5,3)-a(3,3)*a(5,2)))+a(1,2)* &
|
|
(a(2,1)*(a(3,3)*a(5,4)-a(3,4)*a(5,3))-a(2,3)*(a(3,1)*a(5,4)-a(3,4)*a(5,1))+a(2,4)*(a(3,1)*a(5,3)-a(3,3)*a(5,1)))-a(1,3)* &
|
|
(a(2,1)*(a(3,2)*a(5,4)-a(3,4)*a(5,2))-a(2,2)*(a(3,1)*a(5,4)-a(3,4)*a(5,1))+a(2,4)*(a(3,1)*a(5,2)-a(3,2)*a(5,1)))+a(1,4)* &
|
|
(a(2,1)*(a(3,2)*a(5,3)-a(3,3)*a(5,2))-a(2,2)*(a(3,1)*a(5,3)-a(3,3)*a(5,1))+a(2,3)*(a(3,1)*a(5,2)-a(3,2)*a(5,1))))
|
|
|
|
b(1,5) = &
|
|
(a(1,2)*(a(2,3)*(a(3,4)*a(4,5)-a(3,5)*a(4,4))-a(2,4)*(a(3,3)*a(4,5)-a(3,5)*a(4,3))+a(2,5)*(a(3,3)*a(4,4)-a(3,4)*a(4,3)))-a(1,3)* &
|
|
(a(2,2)*(a(3,4)*a(4,5)-a(3,5)*a(4,4))-a(2,4)*(a(3,2)*a(4,5)-a(3,5)*a(4,2))+a(2,5)*(a(3,2)*a(4,4)-a(3,4)*a(4,2)))+a(1,4)* &
|
|
(a(2,2)*(a(3,3)*a(4,5)-a(3,5)*a(4,3))-a(2,3)*(a(3,2)*a(4,5)-a(3,5)*a(4,2))+a(2,5)*(a(3,2)*a(4,3)-a(3,3)*a(4,2)))-a(1,5)* &
|
|
(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))))
|
|
b(2,5) = &
|
|
(-a(1,1)*(a(2,3)*(a(3,4)*a(4,5)-a(3,5)*a(4,4))-a(2,4)*(a(3,3)*a(4,5)-a(3,5)*a(4,3))+a(2,5)*(a(3,3)*a(4,4)-a(3,4)*a(4,3)))+a(1,3)* &
|
|
(a(2,1)*(a(3,4)*a(4,5)-a(3,5)*a(4,4))-a(2,4)*(a(3,1)*a(4,5)-a(3,5)*a(4,1))+a(2,5)*(a(3,1)*a(4,4)-a(3,4)*a(4,1)))-a(1,4)* &
|
|
(a(2,1)*(a(3,3)*a(4,5)-a(3,5)*a(4,3))-a(2,3)*(a(3,1)*a(4,5)-a(3,5)*a(4,1))+a(2,5)*(a(3,1)*a(4,3)-a(3,3)*a(4,1)))+a(1,5)* &
|
|
(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,5) = &
|
|
(a(1,1)*(a(2,2)*(a(3,4)*a(4,5)-a(3,5)*a(4,4))-a(2,4)*(a(3,2)*a(4,5)-a(3,5)*a(4,2))+a(2,5)*(a(3,2)*a(4,4)-a(3,4)*a(4,2)))-a(1,2)* &
|
|
(a(2,1)*(a(3,4)*a(4,5)-a(3,5)*a(4,4))-a(2,4)*(a(3,1)*a(4,5)-a(3,5)*a(4,1))+a(2,5)*(a(3,1)*a(4,4)-a(3,4)*a(4,1)))+a(1,4)* &
|
|
(a(2,1)*(a(3,2)*a(4,5)-a(3,5)*a(4,2))-a(2,2)*(a(3,1)*a(4,5)-a(3,5)*a(4,1))+a(2,5)*(a(3,1)*a(4,2)-a(3,2)*a(4,1)))-a(1,5)* &
|
|
(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,5) = &
|
|
(-a(1,1)*(a(2,2)*(a(3,3)*a(4,5)-a(3,5)*a(4,3))-a(2,3)*(a(3,2)*a(4,5)-a(3,5)*a(4,2))+a(2,5)*(a(3,2)*a(4,3)-a(3,3)*a(4,2)))+a(1,2)* &
|
|
(a(2,1)*(a(3,3)*a(4,5)-a(3,5)*a(4,3))-a(2,3)*(a(3,1)*a(4,5)-a(3,5)*a(4,1))+a(2,5)*(a(3,1)*a(4,3)-a(3,3)*a(4,1)))-a(1,3)* &
|
|
(a(2,1)*(a(3,2)*a(4,5)-a(3,5)*a(4,2))-a(2,2)*(a(3,1)*a(4,5)-a(3,5)*a(4,1))+a(2,5)*(a(3,1)*a(4,2)-a(3,2)*a(4,1)))+a(1,5)* &
|
|
(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(5,5) = &
|
|
(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))))
|
|
|
|
end
|
|
#+end_src
|
|
|
|
#+CALL: generate_c_interface(table=qmckl_adjugate_args,rettyp="qmckl_exit_code",fname="qmckl_adjugate")
|
|
|
|
#+RESULTS:
|
|
#+begin_src f90 :tangle (eval f) :comments org :exports none
|
|
integer(c_int32_t) function qmckl_adjugate &
|
|
(context, n, A, lda, B, ldb, 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 :: n
|
|
real (c_double ) , intent(in) :: A(lda,*)
|
|
integer (c_int64_t) , intent(in) , value :: lda
|
|
real (c_double ) , intent(out) :: B(ldb,*)
|
|
integer (c_int64_t) , intent(in) , value :: ldb
|
|
real (c_double ) , intent(inout) :: det_l
|
|
|
|
integer(c_int32_t), external :: qmckl_adjugate_f
|
|
info = qmckl_adjugate_f &
|
|
(context, n, A, lda, B, ldb, det_l)
|
|
|
|
end function qmckl_adjugate
|
|
#+end_src
|
|
|
|
#+CALL: generate_f_interface(table=qmckl_adjugate_args,rettyp="qmckl_exit_code",fname="qmckl_adjugate")
|
|
|
|
#+RESULTS:
|
|
#+begin_src f90 :tangle (eval fh_func) :comments org :exports none
|
|
interface
|
|
integer(c_int32_t) function qmckl_adjugate &
|
|
(context, n, A, lda, B, ldb, 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 :: n
|
|
real (c_double ) , intent(in) :: A(lda,*)
|
|
integer (c_int64_t) , intent(in) , value :: lda
|
|
real (c_double ) , intent(out) :: B(ldb,*)
|
|
integer (c_int64_t) , intent(in) , value :: ldb
|
|
real (c_double ) , intent(inout) :: det_l
|
|
|
|
end function qmckl_adjugate
|
|
end interface
|
|
#+end_src
|
|
|
|
|
|
#+begin_src f90 :tangle (eval f)
|
|
subroutine adjugate_general(context, na, A, LDA, det_l)
|
|
use qmckl
|
|
implicit none
|
|
integer(qmckl_context) intent(in) :: context
|
|
double precision, intent(in) :: A (LDA,na)
|
|
integer*8, intent(in) :: LDA
|
|
double precision, intent(out) :: B (LDB,na)
|
|
integer*8, intent(in) :: LDB
|
|
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
|
|
|
|
We first copy the array ~A~ into array ~B~.
|
|
|
|
#+begin_src f90 :tangle (eval f)
|
|
B(1:na,1:na) = A(1:na,1:na)
|
|
#+end_src
|
|
|
|
Then, we compute the LU factorization $LU=A$
|
|
using the ~dgetrf~ routine.
|
|
|
|
#+begin_src f90 :tangle (eval f)
|
|
call dgetrf(na, na, B, LDB, 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*B(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, B, LDB, 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)
|
|
B(:,:) = B(:,:)*det_l
|
|
|
|
end subroutine adjugate_general
|
|
#+end_src
|
|
|
|
*** Test :noexport:
|
|
|
|
#+begin_src f90 :tangle (eval f_test)
|
|
integer(qmckl_exit_code) function test_qmckl_adjugate(context) bind(C)
|
|
use qmckl
|
|
implicit none
|
|
integer(qmckl_context), intent(in), value :: context
|
|
|
|
double precision, allocatable :: A(:,:), B(:,:)
|
|
integer*8 :: m, n, k, LDA, LDB
|
|
integer*8 :: i,j,l
|
|
double precision :: x, det_l, det_l_ref
|
|
|
|
LDA = 6
|
|
LDB = 6
|
|
|
|
allocate( A(LDA,6), B(LDB,6))
|
|
|
|
A = 0.1d0
|
|
A(1,1) = 1.0d0;
|
|
A(2,2) = 2.0d0;
|
|
A(3,3) = 3.0d0;
|
|
A(4,4) = 4.0d0;
|
|
A(5,5) = 5.0d0;
|
|
A(6,6) = 6.0d0;
|
|
|
|
test_qmckl_adjugate = QMCKL_SUCCESS
|
|
|
|
#+end_src
|
|
|
|
#+begin_src python :results output :output drawer
|
|
import numpy as np
|
|
import numpy.linalg as la
|
|
N=6
|
|
|
|
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" test_qmckl_adjugate = qmckl_adjugate(context, {i+1}_8, A, LDA, B, LDB, 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_example
|
|
,#+begin_src f90 :tangle (eval f_test)
|
|
! N = 1
|
|
test_qmckl_adjugate = qmckl_adjugate(context, 1_8, A, LDA, B, LDB, 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 = x + (B(1,1) - (1.0) )**2
|
|
if (dabs(x / det_l) > 1.d-12) then
|
|
print *, 'N = 1: x = ', x
|
|
test_qmckl_adjugate = 1
|
|
return
|
|
end if
|
|
|
|
! N = 2
|
|
test_qmckl_adjugate = qmckl_adjugate(context, 2_8, A, LDA, B, LDB, 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
|
|
|
|
x = 0.d0
|
|
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
|
|
test_qmckl_adjugate = qmckl_adjugate(context, 3_8, A, LDA, B, LDB, 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
|
|
test_qmckl_adjugate = qmckl_adjugate(context, 4_8, A, LDA, B, LDB, 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
|
|
test_qmckl_adjugate = qmckl_adjugate(context, 5_8, A, LDA, B, LDB, 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
|
|
test_qmckl_adjugate = qmckl_adjugate(context, 6_8, A, LDA, B, LDB, 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
|
|
#+end_example
|
|
|
|
#+begin_src f90 :tangle (eval f_test)
|
|
|
|
deallocate(A,B)
|
|
|
|
end function test_qmckl_adjugate
|
|
#+end_src
|
|
|
|
#+begin_src c :comments link :tangle (eval c_test)
|
|
qmckl_exit_code test_qmckl_adjugate(qmckl_context context);
|
|
assert(QMCKL_SUCCESS == test_qmckl_adjugate(context));
|
|
#+end_src
|
|
|
|
* End of files :noexport:
|
|
|
|
#+begin_src c :comments link :tangle (eval c_test)
|
|
assert (qmckl_context_destroy(context) == QMCKL_SUCCESS);
|
|
return 0;
|
|
}
|
|
|
|
#+end_src
|
|
|
|
|
|
# -*- mode: org -*-
|
|
# vim: syntax=c
|