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 :comments link :tangle (eval c_test) :noweb yes
#include "qmckl.h"
#include "assert.h"
#ifdef HAVE_CONFIG_H
#include "config.h"
#endif
int main() {
  qmckl_context context;
  context = qmckl_context_create();

  #+end_src

* Matrix operations

** ~qmckl_dgemm~

   Matrix multiply: $C_{ij} = \beta C_{ij} + \alpha \sum_{k} A_{ik} \cdot B_{kj}$ using Fortran ~matmul~ function.

   TODO: Add description about the external library dependence.

   #+NAME: qmckl_dgemm_args
   | qmckl_context | context  | in  | Global state                                 |
   | bool          | TransA   | in  | Number of rows of the input matrix           |
   | bool          | TransB   | in  | Number of rows of the input matrix           |
   | int64_t       | m        | in  | Number of rows of the input matrix           |
   | int64_t       | n        | in  | Number of columns of the input matrix        |
   | int64_t       | k        | in  | Number of columns of the input matrix        |
   | double        | alpha    | in  | Number of columns of the input matrix        |
   | double        | A[][lda] | in  | Array containing the $m \times n$ matrix $A$ |
   | int64_t       | lda      | in  | Leading dimension of array ~A~               |
   | double        | B[][ldb] | in  | Array containing the $n \times m$ matrix $B$ |
   | int64_t       | ldb      | in  | Leading dimension of array ~B~               |
   | double        | beta     | in  | Array containing the $n \times m$ matrix $B$ |
   | double        | C[][ldc] | out | Array containing the $n \times m$ matrix $B$ |
   | int64_t       | ldc      | 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
      
*** C header

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

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

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

  integer*8 :: i,j,l, LDA_2, LDB_2

  info = QMCKL_SUCCESS

  if (TransA) then
  allocate(AT(k,m))
    do i = 1, m
      do j = 1, k
        AT(j,i) = A(i,j)
      end do
    end do
  LDA_2 = M
  else
  LDA_2 = LDA
  endif

  if (TransB) then
  allocate(BT(n,k))
    do i = 1, k
      do j = 1, n
        BT(j,i) = B(i,j)
      end do
    end do
  LDB_2 = K
  else
  LDB_2 = LDB
  endif

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

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

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

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

  if (LDA_2 /= m) then
     info = QMCKL_INVALID_ARG_9
     return
  endif

  if (LDB_2 /= k) then
     info = QMCKL_INVALID_ARG_10
     return
  endif

  if (LDC /= m) then
     info = QMCKL_INVALID_ARG_13
     return
  endif

   if (TransA) then

     if (alpha .eq. 1.0d0 .and. beta .eq. 0.0d0) then
       C = matmul(AT,B)
     else
       C = beta*C + alpha*matmul(AT,B)
     endif
   else if (TransB) then
     if (alpha .eq. 1.0d0 .and. beta .eq. 0.0d0) then
       C = matmul(A,BT)
     else
       C = beta*C + alpha*matmul(A,BT)
     endif
   else
     if (alpha .eq. 1.0d0 .and. beta .eq. 0.0d0) then
       C = matmul(A,B)
     else
       C = beta*C + alpha*matmul(A,B)
     endif
   endif
end function qmckl_dgemm_f
    #+end_src
    
*** C interface                                                    :noexport:

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

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

      use, intrinsic :: iso_c_binding
      implicit none

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

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

    end function qmckl_dgemm
    #+END_src


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

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

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

      end function qmckl_dgemm
    end interface
    #+end_src

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

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

  TransA = .False.
  TransB = .False.
  m = 1_8
  k = 4_8
  n = 6_8
  LDA = m
  LDB = k
  LDC = m

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

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

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

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

  if (test_qmckl_dgemm /= QMCKL_SUCCESS) return

  test_qmckl_dgemm = QMCKL_FAILURE

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

  if (dabs(x) <= 1.d-15) 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_invert~
   
   Matrix invert. Given a matrix M, returns a matrix M⁻¹ such that:


   \[
M · M^{-1} = I
\]

   This is a native Fortran implementation hand written (by: A. Scemama)
   only for small matrices (<=5x5).

   TODO: Add description about the external library dependence.

   #+NAME: qmckl_invert_args
   | qmckl_context | context  | in    | Global state                                 |
   | int64_t       | m        | in    | Number of rows of the input matrix           |
   | int64_t       | n        | in    | Number of columns of the input matrix        |
   | int64_t       | lda      | in    | Leading dimension of array ~A~               |
   | double        | A[][lda] | inout | Array containing the $m \times n$ matrix $A$ |
   | double        | det_l    | inout | determinant of A                             |
   
*** Requirements

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

*** C header

    #+CALL: generate_c_header(table=qmckl_invert_args,rettyp="qmckl_exit_code",fname="qmckl_invert")
    
    #+RESULTS:
    #+begin_src c :tangle (eval h_func) :comments org
    qmckl_exit_code qmckl_invert (
	  const qmckl_context context,
	  const int64_t m,
	  const int64_t n,
	  const int64_t lda,
	  double* A,
	  double det_l ); 
    #+end_src

*** Source
    #+begin_src f90 :tangle (eval f)
integer function qmckl_invert_f(context, ma, na, LDA, A, det_l) &
     result(info)
  use qmckl
  implicit none
  integer(qmckl_context)  , intent(in)  :: context
 double precision, intent(inout) :: A (LDA,na)
 integer*8, intent(in)             :: LDA
 integer*8, intent(in)             :: ma
 integer*8, intent(in)             :: na
 double precision, intent(inout) :: det_l
 integer :: i,j

 info = QMCKL_SUCCESS

 select case (na)
  case default
!DIR$ forceinline
    print *," TODO: Implement general invert"
    stop 0
  case (5)
!DIR$ forceinline
    call invert5(a,LDA,na,det_l)
  case (4)
!DIR$ forceinline
    call invert4(a,LDA,na,det_l)
  case (3)
!DIR$ forceinline
    call invert3(a,LDA,na,det_l)
  case (2)
!DIR$ forceinline
    call invert2(a,LDA,na,det_l)
  case (1)
!DIR$ forceinline
    call invert1(a,LDA,na,det_l)
  case (0)
    det_l=1.d0
 end select
end function qmckl_invert_f

subroutine invert1(a,LDA,na,det_l)
 implicit none
 double precision, intent(inout) :: a (LDA,na)
 integer*8, intent(in)             :: LDA
 integer*8, intent(in)             :: na
 double precision, intent(inout) :: det_l

 det_l = a(1,1)
 a(1,1) = 1.d0
end

subroutine invert2(a,LDA,na,det_l)
 implicit none
 double precision :: a (LDA,na)
 integer*8        :: LDA
 integer*8        :: na
 double precision :: det_l
 double precision :: b(2,2)

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

subroutine invert3(a,LDA,na,det_l)
 implicit none
 double precision, intent(inout) :: a (LDA,na)
 integer*8, intent(in)             :: LDA
 integer*8, intent(in)             :: na
 double precision, intent(inout) :: det_l
 double precision :: b(4,3)
 !DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: b
 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))
 do i=1,4
  b(i,1) = a(i,1)
  b(i,2) = a(i,2)
  b(i,3) = a(i,3)
 enddo
 a(1,1) =  b(2,2)*b(3,3) - b(2,3)*b(3,2)
 a(2,1) =  b(2,3)*b(3,1) - b(2,1)*b(3,3)
 a(3,1) =  b(2,1)*b(3,2) - b(2,2)*b(3,1)

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

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

end

subroutine invert4(a,LDA,na,det_l)
 implicit none
 double precision, intent(inout) :: a (LDA,na)
 integer*8, intent(in)             :: LDA
 integer*8, intent(in)             :: na
 double precision, intent(inout) :: det_l
 double precision :: b(4,4)
 !DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: b
 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)))
 do i=1,4
  b(1,i) = a(1,i)
  b(2,i) = a(2,i)
  b(3,i) = a(3,i)
  b(4,i) = a(4,i)
 enddo
 
 a(1,1) =  b(2,2)*(b(3,3)*b(4,4)-b(3,4)*b(4,3))-b(2,3)*(b(3,2)*b(4,4)-b(3,4)*b(4,2))+b(2,4)*(b(3,2)*b(4,3)-b(3,3)*b(4,2))
 a(2,1) = -b(2,1)*(b(3,3)*b(4,4)-b(3,4)*b(4,3))+b(2,3)*(b(3,1)*b(4,4)-b(3,4)*b(4,1))-b(2,4)*(b(3,1)*b(4,3)-b(3,3)*b(4,1))
 a(3,1) =  b(2,1)*(b(3,2)*b(4,4)-b(3,4)*b(4,2))-b(2,2)*(b(3,1)*b(4,4)-b(3,4)*b(4,1))+b(2,4)*(b(3,1)*b(4,2)-b(3,2)*b(4,1))
 a(4,1) = -b(2,1)*(b(3,2)*b(4,3)-b(3,3)*b(4,2))+b(2,2)*(b(3,1)*b(4,3)-b(3,3)*b(4,1))-b(2,3)*(b(3,1)*b(4,2)-b(3,2)*b(4,1))

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

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

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

end

subroutine invert5(a,LDA,na,det_l)
 implicit none
 double precision, intent(inout) :: a (LDA,na)
 integer*8, intent(in)             :: LDA
 integer*8, intent(in)             :: na
 double precision, intent(inout) :: det_l
 double precision :: b(5,5)
 !DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: b
 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))))

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

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

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

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

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

end
    #+end_src

*** C interface                                                    :noexport:
    
    #+CALL: generate_c_interface(table=qmckl_invert_args,rettyp="qmckl_exit_code",fname="qmckl_invert")

    #+RESULTS:
    #+begin_src f90 :tangle (eval f) :comments org :exports none
    integer(c_int32_t) function qmckl_invert &
	(context, m, n, lda, A, 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 :: m
      integer (c_int64_t) , intent(in)  , value :: n
      integer (c_int64_t) , intent(in)  , value :: lda
      real    (c_double ) , intent(inout)        :: A(lda,*)
      real    (c_double ) , intent(inout)        :: det_l

      integer(c_int32_t), external :: qmckl_invert_f
      info = qmckl_invert_f &
	     (context, m, n, lda, A, det_l)

    end function qmckl_invert
    #+end_src

    #+CALL: generate_f_interface(table=qmckl_invert_args,rettyp="qmckl_exit_code",fname="qmckl_invert")

    #+RESULTS:
    #+begin_src f90 :tangle (eval fh_func) :comments org :exports none
    interface
      integer(c_int32_t) function qmckl_invert &
	  (context, m, n, lda, A, 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 :: m
	integer (c_int64_t) , intent(in)  , value :: n
	integer (c_int64_t) , intent(in)  , value :: lda
	real    (c_double ) , intent(inout)        :: A(lda,*)
	real    (c_double ) , intent(inout)        :: det_l

      end function qmckl_invert
    end interface
    #+end_src
    
*** Test                                                           :noexport:
    #+begin_src f90 :tangle (eval f_test)
integer(qmckl_exit_code) function test_qmckl_invert(context) bind(C)
  use qmckl
  implicit none
  integer(qmckl_context), intent(in), value :: context

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

  m = 4_8
  k = 4_8
  LDA = m
  LDB = m
  LDC = m

  allocate( A(LDA,k), C(LDC,k))

  A = 0.10d0
  C = 0.d0
  A(1,1) = 1.0d0;
  A(2,2) = 2.0d0;
  A(3,3) = 3.0d0;
  A(4,4) = 4.0d0;

  ! Exact inverse (Mathematica)
  C(1,1) = 1.0102367161391992d0 
  C(2,2) = 0.5036819224578257d0
  C(3,3) = 0.33511197860555897d0
  C(4,4) = 0.2510382472105688d0
  C(1,2) = -0.047782608144589914d0
  C(1,3) = -0.031305846715420985d0
  C(1,4) = -0.023278706531979707d0
  C(2,3) = -0.014829085286252043d0
  C(2,4) = -0.011026755725674596d0
  C(3,4) = -0.007224426165097149d0 
  C(2,1) = -0.047782608144589914d0
  C(3,1) = -0.031305846715420985d0
  C(4,1) = -0.023278706531979707d0
  C(3,2) = -0.014829085286252043d0
  C(4,2) = -0.011026755725674596d0
  C(4,3) = -0.007224426165097149d0 
  det_l_ref = 23.6697d0

  test_qmckl_invert = qmckl_invert(context, m, k, LDA, A, det_l)

  if (test_qmckl_invert /= QMCKL_SUCCESS) return

  test_qmckl_invert = QMCKL_FAILURE

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

  if (dabs(x) <= 1.d-15 .and. (dabs(det_l_ref - det_l)) <= 1.d-15) then
     test_qmckl_invert = QMCKL_SUCCESS
  endif

  deallocate(A,C)
end function test_qmckl_invert
    #+end_src
    
    #+begin_src c :comments link :tangle (eval c_test)
qmckl_exit_code test_qmckl_invert(qmckl_context context);
assert(QMCKL_SUCCESS == test_qmckl_invert(context));
    #+end_src
    
* End of files                                                     :noexport:

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

   #+end_src


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