(git:ccc2433)
d3_poly.F
Go to the documentation of this file.
1 !--------------------------------------------------------------------------------------------------!
2 ! CP2K: A general program to perform molecular dynamics simulations !
3 ! Copyright 2000-2024 CP2K developers group <https://cp2k.org> !
4 ! !
5 ! SPDX-License-Identifier: GPL-2.0-or-later !
6 !--------------------------------------------------------------------------------------------------!
7 
8 ! **************************************************************************************************
9 !> \brief
10 !> Routines to efficiently handle dense polynomial in 3 variables up to
11 !> a given degree.
12 !> Multiplication, partial evaluation, affine transform (change of reference
13 !> system), differentiation are efficiently implemented.
14 !> some functions accept or return several polynomial together,
15 !> these have to have all the same size, and are stored one after the other
16 !> in an unique 1d array. This gives them an easy handling and even seem to
17 !> be faster than the transposed layout.
18 !> \note
19 !> not all routines have been fully optimized.
20 !> original available also with a BSD style license
21 !> \author Fawzi Mohamed
22 ! **************************************************************************************************
23 MODULE d3_poly
24 
25  USE kinds, ONLY: dp
26 
27 !$ USE OMP_LIB, ONLY: omp_get_max_threads, omp_get_thread_num, omp_get_num_threads
28 
29 #include "../base/base_uses.f90"
30 
31  IMPLICIT NONE
32 
33  PRIVATE
34 
35  PUBLIC :: poly_size1, poly_size2, poly_size3, &
36  grad_size3, &
37  poly_affine_t3, &
38  poly_p_eval2b, &
39  poly_p_eval3b, poly_cp2k2d3, &
40  poly_padd_uneval3b, poly_padd_uneval2b, &
41  poly_affine_t3t, poly_d32cp2k, init_d3_poly_module
42 
43 #ifdef FD_DEBUG
44 #define IF_CHECK(x,y) x
45 #else
46 #define IF_CHECK(x,y) y
47 #endif
48 
49  CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'd3_poly'
50 ! maximum grad for cached values
51  INTEGER, PUBLIC, PARAMETER :: max_grad2 = 5
52  INTEGER, PUBLIC, PARAMETER :: max_grad3 = 3
53  INTEGER, PUBLIC, PARAMETER :: cached_dim1 = max_grad2 + 1
54  INTEGER, PUBLIC, PARAMETER :: cached_dim2 = (max_grad2 + 1)*(max_grad2 + 2)/2
55  INTEGER, PUBLIC, PARAMETER :: cached_dim3 = (max_grad3 + 1)*(max_grad3 + 2)*(max_grad3 + 3)/6
56 
57 ! cached index -> monomial exponents
58  LOGICAL, SAVE :: module_initialized = .false.
59  INTEGER, SAVE, DIMENSION(2, cached_dim2) :: a_mono_exp2
60  INTEGER, SAVE, DIMENSION(3, cached_dim3) :: a_mono_exp3
61  INTEGER, SAVE, DIMENSION(cached_dim3) :: a_reduce_idx3
62  INTEGER, SAVE, DIMENSION(3, cached_dim3) :: a_deriv_idx3
63  INTEGER, SAVE, DIMENSION(cached_dim2, cached_dim2) :: a_mono_mult2
64  INTEGER, SAVE, DIMENSION(cached_dim3, cached_dim3) :: a_mono_mult3
65  INTEGER, SAVE, DIMENSION(4, cached_dim3) :: a_mono_mult3a
66 
67 CONTAINS
68 
69 ! **************************************************************************************************
70 !> \brief initialization of the cache, is called by functions in this module
71 !> that use cached values
72 ! **************************************************************************************************
73  SUBROUTINE init_d3_poly_module()
74  INTEGER :: grad, i, ii, ij, j, nthreads, subg
75  INTEGER, DIMENSION(2) :: monores2
76  INTEGER, DIMENSION(3) :: monores3
77 
78  nthreads = 1
79 !$ nthreads = OMP_GET_NUM_THREADS()
80  IF (nthreads /= 1) cpabort("init_d3_poly_module must not be called within OMP PARALLEL")
81 
82  IF (module_initialized) RETURN
83 
84  ii = 1
85  DO grad = 0, max_grad2
86  DO i = grad, 0, -1
87  a_mono_exp2(1, ii) = i
88  a_mono_exp2(2, ii) = grad - i
89  ii = ii + 1
90  END DO
91  END DO
92  ii = 1
93  DO grad = 0, max_grad3
94  DO i = grad, 0, -1
95  DO j = grad - i, 0, -1
96  a_mono_exp3(1, ii) = i
97  a_mono_exp3(2, ii) = j
98  a_mono_exp3(3, ii) = grad - i - j
99  ii = ii + 1
100  END DO
101  END DO
102  END DO
103  DO ii = 1, cached_dim3
104  subg = a_mono_exp3(2, ii) + a_mono_exp3(3, ii)
105  a_reduce_idx3(ii) = subg*(subg + 1)/2 + a_mono_exp3(3, ii) + 1
106  END DO
107  DO ii = 1, cached_dim3
108  IF (a_mono_exp3(1, ii) > 0) THEN
109  a_deriv_idx3(1, ii) = mono_index3(a_mono_exp3(1, ii) - 1, a_mono_exp3(2, ii), a_mono_exp3(3, ii))
110  ELSE
111  a_deriv_idx3(1, ii) = 0
112  END IF
113  IF (a_mono_exp3(2, ii) > 0) THEN
114  a_deriv_idx3(2, ii) = mono_index3(a_mono_exp3(1, ii), a_mono_exp3(2, ii) - 1, a_mono_exp3(3, ii))
115  ELSE
116  a_deriv_idx3(2, ii) = 0
117  END IF
118  IF (a_mono_exp3(3, ii) > 0) THEN
119  a_deriv_idx3(3, ii) = mono_index3(a_mono_exp3(1, ii), a_mono_exp3(2, ii), a_mono_exp3(3, ii) - 1)
120  ELSE
121  a_deriv_idx3(3, ii) = 0
122  END IF
123  END DO
124  DO ii = 1, cached_dim2
125  DO ij = ii, cached_dim2
126  monores2 = a_mono_exp2(:, ii) + a_mono_exp2(:, ij)
127  a_mono_mult2(ii, ij) = mono_index2(monores2(1), monores2(2)) + 1
128  a_mono_mult2(ij, ii) = a_mono_mult2(ii, ij)
129  END DO
130  END DO
131  DO ii = 1, cached_dim3
132  DO ij = ii, cached_dim3
133  monores3 = a_mono_exp3(:, ii) + a_mono_exp3(:, ij)
134  a_mono_mult3(ii, ij) = mono_index3(monores3(1), monores3(2), monores3(3)) + 1
135  a_mono_mult3(ij, ii) = a_mono_mult3(ii, ij)
136  END DO
137  END DO
138  ii = 1
139  DO i = 1, cached_dim3
140  DO j = 1, 4
141  a_mono_mult3a(j, i) = a_mono_mult3(j, i)
142  ii = ii + 1
143  END DO
144  END DO
145 
146  module_initialized = .true.
147  END SUBROUTINE
148 
149 ! **************************************************************************************************
150 !> \brief size of a polynomial in x up to the given degree
151 !> \param maxgrad ...
152 !> \return ...
153 ! **************************************************************************************************
154  PURE FUNCTION poly_size1(maxgrad) RESULT(res)
155  INTEGER, INTENT(in) :: maxgrad
156  INTEGER :: res
157 
158  res = maxgrad + 1
159  END FUNCTION
160 
161 ! **************************************************************************************************
162 !> \brief size of a polynomial in x,y up to the given degree
163 !> \param maxgrad ...
164 !> \return ...
165 ! **************************************************************************************************
166  PURE FUNCTION poly_size2(maxgrad) RESULT(res)
167  INTEGER, INTENT(in) :: maxgrad
168  INTEGER :: res
169 
170  res = (maxgrad + 1)*(maxgrad + 2)/2
171  END FUNCTION
172 
173 ! **************************************************************************************************
174 !> \brief size of a polynomial in x,y,z up to the given degree
175 !> \param maxgrad ...
176 !> \return ...
177 ! **************************************************************************************************
178  PURE FUNCTION poly_size3(maxgrad) RESULT(res)
179  INTEGER, INTENT(in) :: maxgrad
180  INTEGER :: res
181 
182  res = (maxgrad + 1)*(maxgrad + 2)*(maxgrad + 3)/6
183  END FUNCTION
184 
185 ! **************************************************************************************************
186 !> \brief max grad for a polynom of the given size
187 !> \param n ...
188 !> \return ...
189 ! **************************************************************************************************
190  PURE FUNCTION grad_size1(n) RESULT(res)
191  INTEGER, INTENT(in) :: n
192  INTEGER :: res
193 
194  res = n - 1
195  END FUNCTION
196 
197 ! **************************************************************************************************
198 !> \brief max grad for a polynom of the given size
199 !> \param n ...
200 !> \return ...
201 ! **************************************************************************************************
202  PURE FUNCTION grad_size2(n) RESULT(res)
203  INTEGER, INTENT(in) :: n
204  INTEGER :: res
205 
206  res = int(floor(0.5_dp*(sqrt(1.0_dp + 8.0_dp*real(n, dp)) - 1.0_dp) - 2.e-6_dp))
207  END FUNCTION
208 
209 ! **************************************************************************************************
210 !> \brief max grad for a polynom of the given size
211 !> \param n ...
212 !> \return ...
213 ! **************************************************************************************************
214  PURE FUNCTION grad_size3(n) RESULT(res)
215  INTEGER, INTENT(in) :: n
216  INTEGER :: res
217 
218  INTEGER :: nn
219  REAL(dp) :: g1
220 
221  IF (n < 1) THEN
222  res = -1
223  ELSE
224  nn = n*6
225  g1 = (108.0_dp*nn + 12.0_dp*sqrt(81.0_dp*nn*nn - 12.0_dp))**(1.0_dp/3.0_dp)
226  res = floor(g1/6.0_dp + 2.0_dp/g1 - 1.0_dp - 2.e-6_dp)
227  END IF
228  END FUNCTION
229 
230 ! **************************************************************************************************
231 !> \brief 0-based index of monomial of the given degree
232 !> \param i ...
233 !> \return ...
234 ! **************************************************************************************************
235  PURE FUNCTION mono_index1(i) RESULT(res)
236  INTEGER, INTENT(in) :: i
237  INTEGER :: res
238 
239  res = i
240  END FUNCTION
241 
242 ! **************************************************************************************************
243 !> \brief 0-based index of monomial of the given degree
244 !> \param i ...
245 !> \param j ...
246 !> \return ...
247 ! **************************************************************************************************
248  PURE FUNCTION mono_index2(i, j) RESULT(res)
249  INTEGER, INTENT(in) :: i, j
250  INTEGER :: res
251 
252  INTEGER :: grad
253 
254  grad = i + j
255  res = grad*(grad + 1)/2 + j
256  END FUNCTION
257 
258 ! **************************************************************************************************
259 !> \brief 0-based index of monomial of the given degree
260 !> \param i ...
261 !> \param j ...
262 !> \param k ...
263 !> \return ...
264 ! **************************************************************************************************
265  PURE FUNCTION mono_index3(i, j, k) RESULT(res)
266  INTEGER, INTENT(in) :: i, j, k
267  INTEGER :: res
268 
269  INTEGER :: grad, sgrad
270 
271  sgrad = j + k
272  grad = i + sgrad
273  res = grad*(grad + 1)*(grad + 2)/6 + (sgrad)*(sgrad + 1)/2 + k
274  END FUNCTION
275 
276 ! **************************************************************************************************
277 !> \brief exponents of the monomial at the given 0-based index
278 !> \param ii ...
279 !> \return ...
280 ! **************************************************************************************************
281  PURE FUNCTION mono_exp1(ii) RESULT(res)
282  INTEGER, INTENT(in) :: ii
283  INTEGER :: res
284 
285  res = ii
286  END FUNCTION
287 
288 ! **************************************************************************************************
289 !> \brief exponents of the monomial at the given 0-based index
290 !> \param ii ...
291 !> \return ...
292 ! **************************************************************************************************
293  PURE FUNCTION mono_exp2(ii) RESULT(res)
294  INTEGER, INTENT(in) :: ii
295  INTEGER, DIMENSION(2) :: res
296 
297  INTEGER :: grad
298 
299  grad = int(floor(0.5_dp*(sqrt(9.0_dp + 8.0_dp*ii) - 1.0_dp) - 2.e-6_dp))
300  res(2) = ii - (grad)*(grad + 1)/2
301  res(1) = grad - res(2)
302  END FUNCTION
303 
304 ! **************************************************************************************************
305 !> \brief exponents of the monomial at the given 0-based index
306 !> \param n ...
307 !> \return ...
308 ! **************************************************************************************************
309  PURE FUNCTION mono_exp3(n) RESULT(res)
310  INTEGER, INTENT(in) :: n
311  INTEGER, DIMENSION(3) :: res
312 
313  INTEGER :: grad, grad1, ii, nn
314  REAL(dp) :: g1
315 
316  nn = (n + 1)*6
317  g1 = (108.0_dp*nn + 12.0_dp*sqrt(81.0_dp*nn*nn - 12.0_dp))**(1.0_dp/3.0_dp)
318  grad1 = int(floor(g1/6.0_dp + 2.0_dp/g1 - 1.0_dp - 2.e-6_dp))
319  ii = n - grad1*(grad1 + 1)*(grad1 + 2)/6
320  grad = int(floor(0.5_dp*(sqrt(9.0_dp + 8.0_dp*ii) - 1.0_dp) - 1.e-6_dp))
321  res(3) = ii - grad*(grad + 1)/2
322  res(2) = grad - res(3)
323  res(1) = grad1 - grad
324  END FUNCTION
325 
326 ! **************************************************************************************************
327 !> \brief the index of the result of the multiplication of the two monomials
328 !> \param ii ...
329 !> \param ij ...
330 !> \return ...
331 ! **************************************************************************************************
332  PURE FUNCTION mono_mult1(ii, ij) RESULT(res)
333  INTEGER, INTENT(in) :: ii, ij
334  INTEGER :: res
335 
336  res = ii + ij
337  END FUNCTION
338 
339 ! **************************************************************************************************
340 !> \brief the index of the result of the multiplication of the two monomials
341 !> \param ii ...
342 !> \param ij ...
343 !> \return ...
344 ! **************************************************************************************************
345  PURE FUNCTION mono_mult2(ii, ij) RESULT(res)
346  INTEGER, INTENT(in) :: ii, ij
347  INTEGER :: res
348 
349  INTEGER, DIMENSION(2) :: monores
350 
351  monores = mono_exp2(ii) + mono_exp2(ij)
352  res = mono_index2(monores(1), monores(2))
353  END FUNCTION
354 
355 ! **************************************************************************************************
356 !> \brief the index of the result of the multiplication of the two monomials
357 !> \param ii ...
358 !> \param ij ...
359 !> \return ...
360 ! **************************************************************************************************
361  PURE FUNCTION mono_mult3(ii, ij) RESULT(res)
362  INTEGER, INTENT(in) :: ii, ij
363  INTEGER :: res
364 
365  INTEGER, DIMENSION(3) :: monores
366 
367  monores = mono_exp3(ii) + mono_exp3(ij)
368  res = mono_index3(monores(1), monores(2), monores(3))
369  END FUNCTION
370 
371 ! **************************************************************************************************
372 !> \brief multiplies the polynomials p1 with p2 using pRes to store the result
373 !> \param p1 ...
374 !> \param p2 ...
375 !> \param pRes ...
376 !> \param np1 ...
377 !> \param sumUp ...
378 ! **************************************************************************************************
379  SUBROUTINE poly_mult1(p1, p2, pRes, np1, sumUp)
380  REAL(dp), DIMENSION(:), INTENT(in) :: p1, p2
381  REAL(dp), DIMENSION(:), INTENT(inout) :: pres
382  INTEGER, INTENT(in), OPTIONAL :: np1
383  LOGICAL, INTENT(in), OPTIONAL :: sumup
384 
385  INTEGER :: i, ipoly, ipos, j, mynp1, newgrad, &
386  newsize, respos, resshift_0, size_p1, &
387  size_p2
388  LOGICAL :: mysumup
389 
390  IF (.NOT. module_initialized) cpabort("module d3_poly not initialized")
391  mysumup = .false.
392  mynp1 = 1
393  IF (PRESENT(np1)) mynp1 = np1
394  IF (PRESENT(sumup)) mysumup = sumup
395  size_p1 = SIZE(p1)/mynp1
396  size_p2 = SIZE(p2)
397  newgrad = grad_size1(size_p1) + grad_size1(size_p2)
398  newsize = SIZE(pres)/mynp1
399  cpassert(newsize >= poly_size1(newgrad))
400  IF (.NOT. mysumup) pres = 0
401  ipos = 1
402  resshift_0 = 0
403  DO ipoly = 0, mynp1 - 1
404  DO i = 1, size_p1
405  respos = resshift_0 + i
406  DO j = 1, size_p2
407  pres(respos) = pres(respos) + p1(ipos)*p2(j)
408  respos = respos + 1
409  END DO
410  ipos = ipos + 1
411  END DO
412  resshift_0 = resshift_0 + newsize
413  END DO
414  END SUBROUTINE
415 
416 ! **************************************************************************************************
417 !> \brief multiplies p1 with p2 using pRes to store the result
418 !> \param p1 ...
419 !> \param p2 ...
420 !> \param pRes ...
421 !> \param np1 ...
422 !> \param sumUp ...
423 ! **************************************************************************************************
424  SUBROUTINE poly_mult2(p1, p2, pRes, np1, sumUp)
425  REAL(dp), DIMENSION(:), INTENT(in) :: p1, p2
426  REAL(dp), DIMENSION(:), INTENT(inout) :: pres
427  INTEGER, INTENT(in), OPTIONAL :: np1
428  LOGICAL, INTENT(in), OPTIONAL :: sumup
429 
430  INTEGER :: g1, g1g2, g2, grad1, grad2, i, ipoly, ishift, j, msize_p1, mynp1, newgrad, &
431  newsize, shift1, shift2, shiftres, shiftres_0, size_p1, size_p2, subg2
432  LOGICAL :: mysumup
433 
434  IF (.NOT. module_initialized) cpabort("module d3_poly not initialized")
435  mysumup = .false.
436  IF (PRESENT(sumup)) mysumup = sumup
437  mynp1 = 1
438  IF (PRESENT(np1)) mynp1 = np1
439  size_p1 = SIZE(p1)/mynp1
440  size_p2 = SIZE(p2)
441  grad1 = grad_size2(size_p1)
442  grad2 = grad_size2(size_p2)
443  newgrad = grad1 + grad2
444  newsize = SIZE(pres)/mynp1
445  cpassert(newsize >= poly_size2(newgrad))
446  IF (.NOT. mysumup) pres = 0
447  ishift = 0
448  shiftres = 0
449  DO ipoly = 1, mynp1
450  DO i = 1, min(size_p1, cached_dim2)
451  DO j = 1, min(size_p2, cached_dim2)
452  pres(shiftres + a_mono_mult2(j, i)) = pres(shiftres + a_mono_mult2(j, i)) &
453  + p1(ishift + i)*p2(j)
454  END DO
455  END DO
456  ishift = ishift + size_p1
457  shiftres = shiftres + newsize
458  END DO
459  IF (grad1 > max_grad2 .OR. grad2 > max_grad2) THEN
460  msize_p1 = size_p1
461  shiftres_0 = 0
462  DO ipoly = 0, mynp1 - 1
463  shift1 = ipoly*size_p1
464  DO g1 = 0, grad1
465  ! shift1=g1*(g1+1)/2
466  IF (g1 > max_grad2) THEN
467  subg2 = 0
468  shift2 = 0
469  g1g2 = shiftres_0 - 1
470  ELSE
471  subg2 = max_grad2 + 1
472  shift2 = cached_dim2
473  g1g2 = shiftres_0 + g1*subg2 - 1
474  END IF
475  DO g2 = subg2, grad2
476  ! shift2=g2*(g2+1)/2
477  shiftres = shift1 + shift2 + g1g2 ! shiftRes=(g1+g2)*(g1+g2+1)/2-1+ipoly*newSize
478  DO i = 1, min(g1 + 1, msize_p1 - shift1)
479  DO j = 1, min(g2 + 1, size_p2 - shift2)
480  pres(shiftres + i + j) = pres(shiftres + i + j) + p1(shift1 + i)*p2(shift2 + j)
481  END DO
482  END DO
483  shift2 = shift2 + g2 + 1 !
484  g1g2 = g1g2 + g1 !
485  END DO
486  shift1 = shift1 + g1 + 1 !
487  END DO
488  shiftres_0 = shiftres_0 + newsize - size_p1
489  msize_p1 = msize_p1 + size_p1
490  END DO
491  END IF
492  END SUBROUTINE
493 
494 ! **************************************************************************************************
495 !> \brief multiplies p1 with p2 using pRes to store the result
496 !> \param p1 ...
497 !> \param p2 ...
498 !> \param pRes ...
499 !> \param np1 ...
500 !> \param sumUp ...
501 ! **************************************************************************************************
502  SUBROUTINE poly_mult3(p1, p2, pRes, np1, sumUp)
503  REAL(dp), DIMENSION(:), INTENT(in) :: p1, p2
504  REAL(dp), DIMENSION(:), INTENT(inout) :: pres
505  INTEGER, INTENT(in), OPTIONAL :: np1
506  LOGICAL, INTENT(in), OPTIONAL :: sumup
507 
508  INTEGER :: grad1, grad2, mynp1, newgrad, newsize, &
509  size_p1, size_p2
510  LOGICAL :: mysumup
511 
512  IF (.NOT. module_initialized) cpabort("module d3_poly not initialized")
513  mysumup = .false.
514  IF (PRESENT(sumup)) mysumup = sumup
515  mynp1 = 1
516  IF (PRESENT(np1)) mynp1 = np1
517  size_p1 = SIZE(p1)/mynp1
518  size_p2 = SIZE(p2)
519  grad1 = grad_size3(size_p1)
520  grad2 = grad_size3(size_p2)
521  newgrad = grad1 + grad2
522  newsize = SIZE(pres)/mynp1
523  cpassert(newsize >= poly_size3(newgrad))
524  CALL poly_mult3b(p1, SIZE(p1), grad1, p2, SIZE(p2), grad2, pres, SIZE(pres), mynp1, mysumup)
525  END SUBROUTINE
526 
527 ! **************************************************************************************************
528 !> \brief low level routine of poly_mult3 without checks
529 !> \param p1 ...
530 !> \param size_p1 ...
531 !> \param grad1 ...
532 !> \param p2 ...
533 !> \param size_p2 ...
534 !> \param grad2 ...
535 !> \param pRes ...
536 !> \param size_pRes ...
537 !> \param np1 ...
538 !> \param sumUp ...
539 ! **************************************************************************************************
540  SUBROUTINE poly_mult3b(p1, size_p1, grad1, p2, size_p2, grad2, pRes, size_pRes, np1, sumUp)
541  INTEGER, INTENT(in) :: size_p1
542  REAL(dp), DIMENSION(IF_CHECK(size_p1, *)), &
543  INTENT(in) :: p1
544  INTEGER, INTENT(in) :: grad1, size_p2
545  REAL(dp), DIMENSION(IF_CHECK(size_p2, *)), &
546  INTENT(in) :: p2
547  INTEGER, INTENT(in) :: grad2, size_pres
548  REAL(dp), DIMENSION(IF_CHECK(size_pRes, *)), &
549  INTENT(inout) :: pres
550  INTEGER, INTENT(in) :: np1
551  LOGICAL, INTENT(in) :: sumup
552 
553  INTEGER :: g1, g2, i, i1, i2, ipoly, ishift, j, j1, j2, msize_p1, my_size_p1, newsize, &
554  shift1, shift1i, shift1j, shift2, shift2i, shift2j, shiftres, shiftres_0, shiftresi, &
555  shiftresi_0, shiftresj, subg2, subgrad
556 
557  IF (.NOT. module_initialized) cpabort("module d3_poly not initialized")
558  my_size_p1 = size_p1/np1
559  newsize = size_pres/np1
560 
561  IF (.NOT. sumup) pres(1:size_pres) = 0.0_dp
562  ishift = 0
563  shiftres = 0
564  DO ipoly = 1, np1
565  DO i = 1, min(my_size_p1, cached_dim3)
566  DO j = 1, min(size_p2, cached_dim3)
567  pres(shiftres + a_mono_mult3(j, i)) = pres(shiftres + a_mono_mult3(j, i)) &
568  + p1(ishift + i)*p2(j)
569  END DO
570  END DO
571  ishift = ishift + my_size_p1
572  shiftres = shiftres + newsize
573  END DO
574  IF (grad1 > max_grad3 .OR. grad2 > max_grad3) THEN
575  ! one could remove multiplications even more aggressively...
576  msize_p1 = my_size_p1
577  DO ipoly = 0, np1 - 1
578  shift1 = 1 + ipoly*my_size_p1
579  shiftres_0 = 1 + ipoly*newsize
580  DO g1 = 0, grad1
581  IF (g1 > max_grad3) THEN
582  subg2 = 0
583  shift2 = 1
584  ELSE
585  subg2 = max_grad3 + 1
586  shift2 = subg2*(subg2 + 1)*(subg2 + 2)/6 + 1
587  END IF
588  DO g2 = subg2, grad2
589  shiftres = (g1 + g2)*(g1 + g2 + 1)*(g1 + g2 + 2)/6 + shiftres_0
590  shift1i = shift1
591  shiftresi_0 = shiftres
592  DO i1 = g1, 0, -1
593  IF (shift1i > msize_p1) EXIT
594  shift2i = shift2
595  shiftresi = shiftresi_0
596  subgrad = g1 - i1
597  DO i2 = g2, 0, -1
598  !subGrad=g1+g2-i1-i2
599  !shiftResI=shiftRes+(subGrad)*(subGrad+1)/2
600  !shift2I=shift2+(g2-i2)*(g2-i2+1)/2
601  IF (shift2i > size_p2) EXIT
602  DO j1 = g1 - i1, 0, -1
603  shift1j = shift1i + g1 - i1 - j1
604  IF (shift1j > msize_p1) EXIT
605  DO j2 = g2 - i2, 0, -1
606  shift2j = shift2i + g2 - i2 - j2
607  IF (shift2j > size_p2) EXIT
608  shiftresj = shiftresi + (subgrad - j1 - j2)
609  ! shift1J=mono_index3(i1,j1,g1-i1-j1)+ipoly*my_size_p1+1
610  ! shift2J=mono_index3(i2,j2,g2-i2-j2)+1
611  ! shiftResJ=mono_index3(i1+i2,j1+j2,g1+g2-i1-i2-j1-j2)+ipoly*newSize+1
612  pres(shiftresj) = pres(shiftresj) + p1(shift1j)*p2(shift2j)
613  END DO
614  END DO
615  subgrad = subgrad + 1
616  shift2i = shift2i + (g2 - i2 + 1)
617  shiftresi = shiftresi + subgrad
618  END DO
619  shift1i = shift1i + (g1 - i1 + 1)
620  shiftresi_0 = shiftresi_0 + (g1 - i1 + 1)
621  END DO
622  shift2 = shift2 + (g2 + 1)*(g2 + 2)/2
623  END DO
624  shift1 = shift1 + (g1 + 1)*(g1 + 2)/2
625  END DO
626  msize_p1 = msize_p1 + my_size_p1
627  END DO
628  END IF
629  END SUBROUTINE
630 
631 ! **************************************************************************************************
632 !> \brief low level routine that multiplies with a polynomial of grad 1
633 !> \param p1 ...
634 !> \param size_p1 ...
635 !> \param grad1 ...
636 !> \param p2 ...
637 !> \param pRes ...
638 !> \param size_pRes ...
639 !> \param np1 ...
640 !> \param sumUp ...
641 ! **************************************************************************************************
642  SUBROUTINE poly_mult3ab(p1, size_p1, grad1, p2, pRes, size_pRes, np1, sumUp)
643  INTEGER, INTENT(in) :: size_p1
644  REAL(dp), DIMENSION(IF_CHECK(size_p1, *)), &
645  INTENT(in) :: p1
646  INTEGER, INTENT(in) :: grad1
647  REAL(dp), DIMENSION(IF_CHECK(4, *)), INTENT(in) :: p2
648  INTEGER, INTENT(in) :: size_pres
649  REAL(dp), DIMENSION(IF_CHECK(size_pRes, *)), &
650  INTENT(inout) :: pres
651  INTEGER, INTENT(in) :: np1
652  LOGICAL, INTENT(in) :: sumup
653 
654  INTEGER, PARAMETER :: grad2 = 1
655 
656  INTEGER :: g1, g2, i, i1, i2, ipoly, ishift, j1, j2, msize_p1, my_size_p1, newsize, shift1, &
657  shift1i, shift1j, shift2, shift2i, shift2j, shiftres, shiftres_0, shiftresi, shiftresi_0, &
658  shiftresj, subg2, subgrad
659 
660  IF (.NOT. module_initialized) cpabort("module d3_poly not initialized")
661  my_size_p1 = size_p1/np1
662  newsize = size_pres/np1
663 
664  IF (.NOT. sumup) pres(1:size_pres) = 0.0_dp
665  ishift = 0
666  shiftres = 0
667  DO ipoly = 1, np1
668  DO i = 1, min(my_size_p1, cached_dim3)
669  pres(shiftres + a_mono_mult3a(1, i)) = pres(shiftres + a_mono_mult3a(1, i)) &
670  + p1(ishift + i)*p2(1)
671  pres(shiftres + a_mono_mult3a(2, i)) = pres(shiftres + a_mono_mult3a(2, i)) &
672  + p1(ishift + i)*p2(2)
673  pres(shiftres + a_mono_mult3a(3, i)) = pres(shiftres + a_mono_mult3a(3, i)) &
674  + p1(ishift + i)*p2(3)
675  pres(shiftres + a_mono_mult3a(4, i)) = pres(shiftres + a_mono_mult3a(4, i)) &
676  + p1(ishift + i)*p2(4)
677  END DO
678  ishift = ishift + my_size_p1
679  shiftres = shiftres + newsize
680  END DO
681  IF (grad1 > max_grad3 .OR. grad2 > max_grad3) THEN
682  ! one could remove multiplications even more aggressively...
683  msize_p1 = my_size_p1
684  DO ipoly = 0, np1 - 1
685  shift1 = 1 + ipoly*my_size_p1 + (max_grad3 + 1)*(max_grad3 + 2)*(max_grad3 + 3)/6
686  shiftres_0 = 1 + ipoly*newsize
687  DO g1 = max_grad3 + 1, grad1
688  subg2 = 0
689  shift2 = 1
690  DO g2 = subg2, grad2
691  shiftres = (g1 + g2)*(g1 + g2 + 1)*(g1 + g2 + 2)/6 + shiftres_0
692  shift1i = shift1
693  shiftresi_0 = shiftres
694  DO i1 = g1, 0, -1
695  IF (shift1i > msize_p1) EXIT
696  shift2i = shift2
697  shiftresi = shiftresi_0
698  subgrad = g1 - i1
699  DO i2 = g2, 0, -1
700  !subGrad=g1+g2-i1-i2
701  !shiftResI=shiftRes+(subGrad)*(subGrad+1)/2
702  !shift2I=shift2+(g2-i2)*(g2-i2+1)/2
703  DO j1 = g1 - i1, 0, -1
704  shift1j = shift1i + g1 - i1 - j1
705  IF (shift1j > msize_p1) EXIT
706  DO j2 = g2 - i2, 0, -1
707  shift2j = shift2i + g2 - i2 - j2
708  shiftresj = shiftresi + (subgrad - j1 - j2)
709  ! shift1J=mono_index3(i1,j1,g1-i1-j1)+ipoly*my_size_p1+1
710  ! shift2J=mono_index3(i2,j2,g2-i2-j2)+1
711  ! shiftResJ=mono_index3(i1+i2,j1+j2,g1+g2-i1-i2-j1-j2)+ipoly*newSize+1
712  pres(shiftresj) = pres(shiftresj) + p1(shift1j)*p2(shift2j)
713  END DO
714  END DO
715  subgrad = subgrad + 1
716  shift2i = shift2i + (g2 - i2 + 1)
717  shiftresi = shiftresi + subgrad
718  END DO
719  shift1i = shift1i + (g1 - i1 + 1)
720  shiftresi_0 = shiftresi_0 + (g1 - i1 + 1)
721  END DO
722  shift2 = shift2 + (g2 + 1)*(g2 + 2)/2
723  END DO
724  shift1 = shift1 + (g1 + 1)*(g1 + 2)/2
725  END DO
726  msize_p1 = msize_p1 + my_size_p1
727  END DO
728  END IF
729  END SUBROUTINE
730 
731 ! **************************************************************************************************
732 !> \brief writes out a 1d polynomial in a human readable form
733 !> \param p ...
734 !> \param out_f ...
735 ! **************************************************************************************************
736  SUBROUTINE poly_write1(p, out_f)
737  REAL(dp), DIMENSION(:), INTENT(in) :: p
738  INTEGER, INTENT(in) :: out_f
739 
740  INTEGER :: i
741  LOGICAL :: did_write
742 
743  did_write = .false.
744  DO i = 1, SIZE(p)
745  IF (p(i) /= 0) THEN
746  IF (p(i) >= 0) WRITE (out_f, '("+")', advance='NO')
747  WRITE (out_f, '(G20.10)', advance='NO') p(i)
748  IF (i /= 1) WRITE (out_f, '("*x^",I3)', advance='NO') i - 1
749  did_write = .true.
750  END IF
751  END DO
752  IF (.NOT. did_write) WRITE (out_f, '("0.0")', advance='NO')
753  END SUBROUTINE
754 
755 ! **************************************************************************************************
756 !> \brief writes out a 2d polynomial in a human readable form
757 !> \param p ...
758 !> \param out_f ...
759 ! **************************************************************************************************
760  SUBROUTINE poly_write2(p, out_f)
761  REAL(dp), DIMENSION(:), INTENT(in) :: p
762  INTEGER, INTENT(in) :: out_f
763 
764  INTEGER :: i
765  INTEGER, DIMENSION(2) :: mono_e
766  LOGICAL :: did_write
767 
768  IF (.NOT. module_initialized) cpabort("module d3_poly not initialized")
769  did_write = .false.
770  DO i = 1, SIZE(p)
771  mono_e = mono_exp2(i - 1)
772  IF (p(i) /= 0) THEN
773  IF (p(i) >= 0) WRITE (out_f, '("+")', advance='NO')
774  WRITE (out_f, '(G20.10)', advance='NO') p(i)
775  IF (mono_e(1) /= 0) WRITE (out_f, '("*x^",I3)', advance='NO') mono_e(1)
776  IF (mono_e(2) /= 0) WRITE (out_f, '("*y^",I3)', advance='NO') mono_e(2)
777  did_write = .true.
778  END IF
779  END DO
780  IF (.NOT. did_write) WRITE (out_f, '("0.0")', advance='NO')
781  END SUBROUTINE
782 
783 ! **************************************************************************************************
784 !> \brief writes out the polynomial in a human readable form
785 !> \param p ...
786 !> \param out_f ...
787 ! **************************************************************************************************
788  SUBROUTINE poly_write3(p, out_f)
789  REAL(dp), DIMENSION(:), INTENT(in) :: p
790  INTEGER, INTENT(in) :: out_f
791 
792  INTEGER :: i
793  INTEGER, DIMENSION(3) :: mono_e
794  LOGICAL :: did_write
795 
796  IF (.NOT. module_initialized) cpabort("module d3_poly not initialized")
797  did_write = .false.
798  DO i = 1, SIZE(p)
799  mono_e = mono_exp3(i - 1)
800  IF (p(i) /= 0) THEN
801  IF (p(i) >= 0) WRITE (out_f, '("+")', advance='NO')
802  WRITE (out_f, '(G20.10)', advance='NO') p(i)
803  IF (mono_e(1) /= 0) WRITE (out_f, '("*x^",I3)', advance='NO') mono_e(1)
804  IF (mono_e(2) /= 0) WRITE (out_f, '("*y^",I3)', advance='NO') mono_e(2)
805  IF (mono_e(3) /= 0) WRITE (out_f, '("*z^",I3)', advance='NO') mono_e(3)
806  did_write = .true.
807  END IF
808  END DO
809  IF (.NOT. did_write) WRITE (out_f, '("0.0")', advance='NO')
810  END SUBROUTINE
811 
812 ! **************************************************************************************************
813 !> \brief random poly with coeffiecents that are easy to print exactly,
814 !> of the given maximum size (for testing purposes)
815 !> \param p ...
816 !> \param maxSize ...
817 !> \param minSize ...
818 !> \return ...
819 ! **************************************************************************************************
820  FUNCTION poly_random(p, maxSize, minSize) RESULT(res)
821  REAL(dp), DIMENSION(:), INTENT(out) :: p
822  INTEGER, INTENT(in) :: maxsize
823  INTEGER, INTENT(in), OPTIONAL :: minsize
824  INTEGER :: res
825 
826  INTEGER :: i, myminsize, psize
827  REAL(dp) :: g
828 
829  IF (.NOT. module_initialized) cpabort("module d3_poly not initialized")
830  myminsize = 1
831  IF (PRESENT(minsize)) myminsize = minsize
832  CALL random_number(g)
833  psize = min(maxsize, myminsize + int((maxsize - myminsize + 1)*g))
834  cpassert(SIZE(p) >= psize)
835  CALL random_number(p)
836  DO i = 1, psize
837  p(i) = real(int(p(i)*200.0_dp - 100.0_dp), dp)/100.0_dp
838  END DO
839  DO i = psize + 1, SIZE(p)
840  p(i) = 0.0_dp
841  END DO
842  res = psize
843  END FUNCTION
844 
845 ! **************************************************************************************************
846 !> \brief returns in the polynomials pRes the transpose of the
847 !> affine transformation x -> m*x+b of p
848 !> \param p ...
849 !> \param m ...
850 !> \param b ...
851 !> \param pRes ...
852 !> \param npoly ...
853 ! **************************************************************************************************
854  SUBROUTINE poly_affine_t3t(p, m, b, pRes, npoly)
855  REAL(dp), DIMENSION(:), INTENT(in) :: p
856  REAL(dp), DIMENSION(3, 3), INTENT(in) :: m
857  REAL(dp), DIMENSION(3), INTENT(in) :: b
858  REAL(dp), DIMENSION(:), INTENT(out) :: pres
859  INTEGER, INTENT(in), OPTIONAL :: npoly
860 
861  INTEGER :: grad, i, igrad, ii, ii1, ipoly, j, k, minressize, monodim1, monodim2, monodimatt, &
862  monodimatt2, monofulldim1, monofulldim2, monosize1, monosize2, my_npoly, pcoeff, pidx, &
863  pshift, rescoeff, resshift, rest_size_p, size_p, size_res, start_idx1
864  REAL(dp), ALLOCATABLE, DIMENSION(:) :: monog1, monog2
865  REAL(dp), DIMENSION(4, 3) :: basepoly
866 
867  IF (.NOT. module_initialized) cpabort("module d3_poly not initialized")
868  my_npoly = 1
869  IF (PRESENT(npoly)) my_npoly = npoly
870  basepoly(1, :) = b
871  DO j = 1, 3
872  DO i = 1, 3
873  basepoly(j + 1, i) = m(i, j)
874  END DO
875  END DO
876  size_p = SIZE(pres)/my_npoly
877  size_res = SIZE(p)/my_npoly
878  grad = grad_size3(size_res)
879  minressize = poly_size3(grad)
880  cpassert(size_res == minressize)
881  cpassert(size_p >= minressize)
882  pres = 0
883  IF (size_p == 0) RETURN
884  ii1 = 1
885  ii = 1
886  DO ipoly = 1, my_npoly
887  pres(ii1) = p(ii)
888  ii = ii + size_res
889  ii1 = ii1 + size_p
890  END DO
891  IF (size_p == 1) RETURN
892 
893  ALLOCATE (monog1((grad + 1)*(grad + 2)/2*minressize), &
894  monog2((grad + 1)*(grad + 2)/2*minressize))
895  !monoG1=0
896  !monoG2=0
897  ii1 = 1
898  DO j = 1, 3
899  DO i = 1, 4
900  monog1(ii1) = basepoly(i, j)
901  ii1 = ii1 + 1
902  END DO
903  END DO
904  ii1 = 2
905  igrad = 1
906  monodim1 = 4
907  monosize1 = 3
908  monofulldim1 = monodim1*monosize1
909  rest_size_p = size_p - 1
910  DO
911  k = min(rest_size_p, monosize1)
912  !call dgemm('T','N',monoDim1,my_npoly,k,&
913  ! 1.0_dp,monoG1,monoDim1,p(ii1:),size_p,1.0_dp,pRes,size_res)
914  resshift = 0
915  pshift = ii1
916  DO ipoly = 1, my_npoly
917  pidx = pshift
918  ii = 1
919  DO pcoeff = 1, k
920  DO rescoeff = 1, monodim1
921  pres(pidx) = pres(pidx) + p(resshift + rescoeff)*monog1(ii)
922  ii = ii + 1
923  END DO
924  pidx = pidx + 1
925  END DO
926  resshift = resshift + size_res
927  pshift = pshift + size_p
928  END DO
929 
930  rest_size_p = rest_size_p - k
931  ii1 = ii1 + k
932  IF (rest_size_p <= 0) EXIT
933 
934  monosize2 = igrad + 2 + monosize1
935  monodim2 = monodim1 + monosize2
936  monofulldim2 = monosize2*monodim2
937  monodimatt = monosize1*monodim2
938  CALL poly_mult3ab(if_check(monog1(1:monofulldim1), monog1(1)), monofulldim1, igrad, &
939  if_check(basepoly(:, 1), basepoly(1, 1)), &
940  if_check(monog2(1:monodimatt), monog2(1)), monodimatt, monosize1, .false.)
941  monodimatt2 = monofulldim2 - monodim2
942  start_idx1 = (monosize1 - igrad - 1)*monodim1
943  CALL poly_mult3ab(if_check(monog1(start_idx1 + 1:monofulldim1), monog1(start_idx1 + 1)), &
944  monofulldim1 - start_idx1, igrad, if_check(basepoly(:, 2), basepoly(1, 2)), &
945  if_check(monog2(monodimatt + 1:monodimatt2), monog2(monodimatt + 1)), &
946  monodimatt2 - monodimatt, igrad + 1, .false.)
947  CALL poly_mult3ab(if_check(monog1(monofulldim1 - monodim1 + 1:monofulldim1), monog1(monofulldim1 - monodim1 + 1)), &
948  monodim1, igrad, if_check(basepoly(:, 3), basepoly(1, 3)), &
949  if_check(monog2(monodimatt2 + 1:monofulldim2), monog2(monodimatt2 + 1)), &
950  monofulldim2 - monodimatt2, 1, .false.)
951  igrad = igrad + 1
952 
953  ! even grads
954 
955  k = min(rest_size_p, monosize2)
956  !call dgemm('T','N',monoDim2,my_npoly,k,&
957  ! 1.0_dp,monoG2,monoDim2,p(ii1:),size_p,1.0_dp,pRes,size_res)
958  resshift = 0
959  pshift = ii1
960  DO ipoly = 1, my_npoly
961  pidx = pshift
962  ii = 1
963  DO pcoeff = 1, k
964  DO rescoeff = 1, monodim2
965  pres(pidx) = pres(pidx) + p(resshift + rescoeff)*monog2(ii)
966  ii = ii + 1
967  END DO
968  pidx = pidx + 1
969  END DO
970  resshift = resshift + size_res
971  pshift = pshift + size_p
972  END DO
973 
974  rest_size_p = rest_size_p - k
975  ii1 = ii1 + k
976  IF (rest_size_p <= 0) EXIT
977 
978  monosize1 = igrad + 2 + monosize2
979  monodim1 = monodim2 + monosize1
980  monofulldim1 = monosize1*monodim1
981  monodimatt = monosize2*monodim1
982  CALL poly_mult3ab(if_check(monog2(1:monofulldim2), monog2(1)), monofulldim2, igrad, &
983  if_check(basepoly(:, 1), basepoly(1, 1)), if_check(monog1(1:monodimatt), monog1(1)), &
984  monodimatt, monosize2, .false.)
985  monodimatt2 = monofulldim1 - monodim1
986  start_idx1 = (monosize2 - igrad - 1)*monodim2
987  CALL poly_mult3ab(if_check(monog2(start_idx1 + 1:monofulldim2), monog2(start_idx1 + 1)), &
988  monofulldim2 - start_idx1, igrad, if_check(basepoly(:, 2), basepoly(1, 2)), &
989  if_check(monog1(monodimatt + 1:monodimatt2), monog1(monodimatt + 1)), monodimatt2 - monodimatt, &
990  igrad + 1, .false.)
991  CALL poly_mult3ab(if_check(monog2(monofulldim2 - monodim2 + 1:monofulldim2), monog2(monofulldim2 - monodim2 + 1)), &
992  monodim2, igrad, if_check(basepoly(:, 3), basepoly(1, 3)), &
993  if_check(monog1(monodimatt2 + 1:monofulldim1), monog1(monodimatt2 + 1)), &
994  monofulldim1 - monodimatt2, 1, .false.)
995  igrad = igrad + 1
996 
997  ! ! alternative to unrolling
998  ! monoG1=monoG2
999  ! monoSize1=monoSize2
1000  ! monoDim1=monoDim2
1001  ! monoFullDim1=monoFullDim2
1002  END DO
1003  DEALLOCATE (monog1, monog2)
1004  END SUBROUTINE
1005 
1006 ! **************************************************************************************************
1007 !> \brief returns in the polynomials pRes the affine transformation x -> m*x+b of p
1008 !> \param p ...
1009 !> \param m ...
1010 !> \param b ...
1011 !> \param pRes ...
1012 !> \param npoly ...
1013 ! **************************************************************************************************
1014  SUBROUTINE poly_affine_t3(p, m, b, pRes, npoly)
1015  REAL(dp), DIMENSION(:), INTENT(in) :: p
1016  REAL(dp), DIMENSION(3, 3), INTENT(in) :: m
1017  REAL(dp), DIMENSION(3), INTENT(in) :: b
1018  REAL(dp), DIMENSION(:), INTENT(out) :: pres
1019  INTEGER, INTENT(in), OPTIONAL :: npoly
1020 
1021  INTEGER :: grad, i, igrad, ii, ii1, ipoly, j, k, minressize, monodim1, monodim2, monodimatt, &
1022  monodimatt2, monofulldim1, monofulldim2, monosize1, monosize2, my_npoly, pcoeff, pidx, &
1023  pshift, rescoeff, resshift, rest_size_p, size_p, size_res, start_idx1
1024  REAL(dp), ALLOCATABLE, DIMENSION(:) :: monog1, monog2
1025  REAL(dp), DIMENSION(4, 3) :: basepoly
1026 
1027  IF (.NOT. module_initialized) cpabort("module d3_poly not initialized")
1028  my_npoly = 1
1029  IF (PRESENT(npoly)) my_npoly = npoly
1030  basepoly(1, :) = b
1031  DO j = 1, 3
1032  DO i = 1, 3
1033  basepoly(j + 1, i) = m(i, j)
1034  END DO
1035  END DO
1036  size_p = SIZE(p)/my_npoly
1037  grad = grad_size3(size_p)
1038  size_res = SIZE(pres)/my_npoly
1039  minressize = poly_size3(grad)
1040  cpassert(size_res >= minressize)
1041  pres = 0
1042  IF (size_p == 0) RETURN
1043  ii1 = 1
1044  ii = 1
1045  DO ipoly = 1, my_npoly
1046  pres(ii) = p(ii1)
1047  ii = ii + size_res
1048  ii1 = ii1 + size_p
1049  END DO
1050  IF (size_p == 1) RETURN
1051 
1052  ALLOCATE (monog1((grad + 1)*(grad + 2)/2*minressize), &
1053  monog2((grad + 1)*(grad + 2)/2*minressize))
1054  monog1 = 0
1055  monog2 = 0
1056  ii1 = 1
1057  DO j = 1, 3
1058  DO i = 1, 4
1059  monog1(ii1) = basepoly(i, j)
1060  ii1 = ii1 + 1
1061  END DO
1062  END DO
1063  ii1 = 2
1064  igrad = 1
1065  monodim1 = 4
1066  monosize1 = 3
1067  monofulldim1 = monodim1*monosize1
1068  rest_size_p = size_p - 1
1069  DO
1070  k = min(rest_size_p, monosize1)
1071  !call dgemm('T','N',monoDim1,my_npoly,k,&
1072  ! 1.0_dp,monoG1,monoDim1,p(ii1:),size_p,1.0_dp,pRes,size_res)
1073  resshift = 0
1074  pshift = ii1
1075  DO ipoly = 1, my_npoly
1076  pidx = pshift
1077  ii = 1
1078  DO pcoeff = 1, k
1079  DO rescoeff = 1, monodim1
1080  pres(resshift + rescoeff) = pres(resshift + rescoeff) + p(pidx)*monog1(ii)
1081  ii = ii + 1
1082  END DO
1083  pidx = pidx + 1
1084  END DO
1085  resshift = resshift + size_res
1086  pshift = pshift + size_p
1087  END DO
1088 
1089  rest_size_p = rest_size_p - k
1090  ii1 = ii1 + k
1091  IF (rest_size_p <= 0) EXIT
1092 
1093  monosize2 = igrad + 2 + monosize1
1094  monodim2 = monodim1 + monosize2
1095  monofulldim2 = monosize2*monodim2
1096  monodimatt = monosize1*monodim2
1097  CALL poly_mult3ab(if_check(monog1(1:monofulldim1), monog1(1)), monofulldim1, igrad, &
1098  if_check(basepoly(:, 1), basepoly(1, 1)), &
1099  if_check(monog2(1:monodimatt), monog2(1)), monodimatt, monosize1, .false.)
1100  monodimatt2 = monofulldim2 - monodim2
1101  start_idx1 = (monosize1 - igrad - 1)*monodim1
1102  CALL poly_mult3ab(if_check(monog1(start_idx1 + 1:monofulldim1), monog1(start_idx1 + 1)), &
1103  monofulldim1 - start_idx1, igrad, if_check(basepoly(:, 2), basepoly(1, 2)), &
1104  if_check(monog2(monodimatt + 1:monodimatt2), monog2(monodimatt + 1)), &
1105  monodimatt2 - monodimatt, igrad + 1, .false.)
1106  CALL poly_mult3ab(if_check(monog1(monofulldim1 - monodim1 + 1:monofulldim1), monog1(monofulldim1 - monodim1 + 1)), &
1107  monodim1, igrad, if_check(basepoly(:, 3), basepoly(1, 3)), &
1108  if_check(monog2(monodimatt2 + 1:monofulldim2), monog2(monodimatt2 + 1)), &
1109  monofulldim2 - monodimatt2, 1, .false.)
1110  igrad = igrad + 1
1111 
1112  ! even grads
1113 
1114  k = min(rest_size_p, monosize2)
1115  !call dgemm('T','N',monoDim2,my_npoly,k,&
1116  ! 1.0_dp,monoG2,monoDim2,p(ii1:),size_p,1.0_dp,pRes,size_res)
1117  resshift = 0
1118  pshift = ii1
1119  DO ipoly = 1, my_npoly
1120  pidx = pshift
1121  ii = 1
1122  DO pcoeff = 1, k
1123  DO rescoeff = 1, monodim2
1124  pres(resshift + rescoeff) = pres(resshift + rescoeff) + p(pidx)*monog2(ii)
1125  ii = ii + 1
1126  END DO
1127  pidx = pidx + 1
1128  END DO
1129  resshift = resshift + size_res
1130  pshift = pshift + size_p
1131  END DO
1132 
1133  rest_size_p = rest_size_p - k
1134  ii1 = ii1 + k
1135  IF (rest_size_p <= 0) EXIT
1136 
1137  monosize1 = igrad + 2 + monosize2
1138  monodim1 = monodim2 + monosize1
1139  monofulldim1 = monosize1*monodim1
1140  monodimatt = monosize2*monodim1
1141  CALL poly_mult3ab(if_check(monog2(1:monofulldim2), monog2(1)), monofulldim2, igrad, &
1142  if_check(basepoly(:, 1), basepoly(1, 1)), &
1143  if_check(monog1(1:monodimatt), monog1(1)), monodimatt, monosize2, .false.)
1144  monodimatt2 = monofulldim1 - monodim1
1145  start_idx1 = (monosize2 - igrad - 1)*monodim2
1146  CALL poly_mult3ab(if_check(monog2(start_idx1 + 1:monofulldim2), monog2(start_idx1 + 1)), &
1147  monofulldim2 - start_idx1, igrad, &
1148  if_check(basepoly(:, 2), basepoly(1, 2)), &
1149  if_check(monog1(monodimatt + 1:monodimatt2), monog1(monodimatt + 1)), monodimatt2 - monodimatt, &
1150  igrad + 1, .false.)
1151  CALL poly_mult3ab(if_check(monog2(monofulldim2 - monodim2 + 1:monofulldim2), monog2(monofulldim2 - monodim2 + 1)), &
1152  monodim2, igrad, if_check(basepoly(:, 3), basepoly(1, 3)), &
1153  if_check(monog1(monodimatt2 + 1:monofulldim1), monog1(monodimatt2 + 1)), &
1154  monofulldim1 - monodimatt2, 1, .false.)
1155  igrad = igrad + 1
1156 
1157  ! ! alternative to unrolling
1158  ! monoG1=monoG2
1159  ! monoSize1=monoSize2
1160  ! monoDim1=monoDim2
1161  ! monoFullDim1=monoFullDim2
1162  END DO
1163  DEALLOCATE (monog1, monog2)
1164  END SUBROUTINE
1165 
1166 ! **************************************************************************************************
1167 !> \brief evaluates the 3d polymial at x (the result is a polynomial in two variables)
1168 !> \param p ...
1169 !> \param x ...
1170 !> \param pRes ...
1171 !> \param npoly ...
1172 ! **************************************************************************************************
1173  SUBROUTINE poly_p_eval3(p, x, pRes, npoly)
1174  REAL(dp), DIMENSION(:), INTENT(in) :: p
1175  REAL(dp), INTENT(in) :: x
1176  REAL(dp), DIMENSION(:), INTENT(inout) :: pres
1177  INTEGER, INTENT(in), OPTIONAL :: npoly
1178 
1179  INTEGER :: grad, my_npoly, newsize, size_p
1180  REAL(dp), ALLOCATABLE, DIMENSION(:) :: xi
1181 
1182  IF (.NOT. module_initialized) cpabort("module d3_poly not initialized")
1183  my_npoly = 1
1184  IF (PRESENT(npoly)) my_npoly = npoly
1185  size_p = SIZE(p)/my_npoly
1186  grad = grad_size3(size_p)
1187  newsize = SIZE(pres)/my_npoly
1188  cpassert(newsize >= poly_size2(grad))
1189  pres = 0.0
1190  ALLOCATE (xi(grad + 1))
1191  CALL poly_p_eval3b(p, SIZE(p), x, pres, SIZE(pres), my_npoly, grad, xi)
1192  DEALLOCATE (xi)
1193  END SUBROUTINE
1194 
1195 ! **************************************************************************************************
1196 !> \brief low level routine of poly_p_eval3 without checks
1197 !> \param p ...
1198 !> \param size_p ...
1199 !> \param x ...
1200 !> \param pRes ...
1201 !> \param size_pRes ...
1202 !> \param npoly ...
1203 !> \param grad ...
1204 !> \param xi ...
1205 ! **************************************************************************************************
1206  SUBROUTINE poly_p_eval3b(p, size_p, x, pRes, size_pRes, npoly, grad, xi)
1207  INTEGER, INTENT(in) :: size_p
1208  REAL(dp), DIMENSION(IF_CHECK(size_p, *)), &
1209  INTENT(in) :: p
1210  REAL(dp), INTENT(in) :: x
1211  INTEGER, INTENT(in) :: size_pres
1212  REAL(dp), DIMENSION(IF_CHECK(size_pRes, *)), &
1213  INTENT(inout) :: pres
1214  INTEGER, INTENT(in) :: npoly, grad
1215  REAL(dp), DIMENSION(IF_CHECK(grad+1, *)), &
1216  INTENT(inout) :: xi
1217 
1218  INTEGER :: i, igrad, ii, ii0, insize, ipoly, j, &
1219  msize_p, newsize, pshift, shiftres, &
1220  shiftres_0, subg
1221 
1222  IF (.NOT. module_initialized) cpabort("module d3_poly not initialized")
1223  insize = size_p/npoly
1224  newsize = size_pres/npoly
1225  pres(1:size_pres) = 0.0
1226  xi(1) = 1.0
1227  DO i = 1, grad
1228  xi(i + 1) = xi(i)*x
1229  END DO
1230  shiftres = 0
1231  pshift = 0
1232  DO ipoly = 1, npoly
1233  DO ii = 1, min(insize, cached_dim3)
1234  pres(shiftres + a_reduce_idx3(ii)) = pres(shiftres + a_reduce_idx3(ii)) + p(pshift + ii)*xi(a_mono_exp3(1, ii) + 1)
1235  END DO
1236  shiftres = shiftres + newsize
1237  pshift = pshift + insize
1238  END DO
1239  IF (grad > max_grad3) THEN
1240  ii0 = (max_grad3 + 1)*(max_grad3 + 2)*(max_grad3 + 3)/6 + 1
1241  shiftres_0 = 1
1242  msize_p = insize
1243  DO ipoly = 1, npoly
1244  ii = ii0
1245  grad_do: DO igrad = max_grad3 + 1, grad
1246  !ii=igrad*(igrad+1)*(igrad+2)/6+1
1247  shiftres = shiftres_0
1248  subg = 0
1249  DO i = igrad, 0, -1
1250  !subG=igrad-i
1251  !shiftRes=subG*(subG+3)/2+1
1252  DO j = subg, 0, -1
1253  IF (msize_p < ii) EXIT grad_do
1254  pres(shiftres - j) = pres(shiftres - j) + p(ii)*xi(i + 1)
1255  ii = ii + 1
1256  END DO
1257  shiftres = shiftres + subg + 2
1258  subg = subg + 1
1259  END DO
1260  END DO grad_do
1261  ii0 = ii0 + insize
1262  shiftres_0 = shiftres_0 + newsize
1263  msize_p = msize_p + insize
1264  END DO
1265  END IF
1266  END SUBROUTINE
1267 
1268 ! **************************************************************************************************
1269 !> \brief unevaluates a 2d polymial to a 3d polynomial at x
1270 !> p(a,b,c)=p(a,b,c)+sum(pRes(b,c)*(x*a)^i,i), this is *not* the inverse of poly_p_eval3
1271 !> adds to p
1272 !> \param p ...
1273 !> \param x ...
1274 !> \param pRes ...
1275 !> \param npoly ...
1276 ! **************************************************************************************************
1277  SUBROUTINE poly_padd_uneval3(p, x, pRes, npoly)
1278  REAL(dp), DIMENSION(:), INTENT(inout) :: p
1279  REAL(dp), INTENT(in) :: x
1280  REAL(dp), DIMENSION(:), INTENT(in) :: pres
1281  INTEGER, INTENT(in), OPTIONAL :: npoly
1282 
1283  INTEGER :: grad, my_npoly, newsize, size_p
1284  REAL(dp), ALLOCATABLE, DIMENSION(:) :: xi
1285 
1286  IF (.NOT. module_initialized) cpabort("module d3_poly not initialized")
1287  my_npoly = 1
1288  IF (PRESENT(npoly)) my_npoly = npoly
1289  size_p = SIZE(p)/my_npoly
1290  newsize = SIZE(pres)/my_npoly
1291  grad = grad_size2(newsize)
1292  cpassert(size_p >= poly_size3(grad))
1293  cpassert(newsize == poly_size2(grad))
1294  ALLOCATE (xi(grad + 1))
1295  CALL poly_padd_uneval3b(p, SIZE(p), x, pres, SIZE(pres), my_npoly, grad, xi)
1296  DEALLOCATE (xi)
1297  END SUBROUTINE
1298 
1299 ! **************************************************************************************************
1300 !> \brief low level routine of poly_padd_uneval3 without checks
1301 !> \param p ...
1302 !> \param size_p ...
1303 !> \param x ...
1304 !> \param pRes ...
1305 !> \param size_pRes ...
1306 !> \param npoly ...
1307 !> \param grad ...
1308 !> \param xi ...
1309 !> \note loop should be structured differently (more contiguous pRes access)
1310 ! **************************************************************************************************
1311  SUBROUTINE poly_padd_uneval3b(p, size_p, x, pRes, size_pRes, npoly, grad, xi)
1312  INTEGER, INTENT(in) :: size_p
1313  REAL(dp), DIMENSION(IF_CHECK(size_p, *)), &
1314  INTENT(inout) :: p
1315  REAL(dp), INTENT(in) :: x
1316  INTEGER, INTENT(in) :: size_pres
1317  REAL(dp), DIMENSION(IF_CHECK(size_pRes, *)), &
1318  INTENT(in) :: pres
1319  INTEGER, INTENT(in) :: npoly, grad
1320  REAL(dp), DIMENSION(IF_CHECK(grad+1, *)), &
1321  INTENT(inout) :: xi
1322 
1323  INTEGER :: i, igrad, ii, ii0, insize, ipoly, j, &
1324  msize_p, newsize, pshift, shiftres, &
1325  shiftres_0, subg, upsize
1326 
1327  IF (.NOT. module_initialized) cpabort("module d3_poly not initialized")
1328  insize = size_p/npoly
1329  newsize = size_pres/npoly
1330  upsize = (grad + 1)*(grad + 2)*(grad + 3)/6
1331  xi(1) = 1.0
1332  DO i = 1, grad
1333  xi(i + 1) = xi(i)*x
1334  END DO
1335  shiftres = 0
1336  pshift = 0
1337  DO ipoly = 1, npoly
1338  DO ii = 1, min(upsize, cached_dim3)
1339  p(pshift + ii) = p(pshift + ii) + pres(shiftres + a_reduce_idx3(ii))*xi(a_mono_exp3(1, ii) + 1)
1340  END DO
1341  shiftres = shiftres + newsize
1342  pshift = pshift + insize
1343  END DO
1344  IF (grad > max_grad3) THEN
1345  ii0 = (max_grad3 + 1)*(max_grad3 + 2)*(max_grad3 + 3)/6 + 1
1346  shiftres_0 = 1
1347  msize_p = upsize
1348  DO ipoly = 1, npoly
1349  ii = ii0
1350  grad_do: DO igrad = max_grad3 + 1, grad
1351  !ii=igrad*(igrad+1)*(igrad+2)/6+1
1352  shiftres = shiftres_0
1353  subg = 0
1354  DO i = igrad, 0, -1
1355  !subG=igrad-i
1356  !shiftRes=subG*(subG+3)/2+1
1357  DO j = subg, 0, -1
1358  IF (msize_p < ii) EXIT grad_do
1359  p(ii) = p(ii) + pres(shiftres - j)*xi(i + 1)
1360  ii = ii + 1
1361  END DO
1362  shiftres = shiftres + subg + 2
1363  subg = subg + 1
1364  END DO
1365  END DO grad_do
1366  ii0 = ii0 + insize
1367  shiftres_0 = shiftres_0 + newsize
1368  msize_p = msize_p + insize
1369  END DO
1370  END IF
1371  END SUBROUTINE
1372 
1373 ! **************************************************************************************************
1374 !> \brief evaluates the 2d polynomial at x (the result is a polynomial in one variable)
1375 !> \param p ...
1376 !> \param x ...
1377 !> \param pRes ...
1378 !> \param npoly ...
1379 ! **************************************************************************************************
1380  SUBROUTINE poly_p_eval2(p, x, pRes, npoly)
1381  REAL(dp), DIMENSION(:), INTENT(in) :: p
1382  REAL(dp), INTENT(in) :: x
1383  REAL(dp), DIMENSION(:), INTENT(inout) :: pres
1384  INTEGER, INTENT(in), OPTIONAL :: npoly
1385 
1386  INTEGER :: grad, my_npoly, newsize, size_p
1387  REAL(dp), ALLOCATABLE, DIMENSION(:) :: xi
1388 
1389  IF (.NOT. module_initialized) cpabort("module d3_poly not initialized")
1390  my_npoly = 1
1391  IF (PRESENT(npoly)) my_npoly = npoly
1392  size_p = SIZE(p)/my_npoly
1393  grad = grad_size2(size_p)
1394  newsize = SIZE(pres)/my_npoly
1395  pres = 0.0_dp
1396  cpassert(newsize >= poly_size1(grad))
1397  ALLOCATE (xi(grad + 1))
1398  CALL poly_p_eval2b(p, SIZE(p), x, pres, SIZE(pres), my_npoly, grad, xi)
1399  DEALLOCATE (xi)
1400  END SUBROUTINE
1401 
1402 ! **************************************************************************************************
1403 !> \brief low level routine of poly_p_eval2 without checks
1404 !> \param p ...
1405 !> \param size_p ...
1406 !> \param x ...
1407 !> \param pRes ...
1408 !> \param size_pRes ...
1409 !> \param npoly ...
1410 !> \param grad ...
1411 !> \param xi ...
1412 ! **************************************************************************************************
1413  SUBROUTINE poly_p_eval2b(p, size_p, x, pRes, size_pRes, npoly, grad, xi)
1414  INTEGER, INTENT(in) :: size_p
1415  REAL(dp), DIMENSION(IF_CHECK(size_p, *)), &
1416  INTENT(in) :: p
1417  REAL(dp), INTENT(in) :: x
1418  INTEGER, INTENT(in) :: size_pres
1419  REAL(dp), DIMENSION(IF_CHECK(size_pRes, *)), &
1420  INTENT(inout) :: pres
1421  INTEGER, INTENT(in) :: npoly
1422  INTEGER :: grad
1423  REAL(dp), DIMENSION(IF_CHECK(grad+1, *)) :: xi
1424 
1425  INTEGER :: i, igrad, ii, ii0, ij, insize, ipoly, &
1426  msize_p, newsize, pshift, shiftres
1427 
1428  IF (.NOT. module_initialized) cpabort("module d3_poly not initialized")
1429  insize = size_p/npoly
1430  newsize = size_pres/npoly
1431  pres(1:size_pres) = 0.0_dp
1432  !CPPreconditionNoFail(newSize>grad,cp_failure_level,routineP)
1433  xi(1) = 1.0_dp
1434  DO i = 1, grad
1435  xi(i + 1) = xi(i)*x
1436  END DO
1437  shiftres = 1
1438  pshift = 0
1439  DO ipoly = 1, npoly
1440  DO ii = 1, min(insize, cached_dim2)
1441  pres(shiftres + a_mono_exp2(2, ii)) = pres(shiftres + a_mono_exp2(2, ii)) + p(pshift + ii)*xi(a_mono_exp2(1, ii) + 1)
1442  END DO
1443  shiftres = shiftres + newsize
1444  pshift = pshift + insize
1445  END DO
1446  IF (grad > max_grad2) THEN
1447  ii0 = (max_grad2 + 1)*(max_grad2 + 2)/2 + 1
1448  shiftres = 1
1449  msize_p = insize
1450  DO ipoly = 1, npoly
1451  ii = ii0
1452  grad_do2: DO igrad = max_grad2 + 1, grad
1453  !ii=igrad*(igrad+1)/2+1
1454  ij = shiftres
1455  DO i = igrad, 0, -1
1456  IF (msize_p < ii) EXIT grad_do2
1457  ! ij=igrad-i
1458  pres(ij) = pres(ij) + p(ii)*xi(i + 1)
1459  ii = ii + 1
1460  ij = ij + 1
1461  END DO
1462  END DO grad_do2
1463  msize_p = msize_p + insize
1464  shiftres = shiftres + newsize
1465  ii0 = ii0 + insize
1466  END DO
1467  END IF
1468  END SUBROUTINE
1469 
1470 ! **************************************************************************************************
1471 !> \brief unevaluates a 1d polynomial to 2d at x
1472 !> p(a,b)=sum(pRes(b)*(x*a)^i,i), this is *not* the inverse of poly_p_eval2
1473 !> overwrites p
1474 !> \param p ...
1475 !> \param x ...
1476 !> \param pRes ...
1477 !> \param npoly ...
1478 ! **************************************************************************************************
1479  SUBROUTINE poly_padd_uneval2(p, x, pRes, npoly)
1480  REAL(dp), DIMENSION(:), INTENT(inout) :: p
1481  REAL(dp), INTENT(in) :: x
1482  REAL(dp), DIMENSION(:), INTENT(in) :: pres
1483  INTEGER, INTENT(in), OPTIONAL :: npoly
1484 
1485  INTEGER :: grad, my_npoly, newsize, size_p
1486  REAL(dp), ALLOCATABLE, DIMENSION(:) :: xi
1487 
1488  IF (.NOT. module_initialized) cpabort("module d3_poly not initialized")
1489  my_npoly = 1
1490  IF (PRESENT(npoly)) my_npoly = npoly
1491  size_p = SIZE(p)/my_npoly
1492  newsize = SIZE(pres)/my_npoly
1493  grad = grad_size1(newsize)
1494  cpassert(size_p >= poly_size2(grad))
1495  cpassert(newsize == poly_size1(grad))
1496  ALLOCATE (xi(grad + 1))
1497  CALL poly_padd_uneval2b(p, SIZE(p), x, pres, SIZE(pres), my_npoly, grad, xi)
1498  DEALLOCATE (xi)
1499  END SUBROUTINE
1500 
1501 ! **************************************************************************************************
1502 !> \brief low level routine of poly_p_uneval2 without checks
1503 !> \param p ...
1504 !> \param size_p ...
1505 !> \param x ...
1506 !> \param pRes ...
1507 !> \param size_pRes ...
1508 !> \param npoly ...
1509 !> \param grad ...
1510 !> \param xi ...
1511 ! **************************************************************************************************
1512  SUBROUTINE poly_padd_uneval2b(p, size_p, x, pRes, size_pRes, npoly, grad, xi)
1513  INTEGER, INTENT(in) :: size_p
1514  REAL(dp), DIMENSION(IF_CHECK(size_p, *)), &
1515  INTENT(inout) :: p
1516  REAL(dp), INTENT(in) :: x
1517  INTEGER, INTENT(in) :: size_pres
1518  REAL(dp), DIMENSION(IF_CHECK(size_pRes, *)), &
1519  INTENT(in) :: pres
1520  INTEGER, INTENT(in) :: npoly
1521  INTEGER :: grad
1522  REAL(dp), DIMENSION(IF_CHECK(grad+1, *)) :: xi
1523 
1524  INTEGER :: i, igrad, ii, ii0, ij, insize, ipoly, &
1525  msize_p, newsize, pshift, shiftres, &
1526  upsize
1527 
1528  IF (.NOT. module_initialized) cpabort("module d3_poly not initialized")
1529  insize = size_p/npoly
1530  upsize = (grad + 1)*(grad + 2)/2
1531  newsize = size_pres/npoly
1532  !CPPreconditionNoFail(newSize>grad,cp_failure_level,routineP)
1533  xi(1) = 1.0_dp
1534  DO i = 1, grad
1535  xi(i + 1) = xi(i)*x
1536  END DO
1537  shiftres = 1
1538  pshift = 0
1539  DO ipoly = 1, npoly
1540  DO ii = 1, min(upsize, cached_dim2)
1541  p(pshift + ii) = p(pshift + ii) + pres(shiftres + a_mono_exp2(2, ii))*xi(a_mono_exp2(1, ii) + 1)
1542  END DO
1543  shiftres = shiftres + newsize
1544  pshift = pshift + insize
1545  END DO
1546  IF (grad > max_grad2) THEN
1547  ii0 = (max_grad2 + 1)*(max_grad2 + 2)/2 + 1
1548  shiftres = 1
1549  msize_p = upsize
1550  DO ipoly = 1, npoly
1551  ii = ii0
1552  grad_do2: DO igrad = max_grad2 + 1, grad
1553  !ii=igrad*(igrad+1)/2+1
1554  ij = shiftres
1555  DO i = igrad, 0, -1
1556  IF (msize_p < ii) EXIT grad_do2
1557  ! ij=igrad-i
1558  p(ii) = p(ii) + pres(ij)*xi(i + 1)
1559  ii = ii + 1
1560  ij = ij + 1
1561  END DO
1562  END DO grad_do2
1563  msize_p = msize_p + insize
1564  shiftres = shiftres + newsize
1565  ii0 = ii0 + insize
1566  END DO
1567  END IF
1568  END SUBROUTINE
1569 
1570 ! **************************************************************************************************
1571 !> \brief evaluates the 1d polynomial at the given place, results are stored contiguosly
1572 !> \param p ...
1573 !> \param x ...
1574 !> \param pRes ...
1575 !> \param npoly ...
1576 ! **************************************************************************************************
1577  SUBROUTINE poly_eval1(p, x, pRes, npoly)
1578  REAL(dp), DIMENSION(:), INTENT(in) :: p
1579  REAL(dp), INTENT(in) :: x
1580  REAL(dp), DIMENSION(:), INTENT(inout) :: pres
1581  INTEGER, INTENT(in), OPTIONAL :: npoly
1582 
1583  INTEGER :: i, ipoly, my_npoly, pshift, size_p
1584  REAL(dp) :: vv, xx
1585 
1586  my_npoly = 1
1587  IF (PRESENT(npoly)) my_npoly = npoly
1588  size_p = SIZE(p)/my_npoly
1589  cpassert(SIZE(pres) >= my_npoly)
1590  pshift = 0
1591  DO ipoly = 1, my_npoly
1592  xx = 1.0_dp
1593  vv = 0.0_dp
1594  DO i = 1, size_p
1595  vv = vv + p(pshift + i)*xx
1596  xx = xx*x
1597  END DO
1598  pres(ipoly) = vv
1599  pshift = pshift + size_p
1600  END DO
1601  END SUBROUTINE
1602 
1603 ! **************************************************************************************************
1604 !> \brief evaluates the 2d polynomial at the given place, results are stored contiguosly
1605 !> \param p ...
1606 !> \param x ...
1607 !> \param y ...
1608 !> \param pRes ...
1609 !> \param npoly ...
1610 ! **************************************************************************************************
1611  SUBROUTINE poly_eval2(p, x, y, pRes, npoly)
1612  REAL(dp), DIMENSION(:), INTENT(in) :: p
1613  REAL(dp), INTENT(in) :: x, y
1614  REAL(dp), DIMENSION(:), INTENT(inout) :: pres
1615  INTEGER, INTENT(in), OPTIONAL :: npoly
1616 
1617  INTEGER :: grad, i, igrad, ii, ipoly, j, msize_p, &
1618  my_npoly, pshift, size_p
1619  REAL(dp) :: v
1620  REAL(dp), ALLOCATABLE, DIMENSION(:) :: xi, yi
1621 
1622  IF (.NOT. module_initialized) cpabort("module d3_poly not initialized")
1623  my_npoly = 1
1624  IF (PRESENT(npoly)) my_npoly = npoly
1625  size_p = SIZE(p)/my_npoly
1626  grad = grad_size2(size_p)
1627  cpassert(SIZE(pres) >= my_npoly)
1628  ALLOCATE (xi(grad + 1), yi(grad + 1))
1629  xi(1) = 1.0_dp
1630  DO i = 1, grad
1631  xi(i + 1) = xi(i)*x
1632  END DO
1633  yi(1) = 1.0_dp
1634  DO i = 1, grad
1635  yi(i + 1) = yi(i)*y
1636  END DO
1637  pshift = 0
1638  DO ipoly = 1, my_npoly
1639  v = 0.0_dp
1640  DO ii = 1, min(size_p, cached_dim2)
1641  v = v + p(pshift + ii)*xi(a_mono_exp2(1, ii) + 1)*yi(a_mono_exp2(2, ii) + 1)
1642  END DO
1643  pres(ipoly) = v
1644  pshift = pshift + size_p
1645  END DO
1646  IF (grad > max_grad2) THEN
1647  pshift = (max_grad2 + 1)*(max_grad2 + 2)/2 + 1
1648  msize_p = size_p
1649  DO ipoly = 1, my_npoly
1650  ii = pshift
1651  v = 0.0_dp
1652  grad_do4: DO igrad = max_grad2 + 1, grad
1653  ! ii=igrad*(igrad+1)*(igrad+2)/6+1
1654  j = 1
1655  DO i = igrad, 0, -1
1656  IF (msize_p < ii) EXIT grad_do4
1657  v = v + p(ii)*xi(i + 1)*yi(j)
1658  j = j + 1
1659  ii = ii + 1
1660  END DO
1661  END DO grad_do4
1662  pres(ipoly) = pres(ipoly) + v
1663  pshift = pshift + size_p
1664  msize_p = msize_p + size_p
1665  END DO
1666  END IF
1667  END SUBROUTINE
1668 
1669 ! **************************************************************************************************
1670 !> \brief evaluates the 3d polynomial at the given place, results are stored contiguosly
1671 !> \param p ...
1672 !> \param x ...
1673 !> \param y ...
1674 !> \param z ...
1675 !> \param pRes ...
1676 !> \param npoly ...
1677 ! **************************************************************************************************
1678  SUBROUTINE poly_eval3(p, x, y, z, pRes, npoly)
1679  REAL(dp), DIMENSION(:), INTENT(in) :: p
1680  REAL(dp), INTENT(in) :: x, y, z
1681  REAL(dp), DIMENSION(:), INTENT(inout) :: pres
1682  INTEGER, INTENT(in), OPTIONAL :: npoly
1683 
1684  INTEGER :: grad, i, igrad, ii, ipoly, j, k, &
1685  msize_p, my_npoly, pshift, size_p
1686  REAL(dp) :: v
1687  REAL(dp), ALLOCATABLE, DIMENSION(:) :: xi, yi, zi
1688 
1689  IF (.NOT. module_initialized) cpabort("module d3_poly not initialized")
1690  my_npoly = 1
1691  IF (PRESENT(npoly)) my_npoly = npoly
1692  size_p = SIZE(p)/my_npoly
1693  grad = grad_size3(size_p)
1694  cpassert(SIZE(pres) >= my_npoly)
1695  ALLOCATE (xi(grad + 1), yi(grad + 1), zi(grad + 1))
1696  xi(1) = 1.0_dp
1697  DO i = 1, grad
1698  xi(i + 1) = xi(i)*x
1699  END DO
1700  yi(1) = 1.0_dp
1701  DO i = 1, grad
1702  yi(i + 1) = yi(i)*y
1703  END DO
1704  zi(1) = 1.0_dp
1705  DO i = 1, grad
1706  zi(i + 1) = zi(i)*z
1707  END DO
1708  pshift = 0
1709  DO ipoly = 1, my_npoly
1710  v = 0.0_dp
1711  DO ii = 1, min(size_p, cached_dim3)
1712  v = v + p(pshift + ii)*xi(a_mono_exp3(1, ii) + 1)*yi(a_mono_exp3(2, ii) + 1) &
1713  *zi(a_mono_exp3(3, ii) + 1)
1714  END DO
1715  pres(ipoly) = v
1716  pshift = pshift + size_p
1717  END DO
1718  IF (grad > max_grad3) THEN
1719  pshift = (max_grad3 + 1)*(max_grad3 + 2)*(max_grad3 + 3)/6 + 1
1720  msize_p = size_p
1721  DO ipoly = 1, my_npoly
1722  ii = pshift
1723  v = 0.0_dp
1724  grad_do3: DO igrad = max_grad3 + 1, grad
1725  ! ii=igrad*(igrad+1)*(igrad+2)/6+1
1726  DO i = igrad, 0, -1
1727  k = 1
1728  DO j = igrad - i, 0, -1
1729  ii = (ipoly - 1)*size_p + mono_index3(i, j, igrad - i - j) + 1
1730  IF (msize_p < ii) EXIT grad_do3
1731  v = v + p(ii)*xi(i + 1)*yi(j + 1)*zi(k)
1732  k = k + 1
1733  ii = ii + 1
1734  END DO
1735  END DO
1736  END DO grad_do3
1737  pres(ipoly) = pres(ipoly) + v
1738  pshift = pshift + size_p
1739  msize_p = msize_p + size_p
1740  END DO
1741  END IF
1742  END SUBROUTINE
1743 
1744 ! **************************************************************************************************
1745 !> \brief returns an array with all dp/dx, the all dp/dy, and finally all dp/dz
1746 !> \param p ...
1747 !> \param pRes ...
1748 !> \param npoly ...
1749 !> \param sumUp ...
1750 ! **************************************************************************************************
1751  SUBROUTINE poly_derive3(p, pRes, npoly, sumUp)
1752  REAL(dp), DIMENSION(:), INTENT(in) :: p
1753  REAL(dp), DIMENSION(:), INTENT(inout) :: pres
1754  INTEGER, INTENT(in), OPTIONAL :: npoly
1755  LOGICAL, INTENT(in), OPTIONAL :: sumup
1756 
1757  INTEGER :: grad, i, igrad, ii, ii2, ipoly, j, k, msize_p, my_npoly, newsize, pshift, size_p, &
1758  xderivshift, yderivshift, yshift, zderivshift, zshift
1759  LOGICAL :: my_sumup
1760 
1761  IF (.NOT. module_initialized) cpabort("module d3_poly not initialized")
1762  my_npoly = 1
1763  IF (PRESENT(npoly)) my_npoly = npoly
1764  my_sumup = .false.
1765  IF (PRESENT(sumup)) my_sumup = sumup
1766  size_p = SIZE(p)/my_npoly
1767  newsize = SIZE(pres)/(3*my_npoly)
1768  grad = grad_size3(size_p)
1769  cpassert(newsize >= poly_size3(grad))
1770  IF (.NOT. my_sumup) pres = 0
1771  xderivshift = 1
1772  yderivshift = my_npoly*newsize + 1
1773  zderivshift = 2*yderivshift - 1
1774  pshift = 0
1775  DO ipoly = 1, my_npoly
1776  ! split derivs to have less streams to follow (3 vs 5)?
1777  DO ii = 1, min(cached_dim3, size_p)
1778  pres(xderivshift + a_deriv_idx3(1, ii)) = pres(xderivshift + a_deriv_idx3(1, ii)) &
1779  + p(pshift + ii)*a_mono_exp3(1, ii)
1780  pres(yderivshift + a_deriv_idx3(2, ii)) = pres(yderivshift + a_deriv_idx3(2, ii)) &
1781  + p(pshift + ii)*a_mono_exp3(2, ii)
1782  pres(zderivshift + a_deriv_idx3(3, ii)) = pres(zderivshift + a_deriv_idx3(3, ii)) &
1783  + p(pshift + ii)*a_mono_exp3(3, ii)
1784  END DO
1785  xderivshift = xderivshift + newsize
1786  yderivshift = yderivshift + newsize
1787  zderivshift = zderivshift + newsize
1788  pshift = pshift + size_p
1789  END DO
1790  IF (grad > max_grad3) THEN
1791  xderivshift = 0
1792  yderivshift = my_npoly*newsize
1793  zderivshift = 2*yderivshift - 1
1794  msize_p = size_p
1795  xderivshift = max_grad3*(max_grad3 + 1)*(max_grad3 + 2)/6 + 1
1796  pshift = xderivshift + (max_grad3 + 1)*(max_grad3 + 2)/2
1797  DO ipoly = 1, my_npoly
1798  ii = pshift
1799  ii2 = xderivshift
1800  grad_do5: DO igrad = max_grad3 + 1, grad
1801  yshift = yderivshift
1802  zshift = zderivshift
1803  DO i = igrad, 0, -1
1804  k = 0
1805  DO j = igrad - i, 0, -1
1806  IF (ii > msize_p) EXIT grad_do5
1807  ! remove ifs?
1808  IF (i > 0) pres(ii2) = pres(ii2) + p(ii)*i
1809  IF (j > 0) pres(yshift + ii2) = pres(yshift + ii2) + p(ii)*j
1810  IF (k > 0) pres(zshift + ii2) = pres(zshift + ii2) + k*p(ii)
1811  ii = ii + 1
1812  ii2 = ii2 + 1
1813  k = k + 1
1814  END DO
1815  yshift = yshift - 1
1816  zshift = zshift - 1
1817  END DO
1818  ii2 = ii2 - igrad - 1
1819  END DO grad_do5
1820  pshift = pshift + size_p
1821  xderivshift = xderivshift + newsize
1822  msize_p = msize_p + size_p
1823  END DO
1824  END IF
1825  END SUBROUTINE
1826 
1827 ! **************************************************************************************************
1828 !> \brief subroutine that converts from the cp2k poly format to the d3 poly format
1829 !> \param poly_cp2k ...
1830 !> \param grad ...
1831 !> \param poly_d3 ...
1832 ! **************************************************************************************************
1833  SUBROUTINE poly_cp2k2d3(poly_cp2k, grad, poly_d3)
1834  REAL(dp), DIMENSION(:), INTENT(in) :: poly_cp2k
1835  INTEGER, INTENT(in) :: grad
1836  REAL(dp), DIMENSION(:), INTENT(out) :: poly_d3
1837 
1838  INTEGER :: cp_ii, i, j, k, mgrad, mgrad2, sgrad, &
1839  sgrad2, sgrad2k, sgrad3, sgrad3k, &
1840  shifti, shiftj, shiftk, size_p
1841 
1842  size_p = (grad + 1)*(grad + 2)*(grad + 3)/6
1843  cpassert(SIZE(poly_cp2k) >= size_p)
1844  cpassert(SIZE(poly_d3) >= size_p)
1845  cp_ii = 0
1846  sgrad2k = 0
1847  sgrad3k = 0
1848  DO k = 0, grad
1849  shiftk = k + 1
1850  sgrad2k = sgrad2k + k
1851  sgrad3k = sgrad3k + sgrad2k
1852  sgrad2 = sgrad2k
1853  sgrad3 = sgrad3k
1854  DO j = 0, grad - k
1855  sgrad = j + k
1856  mgrad2 = sgrad2
1857  shiftj = mgrad2 + shiftk
1858  mgrad = sgrad
1859  shifti = shiftj + sgrad3
1860  DO i = 0, grad - j - k
1861  cp_ii = cp_ii + 1
1862  poly_d3(shifti) = poly_cp2k(cp_ii)
1863  mgrad = mgrad + 1
1864  mgrad2 = mgrad2 + mgrad
1865  shifti = shifti + mgrad2
1866  END DO
1867  sgrad2 = sgrad2 + sgrad + 1
1868  sgrad3 = sgrad3 + sgrad2
1869  END DO
1870  END DO
1871  IF (SIZE(poly_d3) >= size_p) THEN
1872  poly_d3(size_p + 1:) = 0.0_dp
1873  END IF
1874  END SUBROUTINE
1875 
1876 ! **************************************************************************************************
1877 !> \brief subroutine that converts from the d3 poly format to the cp2k poly format
1878 !> \param poly_cp2k ...
1879 !> \param grad ...
1880 !> \param poly_d3 ...
1881 ! **************************************************************************************************
1882  SUBROUTINE poly_d32cp2k(poly_cp2k, grad, poly_d3)
1883  REAL(dp), DIMENSION(:), INTENT(out) :: poly_cp2k
1884  INTEGER, INTENT(in) :: grad
1885  REAL(dp), DIMENSION(:), INTENT(in) :: poly_d3
1886 
1887  INTEGER :: cp_ii, i, j, k, mgrad, mgrad2, sgrad, &
1888  sgrad2, sgrad2k, sgrad3, sgrad3k, &
1889  shifti, shiftj, shiftk, size_p
1890 
1891  size_p = (grad + 1)*(grad + 2)*(grad + 3)/6
1892  cpassert(SIZE(poly_cp2k) >= size_p)
1893  cpassert(SIZE(poly_d3) >= size_p)
1894  cp_ii = 0
1895  sgrad2k = 0
1896  sgrad3k = 0
1897  DO k = 0, grad
1898  shiftk = k + 1
1899  sgrad2k = sgrad2k + k
1900  sgrad3k = sgrad3k + sgrad2k
1901  sgrad2 = sgrad2k
1902  sgrad3 = sgrad3k
1903  DO j = 0, grad - k
1904  sgrad = j + k
1905  mgrad2 = sgrad2
1906  shiftj = mgrad2 + shiftk
1907  mgrad = sgrad
1908  shifti = shiftj + sgrad3
1909  DO i = 0, grad - j - k
1910  cp_ii = cp_ii + 1
1911  poly_cp2k(cp_ii) = poly_d3(shifti)
1912  mgrad = mgrad + 1
1913  mgrad2 = mgrad2 + mgrad
1914  shifti = shifti + mgrad2
1915  END DO
1916  sgrad2 = sgrad2 + sgrad + 1
1917  sgrad3 = sgrad3 + sgrad2
1918  END DO
1919  END DO
1920  IF (SIZE(poly_d3) >= size_p) THEN
1921  poly_cp2k(size_p + 1:) = 0.0_dp
1922  END IF
1923  END SUBROUTINE
1924 
1925 END MODULE d3_poly
d3_poly
Routines to efficiently handle dense polynomial in 3 variables up to a given degree....
Definition: d3_poly.F:23
d3_poly::cached_dim2
integer, parameter, public cached_dim2
Definition: d3_poly.F:54
d3_poly::poly_cp2k2d3
subroutine, public poly_cp2k2d3(poly_cp2k, grad, poly_d3)
subroutine that converts from the cp2k poly format to the d3 poly format
Definition: d3_poly.F:1834
d3_poly::poly_p_eval2b
subroutine, public poly_p_eval2b(p, size_p, x, pRes, size_pRes, npoly, grad, xi)
low level routine of poly_p_eval2 without checks
Definition: d3_poly.F:1414
d3_poly::poly_size1
pure integer function, public poly_size1(maxgrad)
size of a polynomial in x up to the given degree
Definition: d3_poly.F:155
d3_poly::poly_padd_uneval2b
subroutine, public poly_padd_uneval2b(p, size_p, x, pRes, size_pRes, npoly, grad, xi)
low level routine of poly_p_uneval2 without checks
Definition: d3_poly.F:1513
d3_poly::init_d3_poly_module
subroutine, public init_d3_poly_module()
initialization of the cache, is called by functions in this module that use cached values
Definition: d3_poly.F:74
d3_poly::poly_d32cp2k
subroutine, public poly_d32cp2k(poly_cp2k, grad, poly_d3)
subroutine that converts from the d3 poly format to the cp2k poly format
Definition: d3_poly.F:1883
d3_poly::poly_affine_t3
subroutine, public poly_affine_t3(p, m, b, pRes, npoly)
returns in the polynomials pRes the affine transformation x -> m*x+b of p
Definition: d3_poly.F:1015
d3_poly::poly_padd_uneval3b
subroutine, public poly_padd_uneval3b(p, size_p, x, pRes, size_pRes, npoly, grad, xi)
low level routine of poly_padd_uneval3 without checks
Definition: d3_poly.F:1312
d3_poly::poly_size2
pure integer function, public poly_size2(maxgrad)
size of a polynomial in x,y up to the given degree
Definition: d3_poly.F:167
d3_poly::poly_affine_t3t
subroutine, public poly_affine_t3t(p, m, b, pRes, npoly)
returns in the polynomials pRes the transpose of the affine transformation x -> m*x+b of p
Definition: d3_poly.F:855
d3_poly::cached_dim3
integer, parameter, public cached_dim3
Definition: d3_poly.F:55
d3_poly::cached_dim1
integer, parameter, public cached_dim1
Definition: d3_poly.F:53
d3_poly::grad_size3
pure integer function, public grad_size3(n)
max grad for a polynom of the given size
Definition: d3_poly.F:215
d3_poly::poly_p_eval3b
subroutine, public poly_p_eval3b(p, size_p, x, pRes, size_pRes, npoly, grad, xi)
low level routine of poly_p_eval3 without checks
Definition: d3_poly.F:1207
d3_poly::poly_size3
pure integer function, public poly_size3(maxgrad)
size of a polynomial in x,y,z up to the given degree
Definition: d3_poly.F:179
d3_poly::max_grad2
integer, parameter, public max_grad2
Definition: d3_poly.F:51
d3_poly::max_grad3
integer, parameter, public max_grad3
Definition: d3_poly.F:52
kinds
Defines the basic variable types.
Definition: kinds.F:23
kinds::dp
integer, parameter, public dp
Definition: kinds.F:34