1
0
mirror of https://github.com/TREX-CoE/qmckl.git synced 2024-06-18 11:15:38 +02:00

Included auto-generation of Fortran in distances

This commit is contained in:
Anthony Scemama 2021-03-20 16:56:22 +01:00
parent b1f1843e01
commit be44127dba
6 changed files with 399 additions and 131 deletions

View File

@ -87,17 +87,17 @@ endif
.PHONY: clean
.SECONDARY: # Needed to keep the produced C and Fortran files
libqmckl.so: Makefile.generated
libqmckl.so: ../include/qmckl.h
$(MAKE) -f Makefile.generated
../include/qmckl.h: libqmckl.so
../include/qmckl.h: Makefile.generated
../tools/build_qmckl_h.sh
test: Makefile.generated ../include/qmckl.h
test: libqmckl.so
$(MAKE) -f Makefile.generated test
doc: $(ORG_SOURCE_FILES)
$(QMCKL_ROOT)/tools/build_doc.sh
$(QMCKL_ROOT)/tools/create_doc.sh
clean:
$(RM) test_qmckl_* test_qmckl.c test_qmckl \

View File

@ -195,20 +195,6 @@ qmckl_context qmckl_context_create() {
const qmckl_context context = (qmckl_context) ctx;
assert ( qmckl_context_check(context) != QMCKL_NULL_CONTEXT );
/*
qmckl_memory_struct* alloc = (qmckl_memory_struct*)
malloc(sizeof(qmckl_memory_struct));
if (alloc == NULL) {
qmckl_unlock(context);
return QMCKL_NULL_CONTEXT;
}
memset(alloc,0,sizeof(qmckl_memory_struct));
ctx->alloc = alloc;
*/
return context;
}
#+end_src
@ -880,48 +866,48 @@ qmckl_exit_code qmckl_context_update_precision(const qmckl_context context, cons
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT)
return QMCKL_INVALID_CONTEXT;
if (precision < 2) {
return qmckl_failwith(context,
QMCKL_INVALID_ARG_2,
"qmckl_context_update_precision",
"precision < 2");
}
if (precision > 53) {
return qmckl_failwith(context,
QMCKL_INVALID_ARG_2,
"qmckl_context_update_precision",
"precision > 53");
}
qmckl_context_struct* ctx = (qmckl_context_struct*) context;
/* This should be always true */
assert (ctx != NULL);
qmckl_lock(context);
if (ctx->fp == NULL) {
ctx->fp = (qmckl_precision_struct*)
qmckl_malloc(context, sizeof(qmckl_precision_struct));
if (ctx->fp == NULL) {
if (precision < 2) {
return qmckl_failwith(context,
QMCKL_ALLOCATION_FAILED,
"qmckl_context_update_precision",
"ctx->fp");
QMCKL_INVALID_ARG_2,
"qmckl_context_update_precision",
"precision < 2");
}
ctx->fp->precision = QMCKL_DEFAULT_PRECISION;
ctx->fp->range = QMCKL_DEFAULT_RANGE;
}
if (precision > 53) {
return qmckl_failwith(context,
QMCKL_INVALID_ARG_2,
"qmckl_context_update_precision",
"precision > 53");
}
ctx->fp->precision = precision;
qmckl_context_struct* ctx = (qmckl_context_struct*) context;
qmckl_unlock(context);
/* This should be always true */
assert (ctx != NULL);
return QMCKL_SUCCESS;
qmckl_lock(context);
if (ctx->fp == NULL) {
ctx->fp = (qmckl_precision_struct*)
qmckl_malloc(context, sizeof(qmckl_precision_struct));
if (ctx->fp == NULL) {
return qmckl_failwith(context,
QMCKL_ALLOCATION_FAILED,
"qmckl_context_update_precision",
"ctx->fp");
}
ctx->fp->precision = QMCKL_DEFAULT_PRECISION;
ctx->fp->range = QMCKL_DEFAULT_RANGE;
}
ctx->fp->precision = precision;
qmckl_unlock(context);
return QMCKL_SUCCESS;
}
#+end_src
@ -1296,6 +1282,7 @@ qmckl_context_update_ao_basis(qmckl_context context , const char type
/* Memory allocations */
assert (basis->shell_center == NULL);
basis->shell_center = (int64_t*) qmckl_malloc (context, shell_num * sizeof(int64_t));
if (basis->shell_center == NULL) {
qmckl_free(context, basis);
@ -1303,6 +1290,7 @@ qmckl_context_update_ao_basis(qmckl_context context , const char type
return QMCKL_FAILURE;
}
assert (basis->shell_ang_mom == NULL);
basis->shell_ang_mom = (int32_t*) qmckl_malloc (context, shell_num * sizeof(int32_t));
if (basis->shell_ang_mom == NULL) {
qmckl_free(context, basis->shell_center);
@ -1312,6 +1300,7 @@ qmckl_context_update_ao_basis(qmckl_context context , const char type
return QMCKL_FAILURE;
}
assert (basis->shell_prim_num == NULL);
basis->shell_prim_num= (int64_t*) qmckl_malloc (context, shell_num * sizeof(int64_t));
if (basis->shell_prim_num == NULL) {
qmckl_free(context, basis->shell_ang_mom);
@ -1323,7 +1312,8 @@ qmckl_context_update_ao_basis(qmckl_context context , const char type
return QMCKL_FAILURE;
}
basis->shell_factor = (double *) qmckl_malloc (context, shell_num * sizeof(double ));
assert (basis->shell_factor == NULL);
basis->shell_factor = (double *) qmckl_malloc (context, shell_num * sizeof(double));
if (basis->shell_factor == NULL) {
qmckl_free(context, basis->shell_prim_num);
basis->shell_prim_num = NULL;
@ -1336,7 +1326,8 @@ qmckl_context_update_ao_basis(qmckl_context context , const char type
return QMCKL_FAILURE;
}
basis->exponent = (double *) qmckl_malloc (context, prim_num * sizeof(double ));
assert (basis->exponent == NULL);
basis->exponent = (double *) qmckl_malloc (context, prim_num * sizeof(double));
if (basis->exponent == NULL) {
qmckl_free(context, basis->shell_factor);
basis->shell_factor = NULL;
@ -1351,7 +1342,8 @@ qmckl_context_update_ao_basis(qmckl_context context , const char type
return QMCKL_FAILURE;
}
basis->coefficient = (double *) qmckl_malloc (context, prim_num * sizeof(double ));
assert (basis->coefficient == NULL);
basis->coefficient = (double *) qmckl_malloc (context, prim_num * sizeof(double));
if (basis->coefficient == NULL) {
qmckl_free(context, basis->exponent);
basis->exponent = NULL;

View File

@ -1,4 +1,4 @@
#+TITLE: Distances
#+TITLE: Inter-particle distances
#+SETUPFILE: ../docs/theme.setup
Functions for the computation of distances between particles.
@ -10,6 +10,11 @@ Functions for the computation of distances between particles.
(file-name-nondirectory (substring buffer-file-name 0 -4))
#+end_src
#+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 "munit.h"
@ -21,51 +26,65 @@ MunitResult test_<<filename()>>() {
* Squared distance
~qmckl_distance_sq~ computes the matrix of the squared distances
** ~qmckl_distance_sq~
:PROPERTIES:
:Name: qmckl_distance_sq
:CRetType: qmckl_exit_code
:FRetType: integer
:END:
~qmckl_distance_sq~ computes the matrix of the squared distances
between all pairs of points in two sets, one point within each set:
\[
C_{ij} = \sum_{k=1}^3 (A_{k,i}-B_{k,j})^2
\]
| ~context~ | input | Global state |
| ~transa~ | input | Array ~A~ is ~N~: Normal, ~T~: Transposed |
| ~transb~ | input | Array ~B~ is ~N~: Normal, ~T~: Transposed |
| ~m~ | input | Number of points in the first set |
| ~n~ | input | Number of points in the second set |
| ~A(lda,3)~ | input | Array containing the $m \times 3$ matrix $A$ |
| ~lda~ | input | Leading dimension of array ~A~ |
| ~B(ldb,3)~ | input | Array containing the $n \times 3$ matrix $B$ |
| ~ldb~ | input | Leading dimension of array ~B~ |
| ~C(ldc,n)~ | output | Array containing the $m \times n$ matrix $C$ |
| ~ldc~ | input | Leading dimension of array ~C~ |
#+NAME: qmckl_distance_sq_args
| qmckl_context | context | in | Global state |
| char | transa | in | Array ~A~ is ~'N'~: Normal, ~'T'~: Transposed |
| char | transb | in | Array ~B~ is ~'N'~: Normal, ~'T'~: Transposed |
| int64_t | m | in | Number of points in the first set |
| int64_t | n | in | Number of points in the second set |
| double | A[3][lda] | in | Array containing the $m \times 3$ matrix $A$ |
| int64_t | lda | in | Leading dimension of array ~A~ |
| double | B[3][ldb] | in | Array containing the $n \times 3$ matrix $B$ |
| int64_t | ldb | in | Leading dimension of array ~B~ |
| double | C[n][ldc] | out | Array containing the $m \times n$ matrix $C$ |
| int64_t | ldc | in | Leading dimension of array ~C~ |
*** Requirements
- ~context~ is not 0
- ~m~ > 0
- ~n~ > 0
- ~lda~ >= 3 if ~transa~ is ~N~
- ~lda~ >= m if ~transa~ is ~T~
- ~ldb~ >= 3 if ~transb~ is ~N~
- ~ldb~ >= n if ~transb~ is ~T~
- ~ldc~ >= m
- ~context~ is not ~QMCKL_NULL_CONTEXT~
- ~m > 0~
- ~n > 0~
- ~lda >= 3~ if ~transa == 'N'~
- ~lda >= m~ if ~transa == 'T'~
- ~ldb >= 3~ if ~transb == 'N'~
- ~ldb >= n~ if ~transb == 'T'~
- ~ldc >= m~
- ~A~ is allocated with at least $3 \times m \times 8$ bytes
- ~B~ is allocated with at least $3 \times n \times 8$ bytes
- ~C~ is allocated with at least $m \times n \times 8$ bytes
*** C header
*** Performance
#+CALL: generate_c_header(table=qmckl_distance_sq_args,rettyp=get_value("CRetType"),fname=get_value("Name"))
This function might be more efficient when ~A~ and ~B~ are
transposed.
#+begin_src c :comments org :tangle (eval h)
qmckl_exit_code qmckl_distance_sq(const qmckl_context context,
const char transa, const char transb,
const int64_t m, const int64_t n,
const double *A, const int64_t lda,
const double *B, const int64_t ldb,
const double *C, const int64_t ldc);
#+RESULTS:
#+begin_src c :tangle (eval h) :comments org
qmckl_exit_code qmckl_distance_sq (
const qmckl_context context,
const char transa,
const char transb,
const int64_t m,
const int64_t n,
const double* A,
const int64_t lda,
const double* B,
const int64_t ldb,
double* C,
const int64_t ldc );
#+end_src
*** Source
@ -147,7 +166,7 @@ integer function qmckl_distance_sq_f(context, transa, transb, m, n, A, LDA, B, L
select case (transab)
case(0)
do j=1,n
@ -193,49 +212,73 @@ integer function qmckl_distance_sq_f(context, transa, transb, m, n, A, LDA, B, L
end do
end select
end function qmckl_distance_sq_f
#+end_src
*** Performance
This function might be more efficient when ~A~ and ~B~ are
transposed.
*** C interface :noexport:
#+begin_src f90 :tangle (eval f)
integer(c_int32_t) function qmckl_distance_sq(context, transa, transb, m, n, A, LDA, B, LDB, 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, transb
integer (c_int64_t) , intent(in) , value :: m, n
integer (c_int64_t) , intent(in) , value :: lda
real (c_double) , intent(in) :: A(lda,3)
integer (c_int64_t) , intent(in) , value :: ldb
real (c_double) , intent(in) :: B(ldb,3)
integer (c_int64_t) , intent(in) , value :: ldc
real (c_double) , intent(out) :: C(ldc,n)
integer, external :: qmckl_distance_sq_f
info = qmckl_distance_sq_f(context, transa, transb, m, n, A, LDA, B, LDB, C, LDC)
end function qmckl_distance_sq
#+end_src
#+CALL: generate_c_interface(table=qmckl_distance_sq_args,rettyp=get_value("FRetType"),fname=get_value("Name"))
#+begin_src f90 :tangle (eval fh)
interface
integer(c_int32_t) function qmckl_distance_sq(context, transa, transb, m, n, A, LDA, B, LDB, C, LDC) &
bind(C)
#+RESULTS:
#+begin_src f90 :tangle (eval f) :comments org :exports none
integer(c_int32_t) function qmckl_distance_sq &
(context, transa, transb, m, n, A, lda, B, ldb, C, ldc) &
bind(C) result(info)
use, intrinsic :: iso_c_binding
implicit none
integer (c_int64_t) , intent(in) :: context
character , intent(in) :: transa
character , intent(in) :: transb
integer (c_int64_t) , intent(in) :: m
integer (c_int64_t) , intent(in) :: n
real (c_double ) , intent(in) :: A(lda,3)
integer (c_int64_t) , intent(in) :: lda
real (c_double ) , intent(in) :: B(ldb,3)
integer (c_int64_t) , intent(in) :: ldb
real (c_double ) , intent(out) :: C(ldc,n)
integer (c_int64_t) , intent(in) :: ldc
integer(c_int32_t), external :: qmckl_distance_sq_f
info = qmckl_distance_sq_f &
(context, transa, transb, m, n, A, lda, B, ldb, C, ldc)
end function qmckl_distance_sq
#+end_src
#+CALL: generate_f_interface(table=qmckl_distance_sq_args,rettyp=get_value("FRetType"),fname=get_value("Name"))
#+RESULTS:
#+begin_src f90 :tangle (eval fh) :comments org :exports none
interface
integer(c_int32_t) function qmckl_distance_sq &
(context, transa, transb, m, n, A, lda, B, ldb, C, ldc) &
bind(C)
use, intrinsic :: iso_c_binding
implicit none
integer (c_int64_t) , intent(in) , value :: context
character (c_char) , intent(in) , value :: transa, transb
integer (c_int64_t) , intent(in) , value :: m, n
integer (c_int64_t) , intent(in) , value :: lda
integer (c_int64_t) , intent(in) , value :: ldb
integer (c_int64_t) , intent(in) , value :: ldc
real (c_double) , intent(in) :: A(lda,3)
real (c_double) , intent(in) :: B(ldb,3)
real (c_double) , intent(out) :: C(ldc,n)
integer (c_int64_t) , intent(in) :: context
character , intent(in) :: transa
character , intent(in) :: transb
integer (c_int64_t) , intent(in) :: m
integer (c_int64_t) , intent(in) :: n
real (c_double ) , intent(in) :: A(lda,3)
integer (c_int64_t) , intent(in) :: lda
real (c_double ) , intent(in) :: B(ldb,3)
integer (c_int64_t) , intent(in) :: ldb
real (c_double ) , intent(out) :: C(ldc,n)
integer (c_int64_t) , intent(in) :: ldc
end function qmckl_distance_sq
end interface
#+end_src
end interface
#+end_src
*** Test :noexport:
#+begin_src f90 :tangle (eval f_test)
@ -293,7 +336,7 @@ integer(c_int32_t) function test_qmckl_distance_sq(context) bind(C)
if ( dabs(1.d0 - C(i,j)/x) > 1.d-14 ) return
end do
end do
test_qmckl_distance_sq = &
qmckl_distance_sq(context, 'n', 'T', m, n, A, LDA, B, LDB, C, LDC)
@ -341,9 +384,9 @@ integer(c_int32_t) function test_qmckl_distance_sq(context) bind(C)
if ( dabs(1.d0 - C(i,j)/x) > 1.d-14 ) return
end do
end do
test_qmckl_distance_sq = 0
deallocate(A,B,C)
end function test_qmckl_distance_sq
#+end_src

View File

@ -49,11 +49,14 @@ void* qmckl_malloc(qmckl_context context,
const size_t size);
#+end_src
In this implementation, we use ~calloc~ because it initializes the
memory block to zero, so structs will have ~NULL~-initialized pointers.
# Source
#+begin_src c :tangle (eval c)
void* qmckl_malloc(qmckl_context context, const size_t size) {
void * pointer = malloc( (size_t) size );
void * pointer = calloc(size, (size_t) 1);
if (qmckl_context_check(context) != QMCKL_NULL_CONTEXT) {
qmckl_exit_code rc;

View File

@ -159,13 +159,13 @@ endif
.PHONY: clean
.SECONDARY: # Needed to keep the produced C and Fortran files
libqmckl.so: Makefile.generated
libqmckl.so: ../include/qmckl.h
$(MAKE) -f Makefile.generated
../include/qmckl.h: libqmckl.so
../include/qmckl.h: Makefile.generated
../tools/build_qmckl_h.sh
test: Makefile.generated ../include/qmckl.h
test: libqmckl.so
$(MAKE) -f Makefile.generated test
doc: $(ORG_SOURCE_FILES)

230
tools/lib.org Normal file
View File

@ -0,0 +1,230 @@
# -*- mode: org -*-
* Function to get the value of a property.
#+NAME: get_value
#+begin_src elisp :var key="Type"
(setq x (org-property-values key))
(pop x)
#+end_src
#+RESULTS: get_value
* Table of function arguments
#+NAME: test
| qmckl_context | context | in | Global state |
| char | transa | in | Array ~A~ is ~'N'~: Normal, ~'T'~: Transposed |
| char | transb | in | Array ~B~ is ~'N'~: Normal, ~'T'~: Transposed |
| int64_t | m | in | Number of points in the first set |
| int64_t | n | in | Number of points in the second set |
| double | A[3][lda] | in | Array containing the $m \times 3$ matrix $A$ |
| int64_t | lda | in | Leading dimension of array ~A~ |
| double | B[3][ldb] | in | Array containing the $n \times 3$ matrix $B$ |
| int64_t | ldb | in | Leading dimension of array ~B~ |
| double | C[n][ldc] | out | Array containing the $m \times n$ matrix $C$ |
| int64_t | ldc | in | Leading dimension of array ~C~ |
** Fortran-C type conversions
#+NAME:f_of_c
#+BEGIN_SRC python :var table=test :var rettyp="integer" :var fname=[] :results value :noweb yes :wrap "src f90 :tangle (eval f) :comments org :exports none"
f_of_c_d = { '' : ''
, 'qmckl_context' : 'integer (c_int64_t)'
, 'int32_t' : 'integer (c_int32_t)'
, 'int64_t' : 'integer (c_int64_t)'
, 'float' : 'real (c_float )'
, 'double' : 'real (c_double )'
, 'char' : 'character'
}
#+END_SRC
#+NAME:c_of_f
#+BEGIN_SRC python :var table=test :var rettyp="integer" :var fname=[] :results value :noweb yes :wrap "src f90 :tangle (eval f) :comments org :exports none"
ctypeid_d = { '' : ''
, 'integer' : 'integer(c_int32_t)'
, 'integer*8' : 'integer(c_int64_t)'
, 'real' : 'real(c_float)'
, 'real*8' : 'real(c_double)'
, 'character' : 'character(c_char)'
}
#+END_SRC
** Parse the table
#+NAME: parse_table
#+BEGIN_SRC python :results none :noweb yes :exports none
def parse_table(table):
result = []
for line in table:
d = { "c_type" : line[0],
"inout" : line[2].lower(),
"name" : line[1],
"comment" : line[3] }
# Handle inout
if d["inout"] in ["input", "in"]:
d["inout"] == "in"
elif d["inout"] in ["output", "out"]:
d["inout"] == "out"
elif d["inout"] in ["input/output", "inout"]:
d["inout"] == "inout"
# Find dimensions
dims = d["name"].split('[')
d["rank"] = len(dims) - 1
if d["rank"] == 0:
d["dims"] = []
else:
d["name"] = d["name"].split('[')[0].strip()
d["dims"] = [ x.replace(']','').strip() for x in dims[1:] ]
result.append(d)
return result
#+END_SRC
** Generates a C header
#+NAME: generate_c_header
#+BEGIN_SRC python :var table=[] :var rettyp=[] :var fname=[] :results drawer :noweb yes :wrap "src c :tangle (eval h) :comments org"
<<parse_table>>
results = []
for d in parse_table(table):
name = d["name"]
c_type = d["c_type"]
# Add star for arrays
if d["rank"] > 0:
c_type += "*"
# Only inputs are const
if d["inout"] == "in":
const = "const"
else:
const = " "
results += [ f" {const} {c_type} {name}" ]
results=',\n'.join(results)
template = f"""{rettyp} {fname} (
{results} ); """
return template
#+END_SRC
#+RESULTS: generate_c_header
#+begin_src c :tangle (eval h) :comments org
[] [] (
);
#+end_src
** Generates a C interface to the Fortran function
#+NAME: generate_c_interface
#+BEGIN_SRC python :var table=[] :var rettyp="integer" :var fname=[] :results value :noweb yes :wrap "src f90 :tangle (eval f) :comments org :exports none"
<<c_of_f>>
<<f_of_c>>
<<parse_table>>
d = parse_table(table)
args = ", ".join([ x["name"] for x in d ])
rettyp_c = ctypeid_d[rettyp.lower()]
results = [ f"{rettyp_c} function {fname} &"
, f" ({args}) &"
, " bind(C) result(info)"
, ""
, " use, intrinsic :: iso_c_binding"
, " implicit none"
, ""
]
for d in parse_table(table):
f_type = f_of_c_d[d["c_type"]]
inout = "intent("+d["inout"]+")"
name = d["name"]
# Input scalars are passed by value
if d["rank"] == 0 and inout == "in":
value = ", value"
else:
value = " "
# Append dimensions to the name
if d["rank"] == 0:
dims = ""
else:
d["dims"].reverse()
dims = "(" + ",".join(d["dims"]) + ")"
results += [ f" {f_type:20}, {inout:12}{value} :: {name}{dims}" ]
results += [ ""
, f" {rettyp_c}, external :: {fname}_f"
, f" info = {fname}_f &"
, f" ({args})"
, ""
, f"end function {fname}"
]
results='\n'.join(results)
return results
#+END_SRC
** Generates a Fortran interface to the C function
#+NAME: generate_f_interface
#+BEGIN_SRC python :var table=test :var rettyp="integer" :var fname=[] :results value :noweb yes :wrap "src f90 :tangle (eval fh) :comments org :exports none"
<<c_of_f>>
<<f_of_c>>
<<parse_table>>
d = parse_table(table)
args = ", ".join([ x["name"] for x in d ])
rettyp_c = ctypeid_d[rettyp.lower()]
results = [ f"interface"
, f" {rettyp_c} function {fname} &"
, f" ({args}) &"
, " bind(C)"
, " use, intrinsic :: iso_c_binding"
, " implicit none"
, ""
]
for d in parse_table(table):
f_type = f_of_c_d[d["c_type"]]
inout = "intent("+d["inout"]+")"
name = d["name"]
# Input scalars are passed by value
if d["rank"] == 0 and inout == "in":
value = ", value"
else:
value = " "
# Append dimensions to the name
if d["rank"] == 0:
dims = ""
else:
d["dims"].reverse()
dims = "(" + ",".join(d["dims"]) + ")"
results += [ f" {f_type:20}, {inout:12}{value} :: {name}{dims}" ]
results += [ ""
, f" end function {fname}"
, f"end interface"
]
results='\n'.join(results)
return results
#+END_SRC
#+RESULTS: generate_c_interface
#+begin_src f90 :tangle (eval f) :comments org :exports none
#+end_src