TREX/qmckl
1
0
Fork 0
mirror of https://github.com/TREX-CoE/qmckl.git synced 2024-07-18 00:43:51 +02:00
Code Issues Projects Releases Wiki Activity

Tests in Fortran

Browse Source
This commit is contained in:
Anthony Scemama 2020-10-26 19:30:50 +01:00
parent 4c7b2213f4
commit d91cb765e3
4 changed files with 237 additions and 153 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

10
src/Makefile
View File

@ -6,6 +6,14 @@ FFLAGS=-fPIC -g -fcheck=all -Waliasing -Wampersand -Wconversion -Wsurprising -Wi
LIBS=-lgfortran -lm LIBS=-lgfortran -lm
#CC=icc
#CFLAGS=-fPIC -g
#
#FC=ifort
#FFLAGS=-fPIC -g
#
#LIBS=-lm -lifcore -lirc
export CC CFLAGS FC FFLAGS LIBS export CC CFLAGS FC FFLAGS LIBS
@ -25,7 +33,7 @@ doc:$(ORG_SOURCE_FILES)
./create_doc.sh $(ORG_SOURCE_FILES) ./create_doc.sh $(ORG_SOURCE_FILES)
clean: clean:
rm -f qmckl.h test_qmckl_* test_qmckl.c qmckl_*.f90 qmckl_*.c qmckl_*.o qmckl_*.h Makefile.generated libqmckl.so *.html rm -f qmckl.h test_qmckl_* test_qmckl.c qmckl_*.f90 qmckl_*.c qmckl_*.o qmckl_*.h Makefile.generated libqmckl.so *.html *.fh
Makefile.generated: $(ORG_SOURCE_FILES) Makefile create_makefile.sh Makefile.generated: $(ORG_SOURCE_FILES) Makefile create_makefile.sh
./create_makefile.sh $(ORG_SOURCE_FILES) ./create_makefile.sh $(ORG_SOURCE_FILES)

7
src/README.org
View File

@ -72,8 +72,12 @@ rm ${nb}.md
** Writing in Fortran ** Writing in Fortran
The Fortran source files should provide a C interface using The Fortran source files should provide a C interface using
iso-c-binding. The name of the Fortran source files should end =iso_c_binding=. The name of the Fortran source files should end
with =_f.f90= to be properly handled by the Makefile. with =_f.f90= to be properly handled by the Makefile.
The names of the functions defined in fortran should be the same as
those exposed in the API suffixed by =_f=.
Fortran interface files should also be written in a file with a
=.fh= extension.
** Coding style ** Coding style
# TODO: decide on a coding style # TODO: decide on a coding style
@ -124,7 +128,6 @@ rm ${nb}.md
To facilitate the use in other languages than C, we provide some To facilitate the use in other languages than C, we provide some
bindings in other languages in other repositories. bindings in other languages in other repositories.
** Global state ** Global state
Global variables should be avoided in the library, because it is Global variables should be avoided in the library, because it is

245
src/qmckl_ao.org
View File

@ -12,10 +12,11 @@
This files contains all the routines for the computation of the This files contains all the routines for the computation of the
values, gradients and Laplacian of the atomic basis functions. values, gradients and Laplacian of the atomic basis functions.
3 files are produced: 4 files are produced:
- a header file : =qmckl_ao.h= - a header file : =qmckl_ao.h=
- a source file : =qmckl_ao.f90= - a source file : =qmckl_ao.f90=
- a test file : =test_qmckl_ao.c= - a C test file : =test_qmckl_ao.c=
- a Fortran test file : =test_qmckl_ao_f.f90=
*** Header :noexport: *** Header :noexport:
#+BEGIN_SRC C :comments link :tangle qmckl_ao.h #+BEGIN_SRC C :comments link :tangle qmckl_ao.h
@ -25,11 +26,6 @@ values, gradients and Laplacian of the atomic basis functions.
#include "qmckl_distance.h" #include "qmckl_distance.h"
#+END_SRC #+END_SRC
*** Source :noexport:
#+BEGIN_SRC f90 :comments link :tangle qmckl_ao.f90
#+END_SRC
*** Test :noexport: *** Test :noexport:
#+BEGIN_SRC C :comments link :tangle test_qmckl_ao.c #+BEGIN_SRC C :comments link :tangle test_qmckl_ao.c
#include <math.h> #include <math.h>
@ -131,40 +127,64 @@ integer(c_int32_t) function qmckl_ao_powers(context, n, X, LMAX, P, ldp) &
end function qmckl_ao_powers end function qmckl_ao_powers
#+END_SRC #+END_SRC
#+BEGIN_SRC f90 :comments link :tangle qmckl_ao.fh
interface
integer(c_int32_t) function qmckl_ao_powers(context, n, X, LMAX, P, ldp) bind(C)
use, intrinsic :: iso_c_binding
integer (c_int64_t) , intent(in) , value :: context
integer (c_int64_t) , intent(in) , value :: n
integer (c_int64_t) , intent(in) , value :: ldp
real (c_double) , intent(in) :: X(n)
integer (c_int32_t) , intent(in) :: LMAX(n)
real (c_double) , intent(out) :: P(ldp,n)
end function qmckl_ao_powers
end interface
#+END_SRC
*** Test :noexport: *** Test :noexport:
#+BEGIN_SRC C :comments link :tangle test_qmckl_ao.c #+BEGIN_SRC f90 :comments link :tangle test_qmckl_ao_f.f90
{ integer(c_int32_t) function test_qmckl_ao_powers(context) bind(C)
int64_t n, LDP ; use, intrinsic :: iso_c_binding
int32_t *LMAX ; implicit none
double *X, *P ; include 'qmckl_ao.fh'
int i, j;
integer(c_int64_t), intent(in), value :: context
integer*8 :: n, LDP
integer, allocatable :: LMAX(:)
double precision, allocatable :: X(:), P(:,:)
integer*8 :: i,j
n = 100; n = 100;
LDP = 10; LDP = 10;
X = (double*) qmckl_malloc (context, n*sizeof(double)); allocate(X(n), P(LDP,n), LMAX(n))
LMAX = (int32_t*) qmckl_malloc (context, n*sizeof(int32_t));
P = (double*) qmckl_malloc (context, LDP*n*sizeof(double));
for (j=0 ; j<n ; j++) { do j=1,n
X[j] = -5. + 0.1 * (double) (j); X(j) = -5.d0 + 0.1d0 * dble(j)
LMAX[j] = 1 + (j % 9); LMAX(j) = 1 + int(mod(j, 9),4)
} end do
munit_assert_int64(QMCKL_SUCCESS, ==, test_qmckl_ao_powers = qmckl_ao_powers(context, n, X, LMAX, P, LDP)
qmckl_ao_powers(context, n, X, LMAX, P, LDP) ); if (test_qmckl_ao_powers /= 0) return
for (j=0 ; j<n ; j++) { test_qmckl_ao_powers = -1
for (i=0 ; i<LMAX[j] ; i++) {
munit_assert_double_equal( P[i+j*LDP], pow(X[j],i+1), 10 );
}
}
qmckl_free(X);
qmckl_free(P);
qmckl_free(LMAX);
}
#+END_SRC do j=1,n
do i=1,LMAX(j)
if ( dabs(1.d0 - P(i,j) / (X(j)**i)) > 1.d-14 ) return
end do
end do
test_qmckl_ao_powers = 0
deallocate(X,P,LMAX)
end function test_qmckl_ao_powers
#+END_SRC
#+BEGIN_SRC C :comments link :tangle test_qmckl_ao.c
int test_qmckl_ao_powers(qmckl_context context);
munit_assert_int(0, ==, test_qmckl_ao_powers(context));
#+END_SRC
** =qmckl_ao_polynomial_vgl= ** =qmckl_ao_polynomial_vgl=
@ -312,76 +332,115 @@ integer(c_int32_t) function qmckl_ao_polynomial_vgl(context, X, R, lmax, n, L, l
end function qmckl_ao_polynomial_vgl end function qmckl_ao_polynomial_vgl
#+END_SRC #+END_SRC
#+BEGIN_SRC f90 :comments link :tangle qmckl_ao.fh
interface
integer(c_int32_t) function qmckl_ao_polynomial_vgl(context, X, R, lmax, n, L, ldl, VGL, ldv) &
bind(C)
use, intrinsic :: iso_c_binding
integer (c_int64_t) , intent(in) , value :: context
integer (c_int32_t) , intent(in) , value :: lmax
integer (c_int64_t) , intent(in) , value :: ldl
integer (c_int64_t) , intent(in) , value :: ldv
real (c_double) , intent(in) :: X(3), R(3)
integer (c_int64_t) , intent(out) :: n
integer (c_int32_t) , intent(out) :: L(ldl,(lmax+1)*(lmax+2)*(lmax+3)/6)
real (c_double) , intent(out) :: VGL(ldv,(lmax+1)*(lmax+2)*(lmax+3)/6)
end function qmckl_ao_polynomial_vgl
end interface
#+END_SRC
*** Test :noexport: *** Test :noexport:
#+BEGIN_SRC C :comments link :tangle test_qmckl_ao.c #+BEGIN_SRC f90 :comments link :tangle test_qmckl_ao_f.f90
{ integer(c_int32_t) function test_qmckl_ao_polynomial_vgl(context) bind(C)
#include <stdio.h> use, intrinsic :: iso_c_binding
double X[3] = { 1.1 , 2.2 , 3.3 }; implicit none
double R[3] = { 0.1 , 1.2 , -2.3 }; include 'qmckl_ao.fh'
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; integer(c_int64_t), intent(in), value :: context
L_mem = (int32_t*) malloc(ldl*100*sizeof(int32_t)); integer :: lmax, d, i
VGL_mem = (double*) malloc(ldv*100*sizeof(double)); integer, allocatable :: L(:,:)
integer*8 :: n, ldl, ldv, j
double precision :: X(3), R(3), Y(3)
double precision, allocatable :: VGL(:,:)
double precision :: w
munit_assert_int64(QMCKL_SUCCESS, ==, X = (/ 1.1 , 2.2 , 3.3 /)
qmckl_ao_polynomial_vgl(context, X, R, lmax, &n, L_mem, ldl, VGL_mem, ldv) ); R = (/ 0.1 , 1.2 , -2.3 /)
Y(:) = X(:) - R(:)
munit_assert_int64( n, ==, d ); lmax = 4;
for (j=0 ; j<n ; j++) { n = 0;
L[j] = &L_mem[j*ldl]; ldl = 3;
VGL[j] = &VGL_mem[j*ldv]; ldv = 100;
}
Y[0] = X[0] - R[0]; d = (lmax+1)*(lmax+2)*(lmax+3)/6
Y[1] = X[1] - R[1];
Y[2] = X[2] - R[2];
for (j=0 ; j<n ; j++) {
munit_assert_int64( L[j][0], >=, 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.; allocate (L(ldl,100), VGL(ldv,100))
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]); test_qmckl_ao_polynomial_vgl = &
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); qmckl_ao_polynomial_vgl(context, X, R, lmax, n, L, ldl, VGL, ldv)
munit_assert_double_equal( VGL[j][4], w, 10 ); if (test_qmckl_ao_polynomial_vgl /= 0) return
}
free(L_mem); test_qmckl_ao_polynomial_vgl = -1
free(VGL_mem);
} if (n /= d) return
do j=1,n
do i=1,3
if (L(i,j) < 0) return
end do
if (dabs(1.d0 - VGL(1,j) / (&
Y(1)**L(1,j) * Y(2)**L(2,j) * Y(3)**L(3,j) &
)) > 1.d-14 ) return
if (L(1,j) < 1) then
if (VGL(2,j) /= 0.d0) return
else
if (dabs(1.d0 - VGL(2,j) / (&
L(1,j) * Y(1)**(L(1,j)-1) * Y(2)**L(2,j) * Y(3)**L(3,j) &
)) > 1.d-14 ) return
end if
if (L(2,j) < 1) then
if (VGL(3,j) /= 0.d0) return
else
if (dabs(1.d0 - VGL(3,j) / (&
L(2,j) * Y(1)**L(1,j) * Y(2)**(L(2,j)-1) * Y(3)**L(3,j) &
)) > 1.d-14 ) return
end if
if (L(3,j) < 1) then
if (VGL(4,j) /= 0.d0) return
else
if (dabs(1.d0 - VGL(4,j) / (&
L(3,j) * Y(1)**L(1,j) * Y(2)**L(2,j) * Y(3)**(L(3,j)-1) &
)) > 1.d-14 ) return
end if
w = 0.d0
if (L(1,j) > 1) then
w = w + L(1,j) * (L(1,j)-1) * Y(1)**(L(1,j)-2) * Y(2)**L(2,j) * Y(3)**L(3,j)
end if
if (L(2,j) > 1) then
w = w + L(2,j) * (L(2,j)-1) * Y(1)**L(1,j) * Y(2)**(L(2,j)-2) * Y(3)**L(3,j)
end if
if (L(3,j) > 1) then
w = w + L(3,j) * (L(3,j)-1) * Y(1)**L(1,j) * Y(2)**L(2,j) * Y(3)**(L(3,j)-2)
end if
if (dabs(1.d0 - VGL(5,j) / w) > 1.d-14 ) return
end do
test_qmckl_ao_polynomial_vgl = 0
deallocate(L,VGL)
end function test_qmckl_ao_polynomial_vgl
#+END_SRC #+END_SRC
#+BEGIN_SRC C :comments link :tangle test_qmckl_ao.c
int test_qmckl_ao_polynomial_vgl(qmckl_context context);
munit_assert_int(0, ==, test_qmckl_ao_polynomial_vgl(context));
#+END_SRC
#+END_SRC
* TODO Gaussian basis functions * TODO Gaussian basis functions

120
src/qmckl_distance.org
View File

@ -14,7 +14,8 @@ Function for the computation of distances between particles.
3 files are produced: 3 files are produced:
- a header file : =qmckl_distance.h= - a header file : =qmckl_distance.h=
- a source file : =qmckl_distance.f90= - a source file : =qmckl_distance.f90=
- a test file : =test_qmckl_distance.c= - a C test file : =test_qmckl_distance.c=
- a Fortran test file : =test_qmckl_distance_f.f90=
*** Header :noexport: *** Header :noexport:
#+BEGIN_SRC C :comments link :tangle qmckl_distance.h #+BEGIN_SRC C :comments link :tangle qmckl_distance.h
@ -23,10 +24,6 @@ Function for the computation of distances between particles.
#include "qmckl_context.h" #include "qmckl_context.h"
#+END_SRC #+END_SRC
*** Source :noexport:
#+BEGIN_SRC f90 :comments link :tangle qmckl_distance.f90
#+END_SRC
*** Test :noexport: *** Test :noexport:
#+BEGIN_SRC C :comments link :tangle test_qmckl_distance.c #+BEGIN_SRC C :comments link :tangle test_qmckl_distance.c
#include <math.h> #include <math.h>
@ -132,10 +129,10 @@ integer function qmckl_distance_sq_f(context, m, n, A, LDA, B, LDB, C, LDC) resu
do j=1,n do j=1,n
do i=1,m do i=1,m
x = A(i,1) - B(j,1) x = A(i,1) - B(j,1)
y = A(i,2) - B(j,2) y = A(i,2) - B(j,2)
z = A(i,3) - B(j,3) z = A(i,3) - B(j,3)
C(i,j) = x*x + y*y + z*z C(i,j) = x*x + y*y + z*z
end do end do
end do end do
@ -162,60 +159,77 @@ integer(c_int32_t) function qmckl_distance_sq(context, m, n, A, LDA, B, LDB, C,
end function qmckl_distance_sq end function qmckl_distance_sq
#+END_SRC #+END_SRC
#+BEGIN_SRC f90 :comments link :tangle qmckl_distance.fh
interface
integer(c_int32_t) function qmckl_distance_sq(context, 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
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)
end function qmckl_distance_sq
end interface
#+END_SRC
*** Test :noexport: *** Test :noexport:
#+BEGIN_SRC f90 :comments link :tangle test_qmckl_distance_f.f90 #+BEGIN_SRC f90 :comments link :tangle test_qmckl_distance_f.f90
integer(c_int32_t) function test_qmckl_distance_sq(context) bind(C) integer(c_int32_t) function test_qmckl_distance_sq(context) bind(C)
use iso_c_binding use, intrinsic :: iso_c_binding
implicit none implicit none
integer(c_int64_t), intent(in), value :: context include 'qmckl_distance.fh'
integer(c_int64_t), intent(in), value :: context
double precision, allocatable :: A(:,:), B(:,:), C(:,:) double precision, allocatable :: A(:,:), B(:,:), C(:,:)
integer*8 :: m, n, LDA, LDB, LDC integer*8 :: m, n, LDA, LDB, LDC
double precision :: x double precision :: x
integer*8 :: i,j integer*8 :: i,j
integer, external :: qmckl_distance_sq_f m = 5
n = 6
LDA = 6
LDB = 10
LDC = 5
m = 5 allocate( A(LDA,3), B(LDB,3), C(LDC,n) )
n = 6
LDA = 6
LDB = 10
LDC = 5
allocate( A(LDA,3), B(LDB,3), C(LDC,n) ) do j=1,3
do i=1,m
A(i,j) = -10.d0 + dble(i+j)
end do
do i=1,n
B(i,j) = -1.d0 + dble(i*j)
end do
end do
do j=1,3 test_qmckl_distance_sq = qmckl_distance_sq(context, m, n, A, LDA, B, LDB, C, LDC)
do i=1,m if (test_qmckl_distance_sq /= 0) return
A(i,j) = -10.d0 + dble(i+j)
end do
do i=1,n
B(i,j) = -1.d0 + dble(i*j)
end do
end do
test_qmckl_distance_sq = qmckl_distance_sq_f(context, m, n, A, LDA, B, LDB, C, LDC) test_qmckl_distance_sq = -1
if (test_qmckl_distance_sq /= 0) return
test_qmckl_distance_sq = -1 do j=1,n
do i=1,m
x = (A(i,1)-B(j,1))**2 + &
(A(i,2)-B(j,2))**2 + &
(A(i,3)-B(j,3))**2
if ( dabs(1.d0 - C(i,j)/x) > 1.d-12 ) return
end do
end do
test_qmckl_distance_sq = 0
do j=1,n deallocate(A,B,C)
do i=1,m end function test_qmckl_distance_sq
x = (A(i,1)-B(j,1))**2 + & #+END_SRC
(A(i,2)-B(j,2))**2 + &
(A(i,3)-B(j,3))**2
if ( dabs(1.d0 - C(i,j)/x) > 1.d-12 ) return
end do
end do
test_qmckl_distance_sq = 0
deallocate(A,B,C) #+BEGIN_SRC C :comments link :tangle test_qmckl_distance.c
end function test_qmckl_distance_sq int test_qmckl_distance_sq(qmckl_context context);
#+END_SRC munit_assert_int(0, ==, test_qmckl_distance_sq(context));
#+END_SRC
#+BEGIN_SRC C :comments link :tangle test_qmckl_distance.c
int test_qmckl_distance_sq(qmckl_context context);
munit_assert_int(0, ==, test_qmckl_distance_sq(context));
#+END_SRC
* End of files * End of files
*** Header :noexport: *** Header :noexport: