TREX/qmckl
1
0
Fork 0
mirror of https://github.com/TREX-CoE/qmckl.git synced 2024-07-02 01:16:18 +02:00
Code Issues Projects Releases Wiki Activity
18f9f96bcc
qmckl/org/qmckl_blas.org

1717 lines
60 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 :tangle (eval h_private_type)
#ifndef QMCKL_BLAS_HPT
#define QMCKL_BLAS_HPT
#+end_src
#+begin_src c :tangle (eval h_private_func)
#ifndef QMCKL_BLAS_HPF
#define QMCKL_BLAS_HPF
#include "qmckl_blas_private_type.h"
#+end_src
#+begin_src c :tangle (eval c)
#ifdef HAVE_CONFIG_H
#include "config.h"
#endif
#ifdef HAVE_STDINT_H
#include <stdint.h>
#elif HAVE_INTTYPES_H
#include <inttypes.h>
#endif
#include <stdlib.h>
#include <string.h>
#include <stdbool.h>
#include <assert.h>
#include <math.h>
#include "qmckl.h"
#include "qmckl_context_private_type.h"
#include "qmckl_memory_private_type.h"
#include "qmckl_memory_private_func.h"
#include "qmckl_blas_private_type.h"
#include "qmckl_blas_private_func.h"
#+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
#include "qmckl_blas_private_type.h"
#include "qmckl_blas_private_func.h"
int main() {
qmckl_context context;
context = qmckl_context_create();
#+end_src
* Data types
** Vector
| Variable | Type | Description |
|----------+-----------+-------------------------|
| ~size~ | ~int64_t~ | Dimension of the vector |
| ~data~ | ~double*~ | Elements |
#+begin_src c :comments org :tangle (eval h_private_type)
typedef struct qmckl_vector {
int64_t size;
double* data;
} qmckl_vector;
#+end_src
#+begin_src c :comments org :tangle (eval h_private_func)
qmckl_vector
qmckl_vector_alloc( qmckl_context context,
const int64_t size);
#+end_src
Allocates a new vector. If the allocation failed the size is zero.
#+begin_src c :comments org :tangle (eval c)
qmckl_vector
qmckl_vector_alloc( qmckl_context context,
const int64_t size)
{
/* Should always be true by contruction */
assert (size > (int64_t) 0);
qmckl_vector result;
result.size = size;
qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero;
mem_info.size = size * sizeof(double);
result.data = (double*) qmckl_malloc (context, mem_info);
if (result.data == NULL) {
result.size = (int64_t) 0;
}
return result;
}
#+end_src
#+begin_src c :comments org :tangle (eval h_private_func)
qmckl_exit_code
qmckl_vector_free( qmckl_context context,
qmckl_vector vector);
#+end_src
#+begin_src c :comments org :tangle (eval c)
qmckl_exit_code
qmckl_vector_free( qmckl_context context,
qmckl_vector vector)
{
/* Always true */
assert (vector.data != NULL);
qmckl_exit_code rc;
rc = qmckl_free(context, vector.data);
if (rc != QMCKL_SUCCESS) {
return rc;
}
vector.size = (int64_t) 0;
vector.data = NULL;
return QMCKL_SUCCESS;
}
#+end_src
** Matrix
| Variable | Type | Description |
|----------+--------------+-----------------------------|
| ~size~ | ~int64_t[2]~ | Dimension of each component |
| ~data~ | ~double*~ | Elements |
The dimensions use Fortran ordering: two elements differing by one
in the first dimension are consecutive in memory.
#+begin_src c :comments org :tangle (eval h_private_type)
typedef struct qmckl_matrix {
int64_t size[2];
double* data;
} qmckl_matrix;
#+end_src
#+begin_src c :comments org :tangle (eval h_private_func)
qmckl_matrix
qmckl_matrix_alloc( qmckl_context context,
const int64_t size1,
const int64_t size2);
#+end_src
Allocates a new matrix. If the allocation failed the sizes are zero.
#+begin_src c :comments org :tangle (eval c)
qmckl_matrix
qmckl_matrix_alloc( qmckl_context context,
const int64_t size1,
const int64_t size2)
{
/* Should always be true by contruction */
assert (size1 * size2 > (int64_t) 0);
qmckl_matrix result;
result.size[0] = size1;
result.size[1] = size2;
qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero;
mem_info.size = size1 * size2 * sizeof(double);
result.data = (double*) qmckl_malloc (context, mem_info);
if (result.data == NULL) {
result.size[0] = (int64_t) 0;
result.size[1] = (int64_t) 0;
}
return result;
}
#+end_src
#+begin_src c :comments org :tangle (eval h_private_func)
qmckl_exit_code
qmckl_matrix_free( qmckl_context context,
qmckl_matrix matrix);
#+end_src
#+begin_src c :comments org :tangle (eval c)
qmckl_exit_code
qmckl_matrix_free( qmckl_context context,
qmckl_matrix matrix)
{
/* Always true */
assert (matrix.data != NULL);
qmckl_exit_code rc;
rc = qmckl_free(context, matrix.data);
if (rc != QMCKL_SUCCESS) {
return rc;
}
matrix.data = NULL;
matrix.size[0] = (int64_t) 0;
matrix.size[1] = (int64_t) 0;
return QMCKL_SUCCESS;
}
#+end_src
** Tensor
| Variable | Type | Description |
|----------+-----------------------------------+-----------------------------|
| ~order~ | ~int64_t~ | Order of the tensor |
| ~size~ | ~int64_t[QMCKL_TENSOR_ORDER_MAX]~ | Dimension of each component |
| ~data~ | ~double*~ | Elements |
The dimensions use Fortran ordering: two elements differing by one
in the first dimension are consecutive in memory.
#+begin_src c :comments org :tangle (eval h_private_type)
#define QMCKL_TENSOR_ORDER_MAX 16
typedef struct qmckl_tensor {
int64_t order;
int64_t size[QMCKL_TENSOR_ORDER_MAX];
double* data;
} qmckl_tensor;
#+end_src
#+begin_src c :comments org :tangle (eval h_private_func)
qmckl_tensor
qmckl_tensor_alloc( qmckl_context context,
const int64_t order,
const int64_t* size);
#+end_src
Allocates memory for a tensor. If the allocation failed, the size
is zero.
#+begin_src c :comments org :tangle (eval c)
qmckl_tensor
qmckl_tensor_alloc( qmckl_context context,
const int64_t order,
const int64_t* size)
{
/* Should always be true by contruction */
assert (order > 0);
assert (order <= QMCKL_TENSOR_ORDER_MAX);
assert (size != NULL);
qmckl_tensor result;
result.order = order;
int64_t prod_size = (int64_t) 1;
for (int64_t i=0 ; i<order ; ++i) {
assert (size[i] > (int64_t) 0);
result.size[i] = size[i];
prod_size *= size[i];
}
qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero;
mem_info.size = prod_size * sizeof(double);
result.data = (double*) qmckl_malloc (context, mem_info);
if (result.data == NULL) {
memset(&result, 0, sizeof(qmckl_tensor));
}
return result;
}
#+end_src
#+begin_src c :comments org :tangle (eval h_private_func)
qmckl_exit_code
qmckl_tensor_free( qmckl_context context,
qmckl_tensor tensor);
#+end_src
#+begin_src c :comments org :tangle (eval c)
qmckl_exit_code
qmckl_tensor_free( qmckl_context context,
qmckl_tensor tensor)
{
/* Always true */
assert (tensor.data != NULL);
qmckl_exit_code rc;
rc = qmckl_free(context, tensor.data);
if (rc != QMCKL_SUCCESS) {
return rc;
}
memset(&tensor, 0, sizeof(qmckl_tensor));
return QMCKL_SUCCESS;
}
#+end_src
** Reshaping
Reshaping occurs in-place and the pointer to the data is copied.
*** Vector -> Matrix
#+begin_src c :comments org :tangle (eval h_private_func)
qmckl_matrix
qmckl_matrix_of_vector(const qmckl_vector vector,
const int64_t size1,
const int64_t size2);
#+end_src
Reshapes a vector into a matrix.
#+begin_src c :comments org :tangle (eval c)
qmckl_matrix
qmckl_matrix_of_vector(const qmckl_vector vector,
const int64_t size1,
const int64_t size2)
{
/* Always true */
assert (size1 * size2 == vector.size);
qmckl_matrix result;
result.size[0] = size1;
result.size[1] = size2;
result.data = vector.data;
return result;
}
#+end_src
*** Vector -> Tensor
#+begin_src c :comments org :tangle (eval h_private_func)
qmckl_tensor
qmckl_tensor_of_vector(const qmckl_vector vector,
const int64_t order,
const int64_t* size);
#+end_src
Reshapes a vector into a tensor.
#+begin_src c :comments org :tangle (eval c)
qmckl_tensor
qmckl_tensor_of_vector(const qmckl_vector vector,
const int64_t order,
const int64_t* size)
{
qmckl_tensor result;
int64_t prod_size = 1;
for (int64_t i=0 ; i<order ; ++i) {
result.size[i] = size[i];
prod_size *= size[i];
}
assert (prod_size == vector.size);
result.data = vector.data;
return result;
}
#+end_src
*** Matrix -> Vector
#+begin_src c :comments org :tangle (eval h_private_func)
qmckl_vector
qmckl_vector_of_matrix(const qmckl_matrix matrix,
const int64_t size);
#+end_src
Reshapes a matrix into a vector.
#+begin_src c :comments org :tangle (eval c)
qmckl_vector
qmckl_vector_of_matrix(const qmckl_matrix matrix,
const int64_t size)
{
/* Always true */
assert (matrix.size[0] * matrix.size[1] == size);
qmckl_vector result;
result.size = size;
result.data = matrix.data;
return result;
}
#+end_src
*** Matrix -> Tensor
#+begin_src c :comments org :tangle (eval h_private_func)
qmckl_tensor
qmckl_tensor_of_matrix(const qmckl_matrix matrix,
const int64_t order,
const int64_t* size);
#+end_src
Reshapes a matrix into a tensor.
#+begin_src c :comments org :tangle (eval c)
qmckl_tensor
qmckl_tensor_of_matrix(const qmckl_matrix matrix,
const int64_t order,
const int64_t* size)
{
qmckl_tensor result;
int64_t prod_size = 1;
for (int64_t i=0 ; i<order ; ++i) {
result.size[i] = size[i];
prod_size *= size[i];
}
assert (prod_size == matrix.size[0] * matrix.size[1]);
result.data = matrix.data;
return result;
}
#+end_src
*** Tensor -> Vector
#+begin_src c :comments org :tangle (eval h_private_func)
qmckl_vector
qmckl_vector_of_tensor(const qmckl_tensor tensor,
const int64_t size);
#+end_src
Reshapes a tensor into a vector.
#+begin_src c :comments org :tangle (eval c)
qmckl_vector
qmckl_vector_of_tensor(const qmckl_tensor tensor,
const int64_t size)
{
/* Always true */
int64_t prod_size = (int64_t) 1;
for (int64_t i=0 ; i<tensor.order ; i++) {
prod_size *= tensor.size[i];
}
assert (prod_size == size);
qmckl_vector result;
result.size = size;
result.data = tensor.data;
return result;
}
#+end_src
*** Tensor -> Matrix
#+begin_src c :comments org :tangle (eval h_private_func)
qmckl_matrix
qmckl_matrix_of_tensor(const qmckl_tensor tensor,
const int64_t size1,
const int64_t size2);
#+end_src
Reshapes a tensor into a vector.
#+begin_src c :comments org :tangle (eval c)
qmckl_matrix
qmckl_matrix_of_tensor(const qmckl_tensor tensor,
const int64_t size1,
const int64_t size2)
{
/* Always true */
int64_t prod_size = (int64_t) 1;
for (int64_t i=0 ; i<tensor.order ; i++) {
prod_size *= tensor.size[i];
}
assert (prod_size == size1 * size2);
qmckl_matrix result;
result.size[0] = size1;
result.size[1] = size2;
result.data = tensor.data;
return result;
}
#+end_src
** Access macros
#+begin_src c :comments org :tangle (eval h_private_func)
#define qmckl_vec(v, i) v.data[i]
#define qmckl_mat(m, i, j) m.data[(i) + (j)*m.size[0]]
#define qmckl_ten3(t, i, j, k) t.data[(i) + m.size[0]*((j) + size[1]*(k))]
#define qmckl_ten4(t, i, j, k, l) t.data[(i) + m.size[0]*((j) + size[1]*((k) + size[2]*(l)))]
#define qmckl_ten5(t, i, j, k, l, m) t.data[(i) + m.size[0]*((j) + size[1]*((k) + size[2]*((l) + size[3]*(m))))]
#+end_src
** Tests
#+begin_src c :comments link :tangle (eval c_test)
{
int64_t m = 3;
int64_t n = 4;
int64_t p = m*n;
qmckl_vector vec = qmckl_vector_alloc(context, p);
for (int64_t i=0 ; i<p ; ++i)
qmckl_vec(vec, i) = (double) i;
for (int64_t i=0 ; i<p ; ++i)
assert( vec.data[i] == (double) i );
qmckl_matrix mat = qmckl_matrix_of_vector(vec, m, n);
assert (mat.size[0] == m);
assert (mat.size[1] == n);
assert (mat.data == vec.data);
for (int64_t j=0 ; j<n ; ++j)
for (int64_t i=0 ; i<m ; ++i)
assert ( qmckl_mat(mat, i, j) == qmckl_vec(vec, i+j*m)) ;
qmckl_vector vec2 = qmckl_vector_of_matrix(mat, p);
assert (vec2.size == p);
assert (vec2.data == vec.data);
for (int64_t i=0 ; i<p ; ++i)
assert ( qmckl_vec(vec2, i) == qmckl_vec(vec, i) ) ;
qmckl_vector_free(context, vec);
}
#+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~ | ~char~ | in | 'T' is transposed |
| ~TransB~ | ~char~ | in | 'T' is transposed |
| ~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 char TransA,
const char 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
character , 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 <= 0) then
info = QMCKL_INVALID_ARG_9
return
endif
if (LDB <= 0) then
info = QMCKL_INVALID_ARG_11
return
endif
if (LDC <= 0) then
info = QMCKL_INVALID_ARG_14
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
character , intent(in) , value :: TransA
character , 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
character , intent(in) , value :: TransA
character , 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
character :: TransA, TransB
double precision :: x, alpha, beta
TransA = 'N'
TransB = 'N'
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_8,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_8,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,C,4_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(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_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
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
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
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, B, LDB, 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(8) :: 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_8
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_8) /= 0_8) then
det_l = -det_l
endif
#+end_src
Then, the inverse of $A$ is computed using ~dgetri~:
#+begin_src f90 :tangle (eval f)
lwork = SIZE(work)
call dgetri(na, B, LDB, ipiv, work, lwork, inf )
#+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_8
LDB = 6_8
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 :tangle (eval h_private_type)
#endif
#+end_src
#+begin_src c :tangle (eval h_private_func)
#endif
#+end_src
#+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