diff --git a/src/Makefile b/src/Makefile index d98f8e6..c4e7fc1 100644 --- a/src/Makefile +++ b/src/Makefile @@ -6,6 +6,14 @@ FFLAGS=-fPIC -g -fcheck=all -Waliasing -Wampersand -Wconversion -Wsurprising -Wi LIBS=-lgfortran -lm +#CC=icc +#CFLAGS=-fPIC -g +# +#FC=ifort +#FFLAGS=-fPIC -g +# +#LIBS=-lm -lifcore -lirc + export CC CFLAGS FC FFLAGS LIBS @@ -25,7 +33,7 @@ doc:$(ORG_SOURCE_FILES) ./create_doc.sh $(ORG_SOURCE_FILES) 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 ./create_makefile.sh $(ORG_SOURCE_FILES) diff --git a/src/README.org b/src/README.org index 698ea34..01aa7ae 100644 --- a/src/README.org +++ b/src/README.org @@ -72,9 +72,13 @@ rm ${nb}.md ** Writing in Fortran 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. - + 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 # 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 bindings in other languages in other repositories. - ** Global state Global variables should be avoided in the library, because it is diff --git a/src/qmckl_ao.org b/src/qmckl_ao.org index 96d1677..5511e74 100644 --- a/src/qmckl_ao.org +++ b/src/qmckl_ao.org @@ -12,10 +12,11 @@ This files contains all the routines for the computation of the values, gradients and Laplacian of the atomic basis functions. -3 files are produced: +4 files are produced: - a header file : =qmckl_ao.h= - 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: #+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" #+END_SRC -*** Source :noexport: - #+BEGIN_SRC f90 :comments link :tangle qmckl_ao.f90 - - #+END_SRC - *** Test :noexport: #+BEGIN_SRC C :comments link :tangle test_qmckl_ao.c #include @@ -131,40 +127,64 @@ integer(c_int32_t) function qmckl_ao_powers(context, n, X, LMAX, P, ldp) & end function qmckl_ao_powers #+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: - #+BEGIN_SRC C :comments link :tangle test_qmckl_ao.c -{ - int64_t n, LDP ; - int32_t *LMAX ; - double *X, *P ; - int i, j; + #+BEGIN_SRC f90 :comments link :tangle test_qmckl_ao_f.f90 +integer(c_int32_t) function test_qmckl_ao_powers(context) bind(C) + use, intrinsic :: iso_c_binding + implicit none + include 'qmckl_ao.fh' + 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; LDP = 10; + + allocate(X(n), P(LDP,n), LMAX(n)) + + do j=1,n + X(j) = -5.d0 + 0.1d0 * dble(j) + LMAX(j) = 1 + int(mod(j, 9),4) + end do + + test_qmckl_ao_powers = qmckl_ao_powers(context, n, X, LMAX, P, LDP) + if (test_qmckl_ao_powers /= 0) return + + test_qmckl_ao_powers = -1 + + 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 - 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)); + test_qmckl_ao_powers = 0 + deallocate(X,P,LMAX) +end function test_qmckl_ao_powers + #+END_SRC - for (j=0 ; j - 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; + #+BEGIN_SRC f90 :comments link :tangle test_qmckl_ao_f.f90 +integer(c_int32_t) function test_qmckl_ao_polynomial_vgl(context) bind(C) + use, intrinsic :: iso_c_binding + implicit none + include 'qmckl_ao.fh' - int d = (lmax+1)*(lmax+2)*(lmax+3)/6; + integer(c_int64_t), intent(in), value :: context + + integer :: lmax, d, i + integer, allocatable :: L(:,:) + integer*8 :: n, ldl, ldv, j + double precision :: X(3), R(3), Y(3) + double precision, allocatable :: VGL(:,:) + double precision :: w - L_mem = (int32_t*) malloc(ldl*100*sizeof(int32_t)); - VGL_mem = (double*) malloc(ldv*100*sizeof(double)); + X = (/ 1.1 , 2.2 , 3.3 /) + R = (/ 0.1 , 1.2 , -2.3 /) + Y(:) = X(:) - R(:) - munit_assert_int64(QMCKL_SUCCESS, ==, - qmckl_ao_polynomial_vgl(context, X, R, lmax, &n, L_mem, ldl, VGL_mem, ldv) ); + lmax = 4; + n = 0; + ldl = 3; + ldv = 100; - 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 ); - } + allocate (L(ldl,100), VGL(ldv,100)) - 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); -} + test_qmckl_ao_polynomial_vgl = & + qmckl_ao_polynomial_vgl(context, X, R, lmax, n, L, ldl, VGL, ldv) + if (test_qmckl_ao_polynomial_vgl /= 0) return + + test_qmckl_ao_polynomial_vgl = -1 + + 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 + #+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 diff --git a/src/qmckl_distance.org b/src/qmckl_distance.org index d69b133..b34b939 100644 --- a/src/qmckl_distance.org +++ b/src/qmckl_distance.org @@ -14,7 +14,8 @@ 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= +- a C test file : =test_qmckl_distance.c= +- a Fortran test file : =test_qmckl_distance_f.f90= *** Header :noexport: #+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" #+END_SRC -*** Source :noexport: - #+BEGIN_SRC f90 :comments link :tangle qmckl_distance.f90 - #+END_SRC - *** Test :noexport: #+BEGIN_SRC C :comments link :tangle test_qmckl_distance.c #include @@ -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 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 + 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 @@ -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_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: - #+BEGIN_SRC f90 :comments link :tangle test_qmckl_distance_f.f90 - integer(c_int32_t) function test_qmckl_distance_sq(context) bind(C) - use iso_c_binding - implicit none - integer(c_int64_t), intent(in), value :: context + #+BEGIN_SRC f90 :comments link :tangle test_qmckl_distance_f.f90 +integer(c_int32_t) function test_qmckl_distance_sq(context) bind(C) + use, intrinsic :: iso_c_binding + implicit none + include 'qmckl_distance.fh' + integer(c_int64_t), intent(in), value :: context - double precision, allocatable :: A(:,:), B(:,:), C(:,:) - integer*8 :: m, n, LDA, LDB, LDC - double precision :: x - integer*8 :: i,j + double precision, allocatable :: A(:,:), B(:,:), C(:,:) + integer*8 :: m, n, LDA, LDB, LDC + double precision :: x + integer*8 :: i,j - integer, external :: qmckl_distance_sq_f + m = 5 + n = 6 + LDA = 6 + LDB = 10 + LDC = 5 - m = 5 - n = 6 - LDA = 6 - LDB = 10 - LDC = 5 + allocate( A(LDA,3), B(LDB,3), C(LDC,n) ) - 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 - 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 + test_qmckl_distance_sq = qmckl_distance_sq(context, m, n, A, LDA, B, LDB, C, LDC) + if (test_qmckl_distance_sq /= 0) return - test_qmckl_distance_sq = qmckl_distance_sq_f(context, m, n, A, LDA, B, LDB, C, LDC) - if (test_qmckl_distance_sq /= 0) return + test_qmckl_distance_sq = -1 - 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 - 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 + deallocate(A,B,C) +end function test_qmckl_distance_sq + #+END_SRC - deallocate(A,B,C) - end function test_qmckl_distance_sq - #+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 + #+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 *** Header :noexport: