From 8147ad22a72b9a4a03652a6f8c962dee9079ea0b Mon Sep 17 00:00:00 2001
From: Anthony Scemama <scemama@irsamc.ups-tlse.fr>
Date: Wed, 28 Oct 2020 19:28:27 +0100
Subject: [PATCH] Optimizations in polynomials

---
 TODO.org         |  5 +++
 src/qmckl_ao.org | 91 +++++++++++++++++++++++++-----------------------
 2 files changed, 53 insertions(+), 43 deletions(-)

diff --git a/TODO.org b/TODO.org
index 2d5626b..250dd79 100644
--- a/TODO.org
+++ b/TODO.org
@@ -10,3 +10,8 @@ qmckl_malloc, where the domain id is something obtained from the
 context.
 
 
+* TRANSA, TRANSB
+* Performance info
+* Benchmark interpolation of basis functions
+* Complex numbers
+
diff --git a/src/qmckl_ao.org b/src/qmckl_ao.org
index 9f2399b..c8e57a4 100644
--- a/src/qmckl_ao.org
+++ b/src/qmckl_ao.org
@@ -85,9 +85,9 @@ MunitResult test_qmckl_ao() {
 *** 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);
+                int64_t n, 
+                double *X, int32_t *LMAX,
+                double *P, int64_t LDP);
     #+END_SRC
     
 *** Source
@@ -235,10 +235,10 @@ munit_assert_int(0, ==, test_qmckl_ao_powers(context));
 *** 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);
+                double *X, double *R,
+                int32_t lmax, int64_t *n,
+                int32_t *L,   int64_t ldl,
+                double *VGL,  int64_t ldv);
     #+END_SRC
     
 *** 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)
   integer, external :: qmckl_ao_powers_f
   double precision  :: xy, yz, xz
-
+  double precision  :: da, db, dc, dd
+  
   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
@@ -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
   l(1:3,n) = 0
   vgl(1,n) = 1.d0
+  dd = 1.d0
   do d=1,lmax
+     da = 0.d0
      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
+        db = 0.d0
+        do b=0,d-a
+           c  = d  - a  - b
+           dc = dd - da - db
+           n = n+1
+           l(1,n) = a
+           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)
-		 yz = pows(b,2) * pows(c,3)
-		 xz = pows(a,1) * pows(c,3)
-
-		 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
+           db = db + 1.d0
+        end do
+        da = da + 1.d0
      end do
+     dd = dd + 1.d0
   end do
 
 end function qmckl_ao_polynomial_vgl_f