ddx 0.6.0
Libary for domain-decomposition methods for polarizable continuum models
ddx_operators.f90
Go to the documentation of this file.
1
11
13module ddx_operators
14! Get the solver-module
15use ddx_solvers
16! Get lpb core-module
17use ddx_lpb_core
18! Get gradient module
19use ddx_gradients
20implicit none
21
22contains
23
25subroutine lx(params, constants, workspace, x, y, ddx_error)
26 implicit none
27 !! Inputs
28 type(ddx_params_type), intent(in) :: params
29 type(ddx_constants_type), intent(in) :: constants
30 real(dp), intent(in) :: x(constants % nbasis, params % nsph)
31 !! Temporaries
32 type(ddx_workspace_type), intent(inout) :: workspace
33 !! Output
34 real(dp), intent(out) :: y(constants % nbasis, params % nsph)
35 type(ddx_error_type), intent(inout) :: ddx_error
36 !! Local variables
37 integer :: isph, jsph, ij, l, ind, iproc
38
39 ! dummy operation on unused interface arguments
40 if (ddx_error % flag .eq. 0) continue
41
42 !! Initialize
43 y = zero
44!
45! incore code:
46!
47 if (params % matvecmem .eq. 1) then
48 !$omp parallel do default(none) shared(params,constants,x,y) &
49 !$omp private(isph,ij,jsph) schedule(dynamic)
50 do isph = 1, params % nsph
51 do ij = constants % inl(isph), constants % inl(isph + 1) - 1
52 jsph = constants % nl(ij)
53 call dgemv('n', constants % nbasis, constants % nbasis, one, &
54 & constants % l(:,:,ij), constants % nbasis, x(:,jsph), 1, &
55 & one, y(:,isph), 1)
56 end do
57 end do
58 else
59 !$omp parallel do default(none) shared(params,constants,workspace,x,y) &
60 !$omp private(isph,iproc) schedule(dynamic)
61 do isph = 1, params % nsph
62 iproc = omp_get_thread_num() + 1
63 ! Compute NEGATIVE action of off-digonal blocks
64 call calcv(params, constants, isph, workspace % tmp_pot(:, iproc), x, &
65 & workspace % tmp_work(:, iproc))
66 call ddintegrate(1, constants % nbasis, params % ngrid, &
67 & constants % vwgrid, constants % vgrid_nbasis, &
68 & workspace % tmp_pot(:, iproc), y(:, isph))
69 ! now, fix the sign.
70 y(:, isph) = - y(:, isph)
71 end do
72 end if
73!
74! if required, add the diagonal.
75!
76 if (constants % dodiag) then
77 do isph = 1, params % nsph
78 do l = 0, params % lmax
79 ind = l*l + l + 1
80 y(ind-l:ind+l, isph) = y(ind-l:ind+l, isph) + &
81 x(ind-l:ind+l, isph) / (constants % vscales(ind)**2)
82 end do
83 end do
84 end if
85end subroutine lx
86
88subroutine lstarx(params, constants, workspace, x, y, ddx_error)
89 implicit none
90 !! Inputs
91 type(ddx_params_type), intent(in) :: params
92 type(ddx_constants_type), intent(in) :: constants
93 real(dp), intent(in) :: x(constants % nbasis, params % nsph)
94 !! Temporaries
95 type(ddx_workspace_type), intent(inout) :: workspace
96 !! Output
97 real(dp), intent(out) :: y(constants % nbasis, params % nsph)
98 type(ddx_error_type), intent(inout) :: ddx_error
99 !! Local variables
100 integer :: isph, jsph, ij, indmat, igrid, l, ind, iproc
101
102 ! dummy operation on unused interface arguments
103 if (ddx_error % flag .eq. 0) continue
104
105 y = zero
106 if (params % matvecmem .eq. 1) then
107 !$omp parallel do default(none) shared(params,constants,x,y) &
108 !$omp private(isph,ij,jsph,indmat) schedule(dynamic)
109 do isph = 1, params % nsph
110 do ij = constants % inl(isph), constants % inl(isph + 1) - 1
111 jsph = constants % nl(ij)
112 indmat = constants % itrnl(ij)
113 call dgemv('t', constants % nbasis, constants % nbasis, one, &
114 & constants % l(:,:,indmat), constants % nbasis, x(:,jsph), 1, &
115 & one, y(:,isph), 1)
116 end do
117 end do
118 else
119 ! Expand x over spherical harmonics
120 !$omp parallel do default(none) shared(params,constants,workspace,x) &
121 !$omp private(isph,igrid) schedule(dynamic)
122 do isph = 1, params % nsph
123 do igrid = 1, params % ngrid
124 workspace % tmp_grid(igrid, isph) = dot_product(x(:, isph), &
125 & constants % vgrid(:constants % nbasis, igrid))
126 end do
127 end do
128 ! Compute (negative) action
129 !$omp parallel do default(none) shared(params,constants,workspace,x,y) &
130 !$omp private(isph,iproc) schedule(dynamic)
131 do isph = 1, params % nsph
132 iproc = omp_get_thread_num() + 1
133 call adjrhs(params, constants, isph, workspace % tmp_grid, &
134 & y(:, isph), workspace % tmp_work(:, iproc))
135 ! fix the sign
136 y(:, isph) = - y(:, isph)
137 end do
138 end if
139 if (constants % dodiag) then
140 ! Loop over harmonics
141 do isph = 1, params % nsph
142 do l = 0, params % lmax
143 ind = l*l + l + 1
144 y(ind-l:ind+l, isph) = y(ind-l:ind+l, isph) + &
145 & x(ind-l:ind+l, isph) / (constants % vscales(ind)**2)
146 end do
147 end do
148 end if
149end subroutine lstarx
150
154subroutine ldm1x(params, constants, workspace, x, y, ddx_error)
155 implicit none
156 !! Inputs
157 type(ddx_params_type), intent(in) :: params
158 type(ddx_constants_type), intent(in) :: constants
159 real(dp), intent(in) :: x(constants % nbasis, params % nsph)
160 !! Temporaries
161 type(ddx_workspace_type), intent(inout) :: workspace
162 !! Output
163 real(dp), intent(out) :: y(constants % nbasis, params % nsph)
164 type(ddx_error_type), intent(inout) :: ddx_error
165 !! Local variables
166 integer :: isph, l, ind
167
168 ! dummy operation on unused interface arguments
169 if ((ddx_error % flag .eq. 0) .or. &
170 & (allocated(workspace % tmp_pot))) continue
171
172 !! Loop over harmonics
173 !$omp parallel do default(none) shared(params,constants,x,y) &
174 !$omp private(isph,l,ind) schedule(dynamic)
175 do isph = 1, params % nsph
176 do l = 0, params % lmax
177 ind = l*l + l + 1
178 y(ind-l:ind+l, isph) = x(ind-l:ind+l, isph) &
179 & *(constants % vscales(ind)**2)
180 end do
181 end do
182end subroutine ldm1x
183
185subroutine dx(params, constants, workspace, x, y, ddx_error)
186 implicit none
187 !! Inputs
188 type(ddx_params_type), intent(in) :: params
189 type(ddx_constants_type), intent(in) :: constants
190 type(ddx_workspace_type), intent(inout) :: workspace
191 real(dp), intent(in) :: x(constants % nbasis, params % nsph)
192 !! Output
193 real(dp), intent(out) :: y(constants % nbasis, params % nsph)
194 type(ddx_error_type), intent(inout) :: ddx_error
195 !! Local parameter
196 integer :: do_diag
197 !! initialize do_diag
198 do_diag = 0
199 if (constants % dodiag) do_diag = 1
200 !! Select implementation
201 if (params % fmm .eq. 0) then
202 call dx_dense(params, constants, workspace, do_diag, x, y, ddx_error)
203 else
204 call dx_fmm(params, constants, workspace, do_diag, x, y, ddx_error)
205 end if
206end subroutine dx
207
209subroutine dx_dense(params, constants, workspace, do_diag, x, y, ddx_error)
210 implicit none
211 !! Inputs
212 type(ddx_params_type), intent(in) :: params
213 type(ddx_constants_type), intent(in) :: constants
214 type(ddx_workspace_type), intent(inout) :: workspace
215 integer, intent(in) :: do_diag
216 real(dp), intent(in) :: x(constants % nbasis, params % nsph)
217 type(ddx_error_type), intent(inout) :: ddx_error
218 !! Output
219 real(dp), intent(out) :: y(constants % nbasis, params % nsph)
220 !! Local variables
221 real(dp) :: c(3), vij(3), sij(3)
222 real(dp) :: vvij, tij, tt, f, rho, ctheta, stheta, cphi, sphi
223 integer :: its, isph, jsph, l, m, ind
224 real(dp), external :: dnrm2
225
226 ! dummy operation on unused interface arguments
227 if (ddx_error % flag .eq. 0) continue
228
229 y = zero
230 do isph = 1, params % nsph
231 ! compute the "potential" from the other spheres
232 ! at the exposed lebedv points of the i-th sphere
233 workspace % tmp_grid(:, 1) = zero
234 do its = 1, params % ngrid
235 if (constants % ui(its,isph).gt.zero) then
236 c = params % csph(:,isph) + params % rsph(isph)* &
237 & constants % cgrid(:,its)
238 do jsph = 1, params % nsph
239 if (jsph.ne.isph) then
240 ! build the geometrical variables
241 vij = c - params % csph(:,jsph)
242 !vvij = sqrt(dot_product(vij,vij))
243 vvij = dnrm2(3, vij, 1)
244 tij = vvij / params % rsph(jsph)
245 sij = vij/vvij
246 ! build the local basis
247 call ylmbas(sij, rho, ctheta, stheta, cphi, sphi, &
248 & params % lmax, constants % vscales, &
249 & workspace % tmp_vylm, workspace % tmp_vplm, &
250 & workspace % tmp_vcos, workspace % tmp_vsin)
251 ! with all the required stuff, finally compute
252 ! the "potential" at the point
253 tt = one/tij
254 do l = 0, params % lmax
255 ind = l*l + l + 1
256 f = fourpi*dble(l)/(two*dble(l) + one)*tt
257 do m = -l, l
258 workspace % tmp_grid(its, 1) = &
259 & workspace % tmp_grid(its, 1) + &
260 & f*x(ind + m,jsph) * &
261 & workspace % tmp_vylm(ind + m, 1)
262 end do
263 tt = tt/tij
264 end do
265 else if (do_diag .eq. 1) then
266 ! add the diagonal contribution
267 do l = 0, params % lmax
268 ind = l*l + l + 1
269 f = (two*dble(l) + one)/fourpi
270 do m = -l, l
271 workspace % tmp_grid(its, 1) = &
272 & workspace % tmp_grid(its, 1) - &
273 & pt5*x(ind + m,isph) * &
274 & constants % vgrid(ind + m,its)/f
275 end do
276 end do
277 end if
278 end do
279 workspace % tmp_grid(its, 1) = constants % ui(its, isph) * &
280 & workspace % tmp_grid(its, 1)
281 end if
282 end do
283 ! now integrate the potential to get its modal representation
284 call ddintegrate(1, constants % nbasis, params % ngrid, &
285 & constants % vwgrid, constants % vgrid_nbasis, &
286 & workspace % tmp_grid, y(:,isph))
287 end do
288end subroutine dx_dense
289
291subroutine dx_fmm(params, constants, workspace, do_diag, x, y, ddx_error)
292 implicit none
293 !! Inputs
294 type(ddx_params_type), intent(in) :: params
295 type(ddx_constants_type), intent(in) :: constants
296 type(ddx_workspace_type), intent(inout) :: workspace
297 integer, intent(in) :: do_diag
298 real(dp), intent(in) :: x(constants % nbasis, params % nsph)
299 type(ddx_error_type), intent(inout) :: ddx_error
300 !! Output
301 real(dp), intent(out) :: y(constants % nbasis, params % nsph)
302 !! Local variables
303 integer :: isph, inode, l, indl, indl1
304
305 ! dummy operation on unused interface arguments
306 if (ddx_error % flag .eq. 0) continue
307
308 !! Scale input harmonics at first
309 workspace % tmp_sph(1, :) = zero
310 indl = 2
311 do l = 1, params % lmax
312 indl1 = (l+1)**2
313 workspace % tmp_sph(indl:indl1, :) = l * x(indl:indl1, :)
314 indl = indl1 + 1
315 end do
316 ! Load input harmonics into tree data
317 if(params % lmax .lt. params % pm) then
318 do isph = 1, params % nsph
319 inode = constants % snode(isph)
320 workspace % tmp_node_m(1:constants % nbasis, inode) = &
321 & workspace % tmp_sph(:, isph)
322 workspace % tmp_node_m(constants % nbasis+1:, inode) = zero
323 end do
324 else
325 indl = (params % pm+1)**2
326 do isph = 1, params % nsph
327 inode = constants % snode(isph)
328 workspace % tmp_node_m(:, inode) = workspace % tmp_sph(:indl, isph)
329 end do
330 end if
331 ! Do FMM operations
332 call tree_m2m_rotation(params, constants, workspace % tmp_node_m)
333 call tree_m2l_rotation(params, constants, workspace % tmp_node_m, &
334 & workspace % tmp_node_l)
335 call tree_l2l_rotation(params, constants, workspace % tmp_node_l)
336 call tree_l2p(params, constants, one, workspace % tmp_node_l, zero, &
337 & workspace % tmp_grid, workspace % tmp_sph_l)
338 call tree_m2p(params, constants, params % lmax, one, workspace % tmp_sph, one, &
339 & workspace % tmp_grid)
340 ! Apply diagonal contribution if needed
341 if(do_diag .eq. 1) then
342 call dgemm('T', 'N', params % ngrid, params % nsph, &
343 & constants % nbasis, -pt5, constants % vgrid2, &
344 & constants % vgrid_nbasis, x, constants % nbasis, one, &
345 & workspace % tmp_grid, params % ngrid)
346 end if
347 ! Multiply by characteristic function
348 workspace % tmp_grid = workspace % tmp_grid * constants % ui
349 ! now integrate the potential to get its modal representation
350 ! output y is overwritten here
351 call dgemm('N', 'N', constants % nbasis, params % nsph, params % ngrid, &
352 & one, constants % vwgrid, constants % vgrid_nbasis, workspace % tmp_grid, &
353 & params % ngrid, zero, y, constants % nbasis)
354end subroutine dx_fmm
355
357subroutine dstarx(params, constants, workspace, x, y, ddx_error)
358 implicit none
359 !! Inputs
360 type(ddx_params_type), intent(in) :: params
361 type(ddx_constants_type), intent(in) :: constants
362 real(dp), intent(in) :: x(constants % nbasis, params % nsph)
363 !! Temporaries
364 type(ddx_workspace_type), intent(inout) :: workspace
365 !! Output
366 real(dp), intent(out) :: y(constants % nbasis, params % nsph)
367 type(ddx_error_type), intent(inout) :: ddx_error
368 !! Local variables
369 integer :: do_diag
370 !! initialize do_diag
371 do_diag = 0
372 if (constants % dodiag) do_diag = 1
373 !! Select implementation
374 if (params % fmm .eq. 0) then
375 call dstarx_dense(params, constants, workspace, do_diag, x, y, ddx_error)
376 else
377 call dstarx_fmm(params, constants, workspace, do_diag, x, y, ddx_error)
378 end if
379end subroutine dstarx
380
382subroutine dstarx_dense(params, constants, workspace, do_diag, x, y, ddx_error)
383 implicit none
384 ! Inputs
385 type(ddx_params_type), intent(in) :: params
386 type(ddx_constants_type), intent(in) :: constants
387 integer, intent(in) :: do_diag
388 real(dp), intent(in) :: x(constants % nbasis, params % nsph)
389 ! Temporaries
390 type(ddx_workspace_type), intent(inout) :: workspace
391 type(ddx_error_type), intent(inout) :: ddx_error
392 ! Output
393 real(dp), intent(out) :: y(constants % nbasis, params % nsph)
394 ! Local variables
395 real(dp) :: vji(3), sji(3)
396 real(dp) :: vvji, tji, tt, f, rho, ctheta, stheta, cphi, sphi
397 integer :: its, isph, jsph, l, m, ind, iproc
398 real(dp), external :: dnrm2
399
400 ! dummy operation on unused interface arguments
401 if (ddx_error % flag .eq. 0)continue
402
403 y = zero
404 !$omp parallel do default(none) shared(do_diag,params,constants, &
405 !$omp workspace,x,y) private(isph,jsph,its,vji,vvji,tji,sji,rho, &
406 !$omp ctheta,stheta,cphi,sphi,tt,l,ind,f,m,iproc) schedule(dynamic)
407 do isph = 1, params % nsph
408 iproc = omp_get_thread_num() + 1
409 do jsph = 1, params % nsph
410 if (jsph.ne.isph) then
411 do its = 1, params % ngrid
412 if (constants % ui(its,jsph).gt.zero) then
413 ! build the geometrical variables
414 vji = params % csph(:,jsph) + params % rsph(jsph) * &
415 & constants % cgrid(:,its) - params % csph(:,isph)
416 !vvji = sqrt(dot_product(vji,vji))
417 vvji = dnrm2(3, vji, 1)
418 tji = vvji/params % rsph(isph)
419 sji = vji/vvji
420 ! build the local basis
421 call ylmbas(sji, rho, ctheta, stheta, cphi, sphi, &
422 & params % lmax, constants % vscales, &
423 & workspace % tmp_vylm(:(params % lmax + 1)**2, iproc), &
424 & workspace % tmp_vplm(:(params % lmax + 1)**2, iproc), &
425 & workspace % tmp_vcos(:(params % lmax + 1), iproc), &
426 & workspace % tmp_vsin(:(params % lmax + 1), iproc))
427 tt = constants % ui(its,jsph) &
428 & *dot_product(constants % vwgrid(:constants % nbasis, its), &
429 & x(:,jsph))/tji
430 do l = 0, params % lmax
431 ind = l*l + l + 1
432 f = dble(l)*tt/ constants % vscales(ind)**2
433 do m = -l, l
434 y(ind+m,isph) = y(ind+m,isph) + &
435 & f*workspace % tmp_vylm(ind+m, iproc)
436 end do
437 tt = tt/tji
438 end do
439 end if
440 end do
441 else if (do_diag .eq. 1) then
442 do its = 1, params % ngrid
443 f = pt5*constants % ui(its,jsph) &
444 & *dot_product(constants % vwgrid(:constants % nbasis, its), &
445 & x(:,jsph))
446 do l = 0, params % lmax
447 ind = l*l + l + 1
448 y(ind-l:ind+l,isph) = y(ind-l:ind+l,isph) - &
449 & f*constants % vgrid(ind-l:ind+l,its)/ &
450 & constants % vscales(ind)**2
451 end do
452 end do
453 end if
454 end do
455 end do
456end subroutine dstarx_dense
457
459subroutine dstarx_fmm(params, constants, workspace, do_diag, x, y, ddx_error)
460 implicit none
461 ! Inputs
462 type(ddx_params_type), intent(in) :: params
463 type(ddx_constants_type), intent(in) :: constants
464 integer, intent(in) :: do_diag
465 real(dp), intent(in) :: x(constants % nbasis, params % nsph)
466 ! Temporaries
467 type(ddx_workspace_type), intent(inout) :: workspace
468 type(ddx_error_type), intent(inout) :: ddx_error
469 ! Output
470 real(dp), intent(out) :: y(constants % nbasis, params % nsph)
471 ! Local variables
472 integer :: isph, inode, l, indl, indl1
473
474 ! dummy operation on unused interface arguments
475 if (ddx_error % flag .eq. 0) continue
476
477 ! Adjoint integration
478 call dgemm('T', 'N', params % ngrid, params % nsph, constants % nbasis, &
479 & one, constants % vwgrid, constants % vgrid_nbasis, x, &
480 & constants % nbasis, zero, workspace % tmp_grid, params % ngrid)
481 ! Multiply by characteristic function
482 workspace % tmp_grid = workspace % tmp_grid * constants % ui
483 ! Do FMM operations adjointly
484 call tree_m2p_adj(params, constants, params % lmax, one, workspace % tmp_grid, &
485 & zero, y)
486 call tree_l2p_adj(params, constants, one, workspace % tmp_grid, zero, &
487 & workspace % tmp_node_l, workspace % tmp_sph_l)
488 call tree_l2l_rotation_adj(params, constants, workspace % tmp_node_l)
489 call tree_m2l_rotation_adj(params, constants, workspace % tmp_node_l, &
490 & workspace % tmp_node_m)
491 call tree_m2m_rotation_adj(params, constants, workspace % tmp_node_m)
492 ! Adjointly move tree multipole harmonics into output
493 if(params % lmax .lt. params % pm) then
494 do isph = 1, params % nsph
495 inode = constants % snode(isph)
496 y(:, isph) = y(:, isph) + &
497 & workspace % tmp_node_m(1:constants % nbasis, inode)
498 end do
499 else
500 indl = (params % pm+1)**2
501 do isph = 1, params % nsph
502 inode = constants % snode(isph)
503 y(1:indl, isph) = y(1:indl, isph) + workspace % tmp_node_m(:, inode)
504 end do
505 end if
506 ! Scale output harmonics at last
507 y(1, :) = zero
508 indl = 2
509 do l = 1, params % lmax
510 indl1 = (l+1)**2
511 y(indl:indl1, :) = l * y(indl:indl1, :)
512 indl = indl1 + 1
513 end do
514 ! Apply diagonal contribution if needed
515 if(do_diag .eq. 1) then
516 call dgemm('N', 'N', constants % nbasis, params % nsph, &
517 & params % ngrid, -pt5, constants % vgrid2, &
518 & constants % vgrid_nbasis, workspace % tmp_grid, params % ngrid, &
519 & one, y, constants % nbasis)
520 end if
521end subroutine dstarx_fmm
522
527subroutine repsx(params, constants, workspace, x, y, ddx_error)
528 implicit none
529 !! Inputs
530 type(ddx_params_type), intent(in) :: params
531 type(ddx_constants_type), intent(in) :: constants
532 real(dp), intent(in) :: x(constants % nbasis, params % nsph)
533 !! Temporaries
534 type(ddx_workspace_type), intent(inout) :: workspace
535 !! Output
536 real(dp), intent(out) :: y(constants % nbasis, params % nsph)
537 type(ddx_error_type), intent(inout) :: ddx_error
538 !! Local variables
539 real(dp) :: fac
540 !! Output `y` is cleaned here
541 call dx(params, constants, workspace, x, y, ddx_error)
542 y = - y
543 if (constants % dodiag) then
544 ! Apply diagonal
545 fac = twopi * (params % eps + one) / (params % eps - one)
546 y = fac*x + y
547 end if
548end subroutine repsx
549
555subroutine repsstarx(params, constants, workspace, x, y, ddx_error)
556 implicit none
557 ! Inputs
558 type(ddx_params_type), intent(in) :: params
559 type(ddx_constants_type), intent(in) :: constants
560 real(dp), intent(in) :: x(constants % nbasis, params % nsph)
561 ! Temporaries
562 type(ddx_workspace_type), intent(inout) :: workspace
563 ! Output
564 real(dp), intent(out) :: y(constants % nbasis, params % nsph)
565 type(ddx_error_type), intent(inout) :: ddx_error
566 ! Local variables
567 real(dp) :: fac
568 ! Output `y` is cleaned here
569 call dstarx(params, constants, workspace, x, y, ddx_error)
570 y = - y
571 ! Apply diagonal
572 if (constants % dodiag) then
573 fac = twopi * (params % eps + one) / (params % eps - one)
574 y = fac*x + y
575 end if
576end subroutine repsstarx
577
587subroutine rinfx(params, constants, workspace, x, y, ddx_error)
588 implicit none
589 ! Inputs
590 type(ddx_params_type), intent(in) :: params
591 type(ddx_constants_type), intent(in) :: constants
592 real(dp), intent(in) :: x(constants % nbasis, params % nsph)
593 ! Temporaries
594 type(ddx_workspace_type), intent(inout) :: workspace
595 ! Output
596 real(dp), intent(out) :: y(constants % nbasis, params % nsph)
597 type(ddx_error_type), intent(inout) :: ddx_error
598 !! note do_diag hardcoded to 1.
599 !! Select implementation
600 if (params % fmm .eq. 0) then
601 call dx_dense(params, constants, workspace, 1, x, y, ddx_error)
602 else
603 call dx_fmm(params, constants, workspace, 1, x, y, ddx_error)
604 end if
605 ! Apply diagonal
606 y = twopi*x - y
607end subroutine rinfx
608
618subroutine rstarinfx(params, constants, workspace, x, y, ddx_error)
619 implicit none
620 ! Inputs
621 type(ddx_params_type), intent(in) :: params
622 type(ddx_constants_type), intent(in) :: constants
623 real(dp), intent(in) :: x(constants % nbasis, params % nsph)
624 ! Temporaries
625 type(ddx_workspace_type), intent(inout) :: workspace
626 ! Output
627 real(dp), intent(out) :: y(constants % nbasis, params % nsph)
628 type(ddx_error_type), intent(inout) :: ddx_error
629 ! Output `y` is cleaned here
630 call dstarx(params, constants, workspace, x, y, ddx_error)
631 ! Apply diagonal
632 y = twopi*x - y
633end subroutine rstarinfx
634
637subroutine prec_repsx(params, constants, workspace, x, y, ddx_error)
638 implicit none
639 ! Inputs
640 type(ddx_params_type), intent(in) :: params
641 type(ddx_constants_type), intent(in) :: constants
642 real(dp), intent(in) :: x(constants % nbasis, params % nsph)
643 ! Temporaries
644 type(ddx_workspace_type), intent(inout) :: workspace
645 ! Output
646 real(dp), intent(out) :: y(constants % nbasis, params % nsph)
647 type(ddx_error_type), intent(inout) :: ddx_error
648 integer :: isph
649
650 ! dummy operation on unused interface arguments
651 if ((ddx_error % flag .eq. 0) .or. &
652 & (allocated(workspace % tmp_pot))) continue
653
654 ! simply do a matrix-vector product with the transposed preconditioner
655 !$omp parallel do default(shared) schedule(static,1) &
656 !$omp private(isph)
657 do isph = 1, params % nsph
658 call dgemm('N', 'N', constants % nbasis, 1, constants % nbasis, one, &
659 & constants % rx_prc(:, :, isph), constants % nbasis, x(:, isph), &
660 & constants % nbasis, zero, y(:, isph), constants % nbasis)
661 end do
662end subroutine prec_repsx
663
665subroutine prec_repsstarx(params, constants, workspace, x, y, ddx_error)
666 implicit none
667 ! Inputs
668 type(ddx_params_type), intent(in) :: params
669 type(ddx_constants_type), intent(in) :: constants
670 real(dp), intent(in) :: x(constants % nbasis, params % nsph)
671 ! Temporaries
672 type(ddx_workspace_type), intent(inout) :: workspace
673 ! Output
674 real(dp), intent(out) :: y(constants % nbasis, params % nsph)
675 type(ddx_error_type), intent(inout) :: ddx_error
676 ! Local variables
677 integer :: isph
678
679 ! dummy operation on unused interface arguments
680 if ((ddx_error % flag .eq. 0) .or. &
681 & (allocated(workspace % tmp_pot))) continue
682
683 ! simply do a matrix-vector product with the transposed preconditioner
684 !$omp parallel do default(shared) schedule(static,1) &
685 !$omp private(isph)
686 do isph = 1, params % nsph
687 call dgemm('T', 'N', constants % nbasis, 1, constants % nbasis, one, &
688 & constants % rx_prc(:, :, isph), constants % nbasis, x(:, isph), &
689 & constants % nbasis, zero, y(:, isph), constants % nbasis)
690 end do
691end subroutine prec_repsstarx
692
694subroutine bstarx(params, constants, workspace, x, y, ddx_error)
695 ! Inputs
696 type(ddx_params_type), intent(in) :: params
697 type(ddx_constants_type), intent(in) :: constants
698 type(ddx_workspace_type), intent(inout) :: workspace
699 real(dp), intent(in) :: x(constants % nbasis, params % nsph)
700 real(dp), intent(out) :: y(constants % nbasis, params % nsph)
701 type(ddx_error_type), intent(inout) :: ddx_error
702 ! Local variables
703 integer :: isph, jsph, ij, indmat, iproc
704
705 ! dummy operation on unused interface arguments
706 if (ddx_error % flag .eq. 0) continue
707
708 y = zero
709 if (params % matvecmem .eq. 1) then
710 !$omp parallel do default(none) shared(params,constants,x,y) &
711 !$omp private(isph,ij,jsph,indmat) schedule(dynamic)
712 do isph = 1, params % nsph
713 do ij = constants % inl(isph), constants % inl(isph + 1) - 1
714 jsph = constants % nl(ij)
715 indmat = constants % itrnl(ij)
716 call dgemv('t', constants % nbasis, constants % nbasis, one, &
717 & constants % b(:,:,indmat), constants % nbasis, x(:,jsph), 1, &
718 & one, y(:,isph), 1)
719 end do
720 end do
721 else
722 !$omp parallel do default(none) shared(params,constants,workspace,x,y) &
723 !$omp private(isph) schedule(static,1)
724 do isph = 1, params % nsph
725 call dgemv('t', constants % nbasis, params % ngrid, one, constants % vgrid, &
726 & constants % vgrid_nbasis, x(:, isph), 1, zero, &
727 & workspace % tmp_grid(:, isph), 1)
728 end do
729 !$omp parallel do default(none) shared(params,constants,workspace,x,y) &
730 !$omp private(isph,iproc) schedule(dynamic)
731 do isph = 1, params % nsph
732 iproc = omp_get_thread_num() + 1
733 call adjrhs_lpb(params, constants, isph, workspace % tmp_grid, &
734 & y(:, isph), workspace % tmp_vylm(:, iproc), workspace % tmp_vplm(:, iproc), &
735 & workspace % tmp_vcos(:, iproc), workspace % tmp_vsin(:, iproc), &
736 & workspace % tmp_bessel(:, iproc))
737 y(:,isph) = - y(:,isph)
738 end do
739 end if
740
741 ! add the diagonal if required
742 if (constants % dodiag) y = y + x
743end subroutine bstarx
744
746subroutine bx(params, constants, workspace, x, y, ddx_error)
747 implicit none
748 type(ddx_params_type), intent(in) :: params
749 type(ddx_constants_type), intent(in) :: constants
750 type(ddx_workspace_type), intent(inout) :: workspace
751 real(dp), dimension(constants % nbasis, params % nsph), intent(in) :: x
752 real(dp), dimension(constants % nbasis, params % nsph), intent(out) :: y
753 type(ddx_error_type), intent(inout) :: ddx_error
754 integer :: isph, jsph, ij, iproc
755
756 ! dummy operation on unused interface arguments
757 if (ddx_error % flag .eq. 0) continue
758
759 y = zero
760 if (params % matvecmem .eq. 1) then
761 !$omp parallel do default(none) shared(params,constants,x,y) &
762 !$omp private(isph,ij,jsph) schedule(dynamic)
763 do isph = 1, params % nsph
764 do ij = constants % inl(isph), constants % inl(isph + 1) - 1
765 jsph = constants % nl(ij)
766 call dgemv('n', constants % nbasis, constants % nbasis, one, &
767 & constants % b(:,:,ij), constants % nbasis, x(:,jsph), 1, &
768 & one, y(:,isph), 1)
769 end do
770 end do
771 else
772 !$omp parallel do default(none) shared(params,constants,workspace,x,y) &
773 !$omp private(isph,iproc) schedule(dynamic)
774 do isph = 1, params % nsph
775 iproc = omp_get_thread_num() + 1
776 call calcv2_lpb(params, constants, isph, workspace % tmp_pot(:, iproc), x)
777 call ddintegrate(1, constants % nbasis, params % ngrid, constants % vwgrid, &
778 & constants % vgrid_nbasis, workspace % tmp_pot(:, iproc), y(:,isph))
779 y(:,isph) = - y(:,isph)
780 end do
781 end if
782 if (constants % dodiag) y = y + x
783end subroutine bx
784
786subroutine tstarx(params, constants, workspace, x, y, ddx_error)
787 implicit none
788 type(ddx_params_type), intent(in) :: params
789 type(ddx_constants_type), intent(in) :: constants
790 type(ddx_workspace_type), intent(inout) :: workspace
791 real(dp), dimension(constants % nbasis, params % nsph, 2), intent(in) :: x
792 real(dp), dimension(constants % nbasis, params % nsph, 2), intent(out) :: y
793 type(ddx_error_type), intent(inout) :: ddx_error
794
795 real(dp), dimension(constants % nbasis, params % nsph, 2) :: temp_vector
796 y = zero
797 temp_vector = zero
798
799 ! Compute AXr
800 ! call LstarXr
801 call lstarx(params, constants, workspace, x(:,:,1), temp_vector(:,:,1), &
802 & ddx_error)
803 ! Remove the scaling factor
804 call convert_ddcosmo(params, constants, 1, temp_vector(:,:,1))
805 y(:,:,1) = y(:,:,1) + temp_vector(:,:,1)
806 ! Compute BXe
807 call bstarx(params, constants, workspace, x(:,:,2), temp_vector(:,:,2), &
808 & ddx_error)
809 y(:,:,2) = y(:,:,2) + temp_vector(:,:,2)
810
811 ! Reset temp_vector to zero
812 temp_vector = zero
813 ! Call CX
814 call cstarx(params, constants, workspace, x, temp_vector, ddx_error)
815 y = y + temp_vector
816end subroutine tstarx
817
820subroutine bx_prec(params, constants, workspace, x, y, ddx_error)
821 implicit none
822 type(ddx_params_type), intent(in) :: params
823 type(ddx_constants_type), intent(in) :: constants
824 type(ddx_workspace_type), intent(inout) :: workspace
825 real(dp), dimension(constants % nbasis, params % nsph), intent(in) :: x
826 real(dp), dimension(constants % nbasis, params % nsph), intent(out) :: y
827 type(ddx_error_type), intent(inout) :: ddx_error
828
829 ! dummy operation on unused interface arguments
830 if ((ddx_error % flag .eq. 0) .or. &
831 & (allocated(workspace % tmp_pot))) continue
832
833 y = x
834end subroutine bx_prec
835
837subroutine prec_tstarx(params, constants, workspace, x, y, ddx_error)
838 implicit none
839 type(ddx_params_type), intent(in) :: params
840 type(ddx_constants_type), intent(in) :: constants
841 type(ddx_workspace_type), intent(inout) :: workspace
842 real(dp), intent(in) :: x(constants % nbasis, params % nsph, 2)
843 real(dp), intent(inout) :: y(constants % nbasis, params % nsph, 2)
844 type(ddx_error_type), intent(inout) :: ddx_error
845 integer :: n_iter
846 real(dp), dimension(params % maxiter) :: x_rel_diff
847 real(dp) :: start_time
848
849 start_time = omp_get_wtime()
850 y(:,:,1) = x(:,:,1)
851 call convert_ddcosmo(params, constants, 1, y(:,:,1))
852 n_iter = params % maxiter
853 call jacobi_diis(params, constants, workspace, constants % inner_tol, &
854 & y(:,:,1), workspace % ddcosmo_guess, n_iter, x_rel_diff, lstarx, &
855 & ldm1x, hnorm, ddx_error)
856 if (ddx_error % flag .ne. 0) then
857 call update_error(ddx_error, 'prec_tstarx: ddCOSMO failed to ' // &
858 & 'converge, exiting')
859 return
860 end if
861 y(:,:,1) = workspace % ddcosmo_guess
862 workspace % s_time = workspace % s_time + omp_get_wtime() - start_time
863
864 start_time = omp_get_wtime()
865 n_iter = params % maxiter
866 call jacobi_diis(params, constants, workspace, constants % inner_tol, &
867 & x(:,:,2), workspace % hsp_guess, n_iter, x_rel_diff, bstarx, &
868 & bx_prec, hnorm, ddx_error)
869 if (ddx_error % flag .ne. 0) then
870 call update_error(ddx_error, 'prec_tstarx: HSP failed to ' // &
871 & 'converge, exiting')
872 return
873 end if
874 y(:,:,2) = workspace % hsp_guess
875 workspace % hsp_adj_time = workspace % hsp_adj_time &
876 & + omp_get_wtime() - start_time
877
878end subroutine prec_tstarx
879
887subroutine prec_tx(params, constants, workspace, x, y, ddx_error)
888 implicit none
889 type(ddx_params_type), intent(in) :: params
890 type(ddx_constants_type), intent(in) :: constants
891 type(ddx_workspace_type), intent(inout) :: workspace
892 real(dp), intent(in) :: x(constants % nbasis, params % nsph, 2)
893 real(dp), intent(inout) :: y(constants % nbasis, params % nsph, 2)
894 type(ddx_error_type), intent(inout) :: ddx_error
895 integer :: n_iter
896 real(dp), dimension(params % maxiter) :: x_rel_diff
897 real(dp) :: start_time
898
899 ! perform A^-1 * Yr
900 start_time = omp_get_wtime()
901 n_iter = params % maxiter
902 call jacobi_diis(params, constants, workspace, constants % inner_tol, &
903 & x(:,:,1), workspace % ddcosmo_guess, n_iter, x_rel_diff, lx, &
904 & ldm1x, hnorm, ddx_error)
905 if (ddx_error % flag .ne. 0) then
906 call update_error(ddx_error, 'prec_tx: ddCOSMO failed to ' // &
907 & 'converge, exiting')
908 return
909 end if
910
911 ! Scale by the factor of (2l+1)/4Pi
912 y(:,:,1) = workspace % ddcosmo_guess
913 call convert_ddcosmo(params, constants, 1, y(:,:,1))
914 workspace % xs_time = workspace % xs_time + omp_get_wtime() - start_time
915
916 ! perform B^-1 * Ye
917 start_time = omp_get_wtime()
918 n_iter = params % maxiter
919 call jacobi_diis(params, constants, workspace, constants % inner_tol, &
920 & x(:,:,2), workspace % hsp_guess, n_iter, x_rel_diff, bx, &
921 & bx_prec, hnorm, ddx_error)
922 y(:,:,2) = workspace % hsp_guess
923 workspace % hsp_time = workspace % hsp_time + omp_get_wtime() - start_time
924
925 if (ddx_error % flag .ne. 0) then
926 call update_error(ddx_error, 'prec_tx: HSP failed to ' // &
927 & 'converge, exiting')
928 return
929 end if
930
931end subroutine prec_tx
932
934subroutine cstarx(params, constants, workspace, x, y, ddx_error)
935 implicit none
936 ! input/output
937 type(ddx_params_type), intent(in) :: params
938 type(ddx_constants_type), intent(in) :: constants
939 type(ddx_workspace_type), intent(inout) :: workspace
940 real(dp), dimension(constants % nbasis, params % nsph, 2), intent(in) :: x
941 real(dp), dimension(constants % nbasis, params % nsph, 2), intent(out) :: y
942 type(ddx_error_type), intent(inout) :: ddx_error
943 ! local
944 complex(dp) :: work_complex(constants % lmax0+1)
945 real(dp) :: work(constants % lmax0+1)
946 integer :: isph, igrid, jsph, ind, l0, ind0, indl, inode, l, m
947 real(dp), dimension(3) :: vij, vtij
948 real(dp) :: val, epsilon_ratio
949 real(dp), allocatable :: scratch(:,:), scratch0(:,:)
950
951 ! dummy operation on unused interface arguments
952 if (ddx_error % flag .eq. 0) continue
953
954 allocate(scratch(constants % nbasis, params % nsph), &
955 & scratch0(constants % nbasis0, params % nsph))
956
957 epsilon_ratio = params % epsp/params % eps
958
959 ! TODO: maybe use ddeval_grid for code consistency
960 scratch = - x(:,:,1) - x(:,:,2)
961 call dgemm('T', 'N', params % ngrid, params % nsph, constants % nbasis, &
962 & one, constants % vwgrid, constants % vgrid_nbasis, scratch, &
963 & constants % nbasis, zero, workspace % tmp_grid, params % ngrid)
964
965 scratch0 = zero
966 if (params % fmm .eq. 0) then
967 do isph = 1, params % nsph
968 do igrid = 1, params % ngrid
969 if (constants % ui(igrid, isph).gt.zero) then
970 val = workspace % tmp_grid(igrid,isph) &
971 & *constants % ui(igrid,isph)
972 ! quadratically scaling loop
973 do jsph = 1, params % nsph
974 vij = params % csph(:,isph) + &
975 & params % rsph(isph)*constants % cgrid(:,igrid) - &
976 & params % csph(:,jsph)
977 vtij = vij * params % kappa
978 call fmm_m2p_bessel_adj_work(vtij, val, &
979 & constants % SK_ri(:, jsph), constants % lmax0, &
980 & constants % vscales, one, scratch0(:, jsph), &
981 & work_complex, work)
982 end do
983 end if
984 end do
985 end do
986 else
987 ! Multiply by characteristic function
988 workspace % tmp_grid = workspace % tmp_grid * constants % ui
989 workspace % tmp_sph = zero
990 ! Do FMM operations adjointly
991 call tree_m2p_bessel_adj(params, constants, constants % lmax0, &
992 & one, workspace % tmp_grid, zero, params % lmax, &
993 & workspace % tmp_sph)
994 call tree_l2p_bessel_adj(params, constants, one, &
995 & workspace % tmp_grid, zero, workspace % tmp_node_l)
996 call tree_l2l_bessel_rotation_adj(params, constants, &
997 & workspace % tmp_node_l)
998 call tree_m2l_bessel_rotation_adj(params, constants, &
999 & workspace % tmp_node_l, workspace % tmp_node_m)
1000 call tree_m2m_bessel_rotation_adj(params, constants, &
1001 & workspace % tmp_node_m)
1002 ! Adjointly move tree multipole harmonics into output
1003 if(constants % lmax0 .lt. params % pm) then
1004 do isph = 1, params % nsph
1005 inode = constants % snode(isph)
1006 scratch0(:, isph) = workspace % tmp_sph(1:constants % nbasis0, isph) + &
1007 & workspace % tmp_node_m(1:constants % nbasis0, inode)
1008 end do
1009 else
1010 indl = (params % pm+1)**2
1011 do isph = 1, params % nsph
1012 inode = constants % snode(isph)
1013 scratch0(1:indl, isph) = &
1014 & workspace % tmp_sph(1:indl, isph) + &
1015 & workspace % tmp_node_m(:, inode)
1016 scratch0(indl+1:, isph) = zero
1017 end do
1018 end if
1019 end if
1020 ! Scale by C_ik
1021 do isph = 1, params % nsph
1022 do l0 = 0, constants % lmax0
1023 ind0 = l0*l0 + l0 + 1
1024 scratch0(ind0-l0:ind0+l0, isph) = scratch0(ind0-l0:ind0+l0, isph) * &
1025 & constants % C_ik(l0, isph)
1026 end do
1027 end do
1028
1029 do jsph = 1, params % nsph
1030 call dgemv('n', constants % nbasis, constants % nbasis0, one, &
1031 & constants % pchi(1,1,jsph), constants % nbasis, &
1032 & scratch0(1,jsph), 1, zero, scratch(1,jsph), 1)
1033 end do
1034
1035 do isph = 1, params % nsph
1036 do l = 0, params % lmax
1037 do m = -l, l
1038 ind = l**2 + l + m + 1
1039 y(ind,isph,1) = - (epsilon_ratio*l*scratch(ind,isph)) &
1040 & /params % rsph(isph)
1041 y(ind,isph,2) = constants % termimat(l,isph)*scratch(ind,isph)
1042 end do
1043 end do
1044 end do
1045
1046end subroutine cstarx
1047
1049subroutine cx(params, constants, workspace, x, y, ddx_error)
1050 implicit none
1051 type(ddx_params_type), intent(in) :: params
1052 type(ddx_constants_type), intent(in) :: constants
1053 type(ddx_workspace_type), intent(inout) :: workspace
1054 real(dp), dimension(constants % nbasis, params % nsph, 2), intent(in) :: x
1055 real(dp), dimension(constants % nbasis, params % nsph, 2), intent(out) :: y
1056 type(ddx_error_type), intent(inout) :: ddx_error
1057
1058 integer :: isph, jsph, igrid, ind, l, m, ind0
1059 real(dp), dimension(3) :: vij, vtij
1060 real(dp) :: val
1061 complex(dp) :: work_complex(constants % lmax0 + 1)
1062 real(dp) :: work(constants % lmax0 + 1)
1063 integer :: indl, inode, info
1064
1065 real(dp), allocatable :: diff_re(:,:), diff0(:,:)
1066
1067 allocate(diff_re(constants % nbasis, params % nsph), &
1068 & diff0(constants % nbasis0, params % nsph), stat=info)
1069 if (info.ne.0) then
1070 call update_error(ddx_error, "Allocation failed in cx")
1071 return
1072 end if
1073
1074 ! diff_re = params % epsp/eps*l1/ri*Xr - i'(ri)/i(ri)*Xe,
1075 diff_re = zero
1076 !$omp parallel do default(none) shared(params,diff_re, &
1077 !$omp constants,x) private(jsph,l,m,ind)
1078 do jsph = 1, params % nsph
1079 do l = 0, params % lmax
1080 do m = -l,l
1081 ind = l**2 + l + m + 1
1082 diff_re(ind,jsph) = (params % epsp/params % eps)* &
1083 & (l/params % rsph(jsph))*x(ind,jsph,1) &
1084 & - constants % termimat(l,jsph)*x(ind,jsph,2)
1085 end do
1086 end do
1087 end do
1088
1089 ! diff0 = Pchi * diff_er, linear scaling
1090 !$omp parallel do default(none) shared(constants,params, &
1091 !$omp diff_re,diff0) private(jsph)
1092 do jsph = 1, params % nsph
1093 call dgemv('t', constants % nbasis, constants % nbasis0, one, &
1094 & constants % pchi(1,1,jsph), constants % nbasis, &
1095 & diff_re(1,jsph), 1, zero, diff0(1,jsph), 1)
1096 end do
1097
1098 ! Multiply diff0 by C_ik inplace
1099 do isph = 1, params % nsph
1100 do l = 0, constants % lmax0
1101 ind0 = l*l+l+1
1102 diff0(ind0-l:ind0+l, isph) = diff0(ind0-l:ind0+l, isph) * &
1103 & constants % C_ik(l, isph)
1104 end do
1105 end do
1106 ! avoiding N^2 storage, this code does not use the cached coefY
1107 y(:,:,1) = zero
1108 if (params % fmm .eq. 0) then
1109 !$omp parallel do default(none) shared(params,constants, &
1110 !$omp diff0,y) private(isph,igrid,val,vij,work, &
1111 !$omp ind0,ind,vtij,work_complex)
1112 do isph = 1, params % nsph
1113 do igrid = 1, params % ngrid
1114 if (constants % ui(igrid,isph).gt.zero) then
1115 val = zero
1116
1117 ! quadratically scaling loop
1118 do jsph = 1, params % nsph
1119 vij = params % csph(:,isph) + &
1120 & params % rsph(isph)*constants % cgrid(:,igrid) - &
1121 & params % csph(:,jsph)
1122 vtij = vij * params % kappa
1123 call fmm_m2p_bessel_work(vtij, constants % lmax0, &
1124 & constants % vscales, constants % SK_ri(:, jsph), &
1125 & one, diff0(:, jsph), one, val, work_complex, work)
1126 end do
1127 do ind = 1, constants % nbasis
1128 y(ind,isph,1) = y(ind,isph,1) + val*&
1129 & constants % ui(igrid,isph)*&
1130 & constants % vwgrid(ind,igrid)
1131 end do
1132 end if
1133 end do
1134 end do
1135 else
1136 ! Load input harmonics into tree data
1137 workspace % tmp_node_m = zero
1138 workspace % tmp_node_l = zero
1139 workspace % tmp_sph = zero
1140 do isph = 1, params % nsph
1141 do l = 0, constants % lmax0
1142 ind0 = l*l+l+1
1143 workspace % tmp_sph(ind0-l:ind0+l, isph) = &
1144 & diff0(ind0-l:ind0+l, isph)
1145 end do
1146 end do
1147 if(constants % lmax0 .lt. params % pm) then
1148 do isph = 1, params % nsph
1149 inode = constants % snode(isph)
1150 workspace % tmp_node_m(1:constants % nbasis0, inode) = &
1151 & workspace % tmp_sph(1:constants % nbasis0, isph)
1152 workspace % tmp_node_m(constants % nbasis0+1:, inode) = zero
1153 end do
1154 else
1155 indl = (params % pm+1)**2
1156 do isph = 1, params % nsph
1157 inode = constants % snode(isph)
1158 workspace % tmp_node_m(:, inode) = workspace % tmp_sph(1:indl, isph)
1159 end do
1160 end if
1161 ! Do FMM operations
1162 call tree_m2m_bessel_rotation(params, constants, workspace % tmp_node_m)
1163 call tree_m2l_bessel_rotation(params, constants, workspace % tmp_node_m, &
1164 & workspace % tmp_node_l)
1165 call tree_l2l_bessel_rotation(params, constants, workspace % tmp_node_l)
1166 workspace % tmp_grid = zero
1167 call tree_l2p_bessel(params, constants, one, workspace % tmp_node_l, zero, &
1168 & workspace % tmp_grid)
1169 call tree_m2p_bessel(params, constants, constants % lmax0, one, &
1170 & params % lmax, workspace % tmp_sph, one, &
1171 & workspace % tmp_grid)
1172
1173 do isph = 1, params % nsph
1174 do igrid = 1, params % ngrid
1175 do ind = 1, constants % nbasis
1176 y(ind,isph,1) = y(ind,isph,1) + workspace % tmp_grid(igrid, isph)*&
1177 & constants % vwgrid(ind, igrid)*&
1178 & constants % ui(igrid,isph)
1179 end do
1180 end do
1181 end do
1182 end if
1183
1184 y(:,:,2) = y(:,:,1)
1185
1186 deallocate(diff_re, diff0, stat=info)
1187 if (info.ne.0) then
1188 call update_error(ddx_error, "Deallocation failed in cx")
1189 return
1190 end if
1191
1192end subroutine cx
1193
1194end module ddx_operators
ddx_gradients
Core routines and parameters specific to gradients.
Definition: ddx_gradients.f90:11
ddx_lpb_core
Core routines and parameters specific to ddLPB.
Definition: ddx_lpb_core.f90:11
ddx_operators
Operators shared among ddX methods.
Definition: ddx_operators.f90:13
ddx_operators::repsx
subroutine repsx(params, constants, workspace, x, y, ddx_error)
Apply operator to spherical harmonics.
Definition: ddx_operators.f90:528
ddx_operators::prec_repsstarx
subroutine prec_repsstarx(params, constants, workspace, x, y, ddx_error)
Apply preconditioner for ddPCM adjoint linear system.
Definition: ddx_operators.f90:666
ddx_operators::prec_repsx
subroutine prec_repsx(params, constants, workspace, x, y, ddx_error)
Apply preconditioner for ddPCM primal linear system.
Definition: ddx_operators.f90:638
ddx_operators::bx_prec
subroutine bx_prec(params, constants, workspace, x, y, ddx_error)
Apply the preconditioner to the primal HSP linear system.
Definition: ddx_operators.f90:821
ddx_operators::lx
subroutine lx(params, constants, workspace, x, y, ddx_error)
Single layer operator matvec.
Definition: ddx_operators.f90:26
ddx_operators::dstarx_dense
subroutine dstarx_dense(params, constants, workspace, do_diag, x, y, ddx_error)
Baseline implementation of adjoint double layer operator.
Definition: ddx_operators.f90:383
ddx_operators::dx
subroutine dx(params, constants, workspace, x, y, ddx_error)
Double layer operator matvec without diagonal blocks.
Definition: ddx_operators.f90:186
ddx_operators::dstarx_fmm
subroutine dstarx_fmm(params, constants, workspace, do_diag, x, y, ddx_error)
FMM-accelerated implementation of adjoint double layer operator.
Definition: ddx_operators.f90:460
ddx_operators::bstarx
subroutine bstarx(params, constants, workspace, x, y, ddx_error)
Adjoint HSP matrix vector product.
Definition: ddx_operators.f90:695
ddx_operators::prec_tx
subroutine prec_tx(params, constants, workspace, x, y, ddx_error)
Apply the preconditioner to the primal ddLPB linear system |Yr| = |A^-1 0 |*|Xr| |Ye| |0 B^-1 | |Xe|.
Definition: ddx_operators.f90:888
ddx_operators::rinfx
subroutine rinfx(params, constants, workspace, x, y, ddx_error)
Apply operator to spherical harmonics.
Definition: ddx_operators.f90:588
ddx_operators::bx
subroutine bx(params, constants, workspace, x, y, ddx_error)
Primal HSP matrix vector product.
Definition: ddx_operators.f90:747
ddx_operators::cstarx
subroutine cstarx(params, constants, workspace, x, y, ddx_error)
ddLPB adjoint matrix-vector product
Definition: ddx_operators.f90:935
ddx_operators::dx_fmm
subroutine dx_fmm(params, constants, workspace, do_diag, x, y, ddx_error)
FMM-accelerated implementation of double layer operator.
Definition: ddx_operators.f90:292
ddx_operators::cx
subroutine cx(params, constants, workspace, x, y, ddx_error)
ddLPB matrix-vector product
Definition: ddx_operators.f90:1050
ddx_operators::dx_dense
subroutine dx_dense(params, constants, workspace, do_diag, x, y, ddx_error)
Baseline implementation of double layer operator.
Definition: ddx_operators.f90:210
ddx_operators::tstarx
subroutine tstarx(params, constants, workspace, x, y, ddx_error)
Adjoint ddLPB matrix-vector product.
Definition: ddx_operators.f90:787
ddx_operators::ldm1x
subroutine ldm1x(params, constants, workspace, x, y, ddx_error)
Diagonal preconditioning for Lx operator.
Definition: ddx_operators.f90:155
ddx_operators::repsstarx
subroutine repsstarx(params, constants, workspace, x, y, ddx_error)
Apply operator to spherical harmonics.
Definition: ddx_operators.f90:556
ddx_operators::dstarx
subroutine dstarx(params, constants, workspace, x, y, ddx_error)
Adjoint double layer operator matvec.
Definition: ddx_operators.f90:358
ddx_operators::lstarx
subroutine lstarx(params, constants, workspace, x, y, ddx_error)
Adjoint single layer operator matvec.
Definition: ddx_operators.f90:89
ddx_operators::rstarinfx
subroutine rstarinfx(params, constants, workspace, x, y, ddx_error)
Apply operator to spherical harmonics.
Definition: ddx_operators.f90:619
ddx_operators::prec_tstarx
subroutine prec_tstarx(params, constants, workspace, x, y, ddx_error)
Apply the preconditioner to the ddLPB adjoint linear system.
Definition: ddx_operators.f90:838
ddx_solvers
Core routines and parameters of ddX software.
Definition: ddx_solvers.f90:13
ddx_solvers::jacobi_diis
subroutine jacobi_diis(params, constants, workspace, tol, rhs, x, niter, x_rel_diff, matvec, dm1vec, norm_func, ddx_error)
Jacobi iterative solver.
Definition: ddx_solvers.f90:86