Line data Source code
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 28736 : 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 28736 : nthreads = 1
79 28736 : !$ nthreads = OMP_GET_NUM_THREADS()
80 28736 : IF (nthreads /= 1) CPABORT("init_d3_poly_module must not be called within OMP PARALLEL")
81 :
82 28736 : IF (module_initialized) RETURN
83 :
84 : ii = 1
85 44051 : DO grad = 0, max_grad2
86 176204 : DO i = grad, 0, -1
87 132153 : a_mono_exp2(1, ii) = i
88 132153 : a_mono_exp2(2, ii) = grad - i
89 169911 : ii = ii + 1
90 : END DO
91 : END DO
92 : ii = 1
93 31465 : DO grad = 0, max_grad3
94 94395 : DO i = grad, 0, -1
95 213962 : DO j = grad - i, 0, -1
96 125860 : a_mono_exp3(1, ii) = i
97 125860 : a_mono_exp3(2, ii) = j
98 125860 : a_mono_exp3(3, ii) = grad - i - j
99 188790 : ii = ii + 1
100 : END DO
101 : END DO
102 : END DO
103 132153 : DO ii = 1, cached_dim3
104 125860 : subG = a_mono_exp3(2, ii) + a_mono_exp3(3, ii)
105 132153 : a_reduce_idx3(ii) = subG*(subG + 1)/2 + a_mono_exp3(3, ii) + 1
106 : END DO
107 132153 : DO ii = 1, cached_dim3
108 125860 : IF (a_mono_exp3(1, ii) > 0) THEN
109 62930 : 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 62930 : a_deriv_idx3(1, ii) = 0
112 : END IF
113 125860 : IF (a_mono_exp3(2, ii) > 0) THEN
114 62930 : 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 62930 : a_deriv_idx3(2, ii) = 0
117 : END IF
118 132153 : IF (a_mono_exp3(3, ii) > 0) THEN
119 62930 : 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 62930 : a_deriv_idx3(3, ii) = 0
122 : END IF
123 : END DO
124 138446 : DO ii = 1, cached_dim2
125 1592129 : DO ij = ii, cached_dim2
126 4361049 : monoRes2 = a_mono_exp2(:, ii) + a_mono_exp2(:, ij)
127 1453683 : a_mono_mult2(ii, ij) = mono_index2(monoRes2(1), monoRes2(2)) + 1
128 1585836 : a_mono_mult2(ij, ii) = a_mono_mult2(ii, ij)
129 : END DO
130 : END DO
131 132153 : DO ii = 1, cached_dim3
132 1453683 : DO ij = ii, cached_dim3
133 5286120 : monoRes3 = a_mono_exp3(:, ii) + a_mono_exp3(:, ij)
134 1321530 : a_mono_mult3(ii, ij) = mono_index3(monoRes3(1), monoRes3(2), monoRes3(3)) + 1
135 1447390 : a_mono_mult3(ij, ii) = a_mono_mult3(ii, ij)
136 : END DO
137 : END DO
138 : ii = 1
139 132153 : DO i = 1, cached_dim3
140 635593 : DO j = 1, 4
141 503440 : a_mono_mult3a(j, i) = a_mono_mult3(j, i)
142 629300 : ii = ii + 1
143 : END DO
144 : END DO
145 :
146 6293 : 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 0 : PURE FUNCTION poly_size1(maxgrad) RESULT(res)
155 : INTEGER, INTENT(in) :: maxgrad
156 : INTEGER :: res
157 :
158 0 : res = maxgrad + 1
159 0 : 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 0 : PURE FUNCTION poly_size2(maxgrad) RESULT(res)
167 : INTEGER, INTENT(in) :: maxgrad
168 : INTEGER :: res
169 :
170 0 : res = (maxgrad + 1)*(maxgrad + 2)/2
171 0 : 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 0 : PURE FUNCTION poly_size3(maxgrad) RESULT(res)
179 : INTEGER, INTENT(in) :: maxgrad
180 : INTEGER :: res
181 :
182 0 : res = (maxgrad + 1)*(maxgrad + 2)*(maxgrad + 3)/6
183 0 : END FUNCTION
184 :
185 : ! **************************************************************************************************
186 : !> \brief max grad for a polynom of the given size
187 : !> \param n ...
188 : !> \return ...
189 : ! **************************************************************************************************
190 0 : PURE FUNCTION grad_size1(n) RESULT(res)
191 : INTEGER, INTENT(in) :: n
192 : INTEGER :: res
193 :
194 0 : res = n - 1
195 0 : END FUNCTION
196 :
197 : ! **************************************************************************************************
198 : !> \brief max grad for a polynom of the given size
199 : !> \param n ...
200 : !> \return ...
201 : ! **************************************************************************************************
202 0 : PURE FUNCTION grad_size2(n) RESULT(res)
203 : INTEGER, INTENT(in) :: n
204 : INTEGER :: res
205 :
206 0 : res = INT(FLOOR(0.5_dp*(SQRT(1.0_dp + 8.0_dp*REAL(n, dp)) - 1.0_dp) - 2.e-6_dp))
207 0 : END FUNCTION
208 :
209 : ! **************************************************************************************************
210 : !> \brief max grad for a polynom of the given size
211 : !> \param n ...
212 : !> \return ...
213 : ! **************************************************************************************************
214 0 : 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 0 : IF (n < 1) THEN
222 : res = -1
223 : ELSE
224 0 : nn = n*6
225 0 : g1 = (108.0_dp*nn + 12.0_dp*SQRT(81.0_dp*nn*nn - 12.0_dp))**(1.0_dp/3.0_dp)
226 0 : res = FLOOR(g1/6.0_dp + 2.0_dp/g1 - 1.0_dp - 2.e-6_dp)
227 : END IF
228 0 : END FUNCTION
229 :
230 : ! **************************************************************************************************
231 : !> \brief 0-based index of monomial of the given degree
232 : !> \param i ...
233 : !> \return ...
234 : ! **************************************************************************************************
235 0 : PURE FUNCTION mono_index1(i) RESULT(res)
236 : INTEGER, INTENT(in) :: i
237 : INTEGER :: res
238 :
239 0 : res = i
240 0 : 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 1453683 : PURE FUNCTION mono_index2(i, j) RESULT(res)
249 : INTEGER, INTENT(in) :: i, j
250 : INTEGER :: res
251 :
252 : INTEGER :: grad
253 :
254 1453683 : grad = i + j
255 1453683 : res = grad*(grad + 1)/2 + j
256 1453683 : 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 1510320 : 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 1510320 : sgrad = j + k
272 1510320 : grad = i + sgrad
273 1510320 : res = grad*(grad + 1)*(grad + 2)/6 + (sgrad)*(sgrad + 1)/2 + k
274 1510320 : END FUNCTION
275 :
276 : ! **************************************************************************************************
277 : !> \brief exponents of the monomial at the given 0-based index
278 : !> \param ii ...
279 : !> \return ...
280 : ! **************************************************************************************************
281 0 : PURE FUNCTION mono_exp1(ii) RESULT(res)
282 : INTEGER, INTENT(in) :: ii
283 : INTEGER :: res
284 :
285 0 : res = ii
286 0 : END FUNCTION
287 :
288 : ! **************************************************************************************************
289 : !> \brief exponents of the monomial at the given 0-based index
290 : !> \param ii ...
291 : !> \return ...
292 : ! **************************************************************************************************
293 0 : PURE FUNCTION mono_exp2(ii) RESULT(res)
294 : INTEGER, INTENT(in) :: ii
295 : INTEGER, DIMENSION(2) :: res
296 :
297 : INTEGER :: grad
298 :
299 0 : grad = INT(FLOOR(0.5_dp*(SQRT(9.0_dp + 8.0_dp*ii) - 1.0_dp) - 2.e-6_dp))
300 0 : res(2) = ii - (grad)*(grad + 1)/2
301 0 : res(1) = grad - res(2)
302 0 : END FUNCTION
303 :
304 : ! **************************************************************************************************
305 : !> \brief exponents of the monomial at the given 0-based index
306 : !> \param n ...
307 : !> \return ...
308 : ! **************************************************************************************************
309 0 : 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 0 : nn = (n + 1)*6
317 0 : g1 = (108.0_dp*nn + 12.0_dp*SQRT(81.0_dp*nn*nn - 12.0_dp))**(1.0_dp/3.0_dp)
318 0 : grad1 = INT(FLOOR(g1/6.0_dp + 2.0_dp/g1 - 1.0_dp - 2.e-6_dp))
319 0 : ii = n - grad1*(grad1 + 1)*(grad1 + 2)/6
320 0 : grad = INT(FLOOR(0.5_dp*(SQRT(9.0_dp + 8.0_dp*ii) - 1.0_dp) - 1.e-6_dp))
321 0 : res(3) = ii - grad*(grad + 1)/2
322 0 : res(2) = grad - res(3)
323 0 : res(1) = grad1 - grad
324 0 : 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 0 : PURE FUNCTION mono_mult1(ii, ij) RESULT(res)
333 : INTEGER, INTENT(in) :: ii, ij
334 : INTEGER :: res
335 :
336 0 : res = ii + ij
337 0 : 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 0 : 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 0 : monoRes = mono_exp2(ii) + mono_exp2(ij)
352 0 : res = mono_index2(monoRes(1), monoRes(2))
353 0 : 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 0 : 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 0 : monoRes = mono_exp3(ii) + mono_exp3(ij)
368 0 : res = mono_index3(monoRes(1), monoRes(2), monoRes(3))
369 0 : 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 0 : 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 0 : IF (.NOT. module_initialized) CPABORT("module d3_poly not initialized")
391 0 : mySumUp = .FALSE.
392 0 : myNp1 = 1
393 0 : IF (PRESENT(np1)) myNp1 = np1
394 0 : IF (PRESENT(sumUp)) mySumUp = sumUp
395 0 : size_p1 = SIZE(p1)/myNp1
396 0 : size_p2 = SIZE(p2)
397 0 : newGrad = grad_size1(size_p1) + grad_size1(size_p2)
398 0 : newSize = SIZE(pRes)/myNp1
399 0 : CPASSERT(newSize >= poly_size1(newGrad))
400 0 : IF (.NOT. mySumUp) pRes = 0
401 : iPos = 1
402 : resShift_0 = 0
403 0 : DO ipoly = 0, myNp1 - 1
404 0 : DO i = 1, size_p1
405 0 : resPos = resShift_0 + i
406 0 : DO j = 1, size_p2
407 0 : pRes(resPos) = pRes(resPos) + p1(iPos)*p2(j)
408 0 : resPos = resPos + 1
409 : END DO
410 0 : iPos = iPos + 1
411 : END DO
412 0 : resShift_0 = resShift_0 + newSize
413 : END DO
414 0 : 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 0 : 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 0 : IF (.NOT. module_initialized) CPABORT("module d3_poly not initialized")
435 0 : mySumUp = .FALSE.
436 0 : IF (PRESENT(sumUp)) mySumUp = sumUp
437 0 : myNp1 = 1
438 0 : IF (PRESENT(np1)) myNp1 = np1
439 0 : size_p1 = SIZE(p1)/myNp1
440 0 : size_p2 = SIZE(p2)
441 0 : grad1 = grad_size2(size_p1)
442 0 : grad2 = grad_size2(size_p2)
443 0 : newGrad = grad1 + grad2
444 0 : newSize = SIZE(pRes)/myNp1
445 0 : CPASSERT(newSize >= poly_size2(newGrad))
446 0 : IF (.NOT. mySumUp) pRes = 0
447 : iShift = 0
448 : shiftRes = 0
449 0 : DO ipoly = 1, myNp1
450 0 : DO i = 1, MIN(size_p1, cached_dim2)
451 0 : DO j = 1, MIN(size_p2, cached_dim2)
452 : pRes(shiftRes + a_mono_mult2(j, i)) = pRes(shiftRes + a_mono_mult2(j, i)) &
453 0 : + p1(iShift + i)*p2(j)
454 : END DO
455 : END DO
456 0 : iShift = iShift + size_p1
457 0 : shiftRes = shiftRes + newSize
458 : END DO
459 0 : IF (grad1 > max_grad2 .OR. grad2 > max_grad2) THEN
460 : msize_p1 = size_p1
461 : shiftRes_0 = 0
462 0 : DO ipoly = 0, myNp1 - 1
463 0 : shift1 = ipoly*size_p1
464 0 : DO g1 = 0, grad1
465 : ! shift1=g1*(g1+1)/2
466 0 : IF (g1 > max_grad2) THEN
467 0 : subG2 = 0
468 0 : shift2 = 0
469 0 : g1g2 = shiftRes_0 - 1
470 : ELSE
471 0 : subG2 = max_grad2 + 1
472 0 : shift2 = cached_dim2
473 0 : g1g2 = shiftRes_0 + g1*subG2 - 1
474 : END IF
475 0 : DO g2 = subG2, grad2
476 : ! shift2=g2*(g2+1)/2
477 0 : shiftRes = shift1 + shift2 + g1g2 ! shiftRes=(g1+g2)*(g1+g2+1)/2-1+ipoly*newSize
478 0 : DO i = 1, MIN(g1 + 1, msize_p1 - shift1)
479 0 : DO j = 1, MIN(g2 + 1, size_p2 - shift2)
480 0 : pRes(shiftRes + i + j) = pRes(shiftRes + i + j) + p1(shift1 + i)*p2(shift2 + j)
481 : END DO
482 : END DO
483 0 : shift2 = shift2 + g2 + 1 !
484 0 : g1g2 = g1g2 + g1 !
485 : END DO
486 0 : shift1 = shift1 + g1 + 1 !
487 : END DO
488 0 : shiftRes_0 = shiftRes_0 + newSize - size_p1
489 0 : msize_p1 = msize_p1 + size_p1
490 : END DO
491 : END IF
492 0 : 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 0 : 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 0 : IF (.NOT. module_initialized) CPABORT("module d3_poly not initialized")
513 0 : mySumUp = .FALSE.
514 0 : IF (PRESENT(sumUp)) mySumUp = sumUp
515 0 : myNp1 = 1
516 0 : IF (PRESENT(np1)) myNp1 = np1
517 0 : size_p1 = SIZE(p1)/myNp1
518 0 : size_p2 = SIZE(p2)
519 0 : grad1 = grad_size3(size_p1)
520 0 : grad2 = grad_size3(size_p2)
521 0 : newGrad = grad1 + grad2
522 0 : newSize = SIZE(pRes)/myNp1
523 0 : CPASSERT(newSize >= poly_size3(newGrad))
524 0 : CALL poly_mult3b(p1, SIZE(p1), grad1, p2, SIZE(p2), grad2, pRes, SIZE(pRes), myNp1, mySumUp)
525 0 : 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 0 : 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 0 : IF (.NOT. module_initialized) CPABORT("module d3_poly not initialized")
558 0 : my_size_p1 = size_p1/np1
559 0 : newSize = size_pRes/np1
560 :
561 0 : IF (.NOT. sumUp) pRes(1:size_pRes) = 0.0_dp
562 : iShift = 0
563 : shiftRes = 0
564 0 : DO ipoly = 1, np1
565 0 : DO i = 1, MIN(my_size_p1, cached_dim3)
566 0 : DO j = 1, MIN(size_p2, cached_dim3)
567 : pRes(shiftRes + a_mono_mult3(j, i)) = pRes(shiftRes + a_mono_mult3(j, i)) &
568 0 : + p1(iShift + i)*p2(j)
569 : END DO
570 : END DO
571 0 : iShift = iShift + my_size_p1
572 0 : shiftRes = shiftRes + newSize
573 : END DO
574 0 : IF (grad1 > max_grad3 .OR. grad2 > max_grad3) THEN
575 : ! one could remove multiplications even more aggressively...
576 : msize_p1 = my_size_p1
577 0 : DO ipoly = 0, np1 - 1
578 0 : shift1 = 1 + ipoly*my_size_p1
579 0 : shiftRes_0 = 1 + ipoly*newSize
580 0 : DO g1 = 0, grad1
581 0 : IF (g1 > max_grad3) THEN
582 : subG2 = 0
583 : shift2 = 1
584 : ELSE
585 0 : subG2 = max_grad3 + 1
586 0 : shift2 = subG2*(subG2 + 1)*(subG2 + 2)/6 + 1
587 : END IF
588 0 : DO g2 = subG2, grad2
589 0 : shiftRes = (g1 + g2)*(g1 + g2 + 1)*(g1 + g2 + 2)/6 + shiftRes_0
590 0 : shift1I = shift1
591 0 : shiftResI_0 = shiftRes
592 0 : DO i1 = g1, 0, -1
593 0 : IF (shift1I > msize_p1) EXIT
594 0 : shift2I = shift2
595 0 : shiftResI = shiftResI_0
596 0 : subGrad = g1 - i1
597 0 : 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 0 : IF (shift2I > size_p2) EXIT
602 0 : DO j1 = g1 - i1, 0, -1
603 0 : shift1J = shift1I + g1 - i1 - j1
604 0 : IF (shift1J > msize_p1) EXIT
605 0 : DO j2 = g2 - i2, 0, -1
606 0 : shift2J = shift2I + g2 - i2 - j2
607 0 : IF (shift2J > size_p2) EXIT
608 0 : 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 0 : pRes(shiftResJ) = pRes(shiftResJ) + p1(shift1J)*p2(shift2J)
613 : END DO
614 : END DO
615 0 : subGrad = subGrad + 1
616 0 : shift2I = shift2I + (g2 - i2 + 1)
617 0 : shiftResI = shiftResI + subGrad
618 : END DO
619 0 : shift1I = shift1I + (g1 - i1 + 1)
620 0 : shiftResI_0 = shiftResI_0 + (g1 - i1 + 1)
621 : END DO
622 0 : shift2 = shift2 + (g2 + 1)*(g2 + 2)/2
623 : END DO
624 0 : shift1 = shift1 + (g1 + 1)*(g1 + 2)/2
625 : END DO
626 0 : msize_p1 = msize_p1 + my_size_p1
627 : END DO
628 : END IF
629 0 : 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 0 : 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 0 : IF (.NOT. module_initialized) CPABORT("module d3_poly not initialized")
661 0 : my_size_p1 = size_p1/np1
662 0 : newSize = size_pRes/np1
663 :
664 0 : IF (.NOT. sumUp) pRes(1:size_pRes) = 0.0_dp
665 : iShift = 0
666 : shiftRes = 0
667 0 : DO ipoly = 1, np1
668 0 : 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 0 : + p1(iShift + i)*p2(1)
671 : pRes(shiftRes + a_mono_mult3a(2, i)) = pRes(shiftRes + a_mono_mult3a(2, i)) &
672 0 : + p1(iShift + i)*p2(2)
673 : pRes(shiftRes + a_mono_mult3a(3, i)) = pRes(shiftRes + a_mono_mult3a(3, i)) &
674 0 : + p1(iShift + i)*p2(3)
675 : pRes(shiftRes + a_mono_mult3a(4, i)) = pRes(shiftRes + a_mono_mult3a(4, i)) &
676 0 : + p1(iShift + i)*p2(4)
677 : END DO
678 0 : iShift = iShift + my_size_p1
679 0 : shiftRes = shiftRes + newSize
680 : END DO
681 0 : IF (grad1 > max_grad3 .OR. grad2 > max_grad3) THEN
682 : ! one could remove multiplications even more aggressively...
683 : msize_p1 = my_size_p1
684 0 : DO ipoly = 0, np1 - 1
685 0 : shift1 = 1 + ipoly*my_size_p1 + (max_grad3 + 1)*(max_grad3 + 2)*(max_grad3 + 3)/6
686 0 : shiftRes_0 = 1 + ipoly*newSize
687 0 : DO g1 = max_grad3 + 1, grad1
688 : subG2 = 0
689 : shift2 = 1
690 0 : DO g2 = subG2, grad2
691 0 : shiftRes = (g1 + g2)*(g1 + g2 + 1)*(g1 + g2 + 2)/6 + shiftRes_0
692 0 : shift1I = shift1
693 0 : shiftResI_0 = shiftRes
694 0 : DO i1 = g1, 0, -1
695 0 : IF (shift1I > msize_p1) EXIT
696 0 : shift2I = shift2
697 0 : shiftResI = shiftResI_0
698 0 : subGrad = g1 - i1
699 0 : 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 0 : DO j1 = g1 - i1, 0, -1
704 0 : shift1J = shift1I + g1 - i1 - j1
705 0 : IF (shift1J > msize_p1) EXIT
706 0 : DO j2 = g2 - i2, 0, -1
707 0 : shift2J = shift2I + g2 - i2 - j2
708 0 : 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 0 : pRes(shiftResJ) = pRes(shiftResJ) + p1(shift1J)*p2(shift2J)
713 : END DO
714 : END DO
715 0 : subGrad = subGrad + 1
716 0 : shift2I = shift2I + (g2 - i2 + 1)
717 0 : shiftResI = shiftResI + subGrad
718 : END DO
719 0 : shift1I = shift1I + (g1 - i1 + 1)
720 0 : shiftResI_0 = shiftResI_0 + (g1 - i1 + 1)
721 : END DO
722 0 : shift2 = shift2 + (g2 + 1)*(g2 + 2)/2
723 : END DO
724 0 : shift1 = shift1 + (g1 + 1)*(g1 + 2)/2
725 : END DO
726 0 : msize_p1 = msize_p1 + my_size_p1
727 : END DO
728 : END IF
729 0 : 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 0 : 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 0 : did_write = .FALSE.
744 0 : DO i = 1, SIZE(p)
745 0 : IF (p(i) /= 0) THEN
746 0 : IF (p(i) >= 0) WRITE (out_f, '("+")', advance='NO')
747 0 : WRITE (out_f, '(G20.10)', advance='NO') p(i)
748 0 : IF (i /= 1) WRITE (out_f, '("*x^",I3)', advance='NO') i - 1
749 : did_write = .TRUE.
750 : END IF
751 : END DO
752 0 : IF (.NOT. did_write) WRITE (out_f, '("0.0")', advance='NO')
753 0 : 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 0 : 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 0 : IF (.NOT. module_initialized) CPABORT("module d3_poly not initialized")
769 0 : did_write = .FALSE.
770 0 : DO i = 1, SIZE(p)
771 0 : mono_e = mono_exp2(i - 1)
772 0 : IF (p(i) /= 0) THEN
773 0 : IF (p(i) >= 0) WRITE (out_f, '("+")', advance='NO')
774 0 : WRITE (out_f, '(G20.10)', advance='NO') p(i)
775 0 : IF (mono_e(1) /= 0) WRITE (out_f, '("*x^",I3)', advance='NO') mono_e(1)
776 0 : 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 0 : IF (.NOT. did_write) WRITE (out_f, '("0.0")', advance='NO')
781 0 : END SUBROUTINE
782 :
783 : ! **************************************************************************************************
784 : !> \brief writes out the polynomial in a human readable form
785 : !> \param p ...
786 : !> \param out_f ...
787 : ! **************************************************************************************************
788 0 : 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 0 : IF (.NOT. module_initialized) CPABORT("module d3_poly not initialized")
797 0 : did_write = .FALSE.
798 0 : DO i = 1, SIZE(p)
799 0 : mono_e = mono_exp3(i - 1)
800 0 : IF (p(i) /= 0) THEN
801 0 : IF (p(i) >= 0) WRITE (out_f, '("+")', advance='NO')
802 0 : WRITE (out_f, '(G20.10)', advance='NO') p(i)
803 0 : IF (mono_e(1) /= 0) WRITE (out_f, '("*x^",I3)', advance='NO') mono_e(1)
804 0 : IF (mono_e(2) /= 0) WRITE (out_f, '("*y^",I3)', advance='NO') mono_e(2)
805 0 : 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 0 : IF (.NOT. did_write) WRITE (out_f, '("0.0")', advance='NO')
810 0 : 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 0 : 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 0 : IF (.NOT. module_initialized) CPABORT("module d3_poly not initialized")
830 0 : myMinSize = 1
831 0 : IF (PRESENT(minSize)) myMinSize = minSize
832 0 : CALL RANDOM_NUMBER(g)
833 0 : pSize = MIN(maxSize, myMinSize + INT((maxSize - myMinSize + 1)*g))
834 0 : CPASSERT(SIZE(p) >= pSize)
835 0 : CALL RANDOM_NUMBER(p)
836 0 : DO i = 1, pSize
837 0 : p(i) = REAL(INT(p(i)*200.0_dp - 100.0_dp), dp)/100.0_dp
838 : END DO
839 0 : DO i = pSize + 1, SIZE(p)
840 0 : p(i) = 0.0_dp
841 : END DO
842 0 : res = pSize
843 0 : 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 0 : 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 0 : REAL(dp), ALLOCATABLE, DIMENSION(:) :: monoG1, monoG2
865 : REAL(dp), DIMENSION(4, 3) :: basepoly
866 :
867 0 : IF (.NOT. module_initialized) CPABORT("module d3_poly not initialized")
868 0 : my_npoly = 1
869 0 : IF (PRESENT(npoly)) my_npoly = npoly
870 0 : basepoly(1, :) = b
871 0 : DO j = 1, 3
872 0 : DO i = 1, 3
873 0 : basepoly(j + 1, i) = m(i, j)
874 : END DO
875 : END DO
876 0 : size_p = SIZE(pRes)/my_npoly
877 0 : size_res = SIZE(p)/my_npoly
878 0 : grad = grad_size3(size_res)
879 0 : minResSize = poly_size3(grad)
880 0 : CPASSERT(size_res == minResSize)
881 0 : CPASSERT(size_p >= minResSize)
882 0 : pRes = 0
883 0 : IF (size_p == 0) RETURN
884 : ii1 = 1
885 : ii = 1
886 0 : DO ipoly = 1, my_npoly
887 0 : pRes(ii1) = p(ii)
888 0 : ii = ii + size_res
889 0 : ii1 = ii1 + size_p
890 : END DO
891 0 : IF (size_p == 1) RETURN
892 :
893 : ALLOCATE (monoG1((grad + 1)*(grad + 2)/2*minResSize), &
894 0 : monoG2((grad + 1)*(grad + 2)/2*minResSize))
895 : !monoG1=0
896 : !monoG2=0
897 0 : ii1 = 1
898 0 : DO j = 1, 3
899 0 : DO i = 1, 4
900 0 : monoG1(ii1) = basepoly(i, j)
901 0 : ii1 = ii1 + 1
902 : END DO
903 : END DO
904 0 : ii1 = 2
905 0 : igrad = 1
906 0 : monoDim1 = 4
907 0 : monoSize1 = 3
908 0 : monoFullDim1 = monoDim1*monoSize1
909 0 : rest_size_p = size_p - 1
910 0 : DO
911 0 : 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 0 : resShift = 0
915 0 : pShift = ii1
916 0 : DO ipoly = 1, my_npoly
917 : pIdx = pShift
918 : ii = 1
919 0 : DO pcoeff = 1, k
920 0 : DO rescoeff = 1, monoDim1
921 0 : pRes(pIdx) = pRes(pIdx) + p(resShift + rescoeff)*monoG1(ii)
922 0 : ii = ii + 1
923 : END DO
924 0 : pIdx = pIdx + 1
925 : END DO
926 0 : resShift = resShift + size_res
927 0 : pShift = pShift + size_p
928 : END DO
929 :
930 0 : rest_size_p = rest_size_p - k
931 0 : ii1 = ii1 + k
932 0 : IF (rest_size_p <= 0) EXIT
933 :
934 0 : monoSize2 = igrad + 2 + monoSize1
935 0 : monoDim2 = monoDim1 + monoSize2
936 0 : monoFullDim2 = monoSize2*monoDim2
937 0 : 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 0 : IF_CHECK(monoG2(1:monoDimAtt), monoG2(1)), monoDimAtt, monoSize1, .FALSE.)
941 0 : monoDimAtt2 = monoFullDim2 - monoDim2
942 0 : 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 0 : 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 0 : monoFullDim2 - monoDimAtt2, 1, .FALSE.)
951 0 : igrad = igrad + 1
952 :
953 : ! even grads
954 :
955 0 : 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 0 : resShift = 0
959 0 : pShift = ii1
960 0 : DO ipoly = 1, my_npoly
961 : pIdx = pShift
962 : ii = 1
963 0 : DO pcoeff = 1, k
964 0 : DO rescoeff = 1, monoDim2
965 0 : pRes(pIdx) = pRes(pIdx) + p(resShift + rescoeff)*monoG2(ii)
966 0 : ii = ii + 1
967 : END DO
968 0 : pIdx = pIdx + 1
969 : END DO
970 0 : resShift = resShift + size_res
971 0 : pShift = pShift + size_p
972 : END DO
973 :
974 0 : rest_size_p = rest_size_p - k
975 0 : ii1 = ii1 + k
976 0 : IF (rest_size_p <= 0) EXIT
977 :
978 0 : monoSize1 = igrad + 2 + monoSize2
979 0 : monoDim1 = monoDim2 + monoSize1
980 0 : monoFullDim1 = monoSize1*monoDim1
981 0 : 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 0 : monoDimAtt, monoSize2, .FALSE.)
985 0 : monoDimAtt2 = monoFullDim1 - monoDim1
986 0 : 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 0 : 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 0 : monoFullDim1 - monoDimAtt2, 1, .FALSE.)
995 0 : igrad = igrad + 1
996 :
997 : ! ! alternative to unrolling
998 : ! monoG1=monoG2
999 : ! monoSize1=monoSize2
1000 : ! monoDim1=monoDim2
1001 : ! monoFullDim1=monoFullDim2
1002 : END DO
1003 0 : 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 0 : 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 0 : REAL(dp), ALLOCATABLE, DIMENSION(:) :: monoG1, monoG2
1025 : REAL(dp), DIMENSION(4, 3) :: basepoly
1026 :
1027 0 : IF (.NOT. module_initialized) CPABORT("module d3_poly not initialized")
1028 0 : my_npoly = 1
1029 0 : IF (PRESENT(npoly)) my_npoly = npoly
1030 0 : basepoly(1, :) = b
1031 0 : DO j = 1, 3
1032 0 : DO i = 1, 3
1033 0 : basepoly(j + 1, i) = m(i, j)
1034 : END DO
1035 : END DO
1036 0 : size_p = SIZE(p)/my_npoly
1037 0 : grad = grad_size3(size_p)
1038 0 : size_res = SIZE(pRes)/my_npoly
1039 0 : minResSize = poly_size3(grad)
1040 0 : CPASSERT(size_res >= minResSize)
1041 0 : pRes = 0
1042 0 : IF (size_p == 0) RETURN
1043 : ii1 = 1
1044 : ii = 1
1045 0 : DO ipoly = 1, my_npoly
1046 0 : pRes(ii) = p(ii1)
1047 0 : ii = ii + size_res
1048 0 : ii1 = ii1 + size_p
1049 : END DO
1050 0 : IF (size_p == 1) RETURN
1051 :
1052 : ALLOCATE (monoG1((grad + 1)*(grad + 2)/2*minResSize), &
1053 0 : monoG2((grad + 1)*(grad + 2)/2*minResSize))
1054 0 : monoG1 = 0
1055 0 : monoG2 = 0
1056 : ii1 = 1
1057 0 : DO j = 1, 3
1058 0 : DO i = 1, 4
1059 0 : monoG1(ii1) = basepoly(i, j)
1060 0 : ii1 = ii1 + 1
1061 : END DO
1062 : END DO
1063 0 : ii1 = 2
1064 0 : igrad = 1
1065 0 : monoDim1 = 4
1066 0 : monoSize1 = 3
1067 0 : monoFullDim1 = monoDim1*monoSize1
1068 0 : rest_size_p = size_p - 1
1069 0 : DO
1070 0 : 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 0 : resShift = 0
1074 0 : pShift = ii1
1075 0 : DO ipoly = 1, my_npoly
1076 : pIdx = pShift
1077 : ii = 1
1078 0 : DO pcoeff = 1, k
1079 0 : DO rescoeff = 1, monoDim1
1080 0 : pRes(resShift + rescoeff) = pRes(resShift + rescoeff) + p(pIdx)*monoG1(ii)
1081 0 : ii = ii + 1
1082 : END DO
1083 0 : pIdx = pIdx + 1
1084 : END DO
1085 0 : resShift = resShift + size_res
1086 0 : pShift = pShift + size_p
1087 : END DO
1088 :
1089 0 : rest_size_p = rest_size_p - k
1090 0 : ii1 = ii1 + k
1091 0 : IF (rest_size_p <= 0) EXIT
1092 :
1093 0 : monoSize2 = igrad + 2 + monoSize1
1094 0 : monoDim2 = monoDim1 + monoSize2
1095 0 : monoFullDim2 = monoSize2*monoDim2
1096 0 : 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 0 : IF_CHECK(monoG2(1:monoDimAtt), monoG2(1)), monoDimAtt, monoSize1, .FALSE.)
1100 0 : monoDimAtt2 = monoFullDim2 - monoDim2
1101 0 : 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 0 : 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 0 : monoFullDim2 - monoDimAtt2, 1, .FALSE.)
1110 0 : igrad = igrad + 1
1111 :
1112 : ! even grads
1113 :
1114 0 : 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 0 : resShift = 0
1118 0 : pShift = ii1
1119 0 : DO ipoly = 1, my_npoly
1120 : pIdx = pShift
1121 : ii = 1
1122 0 : DO pcoeff = 1, k
1123 0 : DO rescoeff = 1, monoDim2
1124 0 : pRes(resShift + rescoeff) = pRes(resShift + rescoeff) + p(pIdx)*monoG2(ii)
1125 0 : ii = ii + 1
1126 : END DO
1127 0 : pIdx = pIdx + 1
1128 : END DO
1129 0 : resShift = resShift + size_res
1130 0 : pShift = pShift + size_p
1131 : END DO
1132 :
1133 0 : rest_size_p = rest_size_p - k
1134 0 : ii1 = ii1 + k
1135 0 : IF (rest_size_p <= 0) EXIT
1136 :
1137 0 : monoSize1 = igrad + 2 + monoSize2
1138 0 : monoDim1 = monoDim2 + monoSize1
1139 0 : monoFullDim1 = monoSize1*monoDim1
1140 0 : 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 0 : IF_CHECK(monoG1(1:monoDimAtt), monoG1(1)), monoDimAtt, monoSize2, .FALSE.)
1144 0 : monoDimAtt2 = monoFullDim1 - monoDim1
1145 0 : 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 0 : 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 0 : monoFullDim1 - monoDimAtt2, 1, .FALSE.)
1155 0 : igrad = igrad + 1
1156 :
1157 : ! ! alternative to unrolling
1158 : ! monoG1=monoG2
1159 : ! monoSize1=monoSize2
1160 : ! monoDim1=monoDim2
1161 : ! monoFullDim1=monoFullDim2
1162 : END DO
1163 0 : 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 0 : 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 0 : REAL(dp), ALLOCATABLE, DIMENSION(:) :: xi
1181 :
1182 0 : IF (.NOT. module_initialized) CPABORT("module d3_poly not initialized")
1183 0 : my_npoly = 1
1184 0 : IF (PRESENT(npoly)) my_npoly = npoly
1185 0 : size_p = SIZE(p)/my_npoly
1186 0 : grad = grad_size3(size_p)
1187 0 : newSize = SIZE(pRes)/my_npoly
1188 0 : CPASSERT(newSize >= poly_size2(grad))
1189 0 : pRes = 0.0
1190 0 : ALLOCATE (xi(grad + 1))
1191 0 : CALL poly_p_eval3b(p, SIZE(p), x, pRes, SIZE(pRes), my_npoly, grad, xi)
1192 0 : DEALLOCATE (xi)
1193 0 : 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 0 : 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 0 : IF (.NOT. module_initialized) CPABORT("module d3_poly not initialized")
1223 0 : inSize = size_p/npoly
1224 0 : newSize = size_pRes/npoly
1225 0 : pRes(1:size_pRes) = 0.0
1226 0 : xi(1) = 1.0
1227 0 : DO i = 1, grad
1228 0 : xi(i + 1) = xi(i)*x
1229 : END DO
1230 : shiftRes = 0
1231 : pShift = 0
1232 0 : DO ipoly = 1, npoly
1233 0 : DO ii = 1, MIN(inSize, cached_dim3)
1234 0 : 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 0 : shiftRes = shiftRes + newSize
1237 0 : pShift = pShift + inSize
1238 : END DO
1239 0 : 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 0 : DO ipoly = 1, npoly
1244 : ii = ii0
1245 0 : 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 0 : DO i = igrad, 0, -1
1250 : !subG=igrad-i
1251 : !shiftRes=subG*(subG+3)/2+1
1252 0 : DO j = subG, 0, -1
1253 0 : IF (msize_p < ii) EXIT grad_do
1254 0 : pRes(shiftRes - j) = pRes(shiftRes - j) + p(ii)*xi(i + 1)
1255 0 : ii = ii + 1
1256 : END DO
1257 0 : shiftRes = shiftRes + subG + 2
1258 0 : subG = subG + 1
1259 : END DO
1260 : END DO grad_do
1261 0 : ii0 = ii0 + inSize
1262 0 : shiftRes_0 = shiftRes_0 + newSize
1263 0 : msize_p = msize_p + inSize
1264 : END DO
1265 : END IF
1266 0 : 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 0 : 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 0 : REAL(dp), ALLOCATABLE, DIMENSION(:) :: xi
1285 :
1286 0 : IF (.NOT. module_initialized) CPABORT("module d3_poly not initialized")
1287 0 : my_npoly = 1
1288 0 : IF (PRESENT(npoly)) my_npoly = npoly
1289 0 : size_p = SIZE(p)/my_npoly
1290 0 : newSize = SIZE(pRes)/my_npoly
1291 0 : grad = grad_size2(newSize)
1292 0 : CPASSERT(size_p >= poly_size3(grad))
1293 0 : CPASSERT(newSize == poly_size2(grad))
1294 0 : ALLOCATE (xi(grad + 1))
1295 0 : CALL poly_padd_uneval3b(p, SIZE(p), x, pRes, SIZE(pRes), my_npoly, grad, xi)
1296 0 : DEALLOCATE (xi)
1297 0 : 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 0 : 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 0 : IF (.NOT. module_initialized) CPABORT("module d3_poly not initialized")
1328 0 : inSize = size_p/npoly
1329 0 : newSize = size_pRes/npoly
1330 0 : upSize = (grad + 1)*(grad + 2)*(grad + 3)/6
1331 0 : xi(1) = 1.0
1332 0 : DO i = 1, grad
1333 0 : xi(i + 1) = xi(i)*x
1334 : END DO
1335 : shiftRes = 0
1336 : pShift = 0
1337 0 : DO ipoly = 1, npoly
1338 0 : DO ii = 1, MIN(upSize, cached_dim3)
1339 0 : p(pShift + ii) = p(pShift + ii) + pRes(shiftRes + a_reduce_idx3(ii))*xi(a_mono_exp3(1, ii) + 1)
1340 : END DO
1341 0 : shiftRes = shiftRes + newSize
1342 0 : pShift = pShift + inSize
1343 : END DO
1344 0 : 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 0 : DO ipoly = 1, npoly
1349 : ii = ii0
1350 0 : 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 0 : DO i = igrad, 0, -1
1355 : !subG=igrad-i
1356 : !shiftRes=subG*(subG+3)/2+1
1357 0 : DO j = subG, 0, -1
1358 0 : IF (msize_p < ii) EXIT grad_do
1359 0 : p(ii) = p(ii) + pRes(shiftRes - j)*xi(i + 1)
1360 0 : ii = ii + 1
1361 : END DO
1362 0 : shiftRes = shiftRes + subG + 2
1363 0 : subG = subG + 1
1364 : END DO
1365 : END DO grad_do
1366 0 : ii0 = ii0 + inSize
1367 0 : shiftRes_0 = shiftRes_0 + newSize
1368 0 : msize_p = msize_p + inSize
1369 : END DO
1370 : END IF
1371 0 : 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 0 : 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 0 : REAL(dp), ALLOCATABLE, DIMENSION(:) :: xi
1388 :
1389 0 : IF (.NOT. module_initialized) CPABORT("module d3_poly not initialized")
1390 0 : my_npoly = 1
1391 0 : IF (PRESENT(npoly)) my_npoly = npoly
1392 0 : size_p = SIZE(p)/my_npoly
1393 0 : grad = grad_size2(size_p)
1394 0 : newSize = SIZE(pRes)/my_npoly
1395 0 : pRes = 0.0_dp
1396 0 : CPASSERT(newSize >= poly_size1(grad))
1397 0 : ALLOCATE (xi(grad + 1))
1398 0 : CALL poly_p_eval2b(p, SIZE(p), x, pRes, SIZE(pRes), my_npoly, grad, xi)
1399 0 : DEALLOCATE (xi)
1400 0 : 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 0 : 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 0 : IF (.NOT. module_initialized) CPABORT("module d3_poly not initialized")
1429 0 : inSize = size_p/npoly
1430 0 : newSize = size_pRes/npoly
1431 0 : pRes(1:size_pRes) = 0.0_dp
1432 : !CPPreconditionNoFail(newSize>grad,cp_failure_level,routineP)
1433 0 : xi(1) = 1.0_dp
1434 0 : DO i = 1, grad
1435 0 : xi(i + 1) = xi(i)*x
1436 : END DO
1437 : shiftRes = 1
1438 : pShift = 0
1439 0 : DO ipoly = 1, npoly
1440 0 : DO ii = 1, MIN(inSize, cached_dim2)
1441 0 : 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 0 : shiftRes = shiftRes + newSize
1444 0 : pShift = pShift + inSize
1445 : END DO
1446 0 : IF (grad > max_grad2) THEN
1447 : ii0 = (max_grad2 + 1)*(max_grad2 + 2)/2 + 1
1448 : shiftRes = 1
1449 : msize_p = inSize
1450 0 : DO ipoly = 1, npoly
1451 : ii = ii0
1452 0 : grad_do2: DO igrad = max_grad2 + 1, grad
1453 : !ii=igrad*(igrad+1)/2+1
1454 : ij = shiftRes
1455 0 : DO i = igrad, 0, -1
1456 0 : IF (msize_p < ii) EXIT grad_do2
1457 : ! ij=igrad-i
1458 0 : pRes(ij) = pRes(ij) + p(ii)*xi(i + 1)
1459 0 : ii = ii + 1
1460 0 : ij = ij + 1
1461 : END DO
1462 : END DO grad_do2
1463 0 : msize_p = msize_p + inSize
1464 0 : shiftRes = shiftRes + newSize
1465 0 : ii0 = ii0 + inSize
1466 : END DO
1467 : END IF
1468 0 : 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 0 : 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 0 : REAL(dp), ALLOCATABLE, DIMENSION(:) :: xi
1487 :
1488 0 : IF (.NOT. module_initialized) CPABORT("module d3_poly not initialized")
1489 0 : my_npoly = 1
1490 0 : IF (PRESENT(npoly)) my_npoly = npoly
1491 0 : size_p = SIZE(p)/my_npoly
1492 0 : newSize = SIZE(pRes)/my_npoly
1493 0 : grad = grad_size1(newSize)
1494 0 : CPASSERT(size_p >= poly_size2(grad))
1495 0 : CPASSERT(newSize == poly_size1(grad))
1496 0 : ALLOCATE (xi(grad + 1))
1497 0 : CALL poly_padd_uneval2b(p, SIZE(p), x, pRes, SIZE(pRes), my_npoly, grad, xi)
1498 0 : DEALLOCATE (xi)
1499 0 : 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 0 : 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 0 : IF (.NOT. module_initialized) CPABORT("module d3_poly not initialized")
1529 0 : inSize = size_p/npoly
1530 0 : upSize = (grad + 1)*(grad + 2)/2
1531 0 : newSize = size_pRes/npoly
1532 : !CPPreconditionNoFail(newSize>grad,cp_failure_level,routineP)
1533 0 : xi(1) = 1.0_dp
1534 0 : DO i = 1, grad
1535 0 : xi(i + 1) = xi(i)*x
1536 : END DO
1537 : shiftRes = 1
1538 : pShift = 0
1539 0 : DO ipoly = 1, npoly
1540 0 : DO ii = 1, MIN(upSize, cached_dim2)
1541 0 : p(pShift + ii) = p(pShift + ii) + pRes(shiftRes + a_mono_exp2(2, ii))*xi(a_mono_exp2(1, ii) + 1)
1542 : END DO
1543 0 : shiftRes = shiftRes + newSize
1544 0 : pShift = pShift + inSize
1545 : END DO
1546 0 : IF (grad > max_grad2) THEN
1547 : ii0 = (max_grad2 + 1)*(max_grad2 + 2)/2 + 1
1548 : shiftRes = 1
1549 : msize_p = upSize
1550 0 : DO ipoly = 1, npoly
1551 : ii = ii0
1552 0 : grad_do2: DO igrad = max_grad2 + 1, grad
1553 : !ii=igrad*(igrad+1)/2+1
1554 : ij = shiftRes
1555 0 : DO i = igrad, 0, -1
1556 0 : IF (msize_p < ii) EXIT grad_do2
1557 : ! ij=igrad-i
1558 0 : p(ii) = p(ii) + pRes(ij)*xi(i + 1)
1559 0 : ii = ii + 1
1560 0 : ij = ij + 1
1561 : END DO
1562 : END DO grad_do2
1563 0 : msize_p = msize_p + inSize
1564 0 : shiftRes = shiftRes + newSize
1565 0 : ii0 = ii0 + inSize
1566 : END DO
1567 : END IF
1568 0 : 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 0 : 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 0 : my_npoly = 1
1587 0 : IF (PRESENT(npoly)) my_npoly = npoly
1588 0 : size_p = SIZE(p)/my_npoly
1589 0 : CPASSERT(SIZE(pRes) >= my_npoly)
1590 0 : pShift = 0
1591 0 : DO ipoly = 1, my_npoly
1592 : xx = 1.0_dp
1593 : vv = 0.0_dp
1594 0 : DO i = 1, size_p
1595 0 : vv = vv + p(pShift + i)*xx
1596 0 : xx = xx*x
1597 : END DO
1598 0 : pRes(ipoly) = vv
1599 0 : pShift = pShift + size_p
1600 : END DO
1601 0 : 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 0 : 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 0 : REAL(dp), ALLOCATABLE, DIMENSION(:) :: xi, yi
1621 :
1622 0 : IF (.NOT. module_initialized) CPABORT("module d3_poly not initialized")
1623 0 : my_npoly = 1
1624 0 : IF (PRESENT(npoly)) my_npoly = npoly
1625 0 : size_p = SIZE(p)/my_npoly
1626 0 : grad = grad_size2(size_p)
1627 0 : CPASSERT(SIZE(pRes) >= my_npoly)
1628 0 : ALLOCATE (xi(grad + 1), yi(grad + 1))
1629 0 : xi(1) = 1.0_dp
1630 0 : DO i = 1, grad
1631 0 : xi(i + 1) = xi(i)*x
1632 : END DO
1633 0 : yi(1) = 1.0_dp
1634 0 : DO i = 1, grad
1635 0 : yi(i + 1) = yi(i)*y
1636 : END DO
1637 : pShift = 0
1638 0 : DO ipoly = 1, my_npoly
1639 0 : v = 0.0_dp
1640 0 : DO ii = 1, MIN(size_p, cached_dim2)
1641 0 : v = v + p(pShift + ii)*xi(a_mono_exp2(1, ii) + 1)*yi(a_mono_exp2(2, ii) + 1)
1642 : END DO
1643 0 : pRes(ipoly) = v
1644 0 : pShift = pShift + size_p
1645 : END DO
1646 0 : IF (grad > max_grad2) THEN
1647 : pShift = (max_grad2 + 1)*(max_grad2 + 2)/2 + 1
1648 : msize_p = size_p
1649 0 : DO ipoly = 1, my_npoly
1650 : ii = pShift
1651 : v = 0.0_dp
1652 0 : grad_do4: DO igrad = max_grad2 + 1, grad
1653 : ! ii=igrad*(igrad+1)*(igrad+2)/6+1
1654 : j = 1
1655 0 : DO i = igrad, 0, -1
1656 0 : IF (msize_p < ii) EXIT grad_do4
1657 0 : v = v + p(ii)*xi(i + 1)*yi(j)
1658 0 : j = j + 1
1659 0 : ii = ii + 1
1660 : END DO
1661 : END DO grad_do4
1662 0 : pRes(ipoly) = pRes(ipoly) + v
1663 0 : pShift = pShift + size_p
1664 0 : msize_p = msize_p + size_p
1665 : END DO
1666 : END IF
1667 0 : 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 0 : 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 0 : REAL(dp), ALLOCATABLE, DIMENSION(:) :: xi, yi, zi
1688 :
1689 0 : IF (.NOT. module_initialized) CPABORT("module d3_poly not initialized")
1690 0 : my_npoly = 1
1691 0 : IF (PRESENT(npoly)) my_npoly = npoly
1692 0 : size_p = SIZE(p)/my_npoly
1693 0 : grad = grad_size3(size_p)
1694 0 : CPASSERT(SIZE(pRes) >= my_npoly)
1695 0 : ALLOCATE (xi(grad + 1), yi(grad + 1), zi(grad + 1))
1696 0 : xi(1) = 1.0_dp
1697 0 : DO i = 1, grad
1698 0 : xi(i + 1) = xi(i)*x
1699 : END DO
1700 0 : yi(1) = 1.0_dp
1701 0 : DO i = 1, grad
1702 0 : yi(i + 1) = yi(i)*y
1703 : END DO
1704 0 : zi(1) = 1.0_dp
1705 0 : DO i = 1, grad
1706 0 : zi(i + 1) = zi(i)*z
1707 : END DO
1708 : pShift = 0
1709 0 : DO ipoly = 1, my_npoly
1710 0 : v = 0.0_dp
1711 0 : 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 0 : *zi(a_mono_exp3(3, ii) + 1)
1714 : END DO
1715 0 : pRes(ipoly) = v
1716 0 : pShift = pShift + size_p
1717 : END DO
1718 0 : IF (grad > max_grad3) THEN
1719 0 : pShift = (max_grad3 + 1)*(max_grad3 + 2)*(max_grad3 + 3)/6 + 1
1720 0 : msize_p = size_p
1721 0 : DO ipoly = 1, my_npoly
1722 0 : ii = pShift
1723 : v = 0.0_dp
1724 0 : grad_do3: DO igrad = max_grad3 + 1, grad
1725 : ! ii=igrad*(igrad+1)*(igrad+2)/6+1
1726 0 : DO i = igrad, 0, -1
1727 0 : k = 1
1728 0 : DO j = igrad - i, 0, -1
1729 0 : ii = (ipoly - 1)*size_p + mono_index3(i, j, igrad - i - j) + 1
1730 0 : IF (msize_p < ii) EXIT grad_do3
1731 0 : v = v + p(ii)*xi(i + 1)*yi(j + 1)*zi(k)
1732 0 : k = k + 1
1733 0 : ii = ii + 1
1734 : END DO
1735 : END DO
1736 : END DO grad_do3
1737 0 : pRes(ipoly) = pRes(ipoly) + v
1738 0 : pShift = pShift + size_p
1739 0 : msize_p = msize_p + size_p
1740 : END DO
1741 : END IF
1742 0 : 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 0 : 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 0 : IF (.NOT. module_initialized) CPABORT("module d3_poly not initialized")
1762 0 : my_npoly = 1
1763 0 : IF (PRESENT(npoly)) my_npoly = npoly
1764 0 : my_sumUp = .FALSE.
1765 0 : IF (PRESENT(sumUp)) my_sumUp = sumUp
1766 0 : size_p = SIZE(p)/my_npoly
1767 0 : newSize = SIZE(pRes)/(3*my_npoly)
1768 0 : grad = grad_size3(size_p)
1769 0 : CPASSERT(newSize >= poly_size3(grad))
1770 0 : IF (.NOT. my_sumUp) pRes = 0
1771 0 : xDerivShift = 1
1772 0 : yDerivShift = my_npoly*newSize + 1
1773 0 : zDerivShift = 2*yDerivShift - 1
1774 0 : pShift = 0
1775 0 : DO ipoly = 1, my_npoly
1776 : ! split derivs to have less streams to follow (3 vs 5)?
1777 0 : DO ii = 1, MIN(cached_dim3, size_p)
1778 : pRes(xDerivShift + a_deriv_idx3(1, ii)) = pRes(xDerivShift + a_deriv_idx3(1, ii)) &
1779 0 : + p(pShift + ii)*a_mono_exp3(1, ii)
1780 : pRes(yDerivShift + a_deriv_idx3(2, ii)) = pRes(yDerivShift + a_deriv_idx3(2, ii)) &
1781 0 : + p(pShift + ii)*a_mono_exp3(2, ii)
1782 : pRes(zDerivShift + a_deriv_idx3(3, ii)) = pRes(zDerivShift + a_deriv_idx3(3, ii)) &
1783 0 : + p(pShift + ii)*a_mono_exp3(3, ii)
1784 : END DO
1785 0 : xDerivShift = xDerivShift + newSize
1786 0 : yDerivShift = yDerivShift + newSize
1787 0 : zDerivShift = zDerivShift + newSize
1788 0 : pShift = pShift + size_p
1789 : END DO
1790 0 : IF (grad > max_grad3) THEN
1791 0 : xDerivShift = 0
1792 0 : yDerivShift = my_npoly*newSize
1793 0 : zDerivShift = 2*yDerivShift - 1
1794 0 : msize_p = size_p
1795 0 : xDerivShift = max_grad3*(max_grad3 + 1)*(max_grad3 + 2)/6 + 1
1796 0 : pShift = xDerivShift + (max_grad3 + 1)*(max_grad3 + 2)/2
1797 0 : DO ipoly = 1, my_npoly
1798 : ii = pShift
1799 : ii2 = xDerivShift
1800 0 : grad_do5: DO igrad = max_grad3 + 1, grad
1801 : yShift = yDerivShift
1802 : zShift = zDerivShift
1803 0 : DO i = igrad, 0, -1
1804 0 : k = 0
1805 0 : DO j = igrad - i, 0, -1
1806 0 : IF (ii > msize_p) EXIT grad_do5
1807 : ! remove ifs?
1808 0 : IF (i > 0) pRes(ii2) = pRes(ii2) + p(ii)*i
1809 0 : IF (j > 0) pRes(yShift + ii2) = pRes(yShift + ii2) + p(ii)*j
1810 0 : IF (k > 0) pRes(zShift + ii2) = pRes(zShift + ii2) + k*p(ii)
1811 0 : ii = ii + 1
1812 0 : ii2 = ii2 + 1
1813 0 : k = k + 1
1814 : END DO
1815 0 : yShift = yShift - 1
1816 0 : zShift = zShift - 1
1817 : END DO
1818 0 : ii2 = ii2 - igrad - 1
1819 : END DO grad_do5
1820 0 : pShift = pShift + size_p
1821 0 : xDerivShift = xDerivShift + newSize
1822 0 : msize_p = msize_p + size_p
1823 : END DO
1824 : END IF
1825 0 : 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 0 : 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 0 : size_p = (grad + 1)*(grad + 2)*(grad + 3)/6
1843 0 : CPASSERT(SIZE(poly_cp2k) >= size_p)
1844 0 : CPASSERT(SIZE(poly_d3) >= size_p)
1845 0 : cp_ii = 0
1846 0 : sgrad2k = 0
1847 0 : sgrad3k = 0
1848 0 : DO k = 0, grad
1849 0 : shiftk = k + 1
1850 0 : sgrad2k = sgrad2k + k
1851 0 : sgrad3k = sgrad3k + sgrad2k
1852 0 : sgrad2 = sgrad2k
1853 0 : sgrad3 = sgrad3k
1854 0 : DO j = 0, grad - k
1855 0 : sgrad = j + k
1856 0 : mgrad2 = sgrad2
1857 0 : shiftj = mgrad2 + shiftk
1858 0 : mgrad = sgrad
1859 0 : shifti = shiftj + sgrad3
1860 0 : DO i = 0, grad - j - k
1861 0 : cp_ii = cp_ii + 1
1862 0 : poly_d3(shifti) = poly_cp2k(cp_ii)
1863 0 : mgrad = mgrad + 1
1864 0 : mgrad2 = mgrad2 + mgrad
1865 0 : shifti = shifti + mgrad2
1866 : END DO
1867 0 : sgrad2 = sgrad2 + sgrad + 1
1868 0 : sgrad3 = sgrad3 + sgrad2
1869 : END DO
1870 : END DO
1871 0 : IF (SIZE(poly_d3) >= size_p) THEN
1872 0 : poly_d3(size_p + 1:) = 0.0_dp
1873 : END IF
1874 0 : 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 0 : 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 0 : size_p = (grad + 1)*(grad + 2)*(grad + 3)/6
1892 0 : CPASSERT(SIZE(poly_cp2k) >= size_p)
1893 0 : CPASSERT(SIZE(poly_d3) >= size_p)
1894 0 : cp_ii = 0
1895 0 : sgrad2k = 0
1896 0 : sgrad3k = 0
1897 0 : DO k = 0, grad
1898 0 : shiftk = k + 1
1899 0 : sgrad2k = sgrad2k + k
1900 0 : sgrad3k = sgrad3k + sgrad2k
1901 0 : sgrad2 = sgrad2k
1902 0 : sgrad3 = sgrad3k
1903 0 : DO j = 0, grad - k
1904 0 : sgrad = j + k
1905 0 : mgrad2 = sgrad2
1906 0 : shiftj = mgrad2 + shiftk
1907 0 : mgrad = sgrad
1908 0 : shifti = shiftj + sgrad3
1909 0 : DO i = 0, grad - j - k
1910 0 : cp_ii = cp_ii + 1
1911 0 : poly_cp2k(cp_ii) = poly_d3(shifti)
1912 0 : mgrad = mgrad + 1
1913 0 : mgrad2 = mgrad2 + mgrad
1914 0 : shifti = shifti + mgrad2
1915 : END DO
1916 0 : sgrad2 = sgrad2 + sgrad + 1
1917 0 : sgrad3 = sgrad3 + sgrad2
1918 : END DO
1919 : END DO
1920 0 : IF (SIZE(poly_d3) >= size_p) THEN
1921 0 : poly_cp2k(size_p + 1:) = 0.0_dp
1922 : END IF
1923 0 : END SUBROUTINE
1924 :
1925 : END MODULE d3_poly
|
