diff --git a/Makefile.am b/Makefile.am index f6d79a7..472ee62 100644 --- a/Makefile.am +++ b/Makefile.am @@ -118,7 +118,7 @@ python-install: $(qmckl_h) $(lib_LTLIBRARIES) $(dist_python_DATA) cp src/qmckl.py . ; \ export QMCKL_INCLUDEDIR="$(prefix)/include" ; \ export QMCKL_LIBDIR="$(prefix)/lib" ; \ - pip install . + python3 -m pip install . python-test: $(test_py) cd $(abs_srcdir)/python/test/ && \ diff --git a/org/qmckl_blas.org b/org/qmckl_blas.org index c21adee..969e943 100644 --- a/org/qmckl_blas.org +++ b/org/qmckl_blas.org @@ -1084,6 +1084,225 @@ qmckl_exit_code test_qmckl_dgemm(qmckl_context context); assert(QMCKL_SUCCESS == test_qmckl_dgemm(context)); #+end_src +** ~qmckl_dgemm_safe~ + + "Size-safe" proxy function with the same functionality as ~qmckl_dgemm~ + but with 3 additional arguments. These arguments ~size_max_M~ (where ~M~ is a matix) + are required primarily for the Python API, where compatibility with + NumPy arrays implies that sizes of the input and output arrays are provided. + + #+NAME: qmckl_dgemm_safe_args + | Variable | Type | In/Out | Description | + |--------------+-----------------+--------+---------------------------------------| + | ~context~ | ~qmckl_context~ | in | Global state | + | ~TransA~ | ~char~ | in | 'T' is transposed | + | ~TransB~ | ~char~ | in | 'T' is transposed | + | ~m~ | ~int64_t~ | in | Number of rows of the input matrix | + | ~n~ | ~int64_t~ | in | Number of columns of the input matrix | + | ~k~ | ~int64_t~ | in | Number of columns of the input matrix | + | ~alpha~ | ~double~ | in | \alpha | + | ~A~ | ~double[][lda]~ | in | Array containing the matrix $A$ | + | ~size_max_A~ | ~int64_t~ | in | Size of the matrix $A$ | + | ~lda~ | ~int64_t~ | in | Leading dimension of array ~A~ | + | ~B~ | ~double[][ldb]~ | in | Array containing the matrix $B$ | + | ~size_max_B~ | ~int64_t~ | in | Size of the matrix $B$ | + | ~ldb~ | ~int64_t~ | in | Leading dimension of array ~B~ | + | ~beta~ | ~double~ | in | \beta | + | ~C~ | ~double[][ldc]~ | out | Array containing the matrix $C$ | + | ~size_max_C~ | ~int64_t~ | in | Size of the matrix $C$ | + | ~ldc~ | ~int64_t~ | in | Leading dimension of array ~C~ | + + 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 + - ~size_max_A >= m * k~ + - ~size_max_B >= k * n~ + - ~size_max_C >= m * n~ + + #+CALL: generate_c_header(table=qmckl_dgemm_safe_args,rettyp="qmckl_exit_code",fname="qmckl_dgemm_safe") + + #+RESULTS: + #+BEGIN_src c :tangle (eval h_func) :comments org + qmckl_exit_code qmckl_dgemm_safe ( + const qmckl_context context, + const char TransA, + const char TransB, + const int64_t m, + const int64_t n, + const int64_t k, + const double alpha, + const double* A, + const int64_t size_max_A, + const int64_t lda, + const double* B, + const int64_t size_max_B, + const int64_t ldb, + const double beta, + double* const C, + const int64_t size_max_C, + const int64_t ldc ); + #+END_src + + #+begin_src f90 :tangle (eval f) :exports none +integer function qmckl_dgemm_safe_f(context, TransA, TransB, & + m, n, k, alpha, A, size_A, LDA, B, size_B, LDB, beta, C, size_C, LDC) & + result(info) + use qmckl + 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,*) + + 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 <= 0) then + info = QMCKL_INVALID_ARG_10 + return + endif + + if (LDB <= 0) then + info = QMCKL_INVALID_ARG_13 + return + endif + + if (LDC <= 0) then + info = QMCKL_INVALID_ARG_17 + return + endif + + if (size_A <= 0) then + info = QMCKL_INVALID_ARG_9 + return + endif + + if (size_B <= 0) then + info = QMCKL_INVALID_ARG_12 + return + endif + + if (size_C <= 0) then + info = QMCKL_INVALID_ARG_16 + return + endif + + 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_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 , intent(in) , value :: TransA + character , 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) & + bind(C) + use, intrinsic :: iso_c_binding + import + implicit none + + integer (c_int64_t) , intent(in) , value :: context + character , intent(in) , value :: TransA + character , 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 + + end function qmckl_dgemm_safe + end interface + #+end_src + ** ~qmckl_matmul~ Matrix multiplication using the =qmckl_matrix= data type: @@ -1403,8 +1622,6 @@ integer function qmckl_adjugate_f(context, na, A, LDA, B, ldb, det_l) & integer*8, intent(in) :: na double precision, intent(inout) :: det_l - integer :: i,j - info = QMCKL_SUCCESS if (context == QMCKL_NULL_CONTEXT) then @@ -2225,6 +2442,159 @@ qmckl_exit_code test_qmckl_adjugate(qmckl_context context); assert(QMCKL_SUCCESS == test_qmckl_adjugate(context)); #+end_src +** ~qmckl_adjugate_safe~ + + "Size-safe" proxy function with the same functionality as ~qmckl_adjugate~ + but with 2 additional arguments. These arguments ~size_max_M~ (where ~M~ is a matix) + are required primarily for the Python API, where compatibility with + NumPy arrays implies that sizes of the input and output arrays are provided. + + + #+NAME: qmckl_adjugate_safe_args + | Variable | Type | In/Out | Description | + |--------------+-----------------+--------+------------------------------------------------| + | ~context~ | ~qmckl_context~ | in | Global state | + | ~n~ | ~int64_t~ | in | Number of rows and columns of the input matrix | + | ~A~ | ~double[][lda]~ | in | Array containing the $n \times n$ matrix $A$ | + | ~size_max_A~ | ~int64_t~ | in | Size of the matrix $A$ | + | ~lda~ | ~int64_t~ | in | Leading dimension of array ~A~ | + | ~B~ | ~double[][ldb]~ | out | Adjugate of $A$ | + | ~size_max_B~ | ~int64_t~ | in | Size of the matrix $B$ | + | ~ldb~ | ~int64_t~ | in | Leading dimension of array ~B~ | + | ~det_l~ | ~double~ | inout | determinant of $A$ | + + Requirements: + + - ~context~ is not ~QMCKL_NULL_CONTEXT~ + - ~n > 0~ + - ~lda >= m~ + - ~A~ is allocated with at least $m \times m \times 8$ bytes + - ~ldb >= m~ + - ~B~ is allocated with at least $m \times m \times 8$ bytes + - ~size_max_A >= m * m~ + - ~size_max_B >= m * m~ + + #+CALL: generate_c_header(table=qmckl_adjugate_safe_args,rettyp="qmckl_exit_code",fname="qmckl_adjugate_safe") + + #+RESULTS: + #+BEGIN_src c :tangle (eval h_func) :comments org + qmckl_exit_code qmckl_adjugate_safe ( + const qmckl_context context, + const int64_t n, + const double* A, + const int64_t size_max_A, + const int64_t lda, + double* const B, + const int64_t size_max_B, + const int64_t ldb, + double* det_l ); + #+END_src + + 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) :exports none +integer function qmckl_adjugate_safe_f(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 + + integer, external :: qmckl_adjugate_f + + info = QMCKL_SUCCESS + + if (size_A < na) then + info = QMCKL_INVALID_ARG_4 + return + endif + + if (size_B <= 0_8) then + info = QMCKL_INVALID_ARG_7 + return + endif + + info = qmckl_adjugate_f(context, na, A, LDA, B, LDB, det_l) + + if (info == QMCKL_INVALID_ARG_4) then + info = QMCKL_INVALID_ARG_5 + return + endif + + if (info == QMCKL_INVALID_ARG_6) then + info = QMCKL_INVALID_ARG_8 + return + endif + +end function qmckl_adjugate_safe_f + #+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) & + 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 :: 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 + + end function qmckl_adjugate_safe + end interface + #+end_src + ** ~qmckl_transpose~ Transposes a matrix: $A^\dagger_{ji} = A_{ij}$. @@ -2324,9 +2694,7 @@ qmckl_transpose (qmckl_context context, assert (qmckl_context_destroy(context) == QMCKL_SUCCESS); return 0; } - #+end_src - # -*- mode: org -*- # vim: syntax=c diff --git a/python/pyproject.toml b/python/pyproject.toml index c66b022..3b9ceaf 100644 --- a/python/pyproject.toml +++ b/python/pyproject.toml @@ -2,6 +2,6 @@ requires = [ "setuptools>=42", "wheel", - "numpy>=1.17.3" + "oldest-supported-numpy" ] build-backend = "setuptools.build_meta" diff --git a/python/setup.py b/python/setup.py index d2abecf..1a58a9e 100644 --- a/python/setup.py +++ b/python/setup.py @@ -71,5 +71,7 @@ setup(name = MODULE_NAME, "Operating System :: MacOS" ], python_requires = ">=3.0", - install_requires = ["numpy>=1.17.3"] + # The ABI incompatibility of numpy is a pain, for now set the + # min numpy version such that we have wheels for CPython 3.5 & 3.6 + install_requires = ["numpy>=1.13.3"] ) diff --git a/python/src/qmckl.i b/python/src/qmckl.i index e8b2d89..5573bdb 100644 --- a/python/src/qmckl.i +++ b/python/src/qmckl.i @@ -20,6 +20,7 @@ %include typemaps.i %apply int *OUTPUT { qmckl_exit_code *exit_code}; +%apply double *OUTPUT { double* det_l}; /* Avoid passing file_name length as an additiona argument */ %apply (char *STRING, int LENGTH) { (const char* file_name, const int64_t size_max) }; @@ -29,6 +30,9 @@ %cstring_bounded_output(char* function_name, 1024); %cstring_bounded_output(char* message, 1024); +%cstring_bounded_output(char* const basis_type, 2); +/* Required for qmckl_last_error function to work */ +%cstring_bounded_output(char* buffer, 1024); /* This block is needed make SWIG convert NumPy arrays to/from from the C pointer and size_max argument. NOTE: `numpy.i` interface file is not part of SWIG but it is included in the numpy distribution (under numpy/tools/swig/numpy.i) @@ -48,6 +52,12 @@ import_array(); /* Include typemaps generated by the process_header.py script */ %include "qmckl_include.i" +/* Some custom array typemaps which are not generated by process_header.py */ +%apply ( double* IN_ARRAY1 , int64_t DIM1 ) { ( const double * A, const int64_t size_max_A) }; +%apply ( double* IN_ARRAY1 , int64_t DIM1 ) { ( const double * B, const int64_t size_max_B) }; +%apply ( double* ARGOUT_ARRAY1 , int64_t DIM1 ) { ( double* const C, const int64_t size_max_C) }; +%apply ( double* ARGOUT_ARRAY1 , int64_t DIM1 ) { ( double* const B, const int64_t size_max_B) }; + /* Handle properly get_point */ diff --git a/python/test/test_api.py b/python/test/test_api.py index ccffd45..58b8476 100644 --- a/python/test/test_api.py +++ b/python/test/test_api.py @@ -33,6 +33,9 @@ assert mo_num == 404 pq.set_electron_coord(ctx, 'T', walk_num, coord) +ao_type = pq.get_ao_basis_type(ctx) +assert 'G' in ao_type + size_max = 5*walk_num*elec_num*mo_num mo_vgl = pq.get_mo_basis_mo_vgl(ctx, size_max)