From abbc12e160e68cb3c02301ab125c6d6b9468cad5 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Fri, 6 Nov 2020 10:58:22 +0100 Subject: [PATCH 1/3] Transposed VGL in ao_polynomials --- src/Makefile | 33 +++++++++++++++++++----- src/qmckl_ao.org | 67 +++++++++++++++++++++++++++--------------------- 2 files changed, 64 insertions(+), 36 deletions(-) diff --git a/src/Makefile b/src/Makefile index 28de44e..295b783 100644 --- a/src/Makefile +++ b/src/Makefile @@ -1,3 +1,8 @@ +COMPILER=GNU +#COMPILER=INTEL +#COMPILER=LLVM + +ifeq($(COMPILER),GNU) CC=gcc -g CFLAGS=-fPIC -fexceptions -Wall -Werror -Wpedantic -Wextra @@ -5,14 +10,28 @@ FC=gfortran -g FFLAGS=-fPIC -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 +endif -#CC=icc -xHost -#CFLAGS=-fPIC -g -O2 -# -#FC=ifort -xHost -#FFLAGS=-fPIC -g -O2 -# -#LIBS=-lm -lifcore -lirc +ifeq($(COMPILER),INTEL) +CC=icc -xHost +CFLAGS=-fPIC -g -O2 + +FC=ifort -xHost +FFLAGS=-fPIC -g -O2 + +LIBS=-lm -lifcore -lirc +endif + +#TODO +ifeq($(COMPILER),LLVM) +CC=clang +CFLAGS=-fPIC -g -O2 + +FC=flang +FFLAGS=fPIC -g -O2 + +LIBS=-lm +endif export CC CFLAGS FC FFLAGS LIBS diff --git a/src/qmckl_ao.org b/src/qmckl_ao.org index 506beb8..fbf97b7 100644 --- a/src/qmckl_ao.org +++ b/src/qmckl_ao.org @@ -205,7 +205,7 @@ munit_assert_int(0, ==, test_qmckl_ao_powers(context)); | =n= | output | Number of computed polynomials | | =L(ldl,n)= | output | Contains a,b,c for all =n= results | | =ldl= | input | Leading dimension of =L= | - | =VGL(ldv,5)= | output | Value, gradients and Laplacian of the polynomials | + | =VGL(ldv,n)= | output | Value, gradients and Laplacian of the polynomials | | =ldv= | input | Leading dimension of array =VGL= | ***** Requirements @@ -214,12 +214,25 @@ munit_assert_int(0, ==, test_qmckl_ao_powers(context)); - =n= > 0 - =lmax= >= 0 - =ldl= >= 3 - - =ldv= >= (=lmax=+1)(=lmax=+2)(=lmax=+3)/6 + - =ldv= >= 5 - =X= is allocated with at least $3 \times 8$ bytes - =R= is allocated with at least $3 \times 8$ bytes + - =n= >= =(lmax+1)(lmax+2)(lmax+3)/6= - =L= is allocated with at least $3 \times n \times 4$ bytes - - =VGL= is allocated with at least $n \times 5 \times 8$ bytes - - On output, =n= should be equal to (=lmax=+1)(=lmax=+2)(=lmax=+3)/6 + - =VGL= is allocated with at least $5 \times n \times 8$ bytes + - On output, =n= should be equal to =(lmax+1)(lmax+2)(lmax+3)/6= + - On output, the powers are given in the following order (l=a+b+c): + - Increase values of =l= + - Within a given value of =l=, alphabetical order of the + string made by a*"x" + b*"y" + c*"z" (in Python notation). + For example, with a=0, b=2 and c=1 the string is "yyz" + +***** Error codes + + | -1 | Null context | + | -2 | Inconsistent =ldl= | + | -3 | Inconsistent =ldv= | + | -4 | Inconsistent =lmax= | ***** Header #+BEGIN_SRC C :tangle qmckl.h @@ -240,7 +253,7 @@ integer function qmckl_ao_polynomial_vgl_f(context, X, R, lmax, n, L, ldl, VGL, 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,5) + real*8 , intent(out) :: VGL(ldv,(lmax+1)*(lmax+2)*(lmax+3)/6) integer*8 , intent(in) :: ldv integer*8 :: i,j @@ -264,7 +277,7 @@ integer function qmckl_ao_polynomial_vgl_f(context, X, R, lmax, n, L, ldl, VGL, return endif - if (ldv < (lmax+1)*(lmax+2)*(lmax+3)/6) then + if (ldv < 5) then info = -3 return endif @@ -289,8 +302,8 @@ integer function qmckl_ao_polynomial_vgl_f(context, X, R, lmax, n, L, ldl, VGL, if (info /= 0) return - vgl(1,1) = 1.d0 - vgl(1,2:5) = 0.d0 + VGL(1,1) = 1.d0 + vgL(2:5,1) = 0.d0 l(1:3,1) = 0 n=1 dd = 1.d0 @@ -310,17 +323,17 @@ integer function qmckl_ao_polynomial_vgl_f(context, X, R, lmax, n, L, ldl, VGL, yz = pows(b,2) * pows(c,3) xz = pows(a,1) * pows(c,3) - vgl(n,1) = xy * pows(c,3) + vgl(1,n) = xy * pows(c,3) xy = dc * xy xz = db * xz yz = da * yz - vgl(n,2) = pows(a-1,1) * yz - vgl(n,3) = pows(b-1,2) * xz - vgl(n,4) = pows(c-1,3) * xy + vgl(2,n) = pows(a-1,1) * yz + vgl(3,n) = pows(b-1,2) * xz + vgl(4,n) = pows(c-1,3) * xy - vgl(n,5) = & + vgl(5,n) = & (da-1.d0) * pows(a-2,1) * yz + & (db-1.d0) * pows(b-2,2) * xz + & (dc-1.d0) * pows(c-2,3) * xy @@ -332,11 +345,6 @@ integer function qmckl_ao_polynomial_vgl_f(context, X, R, lmax, n, L, ldl, VGL, dd = dd + 1.d0 end do - if (n /= (lmax+1)*(lmax+2)*(lmax+3)/6) then - info = -5 - return - endif - info = 0 end function qmckl_ao_polynomial_vgl_f @@ -354,7 +362,7 @@ integer(c_int32_t) function qmckl_ao_polynomial_vgl(context, X, R, lmax, n, L, l 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,5) + 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 @@ -362,6 +370,7 @@ integer(c_int32_t) function qmckl_ao_polynomial_vgl(context, X, R, lmax, n, L, l end function qmckl_ao_polynomial_vgl #+END_SRC +***** Fortran interface :noexport: #+BEGIN_SRC f90 :tangle qmckl_f.f90 interface integer(c_int32_t) function qmckl_ao_polynomial_vgl(context, X, R, lmax, n, L, ldl, VGL, ldv) & @@ -374,7 +383,7 @@ end function qmckl_ao_polynomial_vgl 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,5) + real (c_double) , intent(out) :: VGL(ldv,(lmax+1)*(lmax+2)*(lmax+3)/6) end function qmckl_ao_polynomial_vgl end interface #+END_SRC @@ -407,7 +416,7 @@ integer(c_int32_t) function test_qmckl_ao_polynomial_vgl(context) bind(C) d = (lmax+1)*(lmax+2)*(lmax+3)/6 - allocate (L(ldl,100), VGL(ldv,5)) + allocate (L(ldl,d), VGL(ldv,d)) test_qmckl_ao_polynomial_vgl = & qmckl_ao_polynomial_vgl(context, X, R, lmax, n, L, ldl, VGL, ldv) @@ -423,33 +432,33 @@ integer(c_int32_t) function test_qmckl_ao_polynomial_vgl(context) bind(C) if (L(i,j) < 0) return end do test_qmckl_ao_polynomial_vgl = -12 - if (dabs(1.d0 - VGL(j,1) / (& + if (dabs(1.d0 - VGL(1,j) / (& Y(1)**L(1,j) * Y(2)**L(2,j) * Y(3)**L(3,j) & )) > epsilon ) return test_qmckl_ao_polynomial_vgl = -13 if (L(1,j) < 1) then - if (VGL(j,2) /= 0.d0) return + if (VGL(2,j) /= 0.d0) return else - if (dabs(1.d0 - VGL(j,2) / (& + 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) & )) > epsilon ) return end if test_qmckl_ao_polynomial_vgl = -14 if (L(2,j) < 1) then - if (VGL(j,3) /= 0.d0) return + if (VGL(3,j) /= 0.d0) return else - if (dabs(1.d0 - VGL(j,3) / (& + 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) & )) > epsilon ) return end if test_qmckl_ao_polynomial_vgl = -15 if (L(3,j) < 1) then - if (VGL(j,4) /= 0.d0) return + if (VGL(4,j) /= 0.d0) return else - if (dabs(1.d0 - VGL(j,4) / (& + 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) & )) > epsilon ) return end if @@ -465,7 +474,7 @@ integer(c_int32_t) function test_qmckl_ao_polynomial_vgl(context) bind(C) 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(j,5) / w) > epsilon ) return + if (dabs(1.d0 - VGL(5,j) / w) > epsilon ) return end do test_qmckl_ao_polynomial_vgl = 0 From 3852dad53d56b7ad0c8810c5f899ae49087b8d5c Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Fri, 6 Nov 2020 12:10:20 +0100 Subject: [PATCH 2/3] Reordering of powers in polynomials --- src/Makefile | 6 +++--- src/qmckl.org | 1 + src/qmckl_ao.org | 15 +++++++-------- src/qmckl_context.org | 4 ++-- src/qmckl_distance.org | 1 - 5 files changed, 13 insertions(+), 14 deletions(-) diff --git a/src/Makefile b/src/Makefile index 295b783..7f78a9b 100644 --- a/src/Makefile +++ b/src/Makefile @@ -2,7 +2,7 @@ COMPILER=GNU #COMPILER=INTEL #COMPILER=LLVM -ifeq($(COMPILER),GNU) +ifeq ($(COMPILER),GNU) CC=gcc -g CFLAGS=-fPIC -fexceptions -Wall -Werror -Wpedantic -Wextra @@ -12,7 +12,7 @@ FFLAGS=-fPIC -fcheck=all -Waliasing -Wampersand -Wconversion -Wsurprising -Wint LIBS=-lgfortran -lm endif -ifeq($(COMPILER),INTEL) +ifeq ($(COMPILER),INTEL) CC=icc -xHost CFLAGS=-fPIC -g -O2 @@ -23,7 +23,7 @@ LIBS=-lm -lifcore -lirc endif #TODO -ifeq($(COMPILER),LLVM) +ifeq ($(COMPILER),LLVM) CC=clang CFLAGS=-fPIC -g -O2 diff --git a/src/qmckl.org b/src/qmckl.org index 7620883..ce5a3ce 100644 --- a/src/qmckl.org +++ b/src/qmckl.org @@ -12,6 +12,7 @@ #define QMCKL_H #include #include +#include #+END_SRC #+BEGIN_SRC f90 :tangle qmckl_f.f90 diff --git a/src/qmckl_ao.org b/src/qmckl_ao.org index fbf97b7..4371cd9 100644 --- a/src/qmckl_ao.org +++ b/src/qmckl_ao.org @@ -11,7 +11,6 @@ *** Test :noexport: #+BEGIN_SRC C :tangle test_qmckl_ao.c -#include #include "qmckl.h" #include "munit.h" MunitResult test_qmckl_ao() { @@ -308,10 +307,10 @@ integer function qmckl_ao_polynomial_vgl_f(context, X, R, lmax, n, L, ldl, VGL, n=1 dd = 1.d0 do d=1,lmax - da = 0.d0 - do a=0,d - db = 0.d0 - do b=0,d-a + da = dd + do a=d,0,-1 + db = dd-da + do b=d-a,0,-1 c = d - a - b dc = dd - da - db n = n+1 @@ -338,9 +337,9 @@ integer function qmckl_ao_polynomial_vgl_f(context, X, R, lmax, n, L, ldl, VGL, (db-1.d0) * pows(b-2,2) * xz + & (dc-1.d0) * pows(c-2,3) * xy - db = db + 1.d0 + db = db - 1.d0 end do - da = da + 1.d0 + da = da - 1.d0 end do dd = dd + 1.d0 end do @@ -413,7 +412,7 @@ integer(c_int32_t) function test_qmckl_ao_polynomial_vgl(context) bind(C) n = 0; ldl = 3; ldv = 100; - + d = (lmax+1)*(lmax+2)*(lmax+3)/6 allocate (L(ldl,d), VGL(ldv,d)) diff --git a/src/qmckl_context.org b/src/qmckl_context.org index e407466..5359470 100644 --- a/src/qmckl_context.org +++ b/src/qmckl_context.org @@ -461,7 +461,7 @@ int qmckl_context_get_range(const qmckl_context context) { ***** TODO Tests :noexport: **** =qmckl_context_get_epsilon= - Returns $\epsilon = 2 / \log_{10} 2^{n-1}$ where =n= is the precision + Returns $\epsilon = 2^{1-n}$ where =n= is the precision #+BEGIN_SRC C :comments org :tangle qmckl.h double qmckl_context_get_epsilon(const qmckl_context context); #+END_SRC @@ -470,7 +470,7 @@ double qmckl_context_get_epsilon(const qmckl_context context); #+BEGIN_SRC C :tangle qmckl_context.c double qmckl_context_get_epsilon(const qmckl_context context) { const qmckl_context_struct* ctx = (qmckl_context_struct*) context; - return 1.0 / ((double) ((int64_t) 1 << (ctx->precision-1))); + return pow(2.0,(double) 1-ctx->precision); } #+END_SRC diff --git a/src/qmckl_distance.org b/src/qmckl_distance.org index 06c8015..5eac91d 100644 --- a/src/qmckl_distance.org +++ b/src/qmckl_distance.org @@ -9,7 +9,6 @@ **** Headers :noexport: #+BEGIN_SRC C :comments link :tangle test_qmckl_distance.c -#include #include "qmckl.h" #include "munit.h" MunitResult test_qmckl_distance() { From d50737687c451b88367d25d7af55a2e090deb452 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Sat, 7 Nov 2020 16:32:23 +0100 Subject: [PATCH 3/3] Inline function in polynomials --- src/qmckl_ao.org | 52 +++++++++++++++++++++++++++++++++++------------- 1 file changed, 38 insertions(+), 14 deletions(-) diff --git a/src/qmckl_ao.org b/src/qmckl_ao.org index 4371cd9..51736d9 100644 --- a/src/qmckl_ao.org +++ b/src/qmckl_ao.org @@ -290,23 +290,46 @@ integer function qmckl_ao_polynomial_vgl_f(context, X, R, lmax, n, L, ldl, VGL, 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 + if (lmax == 0) then + VGL(1,1) = 1.d0 + vgL(2:5,1) = 0.d0 + l(1:3,1) = 0 + n=1 + else if (lmax > 0) then + pows(-2:0,1:3) = 1.d0 + do i=1,lmax + pows(i,1) = pows(i-1,1) * Y(1) + pows(i,2) = pows(i-1,2) * Y(2) + pows(i,3) = pows(i-1,3) * Y(3) + end do - VGL(1,1) = 1.d0 - vgL(2:5,1) = 0.d0 - l(1:3,1) = 0 - n=1 - dd = 1.d0 - do d=1,lmax + VGL(1:5,1:4) = 0.d0 + l(1:3,1:4) = 0 + + VGL(1,1) = 1.d0 + vgl(1:5,2:4) = 0.d0 + + l(1,2) = 1 + vgl(1,2) = pows(1,1) + vgL(2,2) = 1.d0 + + l(2,3) = 1 + vgl(1,3) = pows(1,2) + vgL(3,3) = 1.d0 + + l(3,4) = 1 + vgl(1,4) = pows(1,3) + vgL(4,4) = 1.d0 + + n=4 + endif + + ! l>=2 + dd = 2.d0 + do d=2,lmax da = dd do a=d,0,-1 db = dd-da @@ -314,6 +337,7 @@ integer function qmckl_ao_polynomial_vgl_f(context, X, R, lmax, n, L, ldl, VGL, c = d - a - b dc = dd - da - db n = n+1 + l(1,n) = a l(2,n) = b l(3,n) = c