1
0
mirror of https://github.com/TREX-CoE/qmckl.git synced 2025-01-03 10:06:09 +01:00

Added f90 example file

This commit is contained in:
Anthony Scemama 2020-10-22 00:50:07 +02:00
parent 5f5465eaf9
commit 150518aef0
9 changed files with 291 additions and 69 deletions

1
src/.gitignore vendored
View File

@ -1,5 +1,6 @@
*.o
*.c
*.f90
*.h
*~
*.so

View File

@ -4,8 +4,10 @@ CFLAGS=-fPIC -fexceptions -Wall -Werror -Wpedantic -Wextra -g
FC=gfortran
FFLAGS=-fcheck=all -Waliasing -Wampersand -Wconversion -Wsurprising -Wintrinsics-std -Wno-tabs -Wintrinsic-shadow -Wline-truncation -Wreal-q-constant -Wuninitialized -fbacktrace -ffpe-trap=zero,overflow,underflow -finit-real=nan
LIBS=-lgfortran -lm
export CC CFLAGS FC FFLAGS
export CC CFLAGS FC FFLAGS LIBS
ORG_SOURCE_FILES=$(wildcard qmckl*.org) test_qmckl.org
OBJECT_FILES=$(filter-out $(EXCLUDED_OBJECTS), $(patsubst %.org,%.o,$(ORG_SOURCE_FILES)))

View File

@ -179,11 +179,20 @@ rm ${nb}.md
As QMCkl is a general purpose library, multiple algorithms should
be implemented adapted to different problem sizes.
** Rules for the API
- =stdint= should be used for integers (=int32_t=, =int64_t=)
- integers used for counting should always be =int64_t=
- floats should be by default =double=, unless explicitly mentioned
- pointers are converted to =int64_t= to increase portability
** Documentation
- [[qmckl.org][Main QMCkl header file]]
- [[qmckl_memory.org][Memory management]]
- [[qmckl_context.org][Context]]
- [[qmckldistance.org][Distance]]
** Acknowledgments

View File

@ -53,6 +53,8 @@ FFLAGS=$FFLAGS
OBJECT_FILES=$OBJECTS
TESTS=$TESTS
LIBS=$LIBS
libqmckl.so: \$(OBJECT_FILES)
\$(CC) -shared \$(OBJECT_FILES) -o libqmckl.so
@ -65,7 +67,7 @@ libqmckl.so: \$(OBJECT_FILES)
test_qmckl: test_qmckl.c libqmckl.so \$(TESTS)
echo \$(TESTS)
\$(CC) \$(CFLAGS) -Wl,-rpath,$PWD -L. \
../munit/munit.c \$(TESTS) -lqmckl test_qmckl.c -o test_qmckl
../munit/munit.c \$(TESTS) -lqmckl \$(LIBS) test_qmckl.c -o test_qmckl
test: test_qmckl
./test_qmckl

View File

@ -8,6 +8,8 @@ other C header files. It is the main entry point to the library.
#+BEGIN_SRC C :tangle qmckl.h
#ifndef QMCKL_H
#define QMCKL_H
#include <stdlib.h>
#include <stdint.h>
#+END_SRC
* Constants
@ -21,7 +23,9 @@ other C header files. It is the main entry point to the library.
#define QMCKL_SUCCESS 0
#define QMCKL_FAILURE 1
typedef int qmckl_exit_code;
typedef int32_t qmckl_exit_code;
typedef int64_t qmckl_context ;
#+END_SRC
@ -42,9 +46,11 @@ typedef int qmckl_exit_code;
header files.
#+BEGIN_SRC C :tangle qmckl.h
#include <stdlib.h>
#include "qmckl_memory.h"
#include "qmckl_context.h"
#include "qmckl_distance.h"
#+END_SRC
* End of header

View File

@ -38,19 +38,13 @@ MunitResult test_qmckl_context() {
into a 64-bit signed integer, defined in the =qmckl_context= type.
A value of 0 for the context is equivalent to a NULL pointer.
*** Header
#+BEGIN_SRC C :comments link :tangle qmckl_context.h
/* 64-bit integer */
typedef long long int qmckl_context ;
#+END_SRC
*** Source
#+BEGIN_SRC C :comments link :tangle qmckl_context.c
typedef struct qmckl_context_struct {
struct qmckl_context_struct * prev;
unsigned int tag;
int precision;
int range;
uint32_t tag;
int32_t precision;
int32_t range;
} qmckl_context_struct;
#define VALID_TAG 0xBEEFFACE
@ -75,12 +69,12 @@ typedef struct qmckl_context_struct {
*** Header
#+BEGIN_SRC C :comments link :tangle qmckl_context.h
qmckl_context qmckl_context_check(qmckl_context context) ;
qmckl_context qmckl_context_check(const qmckl_context context) ;
#+END_SRC
*** Source
#+BEGIN_SRC C :comments link :tangle qmckl_context.c
qmckl_context qmckl_context_check(qmckl_context context) {
qmckl_context qmckl_context_check(const qmckl_context context) {
qmckl_context_struct * ctx;
if (context == (qmckl_context) 0) return (qmckl_context) 0;
@ -109,7 +103,7 @@ qmckl_context qmckl_context_create() {
qmckl_context_struct* context;
context = (qmckl_context_struct*) qmckl_malloc (sizeof(qmckl_context_struct));
context = (qmckl_context_struct*) qmckl_malloc ((qmckl_context) 0, sizeof(qmckl_context_struct));
if (context == NULL) {
return (qmckl_context) 0;
}
@ -126,8 +120,8 @@ qmckl_context qmckl_context_create() {
*** Test
#+BEGIN_SRC C :comments link :tangle test_qmckl_context.c
context = qmckl_context_create();
munit_assert_long( context, !=, (qmckl_context) 0);
munit_assert_long( qmckl_context_check(context), ==, context);
munit_assert_int64( context, !=, (qmckl_context) 0);
munit_assert_int64( qmckl_context_check(context), ==, context);
#+END_SRC
** =qmckl_context_copy=
@ -157,7 +151,7 @@ qmckl_context qmckl_context_copy(const qmckl_context context) {
return (qmckl_context) 0;
}
new_context = (qmckl_context_struct*) qmckl_malloc (sizeof(qmckl_context_struct));
new_context = (qmckl_context_struct*) qmckl_malloc (context, sizeof(qmckl_context_struct));
if (new_context == NULL) {
return (qmckl_context) 0;
}
@ -177,9 +171,9 @@ qmckl_context qmckl_context_copy(const qmckl_context context) {
*** Test
#+BEGIN_SRC C :comments link :tangle test_qmckl_context.c
new_context = qmckl_context_copy(context);
munit_assert_long(new_context, !=, (qmckl_context) 0);
munit_assert_long(new_context, !=, context);
munit_assert_long(qmckl_context_check(new_context), ==, new_context);
munit_assert_int64(new_context, !=, (qmckl_context) 0);
munit_assert_int64(new_context, !=, context);
munit_assert_int64(qmckl_context_check(new_context), ==, new_context);
#+END_SRC
** =qmckl_context_previous=
@ -213,10 +207,10 @@ qmckl_context qmckl_context_previous(const qmckl_context context) {
*** Test
#+BEGIN_SRC C :comments link :tangle test_qmckl_context.c
munit_assert_long(qmckl_context_previous(new_context), !=, (qmckl_context) 0);
munit_assert_long(qmckl_context_previous(new_context), ==, context);
munit_assert_long(qmckl_context_previous(context), ==, (qmckl_context) 0);
munit_assert_long(qmckl_context_previous((qmckl_context) 0), ==, (qmckl_context) 0);
munit_assert_int64(qmckl_context_previous(new_context), !=, (qmckl_context) 0);
munit_assert_int64(qmckl_context_previous(new_context), ==, context);
munit_assert_int64(qmckl_context_previous(context), ==, (qmckl_context) 0);
munit_assert_int64(qmckl_context_previous((qmckl_context) 0), ==, (qmckl_context) 0);
#+END_SRC
** =qmckl_context_destroy=
@ -253,12 +247,12 @@ qmckl_exit_code qmckl_context_destroy(qmckl_context context) {
*** Test
#+BEGIN_SRC C :comments link :tangle test_qmckl_context.c
munit_assert_long(qmckl_context_check(new_context), ==, new_context);
munit_assert_long(new_context, !=, (qmckl_context) 0);
munit_assert_int(qmckl_context_destroy(new_context), ==, QMCKL_SUCCESS);
munit_assert_long(qmckl_context_check(new_context), !=, new_context);
munit_assert_long(qmckl_context_check(new_context), ==, (qmckl_context) 0);
munit_assert_long(qmckl_context_destroy((qmckl_context) 0), ==, QMCKL_FAILURE);
munit_assert_int64(qmckl_context_check(new_context), ==, new_context);
munit_assert_int64(new_context, !=, (qmckl_context) 0);
munit_assert_int32(qmckl_context_destroy(new_context), ==, QMCKL_SUCCESS);
munit_assert_int64(qmckl_context_check(new_context), !=, new_context);
munit_assert_int64(qmckl_context_check(new_context), ==, (qmckl_context) 0);
munit_assert_int64(qmckl_context_destroy((qmckl_context) 0), ==, QMCKL_FAILURE);
#+END_SRC
@ -275,11 +269,11 @@ qmckl_exit_code qmckl_context_destroy(qmckl_context context) {
** =qmckl_context_update_precision=
#+BEGIN_SRC C :comments link :tangle qmckl_context.h
qmckl_exit_code qmckl_context_update_precision(const qmckl_context context, int precision);
qmckl_exit_code qmckl_context_update_precision(const qmckl_context context, const int precision);
#+END_SRC
#+BEGIN_SRC C :comments link :tangle qmckl_context.c
qmckl_exit_code qmckl_context_update_precision(const qmckl_context context, int precision) {
qmckl_exit_code qmckl_context_update_precision(const qmckl_context context, const int precision) {
qmckl_context_struct* ctx;
if (precision < 2) return QMCKL_FAILURE;
@ -295,11 +289,11 @@ qmckl_exit_code qmckl_context_update_precision(const qmckl_context context, int
** =qmckl_context_update_range=
#+BEGIN_SRC C :comments link :tangle qmckl_context.h
qmckl_exit_code qmckl_context_update_range(const qmckl_context context, int range);
qmckl_exit_code qmckl_context_update_range(const qmckl_context context, const int range);
#+END_SRC
#+BEGIN_SRC C :comments link :tangle qmckl_context.c
qmckl_exit_code qmckl_context_update_range(const qmckl_context context, int range) {
qmckl_exit_code qmckl_context_update_range(const qmckl_context context, const int range) {
qmckl_context_struct* ctx;
if (range < 2) return QMCKL_FAILURE;
@ -318,7 +312,7 @@ qmckl_exit_code qmckl_context_update_range(const qmckl_context context, int rang
** =qmckl_context_set_precision=
#+BEGIN_SRC C :comments link :tangle qmckl_context.h
qmckl_context qmckl_context_set_precision(const qmckl_context context, int precision);
qmckl_context qmckl_context_set_precision(const qmckl_context context, const int precision);
#+END_SRC
#+BEGIN_SRC C :comments link :tangle qmckl_context.c
@ -336,11 +330,11 @@ qmckl_context qmckl_context_set_precision(const qmckl_context context, const int
** =qmckl_context_set_range=
#+BEGIN_SRC C :comments link :tangle qmckl_context.h
qmckl_context qmckl_context_set_range(const qmckl_context context, int range);
qmckl_context qmckl_context_set_range(const qmckl_context context, const int range);
#+END_SRC
#+BEGIN_SRC C :comments link :tangle qmckl_context.c
qmckl_context qmckl_context_set_range(const qmckl_context context, int range) {
qmckl_context qmckl_context_set_range(const qmckl_context context, const int range) {
qmckl_context new_context;
new_context = qmckl_context_copy(context);

201
src/qmckl_distance.org Normal file
View File

@ -0,0 +1,201 @@
# -*- mode: org -*-
# vim: syntax=c
#+TITLE: Computation of distances
Function for the computation of distances between particles.
3 files are produced:
- a header file : =qmckl_distance.h=
- a source file : =qmckl_distance.f90=
- a test file : =test_qmckl_distance.c=
*** Header
#+BEGIN_SRC C :comments link :tangle qmckl_distance.h
#ifndef QMCKL_DISTANCE_H
#define QMCKL_DISTANCE_H
#include "qmckl_context.h"
#+END_SRC
*** Source
#+BEGIN_SRC f90 :comments link :tangle qmckl_distance.f90
#+END_SRC
*** Test
#+BEGIN_SRC C :comments link :tangle test_qmckl_distance.c
#include <math.h>
#include "qmckl.h"
#include "munit.h"
MunitResult test_qmckl_distance() {
qmckl_context context;
int64_t m, n, LDA, LDB, LDC;
double *A, *B, *C ;
int i, j;
context = qmckl_context_create();
m = 5;
n = 6;
LDA = 6;
LDB = 10;
LDC = 5;
A = (double*) qmckl_malloc (context, LDA*4*sizeof(double));
B = (double*) qmckl_malloc (context, LDB*3*sizeof(double));
C = (double*) qmckl_malloc (context, LDC*n*sizeof(double));
for (j=0 ; j<3 ; j++) {
for (i=0 ; i<m ; i++) {
A[i+j*LDA] = -10. + (double) (i+j);
}
}
for (j=0 ; j<3 ; j++) {
for (i=0 ; i<n ; i++) {
B[i+j*LDB] = -1. + (double) (i*j);
}
}
#+END_SRC
* Squared distance
** =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^2} = \sum_{k=1}^3 (A_{i,k}-B_{j,k})^2
\]
*** Arguments
| =context= | input | Global state |
| =m= | input | Number of points in the first set |
| =n= | input | Number of points in the second set |
| =LDA= | input | Leading dimension of array =A= |
| =A= | input | Array containing the $3 \times m$ matrix $A$ |
| =LDB= | input | Leading dimension of array =B= |
| =B= | input | Array containing the $3 \times n$ matrix $B$ |
| =LDC= | input | Leading dimension of array =C= |
| =C= | output | Array containing the $m \times n$ matrix $C$ |
| =info= | output | exit status is zero upon success |
*** Requirements
- =context= is not 0
- =m= > 0
- =n= > 0
- =LDA= >= m
- =LDB= >= n
- =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
*** Header
#+BEGIN_SRC C :comments link :tangle qmckl_distance.h
qmckl_exit_code qmckl_distance_sq(qmckl_context context,
int64_t m, int64_t n,
double *A, int64_t LDA,
double *B, int64_t LDB,
double *C, int64_t LDC);
#+END_SRC
*** Source
#+BEGIN_SRC f90 :comments link :tangle qmckl_distance.f90
integer(c_int32_t) function qmckl_distance_sq(context, 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
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 (c_int64_t) :: i,j
real (c_double) :: x, y, z
info = 0
if (context == 0_8) then
info = -1
return
endif
if (m <= 0_8) then
info = -2
return
endif
if (n <= 0_8) then
info = -3
return
endif
if (LDA < m) then
info = -4
return
endif
if (LDB < n) then
info = -5
return
endif
if (LDC < m) then
info = -6
return
endif
do j=1,n
do i=1,m
x = A(i,1) - B(j,1)
y = A(i,2) - B(j,2)
z = A(i,3) - B(j,3)
C(i,j) = x*x + y*y + z*z
end do
end do
end function qmckl_distance_sq
#+END_SRC
*** Test
#+BEGIN_SRC C :comments link :tangle test_qmckl_distance.c
munit_assert_int64(QMCKL_SUCCESS, ==,
qmckl_distance_sq(context, m, n, A, LDA, B, LDB, C, LDC) );
for (j=0 ; j<n ; j++) {
for (i=0 ; i<m ; i++) {
munit_assert_double_equal(C[i+j*LDC],
pow(A[i ]-B[j ],2) +
pow(A[i+ LDA]-B[j+ LDB],2) +
pow(A[i+2*LDA]-B[j+2*LDB],2) ,
14 );
}
}
#+END_SRC
* End of files
*** Header
#+BEGIN_SRC C :comments link :tangle qmckl_distance.h
#endif
#+END_SRC
*** Test
#+BEGIN_SRC C :comments link :tangle test_qmckl_distance.c
qmckl_free(A);
qmckl_free(B);
qmckl_free(C);
if (qmckl_context_destroy(context) != QMCKL_SUCCESS)
return QMCKL_FAILURE;
return MUNIT_OK;
}
#+END_SRC

View File

@ -31,15 +31,19 @@ MunitResult test_qmckl_memory() {
#+END_SRC
** =qmckl_malloc=
Analogous of =malloc, but passing signed 64-bit integers as argument.=
Analogous of =malloc, but passing a context and a signed 64-bit integers as argument.=
*** Header
#+BEGIN_SRC C :comments link :tangle qmckl_memory.h
void* qmckl_malloc(long long int size);
void* qmckl_malloc(const qmckl_context ctx, const size_t size);
#+END_SRC
*** Source
#+BEGIN_SRC C :comments link :tangle qmckl_memory.c
void* qmckl_malloc(long long int size) {
void* qmckl_malloc(const qmckl_context ctx, const size_t size) {
if (ctx == (qmckl_context) 0) {
/* Avoids unused parameter error */
return malloc( (size_t) size );
}
return malloc( (size_t) size );
}
@ -48,7 +52,7 @@ void* qmckl_malloc(long long int size) {
*** Test
#+BEGIN_SRC C :comments link :tangle test_qmckl_memory.c
int *a;
a = (int*) qmckl_malloc(3*sizeof(int));
a = (int*) qmckl_malloc( (qmckl_context) 1, 3*sizeof(int));
a[0] = 1;
a[1] = 2;
a[2] = 3;

View File

@ -17,6 +17,7 @@ grep BEGIN_SRC *.org | \
#+RESULTS: test-files
| test_qmckl_context.c |
| test_qmckl_distance.c |
| test_qmckl_memory.c |
We generate the function headers
@ -35,6 +36,7 @@ echo "#+END_SRC"
#+NAME: headers
#+BEGIN_SRC C :tangle no
MunitResult test_qmckl_context();
MunitResult test_qmckl_distance();
MunitResult test_qmckl_memory();
#+END_SRC
@ -54,6 +56,7 @@ echo "#+END_SRC"
#+NAME: calls
#+BEGIN_SRC C :tangle no
{ (char*) "test_qmckl_context", test_qmckl_context, NULL,NULL,MUNIT_TEST_OPTION_NONE,NULL},
{ (char*) "test_qmckl_distance", test_qmckl_distance, NULL,NULL,MUNIT_TEST_OPTION_NONE,NULL},
{ (char*) "test_qmckl_memory", test_qmckl_memory, NULL,NULL,MUNIT_TEST_OPTION_NONE,NULL},
#+END_SRC