diff --git a/configure.ac b/configure.ac
index fa3277d..5850b6a 100644
--- a/configure.ac
+++ b/configure.ac
@@ -124,10 +124,10 @@ PKG_CFLAGS="$PKG_CFLAGS $TREXIO_CFLAGS"
 PKG_LIBS="$PKG_LIBS $TREXIO_LIBS"
 
 ## BLAS
-#AX_BLAS([], [AC_MSG_ERROR([BLAS was not found.])])
+AX_BLAS([], [AC_MSG_ERROR([BLAS was not found.])])
 
 ## LAPACK
-#AX_LAPACK([], [AC_MSG_ERROR([LAPACK was not found.])])
+AX_LAPACK([], [AC_MSG_ERROR([LAPACK was not found.])])
 
 # Options.
 
@@ -252,7 +252,7 @@ fi
 #mkl-dynamic-lp64-seq
 
 PKG_LIBS="$PKG_LIBS $LIBS"
-LIBS="$LAPACK_LIBS $BLAS_LIBS $PKG_LIBS"
+LIBS="$BLAS_LIBS $LAPACK_LIBS $BLAS_LIBS $PKG_LIBS"
 CFLAGS="$CFLAGS $PKG_CFLAGS"
 AC_SUBST([PKG_LIBS])
 AC_SUBST([PKG_CFLAGS])
diff --git a/org/qmckl_blas.org b/org/qmckl_blas.org
index e560095..6048b39 100644
--- a/org/qmckl_blas.org
+++ b/org/qmckl_blas.org
@@ -28,20 +28,20 @@ int main() {
    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~               |
+   | 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 matrix $A$       |
+   | int64_t       | lda      | in  | Leading dimension of array ~A~        |
+   | double        | B[][ldb] | in  | Array containing the matrix $B$       |
+   | int64_t       | ldb      | in  | Leading dimension of array ~B~        |
+   | double        | beta     | in  | Array containing the matrix $B$       |
+   | double        | C[][ldc] | out | Array containing the matrix $B$       |
+   | int64_t       | ldc      | in  | Leading dimension of array ~B~        |
 
 *** Requirements
 
@@ -85,67 +85,67 @@ integer function qmckl_dgemm_f(context, TransA, TransB, m, n, k, alpha, A, LDA,
      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(lda,*)
-  integer*8  , intent(in)  :: ldb
-  real*8     , intent(in)  :: B(ldb,*)
-  integer*8  , intent(in)  :: ldc
-  real*8     , intent(out) :: C(ldc,*)
-  real*8, allocatable      :: AT(:,:), BT(:,:), CT(:,:)
-  integer*4 :: qmckl_dgemm_N_N_f
+  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(lda,*)
+  integer*8             , intent(in)  :: ldb
+  real*8                , intent(in)  :: B(ldb,*)
+  integer*8             , intent(in)  :: ldc
+  real*8                , intent(out) :: C(ldc,*)
 
-  integer*8 :: i,j,l, LDA_2, LDB_2
+  real*8, allocatable  :: AT(:,:), BT(:,:), CT(:,:)
+  integer*4            :: qmckl_dgemm_N_N_f
+  integer*8            :: i,j,l, LDA_2, LDB_2
 
   info = QMCKL_SUCCESS
 
   if (TransA) then
-  allocate(AT(m,k))
-    do i = 1, k
-      do j = 1, m
-        AT(j,i) = A(i,j)
-      end do
-    end do
-  LDA_2 = M
+     allocate(AT(m,k))
+     do i = 1, k
+        do j = 1, m
+           AT(j,i) = A(i,j)
+        end do
+     end do
+     LDA_2 = M
   else
-  LDA_2 = LDA
+     LDA_2 = LDA
   endif
-
+  
   if (TransB) then
-  allocate(BT(k,n))
-    do i = 1, n
-      do j = 1, k
-        BT(j,i) = B(i,j)
-      end do
-    end do
-  LDB_2 = K
+     allocate(BT(k,n))
+     do i = 1, n
+        do j = 1, k
+           BT(j,i) = B(i,j)
+        end do
+     end do
+     LDB_2 = K
   else
-  LDB_2 = LDB
+     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
@@ -161,15 +161,16 @@ integer function qmckl_dgemm_f(context, TransA, TransB, m, n, k, alpha, A, LDA,
      return
   endif
 
-   if (TransA) then
+  if (TransA) then
      info = qmckl_dgemm_N_N_f(context, m, n, k, alpha, AT, LDA_2, B, LDB_2, beta, c, LDC)
-   else if (TransB) then
+  else if (TransB) then
      info = qmckl_dgemm_N_N_f(context, m, n, k, alpha, A, LDA_2, BT, LDB_2, beta, c, LDC)
-   else if (TransA .and. TransB) then
+  else if (TransA .and. TransB) then
      info = qmckl_dgemm_N_N_f(context, m, n, k, alpha, AT, LDA_2, BT, LDB_2, beta, c, LDC)
-   else
+  else
      info = qmckl_dgemm_N_N_f(context, m, n, k, alpha, A, LDA_2, B, LDB_2, beta, c, LDC)
-   endif
+  endif
+
 end function qmckl_dgemm_f
 
 integer function qmckl_dgemm_N_N_f(context, m, n, k, alpha, A, LDA, B, LDB, beta, C, LDC) &
@@ -185,51 +186,52 @@ integer function qmckl_dgemm_N_N_f(context, m, n, k, alpha, A, LDA, B, LDB, beta
   real*8     , intent(in)  :: B(ldb,n)
   integer*8  , intent(in)  :: ldc
   real*8     , intent(out) :: C(ldc,n)
-
+  
   integer*8 :: i,j,l, LDA_2, LDB_2
-
+  
   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
+     info = QMCKL_INVALID_ARG_2
      return
   endif
 
   if (n <= 0_8) then
-     info = QMCKL_INVALID_ARG_5
+     info = QMCKL_INVALID_ARG_3
      return
   endif
 
   if (k <= 0_8) then
-     info = QMCKL_INVALID_ARG_6
+     info = QMCKL_INVALID_ARG_4
      return
   endif
 
   if (LDA /= m) then
-     info = QMCKL_INVALID_ARG_9
+     info = QMCKL_INVALID_ARG_7
      return
   endif
 
   if (LDB /= k) then
-     info = QMCKL_INVALID_ARG_10
+     info = QMCKL_INVALID_ARG_8
      return
   endif
 
   if (LDC /= m) then
-     info = QMCKL_INVALID_ARG_13
+     info = QMCKL_INVALID_ARG_11
      return
   endif
 
   if (alpha == 1.0d0 .and. beta == 0.0d0) then
-    C = matmul(A,B)
+     C = matmul(A,B)
   else
-    C = beta*C + alpha*matmul(A,B)
+     C = beta*C + alpha*matmul(A,B)
   endif
+
 end function qmckl_dgemm_N_N_f
     #+end_src
     
@@ -267,52 +269,52 @@ end function qmckl_dgemm_N_N_f
 
     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 &
+  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
+       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
@@ -323,7 +325,7 @@ integer(qmckl_exit_code) function test_qmckl_dgemm(context) bind(C)
   LDC = m
 
   allocate( A(LDA,k), B(LDB,n) , C(LDC,n), D(LDC,n))
-
+  
   A = 0.d0
   B = 0.d0
   C = 0.d0
@@ -335,34 +337,35 @@ integer(qmckl_exit_code) function test_qmckl_dgemm(context) bind(C)
         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
+        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
+  
+  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
 
@@ -371,377 +374,382 @@ qmckl_exit_code test_qmckl_dgemm(qmckl_context context);
 assert(QMCKL_SUCCESS == test_qmckl_dgemm(context));
     #+end_src
 
-** ~qmckl_adjoint~
+** ~qmckl_adjugate~
    
-   Matrix adjoint. Given a matrix M, returns a matrix M⁻¹ such that:
-
+   Given a matrix $\mathbf{A}$, the adjugate matrix
+   $\text{adj}(\mathbf{A})$ is the transpose of the cofactors matrix
+   of $\mathbf{A}$.
 
    \[
-M · M^{-1} = I
-\]
+   \text{adj}(\mathbf{A}) = \text{det}(\mathbf{A}) \, \mathbf{A}^{-1}
+   \]
 
-   This is a native Fortran implementation hand written (by: A. Scemama)
-   only for small matrices (<=5x5).
+   See also: https://en.wikipedia.org/wiki/Adjugate_matrix
 
-   TODO: Add description about the external library dependence.
-
-   #+NAME: qmckl_adjoint_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                             |
+   #+NAME: qmckl_adjugate_args
+   | qmckl_context | context  | in    | Global state                                   |
+   | int64_t       | n        | in    | Number of rows and columns of the input matrix |
+   | int64_t       | lda      | in    | Leading dimension of array ~A~                 |
+   | double        | A[][lda] | inout | Array containing the $n \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
-
+    - ~A~ is allocated with at least $m \times m \times 8$ bytes
+      
 *** C header
 
-    #+CALL: generate_c_header(table=qmckl_adjoint_args,rettyp="qmckl_exit_code",fname="qmckl_adjoint")
+    #+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_adjoint (
+    qmckl_exit_code qmckl_adjugate (
 	  const qmckl_context context,
-	  const int64_t m,
 	  const int64_t n,
 	  const int64_t lda,
 	  double* A,
 	  double det_l ); 
     #+end_src
-
+    
 *** Source
+
+   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_adjoint_f(context, ma, na, LDA, A, det_l) &
+integer function qmckl_adjugate_f(context, 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
+  double precision, intent(inout)       :: A (LDA,na)
+  integer*8, intent(in)                 :: LDA
+  integer*8, intent(in)                 :: na
+  double precision, intent(inout)       :: det_l
 
- info = QMCKL_SUCCESS
-
- select case (na)
+  integer    :: i,j
+  
+  !TODO CHECK ARGUMENTS
+  info = QMCKL_SUCCESS
+  
+  select case (na)
   case default
-    print *," TODO: Implement general adjoint"
-    stop 0
+     call adjugate_general(context, na, LDA, A, det_l)
   case (5)
-    call adjoint5(a,LDA,na,det_l)
+     call adjugate5(a,LDA,na,det_l)
   case (4)
-    call adjoint4(a,LDA,na,det_l)
-    
+     call adjugate4(a,LDA,na,det_l)
   case (3)
-    call adjoint3(a,LDA,na,det_l)
+     call adjugate3(a,LDA,na,det_l)
   case (2)
-    call adjoint2(a,LDA,na,det_l)
+     call adjugate2(a,LDA,na,det_l)
   case (1)
-    call adjoint1(a,LDA,na,det_l)
+     call adjugate1(a,LDA,na,det_l)
   case (0)
-    det_l=1.d0
- end select
-end function qmckl_adjoint_f
+      det_l=1.d0
+  end select
+   
+end function qmckl_adjugate_f
+    #+end_src
 
-subroutine adjoint1(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
+    #+begin_src f90 :tangle (eval f)
+subroutine adjugate1(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
+  
+  call cofactor1(a,LDA,na,det_l)
+end subroutine adjugate1
 
- call cofactor1(a,LDA,na,det_l)
-end
+subroutine adjugate2(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(2,2)
+  
+  call cofactor2(a,LDA,na,det_l)
+  
+  ! Calculate the transpose
+  b(1,1) = a(1,1)
+  b(1,2) = a(2,1)
+  b(2,1) = a(1,2)
+  b(2,2) = a(2,2)
+  a(1,1) = b(1,1)
+  a(1,2) = b(1,2)
+  a(2,1) = b(2,1)
+  a(2,2) = b(2,2)
+end subroutine adjugate2
 
-subroutine adjoint2(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(2,2)
+subroutine adjugate3(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(3,3)
+  
+  call cofactor3(a,LDA,na,det_l)
+  
+  ! Calculate the transpose
+  b(1,1) = a(1,1)
+  b(1,2) = a(2,1)
+  b(1,3) = a(3,1)
+  b(2,1) = a(1,2)
+  b(2,2) = a(2,2)
+  b(2,3) = a(3,2)
+  b(3,1) = a(1,3)
+  b(3,2) = a(2,3)
+  b(3,3) = a(3,3)
+  ! copy
+  a(1,1) = b(1,1)
+  a(2,1) = b(2,1)
+  a(3,1) = b(3,1)
+  a(1,2) = b(1,2)
+  a(2,2) = b(2,2)
+  a(3,2) = b(3,2)
+  a(1,3) = b(1,3)
+  a(2,3) = b(2,3)
+  a(3,3) = b(3,3)
+end subroutine adjugate3
 
- call cofactor2(a,LDA,na,det_l)
+subroutine adjugate4(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)
+  
+  call cofactor4(a,LDA,na,det_l)
+  
+  ! Calculate the transpose
+  b(1,1) = a(1,1)
+  b(1,2) = a(2,1)
+  b(1,3) = a(3,1)
+  b(1,4) = a(4,1)
+  b(2,1) = a(1,2)
+  b(2,2) = a(2,2)
+  b(2,3) = a(3,2)
+  b(2,4) = a(4,2)
+  b(3,1) = a(1,3)
+  b(3,2) = a(2,3)
+  b(3,3) = a(3,3)
+  b(3,4) = a(4,3)
+  b(4,1) = a(1,4)
+  b(4,2) = a(2,4)
+  b(4,3) = a(3,4)
+  b(4,4) = a(4,4)
+  ! copy
+  a(1,1) = b(1,1)
+  a(2,1) = b(2,1)
+  a(3,1) = b(3,1)
+  a(4,1) = b(4,1)
+  a(1,2) = b(1,2)
+  a(2,2) = b(2,2)
+  a(3,2) = b(3,2)
+  a(4,2) = b(4,2)
+  a(1,3) = b(1,3)
+  a(2,3) = b(2,3)
+  a(3,3) = b(3,3)
+  a(4,3) = b(4,3)
+  a(1,4) = b(1,4)
+  a(2,4) = b(2,4)
+  a(3,4) = b(3,4)
+  a(4,4) = b(4,4)
+end subroutine adjugate4
 
- ! Calculate the transpose
- b(1,1) = a(1,1)
- b(1,2) = a(2,1)
- b(2,1) = a(1,2)
- b(2,2) = a(2,2)
- a(1,1) = b(1,1)
- a(1,2) = b(1,2)
- a(2,1) = b(2,1)
- a(2,2) = b(2,2)
-end
-
-subroutine adjoint3(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(3,3)
-
- call cofactor3(a,LDA,na,det_l)
-
- ! Calculate the transpose
- b(1,1) = a(1,1)
- b(1,2) = a(2,1)
- b(1,3) = a(3,1)
- b(2,1) = a(1,2)
- b(2,2) = a(2,2)
- b(2,3) = a(3,2)
- b(3,1) = a(1,3)
- b(3,2) = a(2,3)
- b(3,3) = a(3,3)
- ! copy
- a(1,1) = b(1,1)
- a(2,1) = b(2,1)
- a(3,1) = b(3,1)
- a(1,2) = b(1,2)
- a(2,2) = b(2,2)
- a(3,2) = b(3,2)
- a(1,3) = b(1,3)
- a(2,3) = b(2,3)
- a(3,3) = b(3,3)
-end
-
-subroutine adjoint4(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)
-
- call cofactor4(a,LDA,na,det_l)
-
- ! Calculate the transpose
- b(1,1) = a(1,1)
- b(1,2) = a(2,1)
- b(1,3) = a(3,1)
- b(1,4) = a(4,1)
- b(2,1) = a(1,2)
- b(2,2) = a(2,2)
- b(2,3) = a(3,2)
- b(2,4) = a(4,2)
- b(3,1) = a(1,3)
- b(3,2) = a(2,3)
- b(3,3) = a(3,3)
- b(3,4) = a(4,3)
- b(4,1) = a(1,4)
- b(4,2) = a(2,4)
- b(4,3) = a(3,4)
- b(4,4) = a(4,4)
- ! copy
- a(1,1) = b(1,1)
- a(2,1) = b(2,1)
- a(3,1) = b(3,1)
- a(4,1) = b(4,1)
- a(1,2) = b(1,2)
- a(2,2) = b(2,2)
- a(3,2) = b(3,2)
- a(4,2) = b(4,2)
- a(1,3) = b(1,3)
- a(2,3) = b(2,3)
- a(3,3) = b(3,3)
- a(4,3) = b(4,3)
- a(1,4) = b(1,4)
- a(2,4) = b(2,4)
- a(3,4) = b(3,4)
- a(4,4) = b(4,4)
-end
-
-subroutine adjoint5(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)
-
- call cofactor5(a,LDA,na,det_l)
-
- ! Calculate the transpose
- b(1,1) = a(1,1)
- b(1,2) = a(2,1)
- b(1,3) = a(3,1)
- b(1,4) = a(4,1)
- b(1,5) = a(5,1)
- b(2,1) = a(1,2)
- b(2,2) = a(2,2)
- b(2,3) = a(3,2)
- b(2,4) = a(4,2)
- b(2,5) = a(5,2)
- b(3,1) = a(1,3)
- b(3,2) = a(2,3)
- b(3,3) = a(3,3)
- b(3,4) = a(4,3)
- b(3,5) = a(5,3)
- b(4,1) = a(1,4)
- b(4,2) = a(2,4)
- b(4,3) = a(3,4)
- b(4,4) = a(4,4)
- b(4,5) = a(5,4)
- b(5,1) = a(1,5)
- b(5,2) = a(2,5)
- b(5,3) = a(3,5)
- b(5,4) = a(4,5)
- b(5,5) = a(5,5)
- ! copy
- a(1,1) = b(1,1)
- a(2,1) = b(2,1)
- a(3,1) = b(3,1)
- a(4,1) = b(4,1)
- a(5,1) = b(5,1)
- a(1,2) = b(1,2)
- a(2,2) = b(2,2)
- a(3,2) = b(3,2)
- a(4,2) = b(4,2)
- a(5,2) = b(5,2)
- a(1,3) = b(1,3)
- a(2,3) = b(2,3)
- a(3,3) = b(3,3)
- a(4,3) = b(4,3)
- a(5,3) = b(5,3)
- a(1,4) = b(1,4)
- a(2,4) = b(2,4)
- a(3,4) = b(3,4)
- a(4,4) = b(4,4)
- a(5,4) = b(5,4)
- a(1,5) = b(1,5)
- a(2,5) = b(2,5)
- a(3,5) = b(3,5)
- a(4,5) = b(4,5)
- a(5,5) = b(5,5)
-end
+subroutine adjugate5(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)
+  
+  call cofactor5(a,LDA,na,det_l)
+  
+  ! Calculate the transpose
+  b(1,1) = a(1,1)
+  b(1,2) = a(2,1)
+  b(1,3) = a(3,1)
+  b(1,4) = a(4,1)
+  b(1,5) = a(5,1)
+  b(2,1) = a(1,2)
+  b(2,2) = a(2,2)
+  b(2,3) = a(3,2)
+  b(2,4) = a(4,2)
+  b(2,5) = a(5,2)
+  b(3,1) = a(1,3)
+  b(3,2) = a(2,3)
+  b(3,3) = a(3,3)
+  b(3,4) = a(4,3)
+  b(3,5) = a(5,3)
+  b(4,1) = a(1,4)
+  b(4,2) = a(2,4)
+  b(4,3) = a(3,4)
+  b(4,4) = a(4,4)
+  b(4,5) = a(5,4)
+  b(5,1) = a(1,5)
+  b(5,2) = a(2,5)
+  b(5,3) = a(3,5)
+  b(5,4) = a(4,5)
+  b(5,5) = a(5,5)
+  ! copy
+  a(1,1) = b(1,1)
+  a(2,1) = b(2,1)
+  a(3,1) = b(3,1)
+  a(4,1) = b(4,1)
+  a(5,1) = b(5,1)
+  a(1,2) = b(1,2)
+  a(2,2) = b(2,2)
+  a(3,2) = b(3,2)
+  a(4,2) = b(4,2)
+  a(5,2) = b(5,2)
+  a(1,3) = b(1,3)
+  a(2,3) = b(2,3)
+  a(3,3) = b(3,3)
+  a(4,3) = b(4,3)
+  a(5,3) = b(5,3)
+  a(1,4) = b(1,4)
+  a(2,4) = b(2,4)
+  a(3,4) = b(3,4)
+  a(4,4) = b(4,4)
+  a(5,4) = b(5,4)
+  a(1,5) = b(1,5)
+  a(2,5) = b(2,5)
+  a(3,5) = b(3,5)
+  a(4,5) = b(4,5)
+  a(5,5) = b(5,5)
+end subroutine adjugate5
 
 subroutine cofactor1(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
+  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 cofactor1
 
 subroutine cofactor2(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
+  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 cofactor2
 
 subroutine cofactor3(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)
- 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)
- end do
- 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)
+  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)
+  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))
 
- 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)
+  do i=1,4
+     b(i,1) = a(i,1)
+     b(i,2) = a(i,2)
+     b(i,3) = a(i,3)
+  end do
 
- 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
+  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 cofactor3
 
 subroutine cofactor4(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)
- 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)
- end do
- 
- 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
+  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)
+  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)
+  end do
+  
+  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 cofactor4
 
 subroutine cofactor5(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)
- integer :: i,j
+  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)
+  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)- &
@@ -781,11 +789,11 @@ subroutine cofactor5(a,LDA,na,det_l)
  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)
+    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)
  end do
  
  a(1,1) = &
@@ -921,126 +929,399 @@ subroutine cofactor5(a,LDA,na,det_l)
 end
     #+end_src
 
+    #+begin_src f90 :tangle (eval f)
+subroutine adjugate_general(context, na, LDA, A, det_l)
+  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)                 :: na
+  double precision, intent(inout)       :: det_l
+
+  double precision :: work(LDA*max(na,64))
+  integer          :: inf
+  integer          :: ipiv(LDA)
+  integer          :: lwork
+  integer          :: i, j
+
+    #+end_src
+
+    For larger matrices, we first compute the LU factorization $LU=A$
+    using the ~dgetrf~ routine.
+
+    #+begin_src f90 :tangle (eval f)
+  call dgetrf(na, na, a, LDA, 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
+  do i=1,na
+   j = j+min(abs(ipiv(i)-i),1)
+   det_l = det_l*a(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) /= 0)  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, a, LDA, 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)
+  a(:,:) = a(:,:)*det_l
+
+end subroutine adjugate_general
+    #+end_src
+
 *** C interface                                                    :noexport:
     
-    #+CALL: generate_c_interface(table=qmckl_adjoint_args,rettyp="qmckl_exit_code",fname="qmckl_adjoint")
+    #+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_adjoint &
-	(context, m, n, lda, A, det_l) &
-	bind(C) result(info)
+    integer(c_int32_t) function qmckl_adjugate &
+        (context, 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_adjoint_f
-      info = qmckl_adjoint_f &
-	     (context, m, n, lda, A, det_l)
+      integer(c_int32_t), external :: qmckl_adjugate_f
+      info = qmckl_adjugate_f &
+             (context, n, lda, A, det_l)
 
-    end function qmckl_adjoint
+    end function qmckl_adjugate
     #+end_src
 
-    #+CALL: generate_f_interface(table=qmckl_adjoint_args,rettyp="qmckl_exit_code",fname="qmckl_adjoint")
+    #+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_adjoint &
-	  (context, m, n, lda, A, det_l) &
-	  bind(C)
-	use, intrinsic :: iso_c_binding
-	import
-	implicit none
+      integer(c_int32_t) function qmckl_adjugate &
+          (context, 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
+        integer (c_int64_t) , intent(in)  , value :: context
+        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_adjoint
+      end function qmckl_adjugate
     end interface
     #+end_src
     
 *** Test                                                           :noexport:
+
     #+begin_src f90 :tangle (eval f_test)
-integer(qmckl_exit_code) function test_qmckl_adjoint(context) bind(C)
+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(:,:), C(:,:)
-  integer*8                     :: m, n, k, LDA, LDB, LDC
+  
+  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
 
-  m = 4_8
-  k = 4_8
-  LDA = m
-  LDB = m
-  LDC = m
+  LDA = 6
+  LDB = 6
+  
+  allocate( A(LDA,6), B(LDB,6))
 
-  allocate( A(LDA,k), C(LDC,k))
-
-  A = 0.10d0
-  C = 0.d0
+  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;
 
-  ! 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_adjugate = QMCKL_SUCCESS
+  
+    #+end_src
+    
+    #+begin_src python :results output :output drawer
+import numpy as np
+import numpy.linalg as la
+N=6
 
-  test_qmckl_adjoint = qmckl_adjoint(context, m, k, LDA, A, det_l)
+A = np.zeros( (N,N) )
+A += 0.1
+for i in range(N):
+    A[i][i] = i+1
 
-  if (test_qmckl_adjoint /= QMCKL_SUCCESS) return
+def adj(A):
+    return la.det(A) * la.inv(A)
 
-  test_qmckl_adjoint = QMCKL_FAILURE
+print ("#+begin_src f90 :tangle (eval f_test)")
 
-  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
+for i in range(N):
+    print ("! N = ", i+1)
+    print (f"  B(:,:) = A(:,:)")
+    print (f"  test_qmckl_adjugate = qmckl_adjugate(context, {i+1}_8, LDB, B, 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"")
+  
 
-  if (dabs(x) <= 1.d-15 .and. (dabs(det_l_ref - det_l)) <= 1.d-15) then
-     test_qmckl_adjoint = QMCKL_SUCCESS
-  endif
+print ("#+end_src")
+#    print(adj(A[0:i+1,0:i+1]))
+    #+end_src
 
-  deallocate(A,C)
-end function test_qmckl_adjoint
+    #+RESULTS:
+    #+begin_src f90 :tangle (eval f_test)
+    ! N =  1
+      B(:,:) = A(:,:)
+      test_qmckl_adjugate = qmckl_adjugate(context, 1_8, LDB, B, 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
+      B(:,:) = A(:,:)
+      test_qmckl_adjugate = qmckl_adjugate(context, 2_8, LDB, B, 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
+      B(:,:) = A(:,:)
+      test_qmckl_adjugate = qmckl_adjugate(context, 3_8, LDB, B, 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
+      B(:,:) = A(:,:)
+      test_qmckl_adjugate = qmckl_adjugate(context, 4_8, LDB, B, 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
+      B(:,:) = A(:,:)
+      test_qmckl_adjugate = qmckl_adjugate(context, 5_8, LDB, B, 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
+      B(:,:) = A(:,:)
+      test_qmckl_adjugate = qmckl_adjugate(context, 6_8, LDB, B, 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
+
+
+    #+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_adjoint(qmckl_context context);
-assert(QMCKL_SUCCESS == test_qmckl_adjoint(context));
+qmckl_exit_code test_qmckl_adjugate(qmckl_context context);
+assert(QMCKL_SUCCESS == test_qmckl_adjugate(context));
     #+end_src
     
 * End of files                                                     :noexport:
diff --git a/org/qmckl_determinant.org b/org/qmckl_determinant.org
index 1e8e559..ddef91c 100644
--- a/org/qmckl_determinant.org
+++ b/org/qmckl_determinant.org
@@ -1817,7 +1817,7 @@ integer function qmckl_compute_det_inv_matrix_alpha_f(context, &
   do iwalk = 1, walk_num
     ! Value
     matA(1:alpha_num,1:alpha_num) = det_vgl_alpha(1:alpha_num, 1:alpha_num, 1, iwalk, idet)
-    res = qmckl_adjoint(context, alpha_num, alpha_num, LDA, matA, det_l)
+    res = qmckl_adjugate(context, alpha_num, LDA, matA, det_l)
     det_adj_matrix_alpha(1:alpha_num, 1:alpha_num, iwalk, idet) = matA
     det_inv_matrix_alpha(1:alpha_num, 1:alpha_num, iwalk, idet) = matA/det_l
     det_value_alpha(iwalk, idet) = det_l
@@ -1948,7 +1948,7 @@ integer function qmckl_compute_det_inv_matrix_beta_f(context, &
   do iwalk = 1, walk_num
     ! Value
     matA(1:beta_num,1:beta_num) = det_vgl_beta(1:beta_num, 1:beta_num, 1, iwalk, idet)
-    res = qmckl_adjoint(context, beta_num, beta_num, LDA, matA, det_l)
+    res = qmckl_adjugate(context, beta_num, LDA, matA, det_l)
     det_adj_matrix_beta(1:beta_num, 1:beta_num, iwalk, idet) = matA
     det_inv_matrix_beta(1:beta_num, 1:beta_num, iwalk, idet) = matA/det_l
     det_value_beta(iwalk, idet) = det_l
diff --git a/org/qmckl_tests.org b/org/qmckl_tests.org
index 21ffe54..16741de 100644
--- a/org/qmckl_tests.org
+++ b/org/qmckl_tests.org
@@ -60040,7 +60040,6 @@ double chbrclf_elec_coord[chbrclf_walk_num][chbrclf_elec_num][3] = { {
 
 #+END_src
 
-
 * N2
 
   This test is mainly for the Jastrow factor and was supplied by