TREX/qmckl
1
0
Fork 0
mirror of https://github.com/TREX-CoE/qmckl.git synced 2024-12-22 20:36:01 +01:00
Code Issues Projects Releases Wiki Activity

Optimizations in polynomials

Browse Source
This commit is contained in:
Anthony Scemama 2020-10-28 19:28:27 +01:00
parent abc0e6b5ab
commit 8147ad22a7
2 changed files with 53 additions and 43 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

5
TODO.org
View File

@ -10,3 +10,8 @@ qmckl_malloc, where the domain id is something obtained from the
context. context.
* TRANSA, TRANSB
* Performance info
* Benchmark interpolation of basis functions
* Complex numbers

91
src/qmckl_ao.org
View File

@ -85,9 +85,9 @@ MunitResult test_qmckl_ao() {
*** Header *** Header
#+BEGIN_SRC C :comments link :tangle qmckl_ao.h #+BEGIN_SRC C :comments link :tangle qmckl_ao.h
qmckl_exit_code qmckl_ao_powers(qmckl_context context, qmckl_exit_code qmckl_ao_powers(qmckl_context context,
int64_t n, int64_t n,
double *X, int32_t *LMAX, double *X, int32_t *LMAX,
double *P, int64_t LDP); double *P, int64_t LDP);
#+END_SRC #+END_SRC
*** Source *** Source
@ -235,10 +235,10 @@ munit_assert_int(0, ==, test_qmckl_ao_powers(context));
*** Header *** Header
#+BEGIN_SRC C :comments link :tangle qmckl_ao.h #+BEGIN_SRC C :comments link :tangle qmckl_ao.h
qmckl_exit_code qmckl_ao_polynomial_vgl(qmckl_context context, qmckl_exit_code qmckl_ao_polynomial_vgl(qmckl_context context,
double *X, double *R, double *X, double *R,
int32_t lmax, int64_t *n, int32_t lmax, int64_t *n,
int32_t *L, int64_t ldl, int32_t *L, int64_t ldl,
double *VGL, int64_t ldv); double *VGL, int64_t ldv);
#+END_SRC #+END_SRC
*** Source *** Source
@ -261,27 +261,28 @@ integer function qmckl_ao_polynomial_vgl_f(context, X, R, lmax, n, L, ldl, VGL,
real*8 :: pows(-2:lmax,3) real*8 :: pows(-2:lmax,3)
integer, external :: qmckl_ao_powers_f integer, external :: qmckl_ao_powers_f
double precision :: xy, yz, xz double precision :: xy, yz, xz
double precision :: da, db, dc, dd
info = 0 info = 0
if (context == 0_8) then if (context == 0_8) then
info = -1 info = -1
return return
endif endif
n = (lmax+1)*(lmax+2)*(lmax+3)/6 n = (lmax+1)*(lmax+2)*(lmax+3)/6
if (ldl < 3) then if (ldl < 3) then
info = -2 info = -2
return return
endif endif
if (ldv < 5) then if (ldv < 5) then
info = -3 info = -3
return return
endif endif
do i=1,3 do i=1,3
Y(i) = X(i) - R(i) Y(i) = X(i) - R(i)
end do end do
@ -300,39 +301,43 @@ integer function qmckl_ao_polynomial_vgl_f(context, X, R, lmax, n, L, ldl, VGL,
vgl(1:5,1:n) = 0.d0 vgl(1:5,1:n) = 0.d0
l(1:3,n) = 0 l(1:3,n) = 0
vgl(1,n) = 1.d0 vgl(1,n) = 1.d0
dd = 1.d0
do d=1,lmax do d=1,lmax
da = 0.d0
do a=0,d do a=0,d
do b=0,d db = 0.d0
do c=0,d do b=0,d-a
if (a+b+c == d) then c = d - a - b
n = n+1 dc = dd - da - db
l(1,n) = a n = n+1
l(2,n) = b l(1,n) = a
l(3,n) = c l(2,n) = b
l(3,n) = c
xy = pows(a,1) * pows(b,2)
yz = pows(b,2) * pows(c,3)
xz = pows(a,1) * pows(c,3)
vgl(1,n) = xy * pows(c,3)
xy = dc * xy
xz = db * xz
yz = da * yz
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(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
xy = pows(a,1) * pows(b,2) db = db + 1.d0
yz = pows(b,2) * pows(c,3) end do
xz = pows(a,1) * pows(c,3) da = da + 1.d0
vgl(1,n) = xy * pows(c,3)
xy = dble(c) * xy
xz = dble(b) * xz
yz = dble(a) * yz
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(5,n) = &
dble(a-1) * pows(a-2,1) * yz + &
dble(b-1) * pows(b-2,2) * xz + &
dble(c-1) * pows(c-2,3) * xy
exit
end if
end do
end do
end do end do
dd = dd + 1.d0
end do end do
end function qmckl_ao_polynomial_vgl_f end function qmckl_ao_polynomial_vgl_f