diff --git a/org/qmckl_blas.org b/org/qmckl_blas.org index fd71021..aae3c5a 100644 --- a/org/qmckl_blas.org +++ b/org/qmckl_blas.org @@ -972,31 +972,41 @@ double* qmckl_alloc_double_of_tensor(const qmckl_context context, const int64_t ldb, const double beta, double* const C, - const int64_t ldc ); + const int64_t ldc ); #+end_src #+begin_src f90 :tangle (eval f) :exports none -integer function qmckl_dgemm_f(context, TransA, TransB, & +function qmckl_dgemm(context, TransA, TransB, & m, n, k, alpha, A, LDA, B, LDB, beta, C, LDC) & - result(info) - use qmckl + bind(C) result(info) + use, intrinsic :: iso_c_binding + use qmckl_constants #ifdef HAVE_LIBQMCKLDGEMM use qmckl_dgemm_tiled_module #endif implicit none - integer(qmckl_context), intent(in) :: context - character , intent(in) :: TransA, TransB - integer*8 , intent(in) :: m, n, k - double precision , intent(in) :: alpha, beta - integer*8 , intent(in) :: lda - double precision , intent(in) :: A(lda,*) - integer*8 , intent(in) :: ldb - double precision , intent(in) :: B(ldb,*) - integer*8 , intent(in) :: ldc - double precision , intent(out) :: C(ldc,*) + integer (qmckl_context), intent(in) , value :: context + character(c_char ) , intent(in) , value :: TransA + character(c_char ) , intent(in) , value :: TransB + integer (c_int64_t) , intent(in) , value :: m + integer (c_int64_t) , intent(in) , value :: n + integer (c_int64_t) , intent(in) , value :: k + real (c_double ) , intent(in) , value :: alpha + real (c_double ) , intent(in) :: A(lda,*) + integer (c_int64_t) , intent(in) , value :: lda + real (c_double ) , intent(in) :: B(ldb,*) + integer (c_int64_t) , intent(in) , value :: ldb + real (c_double ) , intent(in) , value :: beta + real (c_double ) , intent(out) :: C(ldc,*) + integer (c_int64_t) , intent(in) , value :: ldc + integer(qmckl_exit_code) :: info + +#ifdef HAVE_LIBQMCKLDGEMM double precision,allocatable,dimension(:,:) :: A1 double precision,allocatable,dimension(:,:) :: B1 double precision,allocatable,dimension(:,:) :: C1 +#endif + integer*8 :: i, j, LDA1, LDB1, LDC1 info = QMCKL_SUCCESS @@ -1040,25 +1050,25 @@ integer function qmckl_dgemm_f(context, TransA, TransB, & ! Copy A to A1 allocate(A1(k,m)) do j=1,m - do i=1,k - A1(i,j) = A(j,i) - end do + do i=1,k + A1(i,j) = A(j,i) + end do end do ! Copy B to B1 allocate(B1(n,k)) do j=1,k - do i=1,n - B1(i,j) = B(j,i) - end do + do i=1,n + B1(i,j) = B(j,i) + end do end do ! Copy C to C1 allocate(C1(n,m)) do j=1,m - do i=1,n - C1(i,j) = C(j,i) - end do + do i=1,n + C1(i,j) = C(j,i) + end do end do LDA1 = size(A1,1) @@ -1070,7 +1080,7 @@ integer function qmckl_dgemm_f(context, TransA, TransB, & do j=1,n do i=1,m -transpose C(i,j) = alpha*C1(j,i) + beta*C(i,j) + transpose C(i,j) = alpha*C1(j,i) + beta*C(i,j) end do end do @@ -1079,62 +1089,28 @@ transpose C(i,j) = alpha*C1(j,i) + beta*C(i,j) call dgemm(transA, transB, int(m,4), int(n,4), int(k,4), & alpha, A, int(LDA,4), B, int(LDB,4), beta, C, int(LDC,4)) #endif - -end function qmckl_dgemm_f + +end function qmckl_dgemm #+end_src *** C interface :noexport: - #+CALL: generate_c_interface(table=qmckl_dgemm_args,rettyp="qmckl_exit_code",fname="qmckl_dgemm") - - #+RESULTS: - #+begin_src f90 :tangle (eval f) :comments org :exports none - integer(c_int32_t) function qmckl_dgemm & - (context, TransA, TransB, m, n, k, alpha, A, lda, B, ldb, beta, C, ldc) & - bind(C) result(info) - - use, intrinsic :: iso_c_binding - implicit none - - integer (c_int64_t) , intent(in) , value :: context - character(c_char) , intent(in) , value :: TransA - character(c_char) , intent(in) , value :: TransB - integer (c_int64_t) , intent(in) , value :: m - integer (c_int64_t) , intent(in) , value :: n - integer (c_int64_t) , intent(in) , value :: k - real (c_double ) , intent(in) , value :: alpha - real (c_double ) , intent(in) :: A(lda,*) - integer (c_int64_t) , intent(in) , value :: lda - real (c_double ) , intent(in) :: B(ldb,*) - integer (c_int64_t) , intent(in) , value :: ldb - real (c_double ) , intent(in) , value :: beta - real (c_double ) , intent(out) :: C(ldc,*) - integer (c_int64_t) , intent(in) , value :: ldc - - integer(c_int32_t), external :: qmckl_dgemm_f - info = qmckl_dgemm_f & - (context, TransA, TransB, m, n, k, alpha, A, lda, B, ldb, beta, C, ldc) - - end function qmckl_dgemm - #+end_src - - #+CALL: generate_f_interface(table=qmckl_dgemm_args,rettyp="qmckl_exit_code",fname="qmckl_dgemm") #+RESULTS: #+begin_src f90 :tangle (eval fh_func) :comments org :exports none interface - integer(c_int32_t) function qmckl_dgemm & + integer(qmckl_exit_code) function qmckl_dgemm & (context, TransA, TransB, m, n, k, alpha, A, lda, B, ldb, beta, C, ldc) & bind(C) use, intrinsic :: iso_c_binding import implicit none - integer (c_int64_t) , intent(in) , value :: context - character(c_char) , intent(in) , value :: TransA - character(c_char) , intent(in) , value :: TransB + integer (qmckl_context), intent(in) , value :: context + character(c_char ) , intent(in) , value :: TransA + character(c_char ) , 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 @@ -1295,25 +1271,31 @@ printf("qmckl_dgemm ok\n"); #+END_src #+begin_src f90 :tangle (eval f) :exports none -integer function qmckl_dgemm_safe_f(context, TransA, TransB, & +function qmckl_dgemm_safe(context, TransA, TransB, & m, n, k, alpha, A, size_A, LDA, B, size_B, LDB, beta, C, size_C, LDC) & - result(info) - use qmckl + result(info) bind(C) + use, intrinsic :: iso_c_binding + use qmckl_constants implicit none - integer(qmckl_context), intent(in) :: context - character , intent(in) :: TransA, TransB - integer*8 , intent(in) :: m, n, k - double precision , intent(in) :: alpha, beta - integer*8 , intent(in) :: lda - integer*8 , intent(in) :: size_A - double precision , intent(in) :: A(lda,*) - integer*8 , intent(in) :: ldb - integer*8 , intent(in) :: size_B - double precision , intent(in) :: B(ldb,*) - integer*8 , intent(in) :: ldc - integer*8 , intent(in) :: size_C - double precision , intent(out) :: C(ldc,*) + integer (qmckl_context), intent(in) , value :: context + character(c_char ) , intent(in) , value :: TransA + character(c_char ) , intent(in) , value :: TransB + integer (c_int64_t) , intent(in) , value :: m + integer (c_int64_t) , intent(in) , value :: n + integer (c_int64_t) , intent(in) , value :: k + real (c_double ) , intent(in) , value :: alpha + real (c_double ) , intent(in) :: A(lda,*) + integer (c_int64_t) , intent(in) , value :: size_A + integer (c_int64_t) , intent(in) , value :: lda + real (c_double ) , intent(in) :: B(ldb,*) + integer (c_int64_t) , intent(in) , value :: size_B + integer (c_int64_t) , intent(in) , value :: ldb + real (c_double ) , intent(in) , value :: beta + real (c_double ) , intent(out) :: C(ldc,*) + integer (c_int64_t) , intent(in) , value :: size_C + integer (c_int64_t) , intent(in) , value :: ldc + integer(qmckl_exit_code) :: info info = QMCKL_SUCCESS if (context == QMCKL_NULL_CONTEXT) then @@ -1369,76 +1351,39 @@ integer function qmckl_dgemm_safe_f(context, TransA, TransB, & call dgemm(transA, transB, int(m,4), int(n,4), int(k,4), & alpha, A, int(LDA,4), B, int(LDB,4), beta, C, int(LDC,4)) -end function qmckl_dgemm_safe_f +end function qmckl_dgemm_safe #+end_src *** C interface :noexport: - #+CALL: generate_c_interface(table=qmckl_dgemm_safe_args,rettyp="qmckl_exit_code",fname="qmckl_dgemm_safe") - - #+RESULTS: - #+begin_src f90 :tangle (eval f) :comments org :exports none - integer(c_int32_t) function qmckl_dgemm_safe & - (context, TransA, TransB, m, n, k, alpha, A, size_A, lda, B, size_B, ldb, beta, C, size_C, ldc) & - bind(C) result(info) - - use, intrinsic :: iso_c_binding - implicit none - - integer (c_int64_t) , intent(in) , value :: context - character(c_char) , intent(in) , value :: TransA - character(c_char) , intent(in) , value :: TransB - integer (c_int64_t) , intent(in) , value :: m - integer (c_int64_t) , intent(in) , value :: n - integer (c_int64_t) , intent(in) , value :: k - real (c_double ) , intent(in) , value :: alpha - real (c_double ) , intent(in) :: A(lda,*) - integer (c_int64_t) , intent(in) , value :: size_A - integer (c_int64_t) , intent(in) , value :: lda - real (c_double ) , intent(in) :: B(ldb,*) - integer (c_int64_t) , intent(in) , value :: size_B - integer (c_int64_t) , intent(in) , value :: ldb - real (c_double ) , intent(in) , value :: beta - real (c_double ) , intent(out) :: C(ldc,*) - integer (c_int64_t) , intent(in) , value :: size_C - integer (c_int64_t) , intent(in) , value :: ldc - - integer(c_int32_t), external :: qmckl_dgemm_safe_f - info = qmckl_dgemm_safe_f & - (context, TransA, TransB, m, n, k, alpha, A, size_A, lda, B, size_B, ldb, beta, C, size_C, ldc) - - end function qmckl_dgemm_safe - #+end_src - - #+CALL: generate_f_interface(table=qmckl_dgemm_safe_args,rettyp="qmckl_exit_code",fname="qmckl_dgemm_safe") #+RESULTS: #+begin_src f90 :tangle (eval fh_func) :comments org :exports none interface - integer(c_int32_t) function qmckl_dgemm_safe & - (context, TransA, TransB, m, n, k, alpha, A, size_A, lda, B, size_B, ldb, beta, C, size_C, ldc) & + integer(qmckl_exit_code) function qmckl_dgemm_safe & + (context, TransA, TransB, m, n, k, alpha, A, size_max_A, lda, B, size_max_B, ldb, beta, C, size_max_C, ldc) & bind(C) use, intrinsic :: iso_c_binding import implicit none - integer (c_int64_t) , intent(in) , value :: context - character(c_char) , intent(in) , value :: TransA - character(c_char) , intent(in) , value :: TransB + integer (qmckl_context), intent(in) , value :: context + character(c_char ) , intent(in) , value :: TransA + character(c_char ) , intent(in) , value :: TransB integer (c_int64_t) , intent(in) , value :: m integer (c_int64_t) , intent(in) , value :: n integer (c_int64_t) , intent(in) , value :: k real (c_double ) , intent(in) , value :: alpha real (c_double ) , intent(in) :: A(lda,*) - integer (c_int64_t) , intent(in) , value :: size_A + integer (c_int64_t) , intent(in) , value :: size_max_A integer (c_int64_t) , intent(in) , value :: lda real (c_double ) , intent(in) :: B(ldb,*) - integer (c_int64_t) , intent(in) , value :: size_B + integer (c_int64_t) , intent(in) , value :: size_max_B integer (c_int64_t) , intent(in) , value :: ldb real (c_double ) , intent(in) , value :: beta real (c_double ) , intent(out) :: C(ldc,*) - integer (c_int64_t) , intent(in) , value :: size_C + integer (c_int64_t) , intent(in) , value :: size_max_C integer (c_int64_t) , intent(in) , value :: ldc end function qmckl_dgemm_safe @@ -1756,17 +1701,20 @@ print(C.T) LAPACK library. #+begin_src f90 :tangle (eval f) :exports none -integer function qmckl_adjugate_f(context, na, A, LDA, B, ldb, det_l) & - result(info) - use qmckl +function qmckl_adjugate(context, n, A, LDA, B, ldb, det_l) & + result(info) bind(C) + use qmckl_constants + use, intrinsic :: iso_c_binding implicit none - integer(qmckl_context) , intent(in) :: context - double precision, intent(in) :: A (LDA,*) - integer*8, intent(in) :: LDA - double precision, intent(out) :: B (LDB,*) - integer*8, intent(in) :: LDB - integer*8, intent(in) :: na - double precision, intent(inout) :: det_l + + integer (qmckl_context), intent(in) , value :: context + integer (c_int64_t) , intent(in) , value :: n + real (c_double ) , intent(in) :: A(lda,*) + integer (c_int64_t) , intent(in) , value :: lda + real (c_double ) , intent(out) :: B(ldb,*) + integer (c_int64_t) , intent(in) , value :: ldb + real (c_double ) , intent(inout) :: det_l + integer(qmckl_exit_code) :: info info = QMCKL_SUCCESS @@ -1775,7 +1723,7 @@ integer function qmckl_adjugate_f(context, na, A, LDA, B, ldb, det_l) & return endif - if (na <= 0_8) then + if (n <= 0_8) then info = QMCKL_INVALID_ARG_2 return endif @@ -1785,28 +1733,28 @@ integer function qmckl_adjugate_f(context, na, A, LDA, B, ldb, det_l) & return endif - if (LDA < na) then + if (LDA < n) then info = QMCKL_INVALID_ARG_4 return endif - select case (na) + select case (n) case (5) - call adjugate5(A,LDA,B,LDB,na,det_l) + call adjugate5(A,LDA,B,LDB,n,det_l) case (4) - call adjugate4(A,LDA,B,LDB,na,det_l) + call adjugate4(A,LDA,B,LDB,n,det_l) case (3) - call adjugate3(A,LDA,B,LDB,na,det_l) + call adjugate3(A,LDA,B,LDB,n,det_l) case (2) - call adjugate2(A,LDA,B,LDB,na,det_l) + call adjugate2(A,LDA,B,LDB,n,det_l) case (1) - det_l = a(1,1) - b(1,1) = 1.d0 + det_l = a(1,1) + b(1,1) = 1.d0 case default - call adjugate_general(context, na, A, LDA, B, LDB, det_l) + call adjugate_general(context, n, A, LDA, B, LDB, det_l) end select -end function qmckl_adjugate_f +end function qmckl_adjugate #+end_src #+begin_src f90 :tangle (eval f) :exports none @@ -2188,45 +2136,19 @@ subroutine cofactor5(A,LDA,B,LDB,na,det_l) end #+end_src - #+CALL: generate_c_interface(table=qmckl_adjugate_args,rettyp="qmckl_exit_code",fname="qmckl_adjugate") - - #+RESULTS: - #+begin_src f90 :tangle (eval f) :comments org :exports none - integer(c_int32_t) function qmckl_adjugate & - (context, n, A, lda, B, ldb, det_l) & - bind(C) result(info) - - use, intrinsic :: iso_c_binding - implicit none - - integer (c_int64_t) , intent(in) , value :: context - integer (c_int64_t) , intent(in) , value :: n - real (c_double ) , intent(in) :: A(lda,*) - integer (c_int64_t) , intent(in) , value :: lda - real (c_double ) , intent(out) :: B(ldb,*) - integer (c_int64_t) , intent(in) , value :: ldb - real (c_double ) , intent(inout) :: det_l - - integer(c_int32_t), external :: qmckl_adjugate_f - info = qmckl_adjugate_f & - (context, n, A, lda, B, ldb, det_l) - - end function qmckl_adjugate - #+end_src - #+CALL: generate_f_interface(table=qmckl_adjugate_args,rettyp="qmckl_exit_code",fname="qmckl_adjugate") #+RESULTS: #+begin_src f90 :tangle (eval fh_func) :comments org :exports none interface - integer(c_int32_t) function qmckl_adjugate & + integer(qmckl_exit_code) function qmckl_adjugate & (context, n, A, lda, B, ldb, det_l) & bind(C) use, intrinsic :: iso_c_binding import implicit none - integer (c_int64_t) , intent(in) , value :: context + integer (qmckl_context), intent(in) , value :: context integer (c_int64_t) , intent(in) , value :: n real (c_double ) , intent(in) :: A(lda,*) integer (c_int64_t) , intent(in) , value :: lda @@ -2315,7 +2237,7 @@ end subroutine adjugate_general #+begin_src f90 :tangle (eval f_test) integer(qmckl_exit_code) function test_qmckl_adjugate(context) bind(C) - use qmckl + use qmckl_constants implicit none integer(qmckl_context), intent(in), value :: context @@ -2642,22 +2564,26 @@ printf("qmckl_adjugate ok\n"); LAPACK library. #+begin_src f90 :tangle (eval f) :exports none -integer function qmckl_adjugate_safe_f(context, & +function qmckl_adjugate_safe(context, & na, A, size_A, LDA, B, size_B, LDB, det_l) & - result(info) - use qmckl - implicit none - integer(qmckl_context) , intent(in) :: context - double precision, intent(in) :: A (LDA,*) - integer*8, intent(in) :: size_A - integer*8, intent(in) :: LDA - double precision, intent(out) :: B (LDB,*) - integer*8, intent(in) :: size_B - integer*8, intent(in) :: LDB - integer*8, intent(in) :: na - double precision, intent(inout) :: det_l + result(info) bind(C) + use, intrinsic :: iso_c_binding + use qmckl_constants + use qmckl, only: qmckl_adjugate - integer, external :: qmckl_adjugate_f + implicit none + + integer (qmckl_context), intent(in) , value :: context + integer (c_int64_t) , intent(in) , value :: na + real (c_double ) , intent(in) :: A(lda,*) + integer (c_int64_t) , intent(in) , value :: size_A + integer (c_int64_t) , intent(in) , value :: lda + real (c_double ) , intent(out) :: B(ldb,*) + integer (c_int64_t) , intent(in) , value :: size_B + integer (c_int64_t) , intent(in) , value :: ldb + real (c_double ) , intent(inout) :: det_l + + integer(qmckl_exit_code) :: info info = QMCKL_SUCCESS @@ -2671,7 +2597,7 @@ integer function qmckl_adjugate_safe_f(context, & return endif - info = qmckl_adjugate_f(context, na, A, LDA, B, LDB, det_l) + info = qmckl_adjugate(context, na, A, LDA, B, LDB, det_l) if (info == QMCKL_INVALID_ARG_4) then info = QMCKL_INVALID_ARG_5 @@ -2683,59 +2609,31 @@ integer function qmckl_adjugate_safe_f(context, & return endif -end function qmckl_adjugate_safe_f +end function qmckl_adjugate_safe #+end_src *** C interface - #+CALL: generate_c_interface(table=qmckl_adjugate_safe_args,rettyp="qmckl_exit_code",fname="qmckl_adjugate_safe") - - #+RESULTS: - #+begin_src f90 :tangle (eval f) :comments org :exports none - integer(c_int32_t) function qmckl_adjugate_safe & - (context, n, A, size_A, lda, B, size_B, ldb, det_l) & - bind(C) result(info) - - use, intrinsic :: iso_c_binding - implicit none - - integer (c_int64_t) , intent(in) , value :: context - integer (c_int64_t) , intent(in) , value :: n - real (c_double ) , intent(in) :: A(lda,*) - integer (c_int64_t) , intent(in) , value :: lda - integer (c_int64_t) , intent(in) , value :: size_A - real (c_double ) , intent(out) :: B(ldb,*) - integer (c_int64_t) , intent(in) , value :: ldb - integer (c_int64_t) , intent(in) , value :: size_B - real (c_double ) , intent(inout) :: det_l - - integer(c_int32_t), external :: qmckl_adjugate_safe_f - info = qmckl_adjugate_safe_f & - (context, n, A, size_A, lda, B, size_B, ldb, det_l) - - end function qmckl_adjugate_safe - #+end_src - #+CALL: generate_f_interface(table=qmckl_adjugate_safe_args,rettyp="qmckl_exit_code",fname="qmckl_adjugate_safe") #+RESULTS: #+begin_src f90 :tangle (eval fh_func) :comments org :exports none interface - integer(c_int32_t) function qmckl_adjugate_safe & - (context, n, A, size_A, lda, B, size_B, ldb, det_l) & + integer(qmckl_exit_code) function qmckl_adjugate_safe & + (context, n, A, size_max_A, lda, B, size_max_B, ldb, det_l) & bind(C) use, intrinsic :: iso_c_binding import implicit none - integer (c_int64_t) , intent(in) , value :: context + integer (qmckl_context), intent(in) , value :: context integer (c_int64_t) , intent(in) , value :: n real (c_double ) , intent(in) :: A(lda,*) + integer (c_int64_t) , intent(in) , value :: size_max_A integer (c_int64_t) , intent(in) , value :: lda - integer (c_int64_t) , intent(in) , value :: size_A real (c_double ) , intent(out) :: B(ldb,*) + integer (c_int64_t) , intent(in) , value :: size_max_B integer (c_int64_t) , intent(in) , value :: ldb - integer (c_int64_t) , intent(in) , value :: size_B real (c_double ) , intent(inout) :: det_l end function qmckl_adjugate_safe