From a41c67b94fe46cea7b457fb4da07d534f9488adf Mon Sep 17 00:00:00 2001 From: v1j4y Date: Tue, 7 Sep 2021 16:36:26 +0200 Subject: [PATCH] Working on qmckl_dgemm. --- org/qmckl_ao.org | 2 + org/qmckl_blas.org | 287 ++++++++++++++++++++++++++++++++++++++++++ org/table_of_contents | 1 + tools/lib.org | 2 + 4 files changed, 292 insertions(+) create mode 100644 org/qmckl_blas.org diff --git a/org/qmckl_ao.org b/org/qmckl_ao.org index f78350c..011c015 100644 --- a/org/qmckl_ao.org +++ b/org/qmckl_ao.org @@ -135,6 +135,7 @@ int main() { For H_2 with the following basis set, + #+NAME: basis #+BEGIN_EXAMPLE HYDROGEN S 5 @@ -157,6 +158,7 @@ D 1 we have: + #+NAME: params #+BEGIN_EXAMPLE type = 'G' shell_num = 12 diff --git a/org/qmckl_blas.org b/org/qmckl_blas.org new file mode 100644 index 0000000..0869a93 --- /dev/null +++ b/org/qmckl_blas.org @@ -0,0 +1,287 @@ +#+TITLE: BLAS functions +#+SETUPFILE: ../tools/theme.setup +#+INCLUDE: ../tools/lib.org + +* Headers :noexport: + #+begin_src elisp :noexport :results none +(org-babel-lob-ingest "../tools/lib.org") +#+end_src + + #+begin_src c :comments link :tangle (eval c_test) :noweb yes +#include "qmckl.h" +#include "assert.h" +#ifdef HAVE_CONFIG_H +#include "config.h" +#endif +int main() { + qmckl_context context; + context = qmckl_context_create(); + + #+end_src + +* Matrix operations + +** ~qmckl_dgemm~ + + Matrix multiply l$C_{ij} = \beta C_{ij} + \alpha \sum_{k} A_{ik} \cdot B_{kj}$ using Fortran ~matmul~ function. + + #+NAME: qmckl_dgemm_args + | qmckl_context | context | in | Global state | + | bool | TransA | in | Number of rows of the input matrix | + | bool | TransB | in | Number of rows of the input matrix | + | int64_t | m | in | Number of rows of the input matrix | + | int64_t | n | in | Number of columns of the input matrix | + | int64_t | k | in | Number of columns of the input matrix | + | double | alpha | in | Number of columns of the input matrix | + | double | A[][lda] | in | Array containing the $m \times n$ matrix $A$ | + | int64_t | lda | in | Leading dimension of array ~A~ | + | double | B[][ldb] | in | Array containing the $n \times m$ matrix $B$ | + | int64_t | ldb | in | Leading dimension of array ~B~ | + | double | beta | in | Array containing the $n \times m$ matrix $B$ | + | double | C[][ldc] | out | Array containing the $n \times m$ matrix $B$ | + | int64_t | ldc | in | Leading dimension of array ~B~ | + +*** Requirements + + - ~context~ is not ~QMCKL_NULL_CONTEXT~ + - ~m > 0~ + - ~n > 0~ + - ~k > 0~ + - ~lda >= m~ + - ~ldb >= n~ + - ~ldc >= n~ + - ~A~ is allocated with at least $m \times k \times 8$ bytes + - ~B~ is allocated with at least $k \times n \times 8$ bytes + - ~C~ is allocated with at least $m \times n \times 8$ bytes + +*** C header + + #+CALL: generate_c_header(table=qmckl_dgemm_args,rettyp="qmckl_exit_code",fname="qmckl_dgemm") + + #+RESULTS: + #+BEGIN_src c :tangle (eval h_func) :comments org + qmckl_exit_code qmckl_dgemm ( + const qmckl_context context, + const bool TransA, + const bool TransB, + const int64_t m, + const int64_t n, + const int64_t k, + const double alpha, + const double* A, + const int64_t lda, + const double* B, + const int64_t ldb, + const double beta, + double* const C, + const int64_t ldc ); + #+END_src + + +*** Source + #+begin_src f90 :tangle (eval f) +integer function qmckl_dgemm_f(context, TransA, TransB, m, n, k, alpha, A, LDA, B, LDB, beta, C, LDC) & + result(info) + use qmckl + implicit none + integer(qmckl_context) , intent(in) :: context + logical*8 , intent(in) :: TransA, TransN + integer*8 , intent(in) :: m, n, k + 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 + + info = QMCKL_SUCCESS + + if (context == QMCKL_NULL_CONTEXT) then + info = QMCKL_INVALID_CONTEXT + return + endif + + if (m <= 0_8) then + info = QMCKL_INVALID_ARG_4 + return + endif + + if (n <= 0_8) then + info = QMCKL_INVALID_ARG_5 + return + endif + + if (k <= 0_8) then + info = QMCKL_INVALID_ARG_6 + return + endif + + if (LDA < m) then + info = QMCKL_INVALID_ARG_9 + return + endif + + if (LDB < n) then + info = QMCKL_INVALID_ARG_10 + return + endif + + if (LDB < n) then + info = QMCKL_INVALID_ARG_13 + return + endif + + C = matmul(A,B) +end function qmckl_dgemm_f + #+end_src + +*** C interface :noexport: + + #+CALL: generate_c_interface(table=qmckl_dgemm_args,rettyp="qmckl_exit_code",fname="qmckl_dgemm") + + #+RESULTS: + #+BEGIN_src f90 :tangle (eval f) :comments org :exports none + integer(c_int32_t) function qmckl_dgemm & + (context, TransA, TransB, m, n, k, alpha, A, lda, B, ldb, beta, C, ldc) & + bind(C) result(info) + + use, intrinsic :: iso_c_binding + implicit none + + integer (c_int64_t) , intent(in) , value :: context + logical*8 , intent(in) , value :: TransA + logical*8 , intent(in) , value :: TransB + integer (c_int64_t) , intent(in) , value :: m + integer (c_int64_t) , intent(in) , value :: n + integer (c_int64_t) , intent(in) , value :: k + real (c_double ) , intent(in) , value :: alpha + real (c_double ) , intent(in) :: A(lda,*) + integer (c_int64_t) , intent(in) , value :: lda + real (c_double ) , intent(in) :: B(ldb,*) + integer (c_int64_t) , intent(in) , value :: ldb + real (c_double ) , intent(in) , value :: beta + real (c_double ) , intent(out) :: C(ldc,*) + integer (c_int64_t) , intent(in) , value :: ldc + + integer(c_int32_t), external :: qmckl_dgemm_f + info = qmckl_dgemm_f & + (context, TransA, TransB, m, n, k, alpha, A, lda, B, ldb, beta, C, ldc) + + end function qmckl_dgemm + #+END_src + + + #+CALL: generate_f_interface(table=qmckl_dgemm_args,rettyp="qmckl_exit_code",fname="qmckl_dgemm") + + #+RESULTS: + #+begin_src f90 :tangle (eval fh_func) :comments org :exports none + interface + integer(c_int32_t) function qmckl_dgemm & + (context, TransA, TransB, m, n, k, alpha, A, lda, B, ldb, beat, C, ldc) & + bind(C) + use, intrinsic :: iso_c_binding + import + implicit none + + integer (c_int64_t) , intent(in) , value :: context + logical*8 (c_bool) , intent(in) , value :: TransA + logical*8 (c_bool) , 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) :: 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) :: beta + real (c_double ) , intent(out) :: C(ldb,*) + integer (c_int64_t) , intent(in) , value :: ldb + + 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 + logical*8 :: TransA, TransB + double precision :: x + + TransA = .False. + TransB = .False. + m = 5 + k = 4 + n = 6 + LDA = m + LDB = k + LDC = m + + allocate( A(LDA,k), B(LDB,n) , C(LDC,n), D(LDC,n)) + + A = 0.d0 + B = 0.d0 + C = 0.d0 + D = 0.d0 + do j=1,m + do i=1,k + A(i,j) = -10.d0 + dble(i+j) + end do + end do + + do j=1,k + do i=1,n + 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,m + do i=l,n + do l=1,k + D(i,j) += D(i,j) + (A(i,k)*B(k,j)) + end do + x = x + (D(i,j) - C(i,j))**2 + end do + end do + + if (dabs(x) <= 1.d-15) then + test_qmckl_dgemm = QMCKL_SUCCESS + endif + + deallocate(A,B) +end function test_qmckl_dgemm + #+end_src + + #+begin_src c :comments link :tangle (eval c_test) +qmckl_exit_code test_qmckl_dgemm(qmckl_context context); +assert(QMCKL_SUCCESS == test_qmckl_dgemm(context)); + #+end_src + +* End of files :noexport: + + #+begin_src c :comments link :tangle (eval c_test) + assert (qmckl_context_destroy(context) == QMCKL_SUCCESS); + return 0; +} + + #+end_src + + +# -*- mode: org -*- +# vim: syntax=c diff --git a/org/table_of_contents b/org/table_of_contents index 50a6159..01bb7e5 100644 --- a/org/table_of_contents +++ b/org/table_of_contents @@ -9,4 +9,5 @@ qmckl_ao.org qmckl_jastrow.org qmckl_distance.org qmckl_utils.org +qmckl_blas.org qmckl_tests.org diff --git a/tools/lib.org b/tools/lib.org index d8e642a..3fe9fe2 100644 --- a/tools/lib.org +++ b/tools/lib.org @@ -40,6 +40,7 @@ f_of_c_d = { '' : '' , 'qmckl_context' : 'integer (c_int64_t)' , 'qmckl_exit_code' : 'integer (c_int32_t)' + , 'bool' : 'logical*8' , 'int32_t' : 'integer (c_int32_t)' , 'int64_t' : 'integer (c_int64_t)' , 'float' : 'real (c_float )' @@ -58,6 +59,7 @@ ctypeid_d = { '' : '' , 'real' : 'real(c_float)' , 'real*8' : 'real(c_double)' , 'character' : 'character(c_char)' + , 'character' : 'character(c_char)' } #+END_SRC