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
c9450bb1a7
qmckl/org/qmckl_blas.org
Anthony Scemama 5ecb1d6326 Faster AOs
2022-03-21 18:32:39 +01:00

2312 lines
78 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_blas_private_type.h"
#include "qmckl_memory_private_func.h"
#include "qmckl_blas_private_func.h"
#+end_src
#+begin_src c :comments link :tangle (eval c_test) :noweb yes
#include "qmckl.h"
#include <stdio.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
* -
:PROPERTIES:
:UNNUMBERED: t
:END:
Basic linear algebra data types and operations are described in this file.
The data types are private, so that HPC implementations can use
whatever data structures they prefer.
These data types are expected to be used internally in QMCkl. They
are not intended to be passed to external codes.
* 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) :exports none
typedef struct qmckl_vector {
int64_t size;
double* restrict 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) :exports none
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) :exports none
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) :exports none
typedef struct qmckl_matrix {
int64_t size[2];
double* restrict 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) :exports none
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) :exports none
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) :exports none
#define QMCKL_TENSOR_ORDER_MAX 16
typedef struct qmckl_tensor {
int64_t order;
int64_t size[QMCKL_TENSOR_ORDER_MAX];
double* restrict 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) :exports none
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) :exports none
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) :exports none
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) :exports none
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);
#+end_src
Reshapes a matrix into a vector.
#+begin_src c :comments org :tangle (eval c) :exports none
qmckl_vector
qmckl_vector_of_matrix(const qmckl_matrix matrix)
{
qmckl_vector result;
result.size = matrix.size[0] * matrix.size[1];
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) :exports none
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);
#+end_src
Reshapes a tensor into a vector.
#+begin_src c :comments org :tangle (eval c) :exports none
qmckl_vector
qmckl_vector_of_tensor(const qmckl_tensor tensor)
{
int64_t prod_size = (int64_t) tensor.size[0];
for (int64_t i=1 ; i<tensor.order ; i++) {
prod_size *= tensor.size[i];
}
qmckl_vector result;
result.size = prod_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) :exports none
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
Macros are provided to ease the access to vectors, matrices and
tensors. Matrices use column-major ordering, as in Fortran.
#+begin_src c :tangle no
double qmckl_vec (qmckl_vector v, int64_t i); // v[i]
double qmckl_mat (qmckl_matrix m, int64_t i, int64_t j) // m[j][i]
double qmckl_ten3(qmckl_tensor t, int64_t i, int64_t j, int64_t k); // t[k][j][i]
double qmckl_ten4(qmckl_tensor t, int64_t i, int64_t j, int64_t k, int64_t l);
double qmckl_ten5(qmckl_tensor t, int64_t i, int64_t j, int64_t k, int64_t l, int64_t m);
...
#+end_src
#+begin_src c :comments org :tangle (eval h_private_func) :exports none
#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) + t.size[0]*((j) + t.size[1]*(k))]
#define qmckl_ten4(t, i, j, k, l) t.data[(i) + t.size[0]*((j) + t.size[1]*((k) + t.size[2]*(l)))]
#define qmckl_ten5(t, i, j, k, l, m) t.data[(i) + t.size[0]*((j) + t.size[1]*((k) + t.size[2]*((l) + t.size[3]*(m))))]
#+end_src
For example:
** Set all elements
*** Vector
#+begin_src c :comments org :tangle (eval h_private_func)
qmckl_vector
qmckl_vector_set(qmckl_vector vector, double value);
#+end_src
#+begin_src c :comments org :tangle (eval c) :exports none
qmckl_vector
qmckl_vector_set(qmckl_vector vector, double value)
{
for (int64_t i=0 ; i<vector.size ; ++i) {
qmckl_vec(vector, i) = value;
}
return vector;
}
#+end_src
*** Matrix
#+begin_src c :comments org :tangle (eval h_private_func)
qmckl_matrix
qmckl_matrix_set(qmckl_matrix matrix, double value);
#+end_src
#+begin_src c :comments org :tangle (eval c) :exports none
qmckl_matrix
qmckl_matrix_set(qmckl_matrix matrix, double value)
{
qmckl_vector vector = qmckl_vector_of_matrix(matrix);
for (int64_t i=0 ; i<vector.size ; ++i) {
qmckl_vec(vector, i) = value;
}
return qmckl_matrix_of_vector(vector, matrix.size[0], matrix.size[1]);
}
#+end_src
*** Tensor
#+begin_src c :comments org :tangle (eval h_private_func)
qmckl_tensor
qmckl_tensor_set(qmckl_tensor tensor, double value);
#+end_src
#+begin_src c :comments org :tangle (eval c) :exports none
qmckl_tensor
qmckl_tensor_set(qmckl_tensor tensor, double value)
{
qmckl_vector vector = qmckl_vector_of_tensor(tensor);
for (int64_t i=0 ; i<vector.size ; ++i) {
qmckl_vec(vector, i) = value;
}
return qmckl_tensor_of_vector(vector, tensor.order, tensor.size);
}
#+end_src
** Copy to/from to ~double*~
#+begin_src c :comments org :tangle (eval h_private_func)
qmckl_exit_code
qmckl_double_of_vector(const qmckl_context context,
const qmckl_vector vector,
double* const target,
const int64_t size_max);
#+end_src
Converts a vector to a ~double*~.
#+begin_src c :comments org :tangle (eval c) :exports none
qmckl_exit_code
qmckl_double_of_vector(const qmckl_context context,
const qmckl_vector vector,
double* const target,
const int64_t size_max)
{
/* Always true by construction */
assert (qmckl_context_check(context) != QMCKL_NULL_CONTEXT);
assert (vector.size > (int64_t) 0);
assert (target != NULL);
assert (size_max > (int64_t) 0);
assert (size_max >= vector.size);
for (int64_t i=0 ; i<vector.size ; ++i) {
target[i] = vector.data[i];
}
return QMCKL_SUCCESS;
}
#+end_src
#+begin_src c :comments org :tangle (eval h_private_func)
qmckl_exit_code
qmckl_double_of_matrix(const qmckl_context context,
const qmckl_matrix matrix,
double* const target,
const int64_t size_max);
#+end_src
Converts a matrix to a ~double*~.
#+begin_src c :comments org :tangle (eval c) :exports none
qmckl_exit_code
qmckl_double_of_matrix(const qmckl_context context,
const qmckl_matrix matrix,
double* const target,
const int64_t size_max)
{
qmckl_vector vector = qmckl_vector_of_matrix(matrix);
return qmckl_double_of_vector(context, vector, target, size_max);
}
#+end_src
#+begin_src c :comments org :tangle (eval h_private_func)
qmckl_exit_code
qmckl_double_of_tensor(const qmckl_context context,
const qmckl_tensor tensor,
double* const target,
const int64_t size_max);
#+end_src
Converts a tensor to a ~double*~.
#+begin_src c :comments org :tangle (eval c) :exports none
qmckl_exit_code
qmckl_double_of_tensor(const qmckl_context context,
const qmckl_tensor tensor,
double* const target,
const int64_t size_max)
{
qmckl_vector vector = qmckl_vector_of_tensor(tensor);
return qmckl_double_of_vector(context, vector, target, size_max);
}
#+end_src
#+begin_src c :comments org :tangle (eval h_private_func)
qmckl_exit_code
qmckl_vector_of_double(const qmckl_context context,
const double* target,
const int64_t size_max,
qmckl_vector* vector);
#+end_src
Converts a ~double*~ to a vector.
#+begin_src c :comments org :tangle (eval c) :exports none
qmckl_exit_code
qmckl_vector_of_double(const qmckl_context context,
const double* target,
const int64_t size_max,
qmckl_vector* vector_out)
{
qmckl_vector vector = *vector_out;
/* Always true by construction */
assert (qmckl_context_check(context) != QMCKL_NULL_CONTEXT);
if (vector.size == 0) {
return qmckl_failwith( context,
QMCKL_INVALID_ARG_4,
"qmckl_double_of_vector",
"Vector not allocated");
}
if (vector.size != size_max) {
return qmckl_failwith( context,
QMCKL_INVALID_ARG_4,
"qmckl_double_of_vector",
"Wrong vector size");
}
for (int64_t i=0 ; i<vector.size ; ++i) {
vector.data[i] = target[i];
}
*vector_out = vector;
return QMCKL_SUCCESS;
}
#+end_src
#+begin_src c :comments org :tangle (eval h_private_func)
qmckl_exit_code
qmckl_matrix_of_double(const qmckl_context context,
const double* target,
const int64_t size_max,
qmckl_matrix* matrix);
#+end_src
Converts a ~double*~ to a matrix.
#+begin_src c :comments org :tangle (eval c) :exports none
qmckl_exit_code
qmckl_matrix_of_double(const qmckl_context context,
const double* target,
const int64_t size_max,
qmckl_matrix* matrix)
{
qmckl_vector vector = qmckl_vector_of_matrix(*matrix);
qmckl_exit_code rc =
qmckl_vector_of_double(context, target, size_max, &vector);
,*matrix = qmckl_matrix_of_vector(vector, matrix->size[0], matrix->size[1]);
return rc;
}
#+end_src
#+begin_src c :comments org :tangle (eval h_private_func)
qmckl_exit_code
qmckl_tensor_of_double(const qmckl_context context,
const double* target,
const int64_t size_max,
qmckl_tensor* tensor);
#+end_src
Converts a ~double*~ to a tensor.
#+begin_src c :comments org :tangle (eval c) :exports none
qmckl_exit_code
qmckl_tensor_of_double(const qmckl_context context,
const double* target,
const int64_t size_max,
qmckl_tensor* tensor)
{
qmckl_vector vector = qmckl_vector_of_tensor(*tensor);
qmckl_exit_code rc =
qmckl_vector_of_double(context, target, size_max, &vector);
*tensor = qmckl_tensor_of_vector(vector, tensor->order, tensor->size);
return rc;
}
#+end_src
** Tests :noexport:
#+begin_src c :comments link :tangle (eval c_test) :exports none
{
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);
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 with a BLAS interface:
\[
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 | \alpha |
| ~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 | \beta |
| ~C~ | ~double[][ldc]~ | out | Array containing the matrix $C$ |
| ~ldc~ | ~int64_t~ | in | Leading dimension of array ~C~ |
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) :exports none
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) :exports none
qmckl_exit_code test_qmckl_dgemm(qmckl_context context);
assert(QMCKL_SUCCESS == test_qmckl_dgemm(context));
#+end_src
** ~qmckl_matmul~
Matrix multiplication using the =qmckl_matrix= data type:
\[
C_{ij} = \beta C_{ij} + \alpha \sum_{k} A_{ik} \cdot B_{kj}
\]
# TODO: Add description about the external library dependence.
#+NAME: qmckl_matmul_args
| Variable | Type | In/Out | Description |
|-----------+-----------------+--------+-------------------|
| ~context~ | ~qmckl_context~ | in | Global state |
| ~TransA~ | ~char~ | in | 'T' is transposed |
| ~TransB~ | ~char~ | in | 'T' is transposed |
| ~alpha~ | ~double~ | in | \alpha |
| ~A~ | ~qmckl_matrix~ | in | Matrix $A$ |
| ~B~ | ~qmckl_matrix~ | in | Matrix $B$ |
| ~beta~ | ~double~ | in | \beta |
| ~C~ | ~qmckl_matrix~ | out | Matrix $C$ |
#+CALL: generate_c_header(table=qmckl_matmul_args,rettyp="qmckl_exit_code",fname="qmckl_matmul")
#+RESULTS:
#+begin_src c :tangle (eval h_private_func) :comments org
qmckl_exit_code
qmckl_matmul (const qmckl_context context,
const char TransA,
const char TransB,
const double alpha,
const qmckl_matrix A,
const qmckl_matrix B,
const double beta,
qmckl_matrix* const C );
#+end_src
#+begin_src c :tangle (eval c) :comments org :exports none
qmckl_exit_code
qmckl_matmul (const qmckl_context context,
const char TransA,
const char TransB,
const double alpha,
const qmckl_matrix A,
const qmckl_matrix B,
const double beta,
qmckl_matrix* const C )
{
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
return QMCKL_INVALID_CONTEXT;
}
qmckl_exit_code rc = QMCKL_SUCCESS;
if (TransA != 'N' && TransA != 'T') {
return qmckl_failwith( context,
QMCKL_INVALID_ARG_2,
"qmckl_matmul",
"TransA should be 'N' or 'T'");
}
if (TransB != 'N' && TransB != 'T') {
return qmckl_failwith( context,
QMCKL_INVALID_ARG_3,
"qmckl_matmul",
"TransB should be 'N' or 'T'");
}
if (A.size[0] < 1) {
return qmckl_failwith( context,
QMCKL_INVALID_ARG_5,
"qmckl_matmul",
"Invalid size for A");
}
if (B.size[0] < 1) {
return qmckl_failwith( context,
QMCKL_INVALID_ARG_6,
"qmckl_matmul",
"Invalid size for B");
}
if (C == NULL) {
return qmckl_failwith( context,
QMCKL_INVALID_ARG_8,
"qmckl_matmul",
"Null pointer");
}
int t = 0;
if (TransA == 'T') t +=1;
if (TransB == 'T') t +=2;
/*
| t | TransA | TransB |
+---+--------+--------+
| 0 | N | N |
| 1 | T | N |
| 2 | N | T |
| 3 | T | T |
,*/
switch (t) {
case 0:
if (A.size[1] != B.size[0]) {
return qmckl_failwith( context,
QMCKL_INVALID_ARG_2,
"qmckl_matmul",
"A and B have incompatible dimensions");
}
C->size[0] = A.size[0];
C->size[1] = B.size[1];
rc = qmckl_dgemm (context, 'N', 'N',
C->size[0], C->size[1], A.size[1],
alpha,
A.data, A.size[0],
B.data, B.size[0],
beta,
C->data, C->size[0]);
break;
case 1:
if (A.size[0] != B.size[0]) {
return qmckl_failwith( context,
QMCKL_INVALID_ARG_2,
"qmckl_matmul",
"A and B have incompatible dimensions");
}
C->size[0] = A.size[1];
C->size[1] = B.size[1];
rc = qmckl_dgemm (context, 'T', 'N',
C->size[0], C->size[1], A.size[0],
alpha,
A.data, A.size[0],
B.data, B.size[0],
beta,
C->data, C->size[0]);
break;
case 2:
if (A.size[1] != B.size[1]) {
return qmckl_failwith( context,
QMCKL_INVALID_ARG_2,
"qmckl_matmul",
"A and B have incompatible dimensions");
}
C->size[0] = A.size[0];
C->size[1] = B.size[0];
rc = qmckl_dgemm (context, 'N', 'T',
C->size[0], C->size[1], A.size[1],
alpha,
A.data, A.size[0],
B.data, B.size[0],
beta,
C->data, C->size[0]);
break;
case 3:
if (A.size[0] != B.size[1]) {
return qmckl_failwith( context,
QMCKL_INVALID_ARG_2,
"qmckl_matmul",
"A and B have incompatible dimensions");
}
C->size[0] = A.size[1];
C->size[1] = B.size[0];
rc = qmckl_dgemm (context, 'T', 'T',
C->size[0], C->size[1], A.size[0],
alpha,
A.data, A.size[0],
B.data, B.size[0],
beta,
C->data, C->size[0]);
break;
}
return rc;
}
#+end_src
*** Test :noexport:
#+begin_src python :exports none :results output
import numpy as np
A = np.array([[ 1., 2., 3., 4. ],
[ 5., 6., 7., 8. ],
[ 9., 10., 11., 12. ]])
B = np.array([[ 1., -2., 3., 4., 5. ],
[ 5., -6., 7., 8., 9. ],
[ 9., 10., 11., 12., 13. ],
[10., 11., 12., 15., 14. ]])
C = 0.5 * A @ B
print(A.T)
print(B.T)
print(C.T)
#+end_src
#+RESULTS:
#+begin_example
[[ 1. 5. 9.]
[ 2. 6. 10.]
[ 3. 7. 11.]
[ 4. 8. 12.]]
[[ 1. 5. 9. 10.]
[-2. -6. 10. 11.]
[ 3. 7. 11. 12.]
[ 4. 8. 12. 15.]
[ 5. 9. 13. 14.]]
[[ 39. 89. 139.]
[ 30. 56. 82.]
[ 49. 115. 181.]
[ 58. 136. 214.]
[ 59. 141. 223.]]
#+end_example
#+begin_src c :comments link :tangle (eval c_test) :exports none
{
double a[12] = { 1., 5., 9.,
2., 6., 10.,
3., 7., 11.,
4., 8., 12. };
double b[20] = { 1., 5., 9., 10.,
-2., -6., 10., 11.,
3., 7., 11., 12.,
4., 8., 12., 15.,
5., 9., 13., 14. };
double c[15] = { 39., 89., 139.,
30., 56., 82.,
49., 115., 181.,
58., 136., 214.,
59., 141., 223. };
double cnew[15];
qmckl_exit_code rc;
qmckl_matrix A = qmckl_matrix_alloc(context, 3, 4);
rc = qmckl_matrix_of_double(context, a, 12, &A);
assert(rc == QMCKL_SUCCESS);
qmckl_matrix B = qmckl_matrix_alloc(context, 4, 5);
rc = qmckl_matrix_of_double(context, b, 20, &B);
assert(rc == QMCKL_SUCCESS);
qmckl_matrix C = qmckl_matrix_alloc(context, 3, 5);
rc = qmckl_matmul(context, 'N', 'N', 0.5, A, B, 0., &C);
assert(rc == QMCKL_SUCCESS);
rc = qmckl_double_of_matrix(context, C, cnew, 15);
assert(rc == QMCKL_SUCCESS);
for (int i=0 ; i<15 ; ++i) {
printf("%f %f\n", cnew[i], c[i]);
assert (c[i] == cnew[i]);
}
}
#+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) :exports none
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
** ~qmckl_transpose~
Transposes a matrix: $A^\dagger_{ji} = A_{ij}$.
| Variable | Type | In/Out | Description |
|-----------+-----------------+--------+-------------------|
| ~context~ | ~qmckl_context~ | in | Global state |
| ~A~ | ~qmckl_matrix~ | in | Input matrix |
| ~At~ | ~qmckl_matrix~ | out | Transposed matrix |
#+begin_src c :tangle (eval h_private_func) :comments org
qmckl_exit_code
qmckl_transpose (qmckl_context context,
const qmckl_matrix A,
qmckl_matrix At );
#+end_src
#+begin_src c :tangle (eval c) :comments org :exports none
qmckl_exit_code
qmckl_transpose (qmckl_context context,
const qmckl_matrix A,
qmckl_matrix At )
{
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
return QMCKL_INVALID_CONTEXT;
}
if (A.size[0] < 1) {
return qmckl_failwith( context,
QMCKL_INVALID_ARG_2,
"qmckl_transpose",
"Invalid size for A");
}
if (At.data == NULL) {
return qmckl_failwith( context,
QMCKL_INVALID_ARG_3,
"qmckl_transpose",
"Output matrix not allocated");
}
if (At.size[0] != A.size[1] || At.size[1] != A.size[0]) {
return qmckl_failwith( context,
QMCKL_INVALID_ARG_3,
"qmckl_transpose",
"Invalid size for At");
}
for (int64_t j=0 ; j<At.size[1] ; ++j)
for (int64_t i=0 ; i<At.size[0] ; ++i)
qmckl_mat(At, i, j) = qmckl_mat(A, j, i);
return QMCKL_SUCCESS;
}
#+end_src
*** Test :noexport:
#+begin_src c :comments link :tangle (eval c_test)
{
qmckl_matrix A;
qmckl_matrix At;
A = qmckl_matrix_alloc(context, 2, 3);
At = qmckl_matrix_alloc(context, 3, 2);
for (int j=0 ; j<3 ; ++j)
for (int i=0 ; i<2 ; ++i)
qmckl_mat(A, i, j) = (double) 10*i+j;
qmckl_exit_code rc = qmckl_transpose(context, A, At);
assert(rc == QMCKL_SUCCESS);
assert(A.size[0] == At.size[1]);
assert(A.size[1] == At.size[0]);
for (int j=0 ; j<3 ; ++j)
for (int i=0 ; i<2 ; ++i)
assert (qmckl_mat(A, i, j) == qmckl_mat(At, j, i));
qmckl_matrix_free(context, &A);
qmckl_matrix_free(context, &At);
}
#+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