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_memory_private_func.h"
#include "qmckl_blas_private_type.h"
#include "qmckl_blas_private_func.h"
  #+end_src

  #+begin_src c :comments link :tangle (eval c_test) :noweb yes
#include "qmckl.h"
#include "assert.h"
#ifdef HAVE_CONFIG_H
#include "config.h"
#endif

#include "qmckl_blas_private_type.h"
#include "qmckl_blas_private_func.h"

int main() {
  qmckl_context context;
  context = qmckl_context_create();

  #+end_src

* Data types
  
** Vector
   
  | Variable | Type      | Description             |
  |----------+-----------+-------------------------|
  | ~size~   | ~int64_t~ | Dimension of the vector |
  | ~data~   | ~double*~ | Elements                |

   #+begin_src c :comments org :tangle (eval h_private_type)
typedef struct qmckl_vector {
  int64_t size;
  double* data;
} qmckl_vector;
   #+end_src


   #+begin_src c :comments org :tangle (eval h_private_func)
qmckl_vector
qmckl_vector_alloc( qmckl_context context, 
                    const int64_t size);
   #+end_src

   Allocates a new vector. If the allocation failed the size is zero.

   #+begin_src c :comments org :tangle (eval c)
qmckl_vector
qmckl_vector_alloc( qmckl_context context, 
                    const int64_t size)
{
 /* Should always be true by contruction */
  assert (size > (int64_t) 0);
  
  qmckl_vector result;
  result.size = size;

  qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero;
  mem_info.size = size * sizeof(double);
  result.data = (double*) qmckl_malloc (context, mem_info);

  if (result.data == NULL) {
    result.size = (int64_t) 0;
  }

  return result;
}
   #+end_src
   
   #+begin_src c :comments org :tangle (eval h_private_func)
qmckl_exit_code
qmckl_vector_free( qmckl_context context, 
                   qmckl_vector vector);
   #+end_src

   #+begin_src c :comments org :tangle (eval c)
qmckl_exit_code
qmckl_vector_free( qmckl_context context, 
                   qmckl_vector vector)
{
 /* Always true */
  assert (vector.data != NULL);
  
  qmckl_exit_code rc;
  
  rc = qmckl_free(context, vector.data);
  if (rc != QMCKL_SUCCESS) {
    return rc;
  }

  vector.size = (int64_t) 0;
  vector.data = NULL;
  return QMCKL_SUCCESS;
}
   #+end_src

** Matrix
   
  | Variable | Type         | Description                 |
  |----------+--------------+-----------------------------|
  | ~size~   | ~int64_t[2]~ | Dimension of each component |
  | ~data~   | ~double*~    | Elements                    |

  The dimensions use Fortran ordering: two elements differing by one
  in the first dimension are consecutive in memory.
  
   #+begin_src c :comments org :tangle (eval h_private_type)
typedef struct qmckl_matrix {
  int64_t size[2];
  double* data;
} qmckl_matrix;
   #+end_src


   #+begin_src c :comments org :tangle (eval h_private_func)
qmckl_matrix
qmckl_matrix_alloc( qmckl_context context, 
                    const int64_t size1,
                    const int64_t size2);
   #+end_src

   Allocates a new matrix. If the allocation failed the sizes are zero.

   #+begin_src c :comments org :tangle (eval c)
qmckl_matrix
qmckl_matrix_alloc( qmckl_context context, 
                    const int64_t size1,
                    const int64_t size2)
{
 /* Should always be true by contruction */
  assert (size1 * size2 > (int64_t) 0);
  
  qmckl_matrix result;

  result.size[0] = size1;
  result.size[1] = size2;

  qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero;
  mem_info.size = size1 * size2 * sizeof(double);
  result.data = (double*) qmckl_malloc (context, mem_info);

  if (result.data == NULL) {
    result.size[0] = (int64_t) 0;
    result.size[1] = (int64_t) 0;
  }

  return result;
}
   #+end_src
   
   #+begin_src c :comments org :tangle (eval h_private_func)
qmckl_exit_code
qmckl_matrix_free( qmckl_context context, 
                   qmckl_matrix matrix);
   #+end_src

   #+begin_src c :comments org :tangle (eval c)
qmckl_exit_code
qmckl_matrix_free( qmckl_context context, 
                   qmckl_matrix matrix)
{
 /* Always true */
  assert (matrix.data != NULL);
  
  qmckl_exit_code rc;
  
  rc = qmckl_free(context, matrix.data);
  if (rc != QMCKL_SUCCESS) {
    return rc;
  }
  matrix.data = NULL;
  matrix.size[0] = (int64_t) 0;
  matrix.size[1] = (int64_t) 0;

  return QMCKL_SUCCESS;
}
   #+end_src

** Tensor
   
  | Variable | Type                              | Description                 |
  |----------+-----------------------------------+-----------------------------|
  | ~order~  | ~int64_t~                         | Order of the tensor         |
  | ~size~   | ~int64_t[QMCKL_TENSOR_ORDER_MAX]~ | Dimension of each component |
  | ~data~   | ~double*~                         | Elements                    |

  The dimensions use Fortran ordering: two elements differing by one
  in the first dimension are consecutive in memory.
  
   #+begin_src c :comments org :tangle (eval h_private_type)
#define QMCKL_TENSOR_ORDER_MAX 16

typedef struct qmckl_tensor {
  int64_t order;
  int64_t size[QMCKL_TENSOR_ORDER_MAX];
  double* data;
} qmckl_tensor;
   #+end_src


   #+begin_src c :comments org :tangle (eval h_private_func)
qmckl_tensor
qmckl_tensor_alloc( qmckl_context context, 
                    const int64_t order,
                    const int64_t* size);
   #+end_src

   Allocates memory for a tensor. If the allocation failed, the size
   is zero.

   #+begin_src c :comments org :tangle (eval c)
qmckl_tensor
qmckl_tensor_alloc( qmckl_context context, 
                    const int64_t  order,
                    const int64_t* size)
{
 /* Should always be true by contruction */
  assert (order > 0);
  assert (order <= QMCKL_TENSOR_ORDER_MAX);
  assert (size  != NULL);
  
  qmckl_tensor result;
  result.order = order;

  int64_t prod_size = (int64_t) 1;
  for (int64_t i=0 ; i<order ; ++i) {
    assert (size[i] > (int64_t) 0);
    result.size[i] = size[i];
    prod_size *= size[i];
  }

  qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero;
  mem_info.size = prod_size * sizeof(double);

  result.data = (double*) qmckl_malloc (context, mem_info);

  if (result.data == NULL) {
    memset(&result, 0, sizeof(qmckl_tensor));
  }

  return result;
}
   #+end_src
   
   #+begin_src c :comments org :tangle (eval h_private_func)
qmckl_exit_code
qmckl_tensor_free( qmckl_context context, 
                   qmckl_tensor tensor);
   #+end_src

   #+begin_src c :comments org :tangle (eval c)
qmckl_exit_code
qmckl_tensor_free( qmckl_context context, 
                   qmckl_tensor tensor)
{
 /* Always true */
  assert (tensor.data != NULL);
  
  qmckl_exit_code rc;
  
  rc = qmckl_free(context, tensor.data);
  if (rc != QMCKL_SUCCESS) {
    return rc;
  }
  
  memset(&tensor, 0, sizeof(qmckl_tensor));

  return QMCKL_SUCCESS;
}
   #+end_src

** Reshaping

   Reshaping occurs in-place and the pointer to the data is copied.
   
*** Vector -> Matrix

    #+begin_src c :comments org :tangle (eval h_private_func)
qmckl_matrix
qmckl_matrix_of_vector(const qmckl_vector vector,
                       const int64_t size1,
                       const int64_t size2);
    #+end_src

    Reshapes a vector into a matrix.

    #+begin_src c :comments org :tangle (eval c)
qmckl_matrix
qmckl_matrix_of_vector(const qmckl_vector vector,
                       const int64_t size1,
                       const int64_t size2)
{ 
 /* Always true */
  assert (size1 * size2 == vector.size);

  qmckl_matrix result;

  result.size[0] = size1;
  result.size[1] = size2;
  result.data    = vector.data;

  return result;
}
    #+end_src

*** Vector -> Tensor

    #+begin_src c :comments org :tangle (eval h_private_func)
qmckl_tensor
qmckl_tensor_of_vector(const qmckl_vector vector,
                       const int64_t order,
                       const int64_t* size);
    #+end_src

    Reshapes a vector into a tensor.

    #+begin_src c :comments org :tangle (eval c)
qmckl_tensor
qmckl_tensor_of_vector(const qmckl_vector vector,
                       const int64_t order,
                       const int64_t* size)
{ 
  qmckl_tensor result;

  int64_t prod_size = 1;
  for (int64_t i=0 ; i<order ; ++i) {
    result.size[i] = size[i];
    prod_size *= size[i];
  }
  assert (prod_size == vector.size);

  result.data = vector.data;

  return result;
}
    #+end_src

*** Matrix -> Vector

    #+begin_src c :comments org :tangle (eval h_private_func)
qmckl_vector
qmckl_vector_of_matrix(const qmckl_matrix matrix,
                       const int64_t size);
    #+end_src

    Reshapes a matrix into a vector.

    #+begin_src c :comments org :tangle (eval c)
qmckl_vector
qmckl_vector_of_matrix(const qmckl_matrix matrix,
                       const int64_t size)
{ 
 /* Always true */
  assert (matrix.size[0] * matrix.size[1] == size);

  qmckl_vector result;

  result.size = size;
  result.data = matrix.data;

  return result;
}
    #+end_src

*** Matrix -> Tensor

    #+begin_src c :comments org :tangle (eval h_private_func)
qmckl_tensor
qmckl_tensor_of_matrix(const qmckl_matrix matrix,
                       const int64_t order,
                       const int64_t* size);
    #+end_src

    Reshapes a matrix into a tensor.

    #+begin_src c :comments org :tangle (eval c)
qmckl_tensor
qmckl_tensor_of_matrix(const qmckl_matrix matrix,
                       const int64_t order,
                       const int64_t* size)
{ 
  qmckl_tensor result;

  int64_t prod_size = 1;
  for (int64_t i=0 ; i<order ; ++i) {
    result.size[i] = size[i];
    prod_size *= size[i];
  }
  assert (prod_size == matrix.size[0] * matrix.size[1]);

  result.data = matrix.data;

  return result;
}
    #+end_src

*** Tensor -> Vector

    #+begin_src c :comments org :tangle (eval h_private_func)
qmckl_vector
qmckl_vector_of_tensor(const qmckl_tensor tensor,
                       const int64_t size);
    #+end_src

    Reshapes a tensor into a vector.

    #+begin_src c :comments org :tangle (eval c)
qmckl_vector
qmckl_vector_of_tensor(const qmckl_tensor tensor,
                       const int64_t size)
{ 
 /* Always true */
  int64_t prod_size = (int64_t) 1;
  for (int64_t i=0 ; i<tensor.order ; i++) {
    prod_size *= tensor.size[i];
  }
  assert (prod_size == size);

  qmckl_vector result;

  result.size = size;
  result.data = tensor.data;

  return result;
}
    #+end_src

*** Tensor -> Matrix

    #+begin_src c :comments org :tangle (eval h_private_func)
qmckl_matrix
qmckl_matrix_of_tensor(const qmckl_tensor tensor,
                       const int64_t size1,
                       const int64_t size2);
    #+end_src

    Reshapes a tensor into a vector.

    #+begin_src c :comments org :tangle (eval c)
qmckl_matrix
qmckl_matrix_of_tensor(const qmckl_tensor tensor,
                       const int64_t size1,
                       const int64_t size2)
{ 
 /* Always true */
  int64_t prod_size = (int64_t) 1;
  for (int64_t i=0 ; i<tensor.order ; i++) {
    prod_size *= tensor.size[i];
  }
  assert (prod_size == size1 * size2);

  qmckl_matrix result;

  result.size[0] = size1;
  result.size[1] = size2;
  result.data = tensor.data;

  return result;
}
    #+end_src

** Access macros

    #+begin_src c :comments org :tangle (eval h_private_func)
#define qmckl_vec(v, i) v.data[i]
#define qmckl_mat(m, i, j) m.data[(i) + (j)*m.size[0]]

#define qmckl_ten3(t, i, j, k) t.data[(i) + m.size[0]*((j) + size[1]*(k))]
#define qmckl_ten4(t, i, j, k, l) t.data[(i) + m.size[0]*((j) + size[1]*((k) + size[2]*(l)))]
#define qmckl_ten5(t, i, j, k, l, m) t.data[(i) + m.size[0]*((j) + size[1]*((k) + size[2]*((l) + size[3]*(m))))]
    #+end_src
  
** Tests

    #+begin_src c :comments link :tangle (eval c_test)
{
  int64_t m = 3;
  int64_t n = 4;
  int64_t p = m*n;
  qmckl_vector vec = qmckl_vector_alloc(context, p);

  for (int64_t i=0 ; i<p ; ++i) 
    qmckl_vec(vec, i) = (double) i;

  for (int64_t i=0 ; i<p ; ++i) 
    assert( vec.data[i] == (double) i );

  qmckl_matrix mat = qmckl_matrix_of_vector(vec, m, n);
  assert (mat.size[0] == m);
  assert (mat.size[1] == n);
  assert (mat.data == vec.data);
  
  for (int64_t j=0 ; j<n ; ++j)
    for (int64_t i=0 ; i<m ; ++i)
      assert ( qmckl_mat(mat, i, j) == qmckl_vec(vec, i+j*m)) ;

  qmckl_vector vec2 = qmckl_vector_of_matrix(mat, p);
  assert (vec2.size == p);
  assert (vec2.data == vec.data);
  for (int64_t i=0 ; i<p ; ++i) 
      assert ( qmckl_vec(vec2, i) == qmckl_vec(vec, i) ) ;

  qmckl_vector_free(context, vec);

}
    #+end_src
* Matrix operations

** ~qmckl_dgemm~

   Matrix multiplication:

   \[
   C_{ij} = \beta C_{ij} + \alpha \sum_{k} A_{ik} \cdot B_{kj}
   \]

#  TODO: Add description about the external library dependence.

   #+NAME: qmckl_dgemm_args
   | Variable  | Type            | In/Out | Description                           |
   |-----------+-----------------+--------+---------------------------------------|
   | ~context~ | ~qmckl_context~ | in     | Global state                          |
   | ~TransA~  | ~char~          | in     | 'T' is transposed                     |
   | ~TransB~  | ~char~          | in     | 'T' is transposed                     |
   | ~m~       | ~int64_t~       | in     | Number of rows of the input matrix    |
   | ~n~       | ~int64_t~       | in     | Number of columns of the input matrix |
   | ~k~       | ~int64_t~       | in     | Number of columns of the input matrix |
   | ~alpha~   | ~double~        | in     | Number of columns of the input matrix |
   | ~A~       | ~double[][lda]~ | in     | Array containing the matrix $A$       |
   | ~lda~     | ~int64_t~       | in     | Leading dimension of array ~A~        |
   | ~B~       | ~double[][ldb]~ | in     | Array containing the matrix $B$       |
   | ~ldb~     | ~int64_t~       | in     | Leading dimension of array ~B~        |
   | ~beta~    | ~double~        | in     | Array containing the matrix $B$       |
   | ~C~       | ~double[][ldc]~ | out    | Array containing the matrix $B$       |
   | ~ldc~     | ~int64_t~       | in     | Leading dimension of array ~B~        |

   Requirements:

    - ~context~ is not ~QMCKL_NULL_CONTEXT~
    - ~m > 0~
    - ~n > 0~
    - ~k > 0~
    - ~lda >= m~
    - ~ldb >= n~
    - ~ldc >= n~
    - ~A~ is allocated with at least $m \times k \times 8$ bytes
    - ~B~ is allocated with at least $k \times n \times 8$ bytes
    - ~C~ is allocated with at least $m \times n \times 8$ bytes

    #+CALL: generate_c_header(table=qmckl_dgemm_args,rettyp="qmckl_exit_code",fname="qmckl_dgemm")

    #+RESULTS:
    #+begin_src c :tangle (eval h_func) :comments org
    qmckl_exit_code qmckl_dgemm (
          const qmckl_context context,
          const char TransA,
          const char TransB,
          const int64_t m,
          const int64_t n,
          const int64_t k,
          const double alpha,
          const double* A,
          const int64_t lda,
          const double* B,
          const int64_t ldb,
          const double beta,
          double* const C,
          const int64_t ldc );
    #+end_src

    #+begin_src f90 :tangle (eval f)
integer function qmckl_dgemm_f(context, TransA, TransB, &
     m, n, k, alpha, A, LDA, B, LDB, beta, C, LDC) &
     result(info)
  use qmckl
  implicit none
  integer(qmckl_context), intent(in)  :: context
  character             , intent(in)  :: TransA, TransB
  integer*8             , intent(in)  :: m, n, k
  double precision      , intent(in)  :: alpha, beta
  integer*8             , intent(in)  :: lda
  double precision      , intent(in)  :: A(lda,*)
  integer*8             , intent(in)  :: ldb
  double precision      , intent(in)  :: B(ldb,*)
  integer*8             , intent(in)  :: ldc
  double precision      , intent(out) :: C(ldc,*)

  info = QMCKL_SUCCESS

  if (context == QMCKL_NULL_CONTEXT) then
     info = QMCKL_INVALID_CONTEXT
     return
  endif

  if (m <= 0_8) then
     info = QMCKL_INVALID_ARG_4
     return
  endif

  if (n <= 0_8) then
     info = QMCKL_INVALID_ARG_5
     return
  endif

  if (k <= 0_8) then
     info = QMCKL_INVALID_ARG_6
     return
  endif

  if (LDA <= 0) then
     info = QMCKL_INVALID_ARG_9
     return
  endif

  if (LDB <= 0) then
     info = QMCKL_INVALID_ARG_11
     return
  endif

  if (LDC <= 0) then
     info = QMCKL_INVALID_ARG_14
     return
  endif

  call dgemm(transA, transB, int(m,4), int(n,4), int(k,4), &
       alpha, A, int(LDA,4), B, int(LDB,4), beta, C, int(LDC,4))

end function qmckl_dgemm_f
    #+end_src

*** C interface                                                    :noexport:

    #+CALL: generate_c_interface(table=qmckl_dgemm_args,rettyp="qmckl_exit_code",fname="qmckl_dgemm")

    #+RESULTS:
    #+begin_src f90 :tangle (eval f) :comments org :exports none
    integer(c_int32_t) function qmckl_dgemm &
        (context, TransA, TransB, m, n, k, alpha, A, lda, B, ldb, beta, C, ldc) &
        bind(C) result(info)

      use, intrinsic :: iso_c_binding
      implicit none

      integer (c_int64_t) , intent(in)  , value :: context
      character           , intent(in)  , value :: TransA
      character           , intent(in)  , value :: TransB
      integer (c_int64_t) , intent(in)  , value :: m
      integer (c_int64_t) , intent(in)  , value :: n
      integer (c_int64_t) , intent(in)  , value :: k
      real    (c_double ) , intent(in)  , value :: alpha
      real    (c_double ) , intent(in)          :: A(lda,*)
      integer (c_int64_t) , intent(in)  , value :: lda
      real    (c_double ) , intent(in)          :: B(ldb,*)
      integer (c_int64_t) , intent(in)  , value :: ldb
      real    (c_double ) , intent(in)  , value :: beta
      real    (c_double ) , intent(out)         :: C(ldc,*)
      integer (c_int64_t) , intent(in)  , value :: ldc

      integer(c_int32_t), external :: qmckl_dgemm_f
      info = qmckl_dgemm_f &
             (context, TransA, TransB, m, n, k, alpha, A, lda, B, ldb, beta, C, ldc)

    end function qmckl_dgemm
    #+end_src


    #+CALL: generate_f_interface(table=qmckl_dgemm_args,rettyp="qmckl_exit_code",fname="qmckl_dgemm")

    #+RESULTS:
    #+begin_src f90 :tangle (eval fh_func) :comments org :exports none
    interface
      integer(c_int32_t) function qmckl_dgemm &
          (context, TransA, TransB, m, n, k, alpha, A, lda, B, ldb, beta, C, ldc) &
          bind(C)
        use, intrinsic :: iso_c_binding
        import
        implicit none

        integer (c_int64_t) , intent(in)  , value :: context
        character           , intent(in)  , value :: TransA
        character           , intent(in)  , value :: TransB
        integer (c_int64_t) , intent(in)  , value :: m
        integer (c_int64_t) , intent(in)  , value :: n
        integer (c_int64_t) , intent(in)  , value :: k
        real    (c_double ) , intent(in)  , value :: alpha
        real    (c_double ) , intent(in)          :: A(lda,*)
        integer (c_int64_t) , intent(in)  , value :: lda
        real    (c_double ) , intent(in)          :: B(ldb,*)
        integer (c_int64_t) , intent(in)  , value :: ldb
        real    (c_double ) , intent(in)  , value :: beta
        real    (c_double ) , intent(out)         :: C(ldc,*)
        integer (c_int64_t) , intent(in)  , value :: ldc

      end function qmckl_dgemm
    end interface
    #+end_src

*** Test                                                           :noexport:

    #+begin_src f90 :tangle (eval f_test)
integer(qmckl_exit_code) function test_qmckl_dgemm(context) bind(C)
  use qmckl
  implicit none
  integer(qmckl_context), intent(in), value :: context

  double precision, allocatable :: A(:,:), B(:,:), C(:,:), D(:,:)
  integer*8                     :: m, n, k, LDA, LDB, LDC
  integer*8                     :: i,j,l
  character                     :: TransA, TransB
  double precision              :: x, alpha, beta

  TransA = 'N'
  TransB = 'N'
  m = 1_8
  k = 4_8
  n = 6_8
  LDA = m
  LDB = k
  LDC = m

  allocate( A(LDA,k), B(LDB,n) , C(LDC,n), D(LDC,n))

  A = 0.d0
  B = 0.d0
  C = 0.d0
  D = 0.d0
  alpha = 1.0d0
  beta  = 0.0d0
  do j=1,k
     do i=1,m
        A(i,j) = -10.d0 + dble(i+j)
     end do
  end do

  do j=1,n
     do i=1,k
        B(i,j) = -10.d0 + dble(i+j)
     end do
  end do

  test_qmckl_dgemm = qmckl_dgemm(context, TransA, TransB, m, n, k, alpha, A, LDA, B, LDB, beta, C, LDC)

  if (test_qmckl_dgemm /= QMCKL_SUCCESS) return

  test_qmckl_dgemm = QMCKL_FAILURE

  x = 0.d0
  do j=1,n
     do i=1,m
        do l=1,k
           D(i,j) = D(i,j) + A(i,l)*B(l,j)
        end do
        x = x + (D(i,j) - C(i,j))**2
     end do
  end do

  if (dabs(x) <= 1.d-12) then
     test_qmckl_dgemm = QMCKL_SUCCESS
  endif

  deallocate(A,B,C,D)

end function test_qmckl_dgemm
    #+end_src

    #+begin_src c :comments link :tangle (eval c_test)
qmckl_exit_code test_qmckl_dgemm(qmckl_context context);
assert(QMCKL_SUCCESS == test_qmckl_dgemm(context));
    #+end_src

** ~qmckl_adjugate~

   Given a matrix $\mathbf{A}$, the adjugate matrix
   $\text{adj}(\mathbf{A})$ is the transpose of the cofactors matrix
   of $\mathbf{A}$.

   \[
   \mathbf{B} = \text{adj}(\mathbf{A}) = \text{det}(\mathbf{A}) \, \mathbf{A}^{-1}
   \]

   See also: https://en.wikipedia.org/wiki/Adjugate_matrix

   #+NAME: qmckl_adjugate_args
   | Variable  | Type            | In/Out | Description                                    |
   |-----------+-----------------+--------+------------------------------------------------|
   | ~context~ | ~qmckl_context~ | in     | Global state                                   |
   | ~n~       | ~int64_t~       | in     | Number of rows and columns of the input matrix |
   | ~A~       | ~double[][lda]~ | in     | Array containing the $n \times n$ matrix $A$   |
   | ~lda~     | ~int64_t~       | in     | Leading dimension of array ~A~                 |
   | ~B~       | ~double[][ldb]~ | out    | Adjugate of $A$                                |
   | ~ldb~     | ~int64_t~       | in     | Leading dimension of array ~B~                 |
   | ~det_l~   | ~double~        | inout  | determinant of $A$                             |

   Requirements:

    - ~context~ is not ~QMCKL_NULL_CONTEXT~
    - ~n > 0~
    - ~lda >= m~
    - ~A~ is allocated with at least $m \times m \times 8$ bytes
    - ~ldb >= m~
    - ~B~ is allocated with at least $m \times m \times 8$ bytes

    #+CALL: generate_c_header(table=qmckl_adjugate_args,rettyp="qmckl_exit_code",fname="qmckl_adjugate")

    #+RESULTS:
    #+begin_src c :tangle (eval h_func) :comments org
    qmckl_exit_code qmckl_adjugate (
          const qmckl_context context,
          const int64_t n,
          const double* A,
          const int64_t lda,
          double* const B,
          const int64_t ldb,
          double* det_l );
    #+end_src

   For small matrices (\le 5\times 5), we use specialized routines
   for performance motivations. For larger sizes, we rely on the
   LAPACK library.

    #+begin_src f90 :tangle (eval f)
integer function qmckl_adjugate_f(context, na, A, LDA, B, ldb, det_l) &
     result(info)
  use qmckl
  implicit none
  integer(qmckl_context)  , intent(in)  :: context
  double precision, intent(in)          :: A (LDA,*)
  integer*8, intent(in)                 :: LDA
  double precision, intent(out)         :: B (LDB,*)
  integer*8, intent(in)                 :: LDB
  integer*8, intent(in)                 :: na
  double precision, intent(inout)       :: det_l

  integer    :: i,j

  info = QMCKL_SUCCESS

  if (context == QMCKL_NULL_CONTEXT) then
     info = QMCKL_INVALID_CONTEXT
     return
  endif

  if (na <= 0_8) then
     info = QMCKL_INVALID_ARG_2
     return
  endif

  if (LDA <= 0_8) then
     info = QMCKL_INVALID_ARG_4
     return
  endif

  if (LDA < na) then
     info = QMCKL_INVALID_ARG_4
     return
  endif

  select case (na)
  case (5)
     call adjugate5(A,LDA,B,LDB,na,det_l)
  case (4)
     call adjugate4(A,LDA,B,LDB,na,det_l)
  case (3)
     call adjugate3(A,LDA,B,LDB,na,det_l)
  case (2)
     call adjugate2(A,LDA,B,LDB,na,det_l)
  case (1)
    det_l = a(1,1)
    b(1,1) = 1.d0
  case default
     call adjugate_general(context, na, A, LDA, B, LDB, det_l)
  end select

end function qmckl_adjugate_f
    #+end_src

    #+begin_src f90 :tangle (eval f) :exports none
subroutine adjugate2(A,LDA,B,LDB,na,det_l)
  implicit none
  double precision, intent(in)    :: A (LDA,na)
  double precision, intent(out)   :: B (LDA,na)
  integer*8, intent(in)           :: LDA, LDB
  integer*8, intent(in)           :: na
  double precision, intent(inout) :: det_l

  double precision :: C(2,2)

  call cofactor2(A,LDA,C,2_8,na,det_l)

  B(1,1) = C(1,1)
  B(2,1) = C(1,2)
  B(1,2) = C(2,1)
  B(2,2) = C(2,2)

end subroutine adjugate2

subroutine adjugate3(a,LDA,B,LDB,na,det_l)
  implicit none
  double precision, intent(in)    :: A (LDA,na)
  double precision, intent(out)   :: B (LDA,na)
  integer*8, intent(in)           :: LDA, LDB
  integer*8, intent(in)           :: na
  double precision, intent(inout) :: det_l

  double precision :: C(4,3)

  call cofactor3(A,LDA,C,4_8,na,det_l)

  B(1,1) = C(1,1)
  B(1,2) = C(2,1)
  B(1,3) = C(3,1)
  B(2,1) = C(1,2)
  B(2,2) = C(2,2)
  B(2,3) = C(3,2)
  B(3,1) = C(1,3)
  B(3,2) = C(2,3)
  B(3,3) = C(3,3)

end subroutine adjugate3

subroutine adjugate4(a,LDA,B,LDB,na,det_l)
  implicit none
  double precision, intent(in)    :: A (LDA,na)
  double precision, intent(out)   :: B (LDA,na)
  integer*8, intent(in)           :: LDA, LDB
  integer*8, intent(in)           :: na
  double precision, intent(inout) :: det_l

  double precision :: C(4,4)

  call cofactor4(A,LDA,C,4_8,na,det_l)
  B(1,1) = C(1,1)
  B(1,2) = C(2,1)
  B(1,3) = C(3,1)
  B(1,4) = C(4,1)
  B(2,1) = C(1,2)
  B(2,2) = C(2,2)
  B(2,3) = C(3,2)
  B(2,4) = C(4,2)
  B(3,1) = C(1,3)
  B(3,2) = C(2,3)
  B(3,3) = C(3,3)
  B(3,4) = C(4,3)
  B(4,1) = C(1,4)
  B(4,2) = C(2,4)
  B(4,3) = C(3,4)
  B(4,4) = C(4,4)

end subroutine adjugate4

subroutine adjugate5(A,LDA,B,LDB,na,det_l)
  implicit none
  double precision, intent(in)    :: A (LDA,na)
  double precision, intent(out)   :: B (LDA,na)
  integer*8, intent(in)           :: LDA, LDB
  integer*8, intent(in)           :: na
  double precision, intent(inout) :: det_l

  double precision  :: C(8,5)

  call cofactor5(A,LDA,C,8_8,na,det_l)

  B(1,1) = C(1,1)
  B(1,2) = C(2,1)
  B(1,3) = C(3,1)
  B(1,4) = C(4,1)
  B(1,5) = C(5,1)
  B(2,1) = C(1,2)
  B(2,2) = C(2,2)
  B(2,3) = C(3,2)
  B(2,4) = C(4,2)
  B(2,5) = C(5,2)
  B(3,1) = C(1,3)
  B(3,2) = C(2,3)
  B(3,3) = C(3,3)
  B(3,4) = C(4,3)
  B(3,5) = C(5,3)
  B(4,1) = C(1,4)
  B(4,2) = C(2,4)
  B(4,3) = C(3,4)
  B(4,4) = C(4,4)
  B(4,5) = C(5,4)
  B(5,1) = C(1,5)
  B(5,2) = C(2,5)
  B(5,3) = C(3,5)
  B(5,4) = C(4,5)
  B(5,5) = C(5,5)

end subroutine adjugate5

subroutine cofactor2(a,LDA,b,LDB,na,det_l)
  implicit none
  double precision, intent(in)    :: A (LDA,na)
  double precision, intent(out)   :: B (LDA,na)
  integer*8, intent(in)           :: LDA, LDB
  integer*8        :: na
  double precision :: det_l

  det_l  =  a(1,1)*a(2,2) - a(1,2)*a(2,1)
  b(1,1) =  a(2,2)
  b(2,1) = -a(2,1)
  b(1,2) = -a(1,2)
  b(2,2) =  a(1,1)
end subroutine cofactor2

subroutine cofactor3(a,LDA,b,LDB,na,det_l)
  implicit none
  double precision, intent(in)    :: A (LDA,na)
  double precision, intent(out)   :: B (LDA,na)
  integer*8, intent(in)           :: LDA, LDB
  integer*8, intent(in)           :: na
  double precision, intent(inout) :: det_l
  integer :: i

  det_l = a(1,1)*(a(2,2)*a(3,3)-a(2,3)*a(3,2)) &
         -a(1,2)*(a(2,1)*a(3,3)-a(2,3)*a(3,1)) &
         +a(1,3)*(a(2,1)*a(3,2)-a(2,2)*a(3,1))

  b(1,1) =  a(2,2)*a(3,3) - a(2,3)*a(3,2)
  b(2,1) =  a(2,3)*a(3,1) - a(2,1)*a(3,3)
  b(3,1) =  a(2,1)*a(3,2) - a(2,2)*a(3,1)

  b(1,2) =  a(1,3)*a(3,2) - a(1,2)*a(3,3)
  b(2,2) =  a(1,1)*a(3,3) - a(1,3)*a(3,1)
  b(3,2) =  a(1,2)*a(3,1) - a(1,1)*a(3,2)

  b(1,3) =  a(1,2)*a(2,3) - a(1,3)*a(2,2)
  b(2,3) =  a(1,3)*a(2,1) - a(1,1)*a(2,3)
  b(3,3) =  a(1,1)*a(2,2) - a(1,2)*a(2,1)

end subroutine cofactor3

subroutine cofactor4(a,LDA,b,LDB,na,det_l)
  implicit none
  double precision, intent(in)    :: A (LDA,na)
  double precision, intent(out)   :: B (LDA,na)
  integer*8, intent(in)           :: LDA, LDB
  integer*8, intent(in)           :: na
  double precision, intent(inout) :: det_l
  integer :: i,j
  det_l =  a(1,1)*(a(2,2)*(a(3,3)*a(4,4)-a(3,4)*a(4,3))  &
                  -a(2,3)*(a(3,2)*a(4,4)-a(3,4)*a(4,2))  &
                  +a(2,4)*(a(3,2)*a(4,3)-a(3,3)*a(4,2))) &
          -a(1,2)*(a(2,1)*(a(3,3)*a(4,4)-a(3,4)*a(4,3))  &
                  -a(2,3)*(a(3,1)*a(4,4)-a(3,4)*a(4,1))  &
                  +a(2,4)*(a(3,1)*a(4,3)-a(3,3)*a(4,1))) &
          +a(1,3)*(a(2,1)*(a(3,2)*a(4,4)-a(3,4)*a(4,2))  &
                  -a(2,2)*(a(3,1)*a(4,4)-a(3,4)*a(4,1))  &
                  +a(2,4)*(a(3,1)*a(4,2)-a(3,2)*a(4,1))) &
          -a(1,4)*(a(2,1)*(a(3,2)*a(4,3)-a(3,3)*a(4,2))  &
                  -a(2,2)*(a(3,1)*a(4,3)-a(3,3)*a(4,1))  &
                  +a(2,3)*(a(3,1)*a(4,2)-a(3,2)*a(4,1)))

  b(1,1) =  a(2,2)*(a(3,3)*a(4,4)-a(3,4)*a(4,3))-a(2,3)*(a(3,2)*a(4,4)-a(3,4)*a(4,2))+a(2,4)*(a(3,2)*a(4,3)-a(3,3)*a(4,2))
  b(2,1) = -a(2,1)*(a(3,3)*a(4,4)-a(3,4)*a(4,3))+a(2,3)*(a(3,1)*a(4,4)-a(3,4)*a(4,1))-a(2,4)*(a(3,1)*a(4,3)-a(3,3)*a(4,1))
  b(3,1) =  a(2,1)*(a(3,2)*a(4,4)-a(3,4)*a(4,2))-a(2,2)*(a(3,1)*a(4,4)-a(3,4)*a(4,1))+a(2,4)*(a(3,1)*a(4,2)-a(3,2)*a(4,1))
  b(4,1) = -a(2,1)*(a(3,2)*a(4,3)-a(3,3)*a(4,2))+a(2,2)*(a(3,1)*a(4,3)-a(3,3)*a(4,1))-a(2,3)*(a(3,1)*a(4,2)-a(3,2)*a(4,1))

  b(1,2) = -a(1,2)*(a(3,3)*a(4,4)-a(3,4)*a(4,3))+a(1,3)*(a(3,2)*a(4,4)-a(3,4)*a(4,2))-a(1,4)*(a(3,2)*a(4,3)-a(3,3)*a(4,2))
  b(2,2) =  a(1,1)*(a(3,3)*a(4,4)-a(3,4)*a(4,3))-a(1,3)*(a(3,1)*a(4,4)-a(3,4)*a(4,1))+a(1,4)*(a(3,1)*a(4,3)-a(3,3)*a(4,1))
  b(3,2) = -a(1,1)*(a(3,2)*a(4,4)-a(3,4)*a(4,2))+a(1,2)*(a(3,1)*a(4,4)-a(3,4)*a(4,1))-a(1,4)*(a(3,1)*a(4,2)-a(3,2)*a(4,1))
  b(4,2) =  a(1,1)*(a(3,2)*a(4,3)-a(3,3)*a(4,2))-a(1,2)*(a(3,1)*a(4,3)-a(3,3)*a(4,1))+a(1,3)*(a(3,1)*a(4,2)-a(3,2)*a(4,1))

  b(1,3) =  a(1,2)*(a(2,3)*a(4,4)-a(2,4)*a(4,3))-a(1,3)*(a(2,2)*a(4,4)-a(2,4)*a(4,2))+a(1,4)*(a(2,2)*a(4,3)-a(2,3)*a(4,2))
  b(2,3) = -a(1,1)*(a(2,3)*a(4,4)-a(2,4)*a(4,3))+a(1,3)*(a(2,1)*a(4,4)-a(2,4)*a(4,1))-a(1,4)*(a(2,1)*a(4,3)-a(2,3)*a(4,1))
  b(3,3) =  a(1,1)*(a(2,2)*a(4,4)-a(2,4)*a(4,2))-a(1,2)*(a(2,1)*a(4,4)-a(2,4)*a(4,1))+a(1,4)*(a(2,1)*a(4,2)-a(2,2)*a(4,1))
  b(4,3) = -a(1,1)*(a(2,2)*a(4,3)-a(2,3)*a(4,2))+a(1,2)*(a(2,1)*a(4,3)-a(2,3)*a(4,1))-a(1,3)*(a(2,1)*a(4,2)-a(2,2)*a(4,1))

  b(1,4) = -a(1,2)*(a(2,3)*a(3,4)-a(2,4)*a(3,3))+a(1,3)*(a(2,2)*a(3,4)-a(2,4)*a(3,2))-a(1,4)*(a(2,2)*a(3,3)-a(2,3)*a(3,2))
  b(2,4) =  a(1,1)*(a(2,3)*a(3,4)-a(2,4)*a(3,3))-a(1,3)*(a(2,1)*a(3,4)-a(2,4)*a(3,1))+a(1,4)*(a(2,1)*a(3,3)-a(2,3)*a(3,1))
  b(3,4) = -a(1,1)*(a(2,2)*a(3,4)-a(2,4)*a(3,2))+a(1,2)*(a(2,1)*a(3,4)-a(2,4)*a(3,1))-a(1,4)*(a(2,1)*a(3,2)-a(2,2)*a(3,1))
  b(4,4) =  a(1,1)*(a(2,2)*a(3,3)-a(2,3)*a(3,2))-a(1,2)*(a(2,1)*a(3,3)-a(2,3)*a(3,1))+a(1,3)*(a(2,1)*a(3,2)-a(2,2)*a(3,1))

end subroutine cofactor4

subroutine cofactor5(A,LDA,B,LDB,na,det_l)
  implicit none
  double precision, intent(in)    :: A (LDA,na)
  double precision, intent(out)   :: B (LDA,na)
  integer*8, intent(in)           :: LDA, LDB
  integer*8, intent(in)           :: na
  double precision, intent(inout) :: det_l
  integer :: i,j

 det_l = a(1,1)*(a(2,2)*(a(3,3)*(a(4,4)*a(5,5)-a(4,5)*a(5,4))-a(3,4)*( &
 a(4,3)*a(5,5)-a(4,5)*a(5,3))+a(3,5)*(a(4,3)*a(5,4)-a(4,4)*a(5,3)))- &
 a(2,3)*(a(3,2)*(a(4,4)*a(5,5)-a(4,5)*a(5,4))-a(3,4)*(a(4,2)*a(5,5)- &
 a(4,5)*a(5,2))+a(3,5)*(a(4,2)*a(5,4)-a(4,4)*a(5,2)))+a(2,4)*(a(3,2)*( &
 a(4,3)*a(5,5)-a(4,5)*a(5,3))-a(3,3)*(a(4,2)*a(5,5)-a(4,5)*a(5,2))+ &
 a(3,5)*(a(4,2)*a(5,3)-a(4,3)*a(5,2)))-a(2,5)*(a(3,2)*(a(4,3)*a(5,4)- &
 a(4,4)*a(5,3))-a(3,3)*(a(4,2)*a(5,4)-a(4,4)*a(5,2))+a(3,4)*(a(4,2)* &
 a(5,3)-a(4,3)*a(5,2))))-a(1,2)*(a(2,1)*(a(3,3)*(a(4,4)*a(5,5)-a(4,5)* &
 a(5,4))-a(3,4)*(a(4,3)*a(5,5)-a(4,5)*a(5,3))+a(3,5)*(a(4,3)*a(5,4)- &
 a(4,4)*a(5,3)))-a(2,3)*(a(3,1)*(a(4,4)*a(5,5)-a(4,5)*a(5,4))-a(3,4)*( &
 a(4,1)*a(5,5)-a(4,5)*a(5,1))+a(3,5)*(a(4,1)*a(5,4)-a(4,4)*a(5,1)))+ &
 a(2,4)*(a(3,1)*(a(4,3)*a(5,5)-a(4,5)*a(5,3))-a(3,3)*(a(4,1)*a(5,5)- &
 a(4,5)*a(5,1))+a(3,5)*(a(4,1)*a(5,3)-a(4,3)*a(5,1)))-a(2,5)*(a(3,1)*( &
 a(4,3)*a(5,4)-a(4,4)*a(5,3))-a(3,3)*(a(4,1)*a(5,4)-a(4,4)*a(5,1))+ &
 a(3,4)*(a(4,1)*a(5,3)-a(4,3)*a(5,1))))+a(1,3)*(a(2,1)*(a(3,2)*(a(4,4)* &
 a(5,5)-a(4,5)*a(5,4))-a(3,4)*(a(4,2)*a(5,5)-a(4,5)*a(5,2))+a(3,5)*( &
 a(4,2)*a(5,4)-a(4,4)*a(5,2)))-a(2,2)*(a(3,1)*(a(4,4)*a(5,5)-a(4,5)* &
 a(5,4))-a(3,4)*(a(4,1)*a(5,5)-a(4,5)*a(5,1))+a(3,5)*(a(4,1)*a(5,4)- &
 a(4,4)*a(5,1)))+a(2,4)*(a(3,1)*(a(4,2)*a(5,5)-a(4,5)*a(5,2))-a(3,2)*( &
 a(4,1)*a(5,5)-a(4,5)*a(5,1))+a(3,5)*(a(4,1)*a(5,2)-a(4,2)*a(5,1)))- &
 a(2,5)*(a(3,1)*(a(4,2)*a(5,4)-a(4,4)*a(5,2))-a(3,2)*(a(4,1)*a(5,4)- &
 a(4,4)*a(5,1))+a(3,4)*(a(4,1)*a(5,2)-a(4,2)*a(5,1))))-a(1,4)*(a(2,1)*( &
 a(3,2)*(a(4,3)*a(5,5)-a(4,5)*a(5,3))-a(3,3)*(a(4,2)*a(5,5)-a(4,5)* &
 a(5,2))+a(3,5)*(a(4,2)*a(5,3)-a(4,3)*a(5,2)))-a(2,2)*(a(3,1)*(a(4,3)* &
 a(5,5)-a(4,5)*a(5,3))-a(3,3)*(a(4,1)*a(5,5)-a(4,5)*a(5,1))+a(3,5)*( &
 a(4,1)*a(5,3)-a(4,3)*a(5,1)))+a(2,3)*(a(3,1)*(a(4,2)*a(5,5)-a(4,5)* &
 a(5,2))-a(3,2)*(a(4,1)*a(5,5)-a(4,5)*a(5,1))+a(3,5)*(a(4,1)*a(5,2)- &
 a(4,2)*a(5,1)))-a(2,5)*(a(3,1)*(a(4,2)*a(5,3)-a(4,3)*a(5,2))-a(3,2)*( &
 a(4,1)*a(5,3)-a(4,3)*a(5,1))+a(3,3)*(a(4,1)*a(5,2)-a(4,2)*a(5,1))))+ &
 a(1,5)*(a(2,1)*(a(3,2)*(a(4,3)*a(5,4)-a(4,4)*a(5,3))-a(3,3)*(a(4,2)* &
 a(5,4)-a(4,4)*a(5,2))+a(3,4)*(a(4,2)*a(5,3)-a(4,3)*a(5,2)))-a(2,2)*( &
 a(3,1)*(a(4,3)*a(5,4)-a(4,4)*a(5,3))-a(3,3)*(a(4,1)*a(5,4)-a(4,4)* &
 a(5,1))+a(3,4)*(a(4,1)*a(5,3)-a(4,3)*a(5,1)))+a(2,3)*(a(3,1)*(a(4,2)* &
 a(5,4)-a(4,4)*a(5,2))-a(3,2)*(a(4,1)*a(5,4)-a(4,4)*a(5,1))+a(3,4)*( &
 a(4,1)*a(5,2)-a(4,2)*a(5,1)))-a(2,4)*(a(3,1)*(a(4,2)*a(5,3)-a(4,3)* &
 a(5,2))-a(3,2)*(a(4,1)*a(5,3)-a(4,3)*a(5,1))+a(3,3)*(a(4,1)*a(5,2)- &
 a(4,2)*a(5,1))))

 b(1,1) = &
 (a(2,2)*(a(3,3)*(a(4,4)*a(5,5)-a(4,5)*a(5,4))-a(3,4)*(a(4,3)*a(5,5)-a(4,5)*a(5,3))+a(3,5)*(a(4,3)*a(5,4)-a(4,4)*a(5,3)))-a(2,3)* &
 (a(3,2)*(a(4,4)*a(5,5)-a(4,5)*a(5,4))-a(3,4)*(a(4,2)*a(5,5)-a(4,5)*a(5,2))+a(3,5)*(a(4,2)*a(5,4)-a(4,4)*a(5,2)))+a(2,4)* &
 (a(3,2)*(a(4,3)*a(5,5)-a(4,5)*a(5,3))-a(3,3)*(a(4,2)*a(5,5)-a(4,5)*a(5,2))+a(3,5)*(a(4,2)*a(5,3)-a(4,3)*a(5,2)))-a(2,5)* &
 (a(3,2)*(a(4,3)*a(5,4)-a(4,4)*a(5,3))-a(3,3)*(a(4,2)*a(5,4)-a(4,4)*a(5,2))+a(3,4)*(a(4,2)*a(5,3)-a(4,3)*a(5,2))))
 b(2,1) = &
 (-a(2,1)*(a(3,3)*(a(4,4)*a(5,5)-a(4,5)*a(5,4))-a(3,4)*(a(4,3)*a(5,5)-a(4,5)*a(5,3))+a(3,5)*(a(4,3)*a(5,4)-a(4,4)*a(5,3)))+a(2,3)* &
 (a(3,1)*(a(4,4)*a(5,5)-a(4,5)*a(5,4))-a(3,4)*(a(4,1)*a(5,5)-a(4,5)*a(5,1))+a(3,5)*(a(4,1)*a(5,4)-a(4,4)*a(5,1)))-a(2,4)* &
 (a(3,1)*(a(4,3)*a(5,5)-a(4,5)*a(5,3))-a(3,3)*(a(4,1)*a(5,5)-a(4,5)*a(5,1))+a(3,5)*(a(4,1)*a(5,3)-a(4,3)*a(5,1)))+a(2,5)* &
 (a(3,1)*(a(4,3)*a(5,4)-a(4,4)*a(5,3))-a(3,3)*(a(4,1)*a(5,4)-a(4,4)*a(5,1))+a(3,4)*(a(4,1)*a(5,3)-a(4,3)*a(5,1))))
 b(3,1) = &
 (a(2,1)*(a(3,2)*(a(4,4)*a(5,5)-a(4,5)*a(5,4))-a(3,4)*(a(4,2)*a(5,5)-a(4,5)*a(5,2))+a(3,5)*(a(4,2)*a(5,4)-a(4,4)*a(5,2)))-a(2,2)* &
 (a(3,1)*(a(4,4)*a(5,5)-a(4,5)*a(5,4))-a(3,4)*(a(4,1)*a(5,5)-a(4,5)*a(5,1))+a(3,5)*(a(4,1)*a(5,4)-a(4,4)*a(5,1)))+a(2,4)* &
 (a(3,1)*(a(4,2)*a(5,5)-a(4,5)*a(5,2))-a(3,2)*(a(4,1)*a(5,5)-a(4,5)*a(5,1))+a(3,5)*(a(4,1)*a(5,2)-a(4,2)*a(5,1)))-a(2,5)* &
 (a(3,1)*(a(4,2)*a(5,4)-a(4,4)*a(5,2))-a(3,2)*(a(4,1)*a(5,4)-a(4,4)*a(5,1))+a(3,4)*(a(4,1)*a(5,2)-a(4,2)*a(5,1))))
 b(4,1) = &
 (-a(2,1)*(a(3,2)*(a(4,3)*a(5,5)-a(4,5)*a(5,3))-a(3,3)*(a(4,2)*a(5,5)-a(4,5)*a(5,2))+a(3,5)*(a(4,2)*a(5,3)-a(4,3)*a(5,2)))+a(2,2)* &
 (a(3,1)*(a(4,3)*a(5,5)-a(4,5)*a(5,3))-a(3,3)*(a(4,1)*a(5,5)-a(4,5)*a(5,1))+a(3,5)*(a(4,1)*a(5,3)-a(4,3)*a(5,1)))-a(2,3)* &
 (a(3,1)*(a(4,2)*a(5,5)-a(4,5)*a(5,2))-a(3,2)*(a(4,1)*a(5,5)-a(4,5)*a(5,1))+a(3,5)*(a(4,1)*a(5,2)-a(4,2)*a(5,1)))+a(2,5)* &
 (a(3,1)*(a(4,2)*a(5,3)-a(4,3)*a(5,2))-a(3,2)*(a(4,1)*a(5,3)-a(4,3)*a(5,1))+a(3,3)*(a(4,1)*a(5,2)-a(4,2)*a(5,1))))
 b(5,1) = &
 (a(2,1)*(a(3,2)*(a(4,3)*a(5,4)-a(4,4)*a(5,3))-a(3,3)*(a(4,2)*a(5,4)-a(4,4)*a(5,2))+a(3,4)*(a(4,2)*a(5,3)-a(4,3)*a(5,2)))-a(2,2)* &
 (a(3,1)*(a(4,3)*a(5,4)-a(4,4)*a(5,3))-a(3,3)*(a(4,1)*a(5,4)-a(4,4)*a(5,1))+a(3,4)*(a(4,1)*a(5,3)-a(4,3)*a(5,1)))+a(2,3)* &
 (a(3,1)*(a(4,2)*a(5,4)-a(4,4)*a(5,2))-a(3,2)*(a(4,1)*a(5,4)-a(4,4)*a(5,1))+a(3,4)*(a(4,1)*a(5,2)-a(4,2)*a(5,1)))-a(2,4)* &
 (a(3,1)*(a(4,2)*a(5,3)-a(4,3)*a(5,2))-a(3,2)*(a(4,1)*a(5,3)-a(4,3)*a(5,1))+a(3,3)*(a(4,1)*a(5,2)-a(4,2)*a(5,1))))

 b(1,2) = &
 (-a(1,2)*(a(3,3)*(a(4,4)*a(5,5)-a(4,5)*a(5,4))-a(3,4)*(a(4,3)*a(5,5)-a(4,5)*a(5,3))+a(3,5)*(a(4,3)*a(5,4)-a(4,4)*a(5,3)))+a(1,3)* &
 (a(3,2)*(a(4,4)*a(5,5)-a(4,5)*a(5,4))-a(3,4)*(a(4,2)*a(5,5)-a(4,5)*a(5,2))+a(3,5)*(a(4,2)*a(5,4)-a(4,4)*a(5,2)))-a(1,4)* &
 (a(3,2)*(a(4,3)*a(5,5)-a(4,5)*a(5,3))-a(3,3)*(a(4,2)*a(5,5)-a(4,5)*a(5,2))+a(3,5)*(a(4,2)*a(5,3)-a(4,3)*a(5,2)))+a(1,5)* &
 (a(3,2)*(a(4,3)*a(5,4)-a(4,4)*a(5,3))-a(3,3)*(a(4,2)*a(5,4)-a(4,4)*a(5,2))+a(3,4)*(a(4,2)*a(5,3)-a(4,3)*a(5,2))))
 b(2,2) = &
 (a(1,1)*(a(3,3)*(a(4,4)*a(5,5)-a(4,5)*a(5,4))-a(3,4)*(a(4,3)*a(5,5)-a(4,5)*a(5,3))+a(3,5)*(a(4,3)*a(5,4)-a(4,4)*a(5,3)))-a(1,3)* &
 (a(3,1)*(a(4,4)*a(5,5)-a(4,5)*a(5,4))-a(3,4)*(a(4,1)*a(5,5)-a(4,5)*a(5,1))+a(3,5)*(a(4,1)*a(5,4)-a(4,4)*a(5,1)))+a(1,4)* &
 (a(3,1)*(a(4,3)*a(5,5)-a(4,5)*a(5,3))-a(3,3)*(a(4,1)*a(5,5)-a(4,5)*a(5,1))+a(3,5)*(a(4,1)*a(5,3)-a(4,3)*a(5,1)))-a(1,5)* &
 (a(3,1)*(a(4,3)*a(5,4)-a(4,4)*a(5,3))-a(3,3)*(a(4,1)*a(5,4)-a(4,4)*a(5,1))+a(3,4)*(a(4,1)*a(5,3)-a(4,3)*a(5,1))))
 b(3,2) = &
 (-a(1,1)*(a(3,2)*(a(4,4)*a(5,5)-a(4,5)*a(5,4))-a(3,4)*(a(4,2)*a(5,5)-a(4,5)*a(5,2))+a(3,5)*(a(4,2)*a(5,4)-a(4,4)*a(5,2)))+a(1,2)* &
 (a(3,1)*(a(4,4)*a(5,5)-a(4,5)*a(5,4))-a(3,4)*(a(4,1)*a(5,5)-a(4,5)*a(5,1))+a(3,5)*(a(4,1)*a(5,4)-a(4,4)*a(5,1)))-a(1,4)* &
 (a(3,1)*(a(4,2)*a(5,5)-a(4,5)*a(5,2))-a(3,2)*(a(4,1)*a(5,5)-a(4,5)*a(5,1))+a(3,5)*(a(4,1)*a(5,2)-a(4,2)*a(5,1)))+a(1,5)* &
 (a(3,1)*(a(4,2)*a(5,4)-a(4,4)*a(5,2))-a(3,2)*(a(4,1)*a(5,4)-a(4,4)*a(5,1))+a(3,4)*(a(4,1)*a(5,2)-a(4,2)*a(5,1))))
 b(4,2) = &
 (a(1,1)*(a(3,2)*(a(4,3)*a(5,5)-a(4,5)*a(5,3))-a(3,3)*(a(4,2)*a(5,5)-a(4,5)*a(5,2))+a(3,5)*(a(4,2)*a(5,3)-a(4,3)*a(5,2)))-a(1,2)* &
 (a(3,1)*(a(4,3)*a(5,5)-a(4,5)*a(5,3))-a(3,3)*(a(4,1)*a(5,5)-a(4,5)*a(5,1))+a(3,5)*(a(4,1)*a(5,3)-a(4,3)*a(5,1)))+a(1,3)* &
 (a(3,1)*(a(4,2)*a(5,5)-a(4,5)*a(5,2))-a(3,2)*(a(4,1)*a(5,5)-a(4,5)*a(5,1))+a(3,5)*(a(4,1)*a(5,2)-a(4,2)*a(5,1)))-a(1,5)* &
 (a(3,1)*(a(4,2)*a(5,3)-a(4,3)*a(5,2))-a(3,2)*(a(4,1)*a(5,3)-a(4,3)*a(5,1))+a(3,3)*(a(4,1)*a(5,2)-a(4,2)*a(5,1))))
 b(5,2) = &
 (-a(1,1)*(a(3,2)*(a(4,3)*a(5,4)-a(4,4)*a(5,3))-a(3,3)*(a(4,2)*a(5,4)-a(4,4)*a(5,2))+a(3,4)*(a(4,2)*a(5,3)-a(4,3)*a(5,2)))+a(1,2)* &
 (a(3,1)*(a(4,3)*a(5,4)-a(4,4)*a(5,3))-a(3,3)*(a(4,1)*a(5,4)-a(4,4)*a(5,1))+a(3,4)*(a(4,1)*a(5,3)-a(4,3)*a(5,1)))-a(1,3)* &
 (a(3,1)*(a(4,2)*a(5,4)-a(4,4)*a(5,2))-a(3,2)*(a(4,1)*a(5,4)-a(4,4)*a(5,1))+a(3,4)*(a(4,1)*a(5,2)-a(4,2)*a(5,1)))+a(1,4)* &
 (a(3,1)*(a(4,2)*a(5,3)-a(4,3)*a(5,2))-a(3,2)*(a(4,1)*a(5,3)-a(4,3)*a(5,1))+a(3,3)*(a(4,1)*a(5,2)-a(4,2)*a(5,1))))

 b(1,3) = &
 (a(1,2)*(a(2,3)*(a(4,4)*a(5,5)-a(4,5)*a(5,4))-a(2,4)*(a(4,3)*a(5,5)-a(4,5)*a(5,3))+a(2,5)*(a(4,3)*a(5,4)-a(4,4)*a(5,3)))-a(1,3)* &
 (a(2,2)*(a(4,4)*a(5,5)-a(4,5)*a(5,4))-a(2,4)*(a(4,2)*a(5,5)-a(4,5)*a(5,2))+a(2,5)*(a(4,2)*a(5,4)-a(4,4)*a(5,2)))+a(1,4)* &
 (a(2,2)*(a(4,3)*a(5,5)-a(4,5)*a(5,3))-a(2,3)*(a(4,2)*a(5,5)-a(4,5)*a(5,2))+a(2,5)*(a(4,2)*a(5,3)-a(4,3)*a(5,2)))-a(1,5)* &
 (a(2,2)*(a(4,3)*a(5,4)-a(4,4)*a(5,3))-a(2,3)*(a(4,2)*a(5,4)-a(4,4)*a(5,2))+a(2,4)*(a(4,2)*a(5,3)-a(4,3)*a(5,2))))
 b(2,3) = &
 (-a(1,1)*(a(2,3)*(a(4,4)*a(5,5)-a(4,5)*a(5,4))-a(2,4)*(a(4,3)*a(5,5)-a(4,5)*a(5,3))+a(2,5)*(a(4,3)*a(5,4)-a(4,4)*a(5,3)))+a(1,3)* &
 (a(2,1)*(a(4,4)*a(5,5)-a(4,5)*a(5,4))-a(2,4)*(a(4,1)*a(5,5)-a(4,5)*a(5,1))+a(2,5)*(a(4,1)*a(5,4)-a(4,4)*a(5,1)))-a(1,4)* &
 (a(2,1)*(a(4,3)*a(5,5)-a(4,5)*a(5,3))-a(2,3)*(a(4,1)*a(5,5)-a(4,5)*a(5,1))+a(2,5)*(a(4,1)*a(5,3)-a(4,3)*a(5,1)))+a(1,5)* &
 (a(2,1)*(a(4,3)*a(5,4)-a(4,4)*a(5,3))-a(2,3)*(a(4,1)*a(5,4)-a(4,4)*a(5,1))+a(2,4)*(a(4,1)*a(5,3)-a(4,3)*a(5,1))))
 b(3,3) = &
 (a(1,1)*(a(2,2)*(a(4,4)*a(5,5)-a(4,5)*a(5,4))-a(2,4)*(a(4,2)*a(5,5)-a(4,5)*a(5,2))+a(2,5)*(a(4,2)*a(5,4)-a(4,4)*a(5,2)))-a(1,2)* &
 (a(2,1)*(a(4,4)*a(5,5)-a(4,5)*a(5,4))-a(2,4)*(a(4,1)*a(5,5)-a(4,5)*a(5,1))+a(2,5)*(a(4,1)*a(5,4)-a(4,4)*a(5,1)))+a(1,4)* &
 (a(2,1)*(a(4,2)*a(5,5)-a(4,5)*a(5,2))-a(2,2)*(a(4,1)*a(5,5)-a(4,5)*a(5,1))+a(2,5)*(a(4,1)*a(5,2)-a(4,2)*a(5,1)))-a(1,5)* &
 (a(2,1)*(a(4,2)*a(5,4)-a(4,4)*a(5,2))-a(2,2)*(a(4,1)*a(5,4)-a(4,4)*a(5,1))+a(2,4)*(a(4,1)*a(5,2)-a(4,2)*a(5,1))))
 b(4,3) = &
 (-a(1,1)*(a(2,2)*(a(4,3)*a(5,5)-a(4,5)*a(5,3))-a(2,3)*(a(4,2)*a(5,5)-a(4,5)*a(5,2))+a(2,5)*(a(4,2)*a(5,3)-a(4,3)*a(5,2)))+a(1,2)* &
 (a(2,1)*(a(4,3)*a(5,5)-a(4,5)*a(5,3))-a(2,3)*(a(4,1)*a(5,5)-a(4,5)*a(5,1))+a(2,5)*(a(4,1)*a(5,3)-a(4,3)*a(5,1)))-a(1,3)* &
 (a(2,1)*(a(4,2)*a(5,5)-a(4,5)*a(5,2))-a(2,2)*(a(4,1)*a(5,5)-a(4,5)*a(5,1))+a(2,5)*(a(4,1)*a(5,2)-a(4,2)*a(5,1)))+a(1,5)* &
 (a(2,1)*(a(4,2)*a(5,3)-a(4,3)*a(5,2))-a(2,2)*(a(4,1)*a(5,3)-a(4,3)*a(5,1))+a(2,3)*(a(4,1)*a(5,2)-a(4,2)*a(5,1))))
 b(5,3) = &
 (a(1,1)*(a(2,2)*(a(4,3)*a(5,4)-a(4,4)*a(5,3))-a(2,3)*(a(4,2)*a(5,4)-a(4,4)*a(5,2))+a(2,4)*(a(4,2)*a(5,3)-a(4,3)*a(5,2)))-a(1,2)* &
 (a(2,1)*(a(4,3)*a(5,4)-a(4,4)*a(5,3))-a(2,3)*(a(4,1)*a(5,4)-a(4,4)*a(5,1))+a(2,4)*(a(4,1)*a(5,3)-a(4,3)*a(5,1)))+a(1,3)* &
 (a(2,1)*(a(4,2)*a(5,4)-a(4,4)*a(5,2))-a(2,2)*(a(4,1)*a(5,4)-a(4,4)*a(5,1))+a(2,4)*(a(4,1)*a(5,2)-a(4,2)*a(5,1)))-a(1,4)* &
 (a(2,1)*(a(4,2)*a(5,3)-a(4,3)*a(5,2))-a(2,2)*(a(4,1)*a(5,3)-a(4,3)*a(5,1))+a(2,3)*(a(4,1)*a(5,2)-a(4,2)*a(5,1))))

 b(1,4) = &
 (-a(1,2)*(a(2,3)*(a(3,4)*a(5,5)-a(3,5)*a(5,4))-a(2,4)*(a(3,3)*a(5,5)-a(3,5)*a(5,3))+a(2,5)*(a(3,3)*a(5,4)-a(3,4)*a(5,3)))+a(1,3)* &
 (a(2,2)*(a(3,4)*a(5,5)-a(3,5)*a(5,4))-a(2,4)*(a(3,2)*a(5,5)-a(3,5)*a(5,2))+a(2,5)*(a(3,2)*a(5,4)-a(3,4)*a(5,2)))-a(1,4)* &
 (a(2,2)*(a(3,3)*a(5,5)-a(3,5)*a(5,3))-a(2,3)*(a(3,2)*a(5,5)-a(3,5)*a(5,2))+a(2,5)*(a(3,2)*a(5,3)-a(3,3)*a(5,2)))+a(1,5)* &
 (a(2,2)*(a(3,3)*a(5,4)-a(3,4)*a(5,3))-a(2,3)*(a(3,2)*a(5,4)-a(3,4)*a(5,2))+a(2,4)*(a(3,2)*a(5,3)-a(3,3)*a(5,2))))
 b(2,4) = &
 (a(1,1)*(a(2,3)*(a(3,4)*a(5,5)-a(3,5)*a(5,4))-a(2,4)*(a(3,3)*a(5,5)-a(3,5)*a(5,3))+a(2,5)*(a(3,3)*a(5,4)-a(3,4)*a(5,3)))-a(1,3)* &
 (a(2,1)*(a(3,4)*a(5,5)-a(3,5)*a(5,4))-a(2,4)*(a(3,1)*a(5,5)-a(3,5)*a(5,1))+a(2,5)*(a(3,1)*a(5,4)-a(3,4)*a(5,1)))+a(1,4)* &
 (a(2,1)*(a(3,3)*a(5,5)-a(3,5)*a(5,3))-a(2,3)*(a(3,1)*a(5,5)-a(3,5)*a(5,1))+a(2,5)*(a(3,1)*a(5,3)-a(3,3)*a(5,1)))-a(1,5)* &
 (a(2,1)*(a(3,3)*a(5,4)-a(3,4)*a(5,3))-a(2,3)*(a(3,1)*a(5,4)-a(3,4)*a(5,1))+a(2,4)*(a(3,1)*a(5,3)-a(3,3)*a(5,1))))
 b(3,4) = &
 (-a(1,1)*(a(2,2)*(a(3,4)*a(5,5)-a(3,5)*a(5,4))-a(2,4)*(a(3,2)*a(5,5)-a(3,5)*a(5,2))+a(2,5)*(a(3,2)*a(5,4)-a(3,4)*a(5,2)))+a(1,2)* &
 (a(2,1)*(a(3,4)*a(5,5)-a(3,5)*a(5,4))-a(2,4)*(a(3,1)*a(5,5)-a(3,5)*a(5,1))+a(2,5)*(a(3,1)*a(5,4)-a(3,4)*a(5,1)))-a(1,4)* &
 (a(2,1)*(a(3,2)*a(5,5)-a(3,5)*a(5,2))-a(2,2)*(a(3,1)*a(5,5)-a(3,5)*a(5,1))+a(2,5)*(a(3,1)*a(5,2)-a(3,2)*a(5,1)))+a(1,5)* &
 (a(2,1)*(a(3,2)*a(5,4)-a(3,4)*a(5,2))-a(2,2)*(a(3,1)*a(5,4)-a(3,4)*a(5,1))+a(2,4)*(a(3,1)*a(5,2)-a(3,2)*a(5,1))))
 b(4,4) = &
 (a(1,1)*(a(2,2)*(a(3,3)*a(5,5)-a(3,5)*a(5,3))-a(2,3)*(a(3,2)*a(5,5)-a(3,5)*a(5,2))+a(2,5)*(a(3,2)*a(5,3)-a(3,3)*a(5,2)))-a(1,2)* &
 (a(2,1)*(a(3,3)*a(5,5)-a(3,5)*a(5,3))-a(2,3)*(a(3,1)*a(5,5)-a(3,5)*a(5,1))+a(2,5)*(a(3,1)*a(5,3)-a(3,3)*a(5,1)))+a(1,3)* &
 (a(2,1)*(a(3,2)*a(5,5)-a(3,5)*a(5,2))-a(2,2)*(a(3,1)*a(5,5)-a(3,5)*a(5,1))+a(2,5)*(a(3,1)*a(5,2)-a(3,2)*a(5,1)))-a(1,5)* &
 (a(2,1)*(a(3,2)*a(5,3)-a(3,3)*a(5,2))-a(2,2)*(a(3,1)*a(5,3)-a(3,3)*a(5,1))+a(2,3)*(a(3,1)*a(5,2)-a(3,2)*a(5,1))))
 b(5,4) = &
 (-a(1,1)*(a(2,2)*(a(3,3)*a(5,4)-a(3,4)*a(5,3))-a(2,3)*(a(3,2)*a(5,4)-a(3,4)*a(5,2))+a(2,4)*(a(3,2)*a(5,3)-a(3,3)*a(5,2)))+a(1,2)* &
 (a(2,1)*(a(3,3)*a(5,4)-a(3,4)*a(5,3))-a(2,3)*(a(3,1)*a(5,4)-a(3,4)*a(5,1))+a(2,4)*(a(3,1)*a(5,3)-a(3,3)*a(5,1)))-a(1,3)* &
 (a(2,1)*(a(3,2)*a(5,4)-a(3,4)*a(5,2))-a(2,2)*(a(3,1)*a(5,4)-a(3,4)*a(5,1))+a(2,4)*(a(3,1)*a(5,2)-a(3,2)*a(5,1)))+a(1,4)* &
 (a(2,1)*(a(3,2)*a(5,3)-a(3,3)*a(5,2))-a(2,2)*(a(3,1)*a(5,3)-a(3,3)*a(5,1))+a(2,3)*(a(3,1)*a(5,2)-a(3,2)*a(5,1))))

 b(1,5) = &
 (a(1,2)*(a(2,3)*(a(3,4)*a(4,5)-a(3,5)*a(4,4))-a(2,4)*(a(3,3)*a(4,5)-a(3,5)*a(4,3))+a(2,5)*(a(3,3)*a(4,4)-a(3,4)*a(4,3)))-a(1,3)* &
 (a(2,2)*(a(3,4)*a(4,5)-a(3,5)*a(4,4))-a(2,4)*(a(3,2)*a(4,5)-a(3,5)*a(4,2))+a(2,5)*(a(3,2)*a(4,4)-a(3,4)*a(4,2)))+a(1,4)* &
 (a(2,2)*(a(3,3)*a(4,5)-a(3,5)*a(4,3))-a(2,3)*(a(3,2)*a(4,5)-a(3,5)*a(4,2))+a(2,5)*(a(3,2)*a(4,3)-a(3,3)*a(4,2)))-a(1,5)* &
 (a(2,2)*(a(3,3)*a(4,4)-a(3,4)*a(4,3))-a(2,3)*(a(3,2)*a(4,4)-a(3,4)*a(4,2))+a(2,4)*(a(3,2)*a(4,3)-a(3,3)*a(4,2))))
 b(2,5) = &
 (-a(1,1)*(a(2,3)*(a(3,4)*a(4,5)-a(3,5)*a(4,4))-a(2,4)*(a(3,3)*a(4,5)-a(3,5)*a(4,3))+a(2,5)*(a(3,3)*a(4,4)-a(3,4)*a(4,3)))+a(1,3)* &
 (a(2,1)*(a(3,4)*a(4,5)-a(3,5)*a(4,4))-a(2,4)*(a(3,1)*a(4,5)-a(3,5)*a(4,1))+a(2,5)*(a(3,1)*a(4,4)-a(3,4)*a(4,1)))-a(1,4)* &
 (a(2,1)*(a(3,3)*a(4,5)-a(3,5)*a(4,3))-a(2,3)*(a(3,1)*a(4,5)-a(3,5)*a(4,1))+a(2,5)*(a(3,1)*a(4,3)-a(3,3)*a(4,1)))+a(1,5)* &
 (a(2,1)*(a(3,3)*a(4,4)-a(3,4)*a(4,3))-a(2,3)*(a(3,1)*a(4,4)-a(3,4)*a(4,1))+a(2,4)*(a(3,1)*a(4,3)-a(3,3)*a(4,1))))
 b(3,5) = &
 (a(1,1)*(a(2,2)*(a(3,4)*a(4,5)-a(3,5)*a(4,4))-a(2,4)*(a(3,2)*a(4,5)-a(3,5)*a(4,2))+a(2,5)*(a(3,2)*a(4,4)-a(3,4)*a(4,2)))-a(1,2)* &
 (a(2,1)*(a(3,4)*a(4,5)-a(3,5)*a(4,4))-a(2,4)*(a(3,1)*a(4,5)-a(3,5)*a(4,1))+a(2,5)*(a(3,1)*a(4,4)-a(3,4)*a(4,1)))+a(1,4)* &
 (a(2,1)*(a(3,2)*a(4,5)-a(3,5)*a(4,2))-a(2,2)*(a(3,1)*a(4,5)-a(3,5)*a(4,1))+a(2,5)*(a(3,1)*a(4,2)-a(3,2)*a(4,1)))-a(1,5)* &
 (a(2,1)*(a(3,2)*a(4,4)-a(3,4)*a(4,2))-a(2,2)*(a(3,1)*a(4,4)-a(3,4)*a(4,1))+a(2,4)*(a(3,1)*a(4,2)-a(3,2)*a(4,1))))
 b(4,5) = &
 (-a(1,1)*(a(2,2)*(a(3,3)*a(4,5)-a(3,5)*a(4,3))-a(2,3)*(a(3,2)*a(4,5)-a(3,5)*a(4,2))+a(2,5)*(a(3,2)*a(4,3)-a(3,3)*a(4,2)))+a(1,2)* &
 (a(2,1)*(a(3,3)*a(4,5)-a(3,5)*a(4,3))-a(2,3)*(a(3,1)*a(4,5)-a(3,5)*a(4,1))+a(2,5)*(a(3,1)*a(4,3)-a(3,3)*a(4,1)))-a(1,3)* &
 (a(2,1)*(a(3,2)*a(4,5)-a(3,5)*a(4,2))-a(2,2)*(a(3,1)*a(4,5)-a(3,5)*a(4,1))+a(2,5)*(a(3,1)*a(4,2)-a(3,2)*a(4,1)))+a(1,5)* &
 (a(2,1)*(a(3,2)*a(4,3)-a(3,3)*a(4,2))-a(2,2)*(a(3,1)*a(4,3)-a(3,3)*a(4,1))+a(2,3)*(a(3,1)*a(4,2)-a(3,2)*a(4,1))))
 b(5,5) = &
 (a(1,1)*(a(2,2)*(a(3,3)*a(4,4)-a(3,4)*a(4,3))-a(2,3)*(a(3,2)*a(4,4)-a(3,4)*a(4,2))+a(2,4)*(a(3,2)*a(4,3)-a(3,3)*a(4,2)))-a(1,2)* &
 (a(2,1)*(a(3,3)*a(4,4)-a(3,4)*a(4,3))-a(2,3)*(a(3,1)*a(4,4)-a(3,4)*a(4,1))+a(2,4)*(a(3,1)*a(4,3)-a(3,3)*a(4,1)))+a(1,3)* &
 (a(2,1)*(a(3,2)*a(4,4)-a(3,4)*a(4,2))-a(2,2)*(a(3,1)*a(4,4)-a(3,4)*a(4,1))+a(2,4)*(a(3,1)*a(4,2)-a(3,2)*a(4,1)))-a(1,4)* &
 (a(2,1)*(a(3,2)*a(4,3)-a(3,3)*a(4,2))-a(2,2)*(a(3,1)*a(4,3)-a(3,3)*a(4,1))+a(2,3)*(a(3,1)*a(4,2)-a(3,2)*a(4,1))))

end
    #+end_src

    #+CALL: generate_c_interface(table=qmckl_adjugate_args,rettyp="qmckl_exit_code",fname="qmckl_adjugate")

    #+RESULTS:
    #+begin_src f90 :tangle (eval f) :comments org :exports none
    integer(c_int32_t) function qmckl_adjugate &
        (context, n, A, lda, B, ldb, det_l) &
        bind(C) result(info)

      use, intrinsic :: iso_c_binding
      implicit none

      integer (c_int64_t) , intent(in)  , value :: context
      integer (c_int64_t) , intent(in)  , value :: n
      real    (c_double ) , intent(in)          :: A(lda,*)
      integer (c_int64_t) , intent(in)  , value :: lda
      real    (c_double ) , intent(out)         :: B(ldb,*)
      integer (c_int64_t) , intent(in)  , value :: ldb
      real    (c_double ) , intent(inout)        :: det_l

      integer(c_int32_t), external :: qmckl_adjugate_f
      info = qmckl_adjugate_f &
             (context, n, A, lda, B, ldb, det_l)

    end function qmckl_adjugate
    #+end_src

    #+CALL: generate_f_interface(table=qmckl_adjugate_args,rettyp="qmckl_exit_code",fname="qmckl_adjugate")

    #+RESULTS:
    #+begin_src f90 :tangle (eval fh_func) :comments org :exports none
    interface
      integer(c_int32_t) function qmckl_adjugate &
          (context, n, A, lda, B, ldb, det_l) &
          bind(C)
        use, intrinsic :: iso_c_binding
        import
        implicit none

        integer (c_int64_t) , intent(in)  , value :: context
        integer (c_int64_t) , intent(in)  , value :: n
        real    (c_double ) , intent(in)          :: A(lda,*)
        integer (c_int64_t) , intent(in)  , value :: lda
        real    (c_double ) , intent(out)         :: B(ldb,*)
        integer (c_int64_t) , intent(in)  , value :: ldb
        real    (c_double ) , intent(inout)        :: det_l

      end function qmckl_adjugate
    end interface
    #+end_src


    #+begin_src f90 :tangle (eval f)
subroutine adjugate_general(context, na, A, LDA, B, LDB, det_l)
  use qmckl
  implicit none
  integer(qmckl_context), intent(in) :: context
  double precision, intent(in)       :: A (LDA,na)
  integer*8, intent(in)              :: LDA
  double precision, intent(out)      :: B (LDB,na)
  integer*8, intent(in)              :: LDB
  integer*8, intent(in)              :: na
  double precision, intent(inout)    :: det_l

  double precision :: work(LDA*max(na,64))
  integer          :: inf
  integer          :: ipiv(LDA)
  integer          :: lwork
  integer(8)       :: i, j

    #+end_src

    We first copy the array ~A~ into array ~B~.

    #+begin_src f90 :tangle (eval f)
  B(1:na,1:na) = A(1:na,1:na)
    #+end_src

    Then, we compute the LU factorization $LU=A$
    using the ~dgetrf~ routine.

    #+begin_src f90 :tangle (eval f)
  call dgetrf(na, na, B, LDB, ipiv, inf )
    #+end_src

    By convention, the determinant of $L$ is equal to one, so the
    determinant of $A$ is equal to the determinant of $U$, which is
    simply computed as the product of its diagonal elements.

    #+begin_src f90 :tangle (eval f)
  det_l = 1.d0
  j=0_8
  do i=1,na
   j = j+min(abs(ipiv(i)-i),1)
   det_l = det_l*B(i,i)
  enddo
    #+end_src

    As ~dgetrf~ returns $PLU=A$ where $P$ is a permutation matrix, the
    sign of the determinant is computed as $-1^m$ where $m$ is the
    number of permutations.

    #+begin_src f90 :tangle (eval f)
  if (iand(j,1_8) /= 0_8)  then
    det_l = -det_l
  endif
    #+end_src

    Then, the inverse of $A$ is computed using ~dgetri~:

    #+begin_src f90 :tangle (eval f)
  lwork = SIZE(work)
  call dgetri(na, B, LDB, ipiv, work, lwork, inf )
    #+end_src

    and the adjugate matrix is computed as the product of the
    determinant with the inverse:

    #+begin_src f90 :tangle (eval f)
  B(:,:) = B(:,:)*det_l

end subroutine adjugate_general
    #+end_src

*** Test                                                           :noexport:

    #+begin_src f90 :tangle (eval f_test)
integer(qmckl_exit_code) function test_qmckl_adjugate(context) bind(C)
  use qmckl
  implicit none
  integer(qmckl_context), intent(in), value :: context

  double precision, allocatable :: A(:,:), B(:,:)
  integer*8                     :: m, n, k, LDA, LDB
  integer*8                     :: i,j,l
  double precision              :: x, det_l, det_l_ref

  LDA = 6_8
  LDB = 6_8

  allocate( A(LDA,6), B(LDB,6))

  A = 0.1d0
  A(1,1) = 1.0d0;
  A(2,2) = 2.0d0;
  A(3,3) = 3.0d0;
  A(4,4) = 4.0d0;
  A(5,5) = 5.0d0;
  A(6,6) = 6.0d0;

  test_qmckl_adjugate = QMCKL_SUCCESS

    #+end_src

    #+begin_src python :results output :output drawer
import numpy as np
import numpy.linalg as la
N=6

A = np.zeros( (N,N) )
A += 0.1
for i in range(N):
    A[i][i] = i+1

def adj(A):
    return la.det(A) * la.inv(A)

print ("#+begin_src f90 :tangle (eval f_test)")

for i in range(N):
    print ("! N = ", i+1)
    print (f"  test_qmckl_adjugate = qmckl_adjugate(context, {i+1}_8, A, LDA, B, LDB, det_l)")
    print (f"  if (test_qmckl_adjugate /= QMCKL_SUCCESS) return")
    print (f"  if (dabs((det_l - ({la.det(A[0:i+1,0:i+1])}d0))/det_l) > 1.d-13) then")
    print (f"      print *, 'N = {i+1}: det = ', det_l, {la.det(A[0:i+1,0:i+1])}d0")
    print (f"      test_qmckl_adjugate = {i+1}")
    print (f"      return")
    print (f"  end if")
    print (f"")
    print (f"  x = 0.d0")
    for j in range(i+1):
        C = adj(A[0:i+1,0:i+1])
        for k in range(i+1):
            print (f"  x = x + (B({j+1},{k+1}) - ({C[j,k]}) )**2")
    print (f"  if (dabs(x / det_l) > 1.d-12) then")
    print (f"      print *, 'N = {i+1}: x = ', x")
    print (f"      test_qmckl_adjugate = {i+1}")
    print (f"      return")
    print (f"  end if")
    print (f"")


print ("#+end_src")
#    print(adj(A[0:i+1,0:i+1]))
    #+end_src

    #+RESULTS:
    #+begin_example
    ,#+begin_src f90 :tangle (eval f_test)
    ! N =  1
      test_qmckl_adjugate = qmckl_adjugate(context, 1_8, A, LDA, B, LDB, det_l)
      if (test_qmckl_adjugate /= QMCKL_SUCCESS) return
      if (dabs((det_l - (1.0d0))/det_l) > 1.d-13) then
          print *, 'N = 1: det = ', det_l, 1.0d0
          test_qmckl_adjugate = 1
          return
      end if

      x = 0.d0
      x = x + (B(1,1) - (1.0) )**2
      if (dabs(x / det_l) > 1.d-12) then
          print *, 'N = 1: x = ', x
          test_qmckl_adjugate = 1
          return
      end if

    ! N =  2
      test_qmckl_adjugate = qmckl_adjugate(context, 2_8, A, LDA, B, LDB, det_l)
      if (test_qmckl_adjugate /= QMCKL_SUCCESS) return
      if (dabs((det_l - (1.99d0))/det_l) > 1.d-13) then
          print *, 'N = 2: det = ', det_l, 1.99d0
          test_qmckl_adjugate = 2
          return
      end if

      x = 0.d0
      x = x + (B(1,1) - (1.9999999999999998) )**2
      x = x + (B(1,2) - (-0.09999999999999999) )**2
      x = x + (B(2,1) - (-0.09999999999999999) )**2
      x = x + (B(2,2) - (0.9999999999999999) )**2
      if (dabs(x / det_l) > 1.d-12) then
          print *, 'N = 2: x = ', x
          test_qmckl_adjugate = 2
          return
      end if

    ! N =  3
      test_qmckl_adjugate = qmckl_adjugate(context, 3_8, A, LDA, B, LDB, det_l)
      if (test_qmckl_adjugate /= QMCKL_SUCCESS) return
      if (dabs((det_l - (5.942000000000001d0))/det_l) > 1.d-13) then
          print *, 'N = 3: det = ', det_l, 5.942000000000001d0
          test_qmckl_adjugate = 3
          return
      end if

      x = 0.d0
      x = x + (B(1,1) - (5.990000000000001) )**2
      x = x + (B(1,2) - (-0.29000000000000004) )**2
      x = x + (B(1,3) - (-0.19000000000000003) )**2
      x = x + (B(2,1) - (-0.29000000000000004) )**2
      x = x + (B(2,2) - (2.9900000000000007) )**2
      x = x + (B(2,3) - (-0.09000000000000001) )**2
      x = x + (B(3,1) - (-0.19000000000000003) )**2
      x = x + (B(3,2) - (-0.09) )**2
      x = x + (B(3,3) - (1.9900000000000002) )**2
      if (dabs(x / det_l) > 1.d-12) then
          print *, 'N = 3: x = ', x
          test_qmckl_adjugate = 3
          return
      end if

    ! N =  4
      test_qmckl_adjugate = qmckl_adjugate(context, 4_8, A, LDA, B, LDB, det_l)
      if (test_qmckl_adjugate /= QMCKL_SUCCESS) return
      if (dabs((det_l - (23.669700000000006d0))/det_l) > 1.d-13) then
          print *, 'N = 4: det = ', det_l, 23.669700000000006d0
          test_qmckl_adjugate = 4
          return
      end if

      x = 0.d0
      x = x + (B(1,1) - (23.91200000000001) )**2
      x = x + (B(1,2) - (-1.1310000000000004) )**2
      x = x + (B(1,3) - (-0.7410000000000001) )**2
      x = x + (B(1,4) - (-0.5510000000000002) )**2
      x = x + (B(2,1) - (-1.1310000000000002) )**2
      x = x + (B(2,2) - (11.922000000000002) )**2
      x = x + (B(2,3) - (-0.351) )**2
      x = x + (B(2,4) - (-0.261) )**2
      x = x + (B(3,1) - (-0.7410000000000002) )**2
      x = x + (B(3,2) - (-0.351) )**2
      x = x + (B(3,3) - (7.932000000000001) )**2
      x = x + (B(3,4) - (-0.17100000000000004) )**2
      x = x + (B(4,1) - (-0.5510000000000002) )**2
      x = x + (B(4,2) - (-0.261) )**2
      x = x + (B(4,3) - (-0.17100000000000004) )**2
      x = x + (B(4,4) - (5.942000000000001) )**2
      if (dabs(x / det_l) > 1.d-12) then
          print *, 'N = 4: x = ', x
          test_qmckl_adjugate = 4
          return
      end if

    ! N =  5
      test_qmckl_adjugate = qmckl_adjugate(context, 5_8, A, LDA, B, LDB, det_l)
      if (test_qmckl_adjugate /= QMCKL_SUCCESS) return
      if (dabs((det_l - (117.91554000000008d0))/det_l) > 1.d-13) then
          print *, 'N = 5: det = ', det_l, 117.91554000000008d0
          test_qmckl_adjugate = 5
          return
      end if

      x = 0.d0
      x = x + (B(1,1) - (119.31770000000006) )**2
      x = x + (B(1,2) - (-5.541900000000004) )**2
      x = x + (B(1,3) - (-3.6309000000000022) )**2
      x = x + (B(1,4) - (-2.6999000000000017) )**2
      x = x + (B(1,5) - (-2.1489000000000016) )**2
      x = x + (B(2,1) - (-5.541900000000004) )**2
      x = x + (B(2,2) - (59.435700000000026) )**2
      x = x + (B(2,3) - (-1.7199000000000007) )**2
      x = x + (B(2,4) - (-1.2789000000000006) )**2
      x = x + (B(2,5) - (-1.0179000000000005) )**2
      x = x + (B(3,1) - (-3.6309000000000027) )**2
      x = x + (B(3,2) - (-1.7199000000000007) )**2
      x = x + (B(3,3) - (39.53370000000002) )**2
      x = x + (B(3,4) - (-0.8379000000000005) )**2
      x = x + (B(3,5) - (-0.6669000000000004) )**2
      x = x + (B(4,1) - (-2.699900000000002) )**2
      x = x + (B(4,2) - (-1.2789000000000006) )**2
      x = x + (B(4,3) - (-0.8379000000000004) )**2
      x = x + (B(4,4) - (29.611700000000017) )**2
      x = x + (B(4,5) - (-0.4959000000000003) )**2
      x = x + (B(5,1) - (-2.1489000000000016) )**2
      x = x + (B(5,2) - (-1.0179000000000005) )**2
      x = x + (B(5,3) - (-0.6669000000000004) )**2
      x = x + (B(5,4) - (-0.4959000000000003) )**2
      x = x + (B(5,5) - (23.669700000000013) )**2
      if (dabs(x / det_l) > 1.d-12) then
          print *, 'N = 5: x = ', x
          test_qmckl_adjugate = 5
          return
      end if

    ! N =  6
      test_qmckl_adjugate = qmckl_adjugate(context, 6_8, A, LDA, B, LDB, det_l)
      if (test_qmckl_adjugate /= QMCKL_SUCCESS) return
      if (dabs((det_l - (705.1783350000001d0))/det_l) > 1.d-13) then
          print *, 'N = 6: det = ', det_l, 705.1783350000001d0
          test_qmckl_adjugate = 6
          return
      end if

      x = 0.d0
      x = x + (B(1,1) - (714.5040400000001) )**2
      x = x + (B(1,2) - (-32.697210000000005) )**2
      x = x + (B(1,3) - (-21.422310000000003) )**2
      x = x + (B(1,4) - (-15.929410000000006) )**2
      x = x + (B(1,5) - (-12.678510000000003) )**2
      x = x + (B(1,6) - (-10.529610000000003) )**2
      x = x + (B(2,1) - (-32.69721) )**2
      x = x + (B(2,2) - (355.65834) )**2
      x = x + (B(2,3) - (-10.147409999999997) )**2
      x = x + (B(2,4) - (-7.54551) )**2
      x = x + (B(2,5) - (-6.005610000000001) )**2
      x = x + (B(2,6) - (-4.987709999999999) )**2
      x = x + (B(3,1) - (-21.422310000000003) )**2
      x = x + (B(3,2) - (-10.147409999999999) )**2
      x = x + (B(3,3) - (236.51663999999997) )**2
      x = x + (B(3,4) - (-4.943610000000001) )**2
      x = x + (B(3,5) - (-3.93471) )**2
      x = x + (B(3,6) - (-3.267810000000001) )**2
      x = x + (B(4,1) - (-15.929410000000003) )**2
      x = x + (B(4,2) - (-7.54551) )**2
      x = x + (B(4,3) - (-4.9436100000000005) )**2
      x = x + (B(4,4) - (177.13894000000002) )**2
      x = x + (B(4,5) - (-2.92581) )**2
      x = x + (B(4,6) - (-2.42991) )**2
      x = x + (B(5,1) - (-12.678510000000001) )**2
      x = x + (B(5,2) - (-6.005609999999999) )**2
      x = x + (B(5,3) - (-3.9347100000000004) )**2
      x = x + (B(5,4) - (-2.92581) )**2
      x = x + (B(5,5) - (141.58524) )**2
      x = x + (B(5,6) - (-1.93401) )**2
      x = x + (B(6,1) - (-10.529610000000003) )**2
      x = x + (B(6,2) - (-4.98771) )**2
      x = x + (B(6,3) - (-3.2678100000000003) )**2
      x = x + (B(6,4) - (-2.42991) )**2
      x = x + (B(6,5) - (-1.9340100000000005) )**2
      x = x + (B(6,6) - (117.91554000000001) )**2
      if (dabs(x / det_l) > 1.d-12) then
          print *, 'N = 6: x = ', x
          test_qmckl_adjugate = 6
          return
      end if

    ,#+end_src
    #+end_example

    #+begin_src f90 :tangle (eval f_test)

  deallocate(A,B)

end function test_qmckl_adjugate
    #+end_src

    #+begin_src c :comments link :tangle (eval c_test)
qmckl_exit_code test_qmckl_adjugate(qmckl_context context);
assert(QMCKL_SUCCESS == test_qmckl_adjugate(context));
    #+end_src

* End of files                                                     :noexport:


  #+begin_src c :tangle (eval h_private_type)
#endif
  #+end_src

  #+begin_src c :tangle (eval h_private_func)
#endif
  #+end_src

   #+begin_src c :comments link :tangle (eval c_test)
  assert (qmckl_context_destroy(context) == QMCKL_SUCCESS);
  return 0;
}

   #+end_src


# -*- mode: org -*-
# vim: syntax=c