(git:374b731)
Loading...
Searching...
No Matches
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! **************************************************************************************************
23MODULE 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
67CONTAINS
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
1925END 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_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::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_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_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::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_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_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_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::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_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::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::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::max_grad3
integer, parameter, public max_grad3
Definition d3_poly.F:52
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
kinds
Defines the basic variable types.
Definition kinds.F:23
kinds::dp
integer, parameter, public dp
Definition kinds.F:34