diff --git a/src/Makefile b/src/Makefile
index d5b363d..d98f8e6 100644
--- a/src/Makefile
+++ b/src/Makefile
@@ -2,7 +2,7 @@ CC=gcc
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
+FFLAGS=-fPIC -g -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
@@ -25,7 +25,7 @@ doc:$(ORG_SOURCE_FILES)
./create_doc.sh $(ORG_SOURCE_FILES)
clean:
- rm -f qmckl.h test_qmckl_* qmckl_*.f90 qmckl_*.c qmckl_*.o qmckl_*.h Makefile.generated libqmckl.so
+ rm -f qmckl.h test_qmckl_* test_qmckl.c qmckl_*.f90 qmckl_*.c qmckl_*.o qmckl_*.h Makefile.generated libqmckl.so *.html
Makefile.generated: $(ORG_SOURCE_FILES) Makefile create_makefile.sh
./create_makefile.sh $(ORG_SOURCE_FILES)
diff --git a/src/README.org b/src/README.org
index 5b2ade4..698ea34 100644
--- a/src/README.org
+++ b/src/README.org
@@ -103,6 +103,8 @@ rm ${nb}.md
If the name of the org-mode file is =xxx.org=, the name of the
produced C files should be =xxx.c= and =xxx.h= and the name of the
produced Fortran files should be =xxx.f90=
+
+ Arrays are in uppercase and scalars are in lowercase.
** Application programming interface
@@ -111,17 +113,17 @@ rm ${nb}.md
that the library will be easily usable in any language.
This implies that only the following data types are allowed in the API:
- - 32-bit and 64-bit floats and arrays
- - 32-bit and 64-bit integers and arrays
+ - 32-bit and 64-bit floats and arrays (=real= and =double=)
+ - 32-bit and 64-bit integers and arrays (=int32_t= and =int64_t=)
- Pointers should be represented as 64-bit integers (even on
32-bit architectures)
- ASCII strings are represented as a pointers to a character arrays
and terminated by a zero character (C convention).
+ # TODO : Link to repositories for bindings
To facilitate the use in other languages than C, we provide some
bindings in other languages in other repositories.
- # TODO : Link to repositories for bindings
** Global state
@@ -178,7 +180,7 @@ rm ${nb}.md
=context= variable.
* Algorithms
-
+
Reducing the scaling of an algorithm usually implies also reducing
its arithmetic complexity (number of flops per byte). Therefore,
for small sizes \(\mathcal{O}(N^3)\) and \(\mathcal{O}(N^2)\) algorithms
@@ -186,7 +188,6 @@ 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=)
@@ -200,6 +201,7 @@ rm ${nb}.md
- [[./qmckl_memory.org][Memory management]]
- [[./qmckl_context.org][Context]]
- [[./qmckl_distance.org][Distance]]
+ - [[./qmckl_ao.org][Atomic orbitals]]
* Acknowledgments
diff --git a/src/qmckl.org b/src/qmckl.org
index 2d7e4f7..c48e771 100644
--- a/src/qmckl.org
+++ b/src/qmckl.org
@@ -58,6 +58,7 @@ typedef int64_t qmckl_context ;
#include "qmckl_context.h"
#include "qmckl_distance.h"
+#include "qmckl_ao.h"
#+END_SRC
* End of header
diff --git a/src/qmckl_ao.org b/src/qmckl_ao.org
new file mode 100644
index 0000000..ba01f8d
--- /dev/null
+++ b/src/qmckl_ao.org
@@ -0,0 +1,400 @@
+# -*- mode: org -*-
+# vim: syntax=c
+#+TITLE: Atomic Orbitals
+
+#+HTML_HEAD:
+#+HTML_HEAD:
+#+HTML_HEAD:
+#+HTML_HEAD:
+#+HTML_HEAD:
+#+HTML_HEAD:
+
+This files contains all the routines for the computation of the
+values, gradients and Laplacian of the atomic basis functions.
+
+3 files are produced:
+- a header file : =qmckl_ao.h=
+- a source file : =qmckl_ao.f90=
+- a test file : =test_qmckl_ao.c=
+
+*** Header
+ #+BEGIN_SRC C :comments link :tangle qmckl_ao.h
+#ifndef QMCKL_AO_H
+#define QMCKL_AO_H
+#include "qmckl_context.h"
+#include "qmckl_distance.h"
+ #+END_SRC
+
+*** Source
+ #+BEGIN_SRC f90 :comments link :tangle qmckl_ao.f90
+
+ #+END_SRC
+
+*** Test
+ #+BEGIN_SRC C :comments link :tangle test_qmckl_ao.c
+#include
+#include "qmckl.h"
+#include "munit.h"
+MunitResult test_qmckl_ao() {
+ qmckl_context context;
+ context = qmckl_context_create();
+ #+END_SRC
+
+
+* Polynomials
+
+ \[ P_l(\mathbf{r},\mathbf{R}_i) = (x-X_i)^a (y-Y_i)^b (z-Z_i)^c \]
+
+** =qmckl_ao_powers=
+
+ Computes all the powers of the =n= input data up to the given
+ maximum value given in input for each of the $n$ points:
+
+ \[ P_{ij} = X_j^i \]
+
+*** Arguments
+
+ | =context= | input | Global state |
+ | =n= | input | Number of values |
+ | =X(n)= | input | Array containing the input values |
+ | =LMAX(n)= | input | Array containing the maximum power for each value |
+ | =P(LDP,n)= | output | Array containing all the powers of $X$ |
+ | =LDP= | input | Leading dimension of array =P= |
+
+*** Requirements
+
+ - =context= is not 0
+ - =n= > 0
+ - =X= is allocated with at least $n \times 8$ bytes
+ - =LMAX= is allocated with at least $n \times 4$ bytes
+ - =P= is allocated with at least $n \times \max_i \text{LMAX}_i \times 8$ bytes
+ - =LDP= >= $\max_i$ =LMAX[i]=
+
+*** Header
+ #+BEGIN_SRC C :comments link :tangle qmckl_ao.h
+qmckl_exit_code qmckl_ao_powers(qmckl_context context,
+ int64_t n,
+ double *X, int32_t *LMAX,
+ double *P, int64_t LDP);
+ #+END_SRC
+
+*** Source
+ #+BEGIN_SRC f90 :comments link :tangle qmckl_ao.f90
+integer function qmckl_ao_powers_f(context, n, X, LMAX, P, ldp) result(info)
+ implicit none
+ integer*8 , intent(in) :: context
+ integer*8 , intent(in) :: n
+ real*8 , intent(in) :: X(n)
+ integer , intent(in) :: LMAX(n)
+ real*8 , intent(out) :: P(ldp,n)
+ integer*8 , intent(in) :: ldp
+
+ integer*8 :: i,j
+
+ info = 0
+
+ if (context == 0_8) then
+ info = -1
+ return
+ endif
+
+ if (LDP < MAXVAL(LMAX)) then
+ info = -2
+ return
+ endif
+
+ do j=1,n
+ P(1,j) = X(j)
+ do i=2,LMAX(j)
+ P(i,j) = P(i-1,j) * X(j)
+ end do
+ end do
+
+end function qmckl_ao_powers_f
+
+
+integer(c_int32_t) function qmckl_ao_powers(context, n, X, LMAX, P, ldp) &
+ 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) :: X(n)
+ integer (c_int32_t) , intent(in) :: LMAX(n)
+ real (c_double) , intent(out) :: P(ldp,n)
+ integer (c_int64_t) , intent(in) , value :: ldp
+
+ integer, external :: qmckl_ao_powers_f
+ info = qmckl_ao_powers_f(context, n, X, LMAX, P, ldp)
+end function qmckl_ao_powers
+ #+END_SRC
+
+*** Test
+ #+BEGIN_SRC C :comments link :tangle test_qmckl_ao.c
+{
+ int64_t n, LDP ;
+ int32_t *LMAX ;
+ double *X, *P ;
+ int i, j;
+
+ n = 100;
+ LDP = 10;
+
+ X = (double*) qmckl_malloc (context, n*sizeof(double));
+ LMAX = (int32_t*) qmckl_malloc (context, n*sizeof(int32_t));
+ P = (double*) qmckl_malloc (context, LDP*n*sizeof(double));
+
+ for (j=0 ; j 0
+ - =X= is allocated with at least $3 \times 8$ bytes
+ - =R= is allocated with at least $3 \times 8$ bytes
+ - =lmax= >= 0
+ - On output, =n= should be equal to (=lmax=+1)(=lmax=+2)(=lmax=+3)/6
+ - =L= is allocated with at least $3 \times n \times 4$ bytes
+ - =ldl= >= 3
+ - =VGL= is allocated with at least $5 \times n \times 8$ bytes
+ - =ldv= >= 5
+
+*** Header
+ #+BEGIN_SRC C :comments link :tangle qmckl_ao.h
+qmckl_exit_code qmckl_ao_polynomial_vgl(qmckl_context context,
+ double *X, double *R,
+ int32_t lmax, int64_t *n,
+ int32_t *L, int64_t ldl,
+ double *VGL, int64_t ldv);
+ #+END_SRC
+
+*** Source
+ #+BEGIN_SRC f90 :comments link :tangle qmckl_ao.f90
+integer function qmckl_ao_polynomial_vgl_f(context, X, R, lmax, n, L, ldl, VGL, ldv) result(info)
+ implicit none
+ integer*8 , intent(in) :: context
+ real*8 , intent(in) :: X(3), R(3)
+ integer , intent(in) :: lmax
+ integer*8 , intent(out) :: n
+ integer , intent(out) :: L(ldl,(lmax+1)*(lmax+2)*(lmax+3)/6)
+ integer*8 , intent(in) :: ldl
+ real*8 , intent(out) :: VGL(ldv,(lmax+1)*(lmax+2)*(lmax+3)/6)
+ integer*8 , intent(in) :: ldv
+
+ integer*8 :: i,j
+ integer :: a,b,c,d
+ real*8 :: Y(3)
+ integer :: lmax_array(3)
+ real*8 :: pows(-2:lmax,3)
+ integer, external :: qmckl_ao_powers_f
+
+ info = 0
+
+ if (context == 0_8) then
+ info = -1
+ return
+ endif
+
+ n = (lmax+1)*(lmax+2)*(lmax+3)/6
+
+ if (ldl < 3) then
+ info = -2
+ return
+ endif
+
+ if (ldv < 5) then
+ info = -3
+ return
+ endif
+
+
+ do i=1,3
+ Y(i) = X(i) - R(i)
+ end do
+ pows(-2:-1,1:3) = 0.d0
+ pows(0,1:3) = 1.d0
+ lmax_array(1:3) = lmax
+ info = qmckl_ao_powers_f(context, 1_8, Y(1), (/lmax/), pows(1,1), size(pows,1,kind=8))
+ if (info /= 0) return
+ info = qmckl_ao_powers_f(context, 1_8, Y(2), (/lmax/), pows(1,2), size(pows,1,kind=8))
+ if (info /= 0) return
+ info = qmckl_ao_powers_f(context, 1_8, Y(3), (/lmax/), pows(1,3), size(pows,1,kind=8))
+ if (info /= 0) return
+
+
+ n=1
+ vgl(1:5,1:n) = 0.d0
+ l(1:3,n) = 0
+ vgl(1,n) = 1.d0
+ do d=1,lmax
+ do a=0,d
+ do b=0,d
+ do c=0,d
+ if (a+b+c == d) then
+ n = n+1
+ l(1,n) = a
+ l(2,n) = b
+ l(3,n) = c
+
+ vgl(1,n) = pows(a,1) * pows(b,2) * pows(c,3)
+
+ vgl(2,n) = dble(a) * pows(a-1,1) * pows(b ,2) * pows(c ,3)
+ vgl(3,n) = dble(b) * pows(a ,1) * pows(b-1,2) * pows(c ,3)
+ vgl(4,n) = dble(c) * pows(a ,1) * pows(b ,2) * pows(c-1,3)
+
+ vgl(5,n) = dble(a) * dble(a-1) * pows(a-2,1) * pows(b ,2) * pows(c ,3) &
+ + dble(b) * dble(b-1) * pows(a ,1) * pows(b-2,2) * pows(c ,3) &
+ + dble(c) * dble(c-1) * pows(a ,1) * pows(b ,2) * pows(c-2,3)
+ exit
+ end if
+ end do
+ end do
+ end do
+ end do
+
+end function qmckl_ao_polynomial_vgl_f
+
+integer(c_int32_t) function qmckl_ao_polynomial_vgl(context, X, R, lmax, n, L, ldl, VGL, ldv) &
+ bind(C) result(info)
+ use, intrinsic :: iso_c_binding
+ implicit none
+ integer (c_int64_t) , intent(in) , value :: context
+ real (c_double) , intent(in) :: X(3), R(3)
+ integer (c_int32_t) , intent(in) , value :: lmax
+ integer (c_int64_t) , intent(out) :: n
+ integer (c_int32_t) , intent(out) :: L(ldl,(lmax+1)*(lmax+2)*(lmax+3)/6)
+ integer (c_int64_t) , intent(in) , value :: ldl
+ real (c_double) , intent(out) :: VGL(ldv,(lmax+1)*(lmax+2)*(lmax+3)/6)
+ integer (c_int64_t) , intent(in) , value :: ldv
+
+ integer, external :: qmckl_ao_polynomial_vgl_f
+ info = qmckl_ao_polynomial_vgl_f(context, X, R, lmax, n, L, ldl, VGL, ldv)
+end function qmckl_ao_polynomial_vgl
+ #+END_SRC
+
+*** Test
+ #+BEGIN_SRC C :comments link :tangle test_qmckl_ao.c
+{
+#include
+ double X[3] = { 1.1 , 2.2 , 3.3 };
+ double R[3] = { 0.1 , 1.2 , -2.3 };
+ double Y[3];
+ int32_t lmax = 4;
+ int64_t n = 0;
+ int64_t ldl = 3;
+ int64_t ldv = 100;
+ int32_t* L_mem;
+ int32_t* L[100];
+ double* VGL_mem;
+ double* VGL[100];
+ int j;
+
+ int d = (lmax+1)*(lmax+2)*(lmax+3)/6;
+
+ L_mem = (int32_t*) malloc(ldl*100*sizeof(int32_t));
+ VGL_mem = (double*) malloc(ldv*100*sizeof(double));
+
+ munit_assert_int64(QMCKL_SUCCESS, ==,
+ qmckl_ao_polynomial_vgl(context, X, R, lmax, &n, L_mem, ldl, VGL_mem, ldv) );
+
+ munit_assert_int64( n, ==, d );
+ for (j=0 ; j=, 0 );
+ munit_assert_int64( L[j][1], >=, 0 );
+ munit_assert_int64( L[j][2], >=, 0 );
+ munit_assert_double_equal( VGL[j][0],
+ pow(Y[0],L[j][0]) * pow(Y[1],L[j][1]) * pow(Y[2],L[j][2]), 10 );
+ if (L[j][0] < 1) {
+ munit_assert_double_equal( VGL[j][1], 0., 10);
+ } else {
+ munit_assert_double_equal( VGL[j][1],
+ L[j][0] * pow(Y[0],L[j][0]-1) * pow(Y[1],L[j][1]) * pow(Y[2],L[j][2]), 10 );
+ }
+ if (L[j][1] < 1) {
+ munit_assert_double_equal( VGL[j][2], 0., 10);
+ } else {
+ munit_assert_double_equal( VGL[j][2],
+ L[j][1] * pow(Y[0],L[j][0]) * pow(Y[1],L[j][1]-1) * pow(Y[2],L[j][2]), 10 );
+ }
+ if (L[j][2] < 1) {
+ munit_assert_double_equal( VGL[j][3], 0., 10);
+ } else {
+ munit_assert_double_equal( VGL[j][3],
+ L[j][2] * pow(Y[0],L[j][0]) * pow(Y[1],L[j][1]) * pow(Y[2],L[j][2]-1), 10 );
+ }
+
+ double w = 0.;
+ if (L[j][0] > 1) w += L[j][0] * (L[j][0]-1) * pow(Y[0],L[j][0]-2) * pow(Y[1],L[j][1]) * pow(Y[2],L[j][2]);
+ if (L[j][1] > 1) w += L[j][1] * (L[j][1]-1) * pow(Y[0],L[j][0]) * pow(Y[1],L[j][1]-2) * pow(Y[2],L[j][2]);
+ if (L[j][2] > 1) w += L[j][2] * (L[j][2]-1) * pow(Y[0],L[j][0]) * pow(Y[1],L[j][1]) * pow(Y[2],L[j][2]-2);
+ munit_assert_double_equal( VGL[j][4], w, 10 );
+ }
+ free(L_mem);
+ free(VGL_mem);
+}
+ #+END_SRC
+
+
+
+* TODO Gaussian basis functions
+
+* TODO Slater basis functions
+
+* End of files
+
+*** Header
+ #+BEGIN_SRC C :comments link :tangle qmckl_ao.h
+#endif
+ #+END_SRC
+
+*** Test
+ #+BEGIN_SRC C :comments link :tangle test_qmckl_ao.c
+ if (qmckl_context_destroy(context) != QMCKL_SUCCESS)
+ return QMCKL_FAILURE;
+ return MUNIT_OK;
+}
+
+ #+END_SRC
diff --git a/src/qmckl_distance.org b/src/qmckl_distance.org
index b5deee5..cba53ca 100644
--- a/src/qmckl_distance.org
+++ b/src/qmckl_distance.org
@@ -77,25 +77,24 @@ MunitResult test_qmckl_distance() {
*** 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 |
+ | =context= | input | Global state |
+ | =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= |
*** Requirements
- =context= is not 0
- =m= > 0
- =n= > 0
- - =LDA= >= m
- - =LDB= >= n
- - =LDC= >= m
+ - =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
@@ -104,28 +103,26 @@ MunitResult test_qmckl_distance() {
#+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);
+ 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
+integer function qmckl_distance_sq_f(context, m, n, A, LDA, B, LDB, C, LDC) result(info)
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*8 , intent(in) :: context
+ integer*8 , intent(in) :: m, n
+ integer*8 , intent(in) :: lda
+ real*8 , intent(in) :: A(lda,3)
+ integer*8 , intent(in) :: ldb
+ real*8 , intent(in) :: B(ldb,3)
+ integer*8 , intent(in) :: ldc
+ real*8 , intent(out) :: C(ldc,n)
- integer (c_int64_t) :: i,j
- real (c_double) :: x, y, z
+ integer*8 :: i,j
+ real*8 :: x, y, z
info = 0
@@ -168,6 +165,24 @@ integer(c_int32_t) function qmckl_distance_sq(context, m, n, A, LDA, B, LDB, C,
end do
end do
+end function qmckl_distance_sq_f
+
+! C interface
+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, external :: qmckl_distance_sq_f
+ info = qmckl_distance_sq_f(context, m, n, A, LDA, B, LDB, C, LDC)
end function qmckl_distance_sq
#+END_SRC
diff --git a/src/test_qmckl.org b/src/test_qmckl.org
index f2ca850..bcd0fce 100644
--- a/src/test_qmckl.org
+++ b/src/test_qmckl.org
@@ -23,6 +23,7 @@ grep BEGIN_SRC *.org | \
#+END_SRC
#+RESULTS: test-files
+| test_qmckl_ao.c |
| test_qmckl_context.c |
| test_qmckl_distance.c |
| test_qmckl_memory.c |
@@ -42,6 +43,7 @@ echo "#+END_SRC"
#+RESULTS:
#+NAME: headers
#+BEGIN_SRC C :tangle no
+MunitResult test_qmckl_ao();
MunitResult test_qmckl_context();
MunitResult test_qmckl_distance();
MunitResult test_qmckl_memory();
@@ -62,6 +64,7 @@ echo "#+END_SRC"
#+RESULTS:
#+NAME: calls
#+BEGIN_SRC C :tangle no
+ { (char*) "test_qmckl_ao", test_qmckl_ao, NULL,NULL,MUNIT_TEST_OPTION_NONE,NULL},
{ (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},