1
0
mirror of https://github.com/TREX-CoE/qmckl.git synced 2024-12-23 04:44:03 +01:00

Add qmckl_adjugate_safe proxy function

This commit is contained in:
q-posev 2022-08-22 17:31:15 +02:00
parent 578bca9d52
commit e8c7a69c5a
2 changed files with 156 additions and 5 deletions

View File

@ -1086,7 +1086,7 @@ assert(QMCKL_SUCCESS == test_qmckl_dgemm(context));
** ~qmckl_dgemm_safe~ ** ~qmckl_dgemm_safe~
"Size-safe" is a proxy function with the same functionality as ~qmckl_dgemm~ "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) 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 are required primarily for the Python API, where compatibility with
NumPy arrays implies that sizes of the input and output arrays are provided. NumPy arrays implies that sizes of the input and output arrays are provided.
@ -1622,8 +1622,6 @@ integer function qmckl_adjugate_f(context, na, A, LDA, B, ldb, det_l) &
integer*8, intent(in) :: na integer*8, intent(in) :: na
double precision, intent(inout) :: det_l double precision, intent(inout) :: det_l
integer :: i,j
info = QMCKL_SUCCESS info = QMCKL_SUCCESS
if (context == QMCKL_NULL_CONTEXT) then if (context == QMCKL_NULL_CONTEXT) then
@ -2444,6 +2442,159 @@ qmckl_exit_code test_qmckl_adjugate(qmckl_context context);
assert(QMCKL_SUCCESS == test_qmckl_adjugate(context)); assert(QMCKL_SUCCESS == test_qmckl_adjugate(context));
#+end_src #+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~ ** ~qmckl_transpose~
Transposes a matrix: $A^\dagger_{ji} = A_{ij}$. Transposes a matrix: $A^\dagger_{ji} = A_{ij}$.
@ -2543,9 +2694,7 @@ qmckl_transpose (qmckl_context context,
assert (qmckl_context_destroy(context) == QMCKL_SUCCESS); assert (qmckl_context_destroy(context) == QMCKL_SUCCESS);
return 0; return 0;
} }
#+end_src #+end_src
# -*- mode: org -*- # -*- mode: org -*-
# vim: syntax=c # vim: syntax=c

View File

@ -20,6 +20,7 @@
%include typemaps.i %include typemaps.i
%apply int *OUTPUT { qmckl_exit_code *exit_code}; %apply int *OUTPUT { qmckl_exit_code *exit_code};
%apply double *OUTPUT { double* det_l};
/* Avoid passing file_name length as an additiona argument */ /* Avoid passing file_name length as an additiona argument */
%apply (char *STRING, int LENGTH) { (const char* file_name, const int64_t size_max) }; %apply (char *STRING, int LENGTH) { (const char* file_name, const int64_t size_max) };
@ -55,6 +56,7 @@ import_array();
%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 * 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* 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 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 */ /* Handle properly get_point */