Added Lapack for general determinants

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])

org/qmckl_blas.org

@@ -35,12 +35,12 @@ int main() {
    | 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$ |
+   | 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 $n \times m$ matrix $B$ |
+   | double        | B[][ldb] | in  | Array containing the matrix $B$       |
    | int64_t       | ldb      | in  | Leading dimension of array ~B~        |
-   | double        | beta     | in  | Array containing the $n \times m$ matrix $B$ |
+   | double        | beta     | in  | Array containing the matrix $B$       |
-   | double        | C[][ldc] | out | Array containing the $n \times m$ matrix $B$ |
+   | double        | C[][ldc] | out | Array containing the matrix $B$       |
    | int64_t       | ldc      | in  | Leading dimension of array ~B~        |
 
 *** Requirements
@@ -85,7 +85,7 @@ 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
+  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
@@ -95,9 +95,9 @@ integer function qmckl_dgemm_f(context, TransA, TransB, m, n, k, alpha, A, LDA,
   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*8            :: i,j,l, LDA_2, LDB_2
 
   info = QMCKL_SUCCESS
@@ -170,6 +170,7 @@ integer function qmckl_dgemm_f(context, TransA, TransB, m, n, k, alpha, A, LDA,
   else
      info = qmckl_dgemm_N_N_f(context, m, n, k, alpha, A, LDA_2, B, LDB_2, beta, c, LDC)
   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) &
@@ -196,32 +197,32 @@ integer function qmckl_dgemm_N_N_f(context, m, n, k, alpha, A, LDA, B, LDB, beta
   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
 
@@ -230,6 +231,7 @@ integer function qmckl_dgemm_N_N_f(context, m, n, k, alpha, A, LDA, B, LDB, beta
   else
      C = beta*C + alpha*matmul(A,B)
   endif
+
 end function qmckl_dgemm_N_N_f
     #+end_src
     
@@ -358,11 +360,12 @@ integer(qmckl_exit_code) function test_qmckl_dgemm(context) bind(C)
      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,45 +374,40 @@ 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)
+   See also: https://en.wikipedia.org/wiki/Adjugate_matrix
-   only for small matrices (<=5x5).
 
-   TODO: Add description about the external library dependence.
+   #+NAME: qmckl_adjugate_args
-
-   #+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 rows and columns 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        | A[][lda] | inout | Array containing the $n \times n$ matrix $A$   |
-   | double        | det_l    | inout | determinant of 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,
@@ -417,42 +415,49 @@ M · M^{-1} = I
     #+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
   
+  !TODO CHECK ARGUMENTS
   info = QMCKL_SUCCESS
   
   select case (na)
   case default
-    print *," TODO: Implement general adjoint"
+     call adjugate_general(context, na, LDA, A, det_l)
-    stop 0
   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
    
-subroutine adjoint1(a,LDA,na,det_l)
+end function qmckl_adjugate_f
+    #+end_src
+
+    #+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
@@ -460,9 +465,9 @@ subroutine adjoint1(a,LDA,na,det_l)
   double precision, intent(inout) :: det_l
   
   call cofactor1(a,LDA,na,det_l)
-end
+end subroutine adjugate1
 
-subroutine adjoint2(a,LDA,na,det_l)
+subroutine adjugate2(a,LDA,na,det_l)
   implicit none
   double precision, intent(inout) :: a (LDA,na)
   integer*8, intent(in)             :: LDA
@@ -481,9 +486,9 @@ subroutine adjoint2(a,LDA,na,det_l)
   a(1,2) = b(1,2)
   a(2,1) = b(2,1)
   a(2,2) = b(2,2)
-end
+end subroutine adjugate2
 
-subroutine adjoint3(a,LDA,na,det_l)
+subroutine adjugate3(a,LDA,na,det_l)
   implicit none
   double precision, intent(inout) :: a (LDA,na)
   integer*8, intent(in)             :: LDA
@@ -513,9 +518,9 @@ subroutine adjoint3(a,LDA,na,det_l)
   a(1,3) = b(1,3)
   a(2,3) = b(2,3)
   a(3,3) = b(3,3)
-end
+end subroutine adjugate3
 
-subroutine adjoint4(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
@@ -559,9 +564,9 @@ subroutine adjoint4(a,LDA,na,det_l)
   a(2,4) = b(2,4)
   a(3,4) = b(3,4)
   a(4,4) = b(4,4)
-end
+end subroutine adjugate4
 
-subroutine adjoint5(a,LDA,na,det_l)
+subroutine adjugate5(a,LDA,na,det_l)
   implicit none
   double precision, intent(inout) :: a (LDA,na)
   integer*8, intent(in)             :: LDA
@@ -623,7 +628,7 @@ subroutine adjoint5(a,LDA,na,det_l)
   a(3,5) = b(3,5)
   a(4,5) = b(4,5)
   a(5,5) = b(5,5)
-end
+end subroutine adjugate5
 
 subroutine cofactor1(a,LDA,na,det_l)
   implicit none
@@ -634,7 +639,7 @@ subroutine cofactor1(a,LDA,na,det_l)
   
   det_l = a(1,1)
   a(1,1) = 1.d0
-end
+end subroutine cofactor1
 
 subroutine cofactor2(a,LDA,na,det_l)
   implicit none
@@ -653,7 +658,7 @@ subroutine cofactor2(a,LDA,na,det_l)
   a(2,1) = -b(2,1)
   a(1,2) = -b(1,2)
   a(2,2) = b(1,1)
-end
+end subroutine cofactor2
 
 subroutine cofactor3(a,LDA,na,det_l)
   implicit none
@@ -666,11 +671,13 @@ subroutine cofactor3(a,LDA,na,det_l)
   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)
@@ -683,7 +690,7 @@ subroutine cofactor3(a,LDA,na,det_l)
   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
+end subroutine cofactor3
 
 subroutine cofactor4(a,LDA,na,det_l)
   implicit none
@@ -732,7 +739,7 @@ subroutine cofactor4(a,LDA,na,det_l)
   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
+end subroutine cofactor4
 
 subroutine cofactor5(a,LDA,na,det_l)
   implicit none
@@ -742,6 +749,7 @@ subroutine cofactor5(a,LDA,na,det_l)
   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)- &
@@ -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 &
+    integer(c_int32_t) function qmckl_adjugate &
-	(context, m, n, lda, A, det_l) &
+        (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
+      integer(c_int32_t), external :: qmckl_adjugate_f
-      info = qmckl_adjoint_f &
+      info = qmckl_adjugate_f &
-	     (context, m, n, lda, A, det_l)
+             (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 &
+      integer(c_int32_t) function qmckl_adjugate &
-	  (context, m, n, lda, A, det_l) &
+          (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
 
-      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(:,:)
+  double precision, allocatable :: A(:,:), B(:,:)
-  integer*8                     :: m, n, k, LDA, LDB, LDC
+  integer*8                     :: m, n, k, LDA, LDB
   integer*8                     :: i,j,l
   double precision              :: x, det_l, det_l_ref
 
-  m = 4_8
+  LDA = 6
-  k = 4_8
+  LDB = 6
-  LDA = m
-  LDB = m
-  LDC = m
   
-  allocate( A(LDA,k), C(LDC,k))
+  allocate( A(LDA,6), B(LDB,6))
 
-  A = 0.10d0
+  A = 0.1d0
-  C = 0.d0
   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)
+  test_qmckl_adjugate = QMCKL_SUCCESS
-  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_adjoint = qmckl_adjoint(context, m, k, LDA, A, det_l)
+    #+end_src
     
-  if (test_qmckl_adjoint /= QMCKL_SUCCESS) return
+    #+begin_src python :results output :output drawer
+import numpy as np
+import numpy.linalg as la
+N=6
 
-  test_qmckl_adjoint = QMCKL_FAILURE
+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"  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"")
+  
+
+print ("#+end_src")
+#    print(adj(A[0:i+1,0:i+1]))
+    #+end_src
+
+    #+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
-  do j=1,m
+      x = x + (B(1,1) - (1.0) )**2
-     do i=1,k
+      if (dabs(x / det_l) > 1.d-12) then
-       x = x + (A(i,j) - (C(i,j) * det_l_ref))**2
+          print *, 'N = 1: x = ', x
-     end do
+          test_qmckl_adjugate = 1
-  end do
+          return
+      end if
 
-  if (dabs(x) <= 1.d-15 .and. (dabs(det_l_ref - det_l)) <= 1.d-15) then
+    ! N =  2
-     test_qmckl_adjoint = QMCKL_SUCCESS
+      B(:,:) = A(:,:)
-  endif
+      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
 
-  deallocate(A,C)
+      x = 0.d0
-end function test_qmckl_adjoint
+      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);
+qmckl_exit_code test_qmckl_adjugate(qmckl_context context);
-assert(QMCKL_SUCCESS == test_qmckl_adjoint(context));
+assert(QMCKL_SUCCESS == test_qmckl_adjugate(context));
     #+end_src
     
 * End of files                                                     :noexport:

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

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