qmckl/org/qmckl_blas.org

#+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 "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* 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* 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* 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) + 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
  
    For example:
    
** 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);

  #include <stdio.h>
  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