Line data Source code
1 : !-----------------------------------------------------------------------
2 : ! $Id$
3 : !===============================================================================
4 : module matrix_solver_wrapper
5 :
6 : use clubb_precision, only: &
7 : core_rknd ! Variable(s)
8 :
9 : use error_code, only: &
10 : clubb_at_least_debug_level, & ! Procedure
11 : err_code, & ! Error indicator
12 : clubb_no_error, & ! Constant
13 : clubb_fatal_error
14 :
15 : use constants_clubb, only: &
16 : fstderr ! Constant(s)
17 :
18 : use model_flags, only: &
19 : lapack , & ! Variable(s)
20 : penta_lu , &
21 : tridiag_lu, &
22 : penta_bicgstab
23 :
24 : implicit none
25 :
26 : public :: band_solve, &
27 : tridiag_solve
28 :
29 : interface band_solve
30 : module procedure band_solve_single_rhs_multiple_lhs
31 : module procedure band_solve_multiple_rhs_lhs
32 : end interface
33 :
34 : interface tridiag_solve
35 : module procedure tridiag_solve_single_rhs_lhs
36 : module procedure tridiag_solve_single_rhs_multiple_lhs
37 : module procedure tridiag_solve_multiple_rhs_lhs
38 : end interface
39 :
40 : contains
41 :
42 :
43 : !-------------------------------------------------------------------
44 : ! Band Solvers Procedures
45 : !-------------------------------------------------------------------
46 :
47 8935056 : subroutine band_solve_single_rhs_multiple_lhs( &
48 : solve_name, penta_solve_method, & ! Intent(in)
49 : ngrdcol, nsup, nsub, ndim, & ! Intent(in)
50 : old_soln, & ! Intent(in)
51 8935056 : lhs, rhs, & ! Intent(inout)
52 8935056 : soln, rcond ) ! Intent(out)
53 :
54 : use lapack_wrap, only: &
55 : lapack_band_solve, & ! Procedure(s)
56 : lapack_band_solvex
57 :
58 : use penta_lu_solvers, only: &
59 : penta_lu_solve ! Procedure(s)
60 :
61 : implicit none
62 :
63 : ! ----------------------- Input Variables -----------------------
64 : character(len=*), intent(in) :: &
65 : solve_name
66 :
67 : integer, intent(in) :: &
68 : penta_solve_method
69 :
70 : integer, intent(in) :: &
71 : ngrdcol, & ! Number of grid columns
72 : nsup, & ! Number of superdiagonals
73 : nsub, & ! Number of subdiagonals
74 : ndim ! The order of the LHS Matrix, i.e. the # of linear equations
75 :
76 : real( kind = core_rknd ), dimension(ngrdcol,ndim), intent(in) :: &
77 : old_soln ! Old soln, used as an initial guess in the bicgstab method
78 :
79 : real( kind = core_rknd ), dimension(nsup+nsub+1,ngrdcol,ndim), intent(inout) :: &
80 : lhs ! Left hand side
81 :
82 : real( kind = core_rknd ), dimension(ngrdcol,ndim), intent(inout) :: &
83 : rhs ! Right hand side(s)
84 :
85 : ! ----------------------- Output Variables -----------------------
86 : real( kind = core_rknd ), dimension(ngrdcol,ndim), intent(out) :: &
87 : soln
88 :
89 : ! ----------------------- Optional Out -----------------------
90 :
91 : ! The estimate of the reciprocal condition number of matrix
92 : ! after equilibration (if done).
93 : real( kind = core_rknd ), optional, dimension(ngrdcol), intent(out) :: &
94 : rcond
95 :
96 : ! ----------------------- Local Variables -----------------------
97 : real( kind = core_rknd ), dimension(nsup+nsub+1,ngrdcol,ndim) :: &
98 17870112 : lhs_copy ! Copy of left hand side
99 :
100 : real( kind = core_rknd ), dimension(ngrdcol,ndim) :: &
101 17870112 : rhs_copy ! Copy of right hand side
102 :
103 : real( kind = core_rknd ), dimension(ngrdcol,ndim) :: &
104 17870112 : dummy_soln
105 :
106 : integer :: i
107 :
108 : ! ----------------------- Begin Code -----------------------
109 :
110 8935056 : if ( present(rcond) ) then
111 :
112 : !$acc update host( lhs, rhs )
113 :
114 : ! Lapack overwrites lhs and rhs, so we'll give it copies of them.
115 0 : lhs_copy = lhs
116 0 : rhs_copy = rhs
117 :
118 : ! Perform LU decomp and solve system (LAPACK with diagnostics)
119 : ! Using dummy_soln, since we only want this routine for diagnostics
120 : call lapack_band_solvex( "xm_wpxp", nsup, nsub, & ! Intent(in)
121 : ndim, 1, ngrdcol, & ! Intent(in)
122 : lhs_copy, rhs_copy, & ! Intent(inout)
123 0 : dummy_soln, rcond ) ! Intent(out)
124 :
125 : !$acc update device( rcond )
126 :
127 : end if
128 :
129 :
130 8935056 : if ( penta_solve_method == lapack ) then
131 :
132 : !$acc update host( lhs, rhs )
133 :
134 : ! Perform LU decomp and solve system (LAPACK)
135 : call lapack_band_solve( "xm_wpxp", nsup, nsub, & ! Intent(in)
136 : ndim, 1, ngrdcol, & ! Intent(in)
137 : lhs, rhs, & ! Intent(inout)
138 8935056 : soln ) ! Intent(out)
139 :
140 : !$acc update device( soln )
141 :
142 0 : else if ( penta_solve_method == penta_lu ) then
143 :
144 : ! Solve the system with a penta-diagonal specific LU decomp
145 : call penta_lu_solve( ndim, ngrdcol, & ! Intent(in)
146 : lhs(:,:,:), rhs(:,:), & ! Intent(in)
147 0 : soln(:,:) ) ! Intent(out)
148 :
149 : else
150 :
151 : ! The solve method should match one of the above
152 0 : if ( clubb_at_least_debug_level( 0 ) ) then
153 0 : write(fstderr,*) "Error in band_solve_single_rhs_multiple_lhs: "
154 0 : write(fstderr,*) " no case for penta_solve_method = ", penta_solve_method
155 0 : err_code = clubb_fatal_error
156 : end if
157 :
158 : end if
159 :
160 8935056 : return
161 :
162 : end subroutine band_solve_single_rhs_multiple_lhs
163 :
164 8935056 : subroutine band_solve_multiple_rhs_lhs( &
165 : solve_name, penta_solve_method, & ! Intent(in)
166 : ngrdcol, nsup, nsub, ndim, nrhs, & ! Intent(in)
167 : old_soln, & ! Intent(in)
168 8935056 : lhs, rhs, & ! Intent(inout)
169 8935056 : soln, rcond ) ! Intent(out)
170 :
171 : use lapack_wrap, only: &
172 : lapack_band_solve, & ! Procedure(s)
173 : lapack_band_solvex
174 :
175 : use penta_lu_solvers, only: &
176 : penta_lu_solve ! Procedure(s)
177 :
178 : implicit none
179 :
180 : ! ----------------------- Input Variables -----------------------
181 : character(len=*), intent(in) :: &
182 : solve_name
183 :
184 : integer, intent(in) :: &
185 : penta_solve_method
186 :
187 : integer, intent(in) :: &
188 : ngrdcol, & ! Number of grid columns
189 : nsup, & ! Number of superdiagonals
190 : nsub, & ! Number of subdiagonals
191 : ndim, & ! The order of the LHS Matrix, i.e. the # of linear equations
192 : nrhs ! Number of RHS's to back substitute for
193 :
194 : real( kind = core_rknd ), dimension(ngrdcol,ndim,nrhs), intent(in) :: &
195 : old_soln ! Old soln, used as an initial guess in the bicgstab method
196 :
197 : real( kind = core_rknd ), dimension(nsup+nsub+1,ngrdcol,ndim), intent(inout) :: &
198 : lhs ! Left hand side
199 :
200 : real( kind = core_rknd ), dimension(ngrdcol,ndim,nrhs), intent(inout) :: &
201 : rhs ! Right hand side(s)
202 :
203 : ! ----------------------- Output Variables -----------------------
204 : real( kind = core_rknd ), dimension(ngrdcol,ndim,nrhs), intent(out) :: &
205 : soln
206 :
207 : ! ----------------------- Optional Out -----------------------
208 :
209 : ! The estimate of the reciprocal condition number of matrix
210 : ! after equilibration (if done).
211 : real( kind = core_rknd ), optional, dimension(ngrdcol), intent(out) :: &
212 : rcond
213 :
214 : ! ----------------------- Local Variables -----------------------
215 : real( kind = core_rknd ), dimension(nsup+nsub+1,ngrdcol,ndim) :: &
216 17870112 : lhs_copy ! Copy of left hand side
217 :
218 : real( kind = core_rknd ), dimension(ngrdcol,ndim,nrhs) :: &
219 17870112 : rhs_copy ! Copy of right hand side
220 :
221 : real( kind = core_rknd ), dimension(ngrdcol,ndim,nrhs) :: &
222 17870112 : dummy_soln
223 :
224 : integer :: i
225 :
226 : ! ----------------------- Begin Code -----------------------
227 :
228 8935056 : if ( present(rcond) ) then
229 :
230 : !$acc update host( lhs, rhs )
231 :
232 : ! Lapack overwrites lhs and rhs, so we'll give it copies of them.
233 0 : lhs_copy = lhs
234 0 : rhs_copy = rhs
235 :
236 : ! Perform LU decomp and solve system (LAPACK with diagnostics)
237 : ! Using dummy_soln, since we only want this routine for diagnostics
238 : call lapack_band_solvex( "xm_wpxp", nsup, nsub, & ! Intent(in)
239 : ndim, nrhs, ngrdcol, & ! Intent(in)
240 : lhs_copy, rhs_copy, & ! Intent(inout)
241 0 : dummy_soln, rcond ) ! Intent(out)
242 :
243 : !$acc update device( rcond )
244 :
245 : end if
246 :
247 :
248 8935056 : if ( penta_solve_method == lapack ) then
249 :
250 : !$acc update host( lhs, rhs )
251 :
252 : ! Perform LU decomp and solve system (LAPACK)
253 : call lapack_band_solve( "xm_wpxp", nsup, nsub, & ! Intent(in)
254 : ndim, nrhs, ngrdcol, & ! Intent(in)
255 : lhs, rhs, & ! Intent(inout)
256 8935056 : soln ) ! Intent(out)
257 :
258 : !$acc update device( soln )
259 :
260 0 : else if ( penta_solve_method == penta_lu ) then
261 :
262 : ! Solve the system with a penta-diagonal specific LU decomp
263 : call penta_lu_solve( ndim, nrhs, ngrdcol, & ! Intent(in)
264 : lhs(:,:,:), rhs(:,:,:), & ! Intent(in)
265 0 : soln(:,:,:) ) ! Intent(out)
266 :
267 : else
268 :
269 : ! The solve method should match one of the above
270 0 : if ( clubb_at_least_debug_level( 0 ) ) then
271 0 : write(fstderr,*) "Error in band_solve_multiple_rhs_lhs: "
272 0 : write(fstderr,*) " no case for penta_solve_method = ", penta_solve_method
273 0 : err_code = clubb_fatal_error
274 : end if
275 :
276 : end if
277 :
278 8935056 : return
279 :
280 : end subroutine band_solve_multiple_rhs_lhs
281 :
282 :
283 :
284 : !-------------------------------------------------------------------
285 : ! Tridiag Solver Procedures
286 : !-------------------------------------------------------------------
287 :
288 0 : subroutine tridiag_solve_single_rhs_lhs( &
289 : solve_name, tridiag_solve_method, & ! Intent(in)
290 : ndim, & ! Intent(in)
291 0 : lhs, rhs, & ! Intent(inout)
292 0 : soln, rcond ) ! Intent(out)
293 :
294 : use lapack_wrap, only: &
295 : lapack_tridiag_solve, & ! Procedure(s)
296 : lapack_tridiag_solvex
297 :
298 : use tridiag_lu_solvers, only: &
299 : tridiag_lu_solve ! Procedure(s)
300 :
301 : implicit none
302 :
303 : ! ----------------------- Input Variables -----------------------
304 : character(len=*), intent(in) :: &
305 : solve_name
306 :
307 : integer, intent(in) :: &
308 : tridiag_solve_method
309 :
310 : integer, intent(in) :: &
311 : ndim ! The order of the LHS Matrix, i.e. the # of linear equations
312 :
313 : real( kind = core_rknd ), dimension(3,ndim), intent(inout) :: &
314 : lhs ! Left hand side
315 :
316 : real( kind = core_rknd ), dimension(ndim), intent(inout) :: &
317 : rhs ! Right hand side(s)
318 :
319 : ! ----------------------- Output Variables -----------------------
320 : real( kind = core_rknd ), dimension(ndim), intent(out) :: &
321 : soln
322 :
323 : ! ----------------------- Optional Out -----------------------
324 :
325 : ! The estimate of the reciprocal condition number of matrix
326 : ! after equilibration (if done).
327 : real( kind = core_rknd ), dimension(1), optional, intent(out) :: &
328 : rcond
329 :
330 : ! ----------------------- Local Variables -----------------------
331 : real( kind = core_rknd ), dimension(3,ndim) :: &
332 0 : lhs_copy ! Copy of left hand side
333 :
334 : real( kind = core_rknd ), dimension(ndim) :: &
335 0 : rhs_copy ! Copy of right hand side
336 :
337 : real( kind = core_rknd ), dimension(ndim) :: &
338 0 : dummy_soln
339 :
340 : ! ----------------------- Begin Code -----------------------
341 :
342 0 : if ( present(rcond) ) then
343 :
344 : ! Lapack overwrites lhs and rhs, so we'll give it copies of them.
345 0 : lhs_copy = lhs
346 0 : rhs_copy = rhs
347 :
348 : ! Perform LU decomp and solve system (LAPACK with diagnostics)
349 : call lapack_tridiag_solvex( solve_name, ndim, 1, 1, & ! Intent(in)
350 : lhs_copy, rhs_copy, & ! Intent(inout)
351 0 : dummy_soln, rcond ) ! Intent(out)
352 : end if
353 :
354 :
355 0 : if ( tridiag_solve_method == lapack ) then
356 :
357 : ! Perform LU decomp and solve system (LAPACK)
358 : call lapack_tridiag_solve( solve_name, ndim, 1, 1, & ! Intent(in)
359 : lhs, rhs, & ! Intent(inout)
360 0 : soln ) ! Intent(out)
361 :
362 0 : else if ( tridiag_solve_method == tridiag_lu ) then
363 :
364 : ! Solve the system with a tri-diagonal specific LU decomp
365 : call tridiag_lu_solve( ndim, & ! Intent(in)
366 : lhs, rhs, & ! Intent(in)
367 0 : soln ) ! Intent(out)
368 :
369 : else
370 :
371 : ! The solve method should match one of the above
372 0 : if ( clubb_at_least_debug_level( 0 ) ) then
373 0 : write(fstderr,*) "Error in tridiag_solve_single_rhs_lhs: "
374 0 : write(fstderr,*) " no case for tridiag_solve_method = ", tridiag_solve_method
375 0 : err_code = clubb_fatal_error
376 : end if
377 :
378 : end if
379 :
380 0 : return
381 :
382 : end subroutine tridiag_solve_single_rhs_lhs
383 :
384 :
385 405213 : subroutine tridiag_solve_single_rhs_multiple_lhs( &
386 : solve_name, tridiag_solve_method, & ! Intent(in)
387 : ngrdcol, ndim, & ! Intent(in)
388 405213 : lhs, rhs, & ! Intent(inout)
389 405213 : soln, rcond ) ! Intent(out)
390 :
391 : use lapack_wrap, only: &
392 : lapack_tridiag_solve, & ! Procedure(s)
393 : lapack_tridiag_solvex
394 :
395 : use tridiag_lu_solvers, only: &
396 : tridiag_lu_solve ! Procedure(s)
397 :
398 : implicit none
399 :
400 : ! ----------------------- Input Variables -----------------------
401 : character(len=*), intent(in) :: &
402 : solve_name
403 :
404 : integer, intent(in) :: &
405 : tridiag_solve_method
406 :
407 : integer, intent(in) :: &
408 : ngrdcol, & ! Number of grid columns
409 : ndim ! The order of the LHS Matrix, i.e. the # of linear equations
410 :
411 : real( kind = core_rknd ), dimension(3,ngrdcol,ndim), intent(inout) :: &
412 : lhs ! Left hand side
413 :
414 : real( kind = core_rknd ), dimension(ngrdcol,ndim), intent(inout) :: &
415 : rhs ! Right hand side(s)
416 :
417 : ! ----------------------- Output Variables -----------------------
418 : real( kind = core_rknd ), dimension(ngrdcol,ndim), intent(out) :: &
419 : soln
420 :
421 : ! ----------------------- Optional Out -----------------------
422 :
423 : ! The estimate of the reciprocal condition number of matrix
424 : ! after equilibration (if done).
425 : real( kind = core_rknd ), optional, dimension(ngrdcol), intent(out) :: &
426 : rcond
427 :
428 : ! ----------------------- Local Variables -----------------------
429 : real( kind = core_rknd ), dimension(3,ngrdcol,ndim) :: &
430 810426 : lhs_copy ! Copy of left hand side
431 :
432 : real( kind = core_rknd ), dimension(ngrdcol,ndim) :: &
433 810426 : rhs_copy ! Copy of right hand side
434 :
435 : real( kind = core_rknd ), dimension(ngrdcol,ndim) :: &
436 810426 : dummy_soln
437 :
438 : integer :: i
439 :
440 : ! ----------------------- Begin Code -----------------------
441 :
442 405213 : if ( present(rcond) ) then
443 :
444 : !$acc update host( lhs, rhs )
445 :
446 : ! Lapack overwrites lhs and rhs, so we'll give it copies of them.
447 0 : lhs_copy = lhs
448 0 : rhs_copy = rhs
449 :
450 : ! Perform LU decomp and solve system (LAPACK with diagnostics)
451 : call lapack_tridiag_solvex( solve_name, ndim, 1, ngrdcol, & ! Intent(in)
452 : lhs_copy, rhs_copy, & ! Intent(inout)
453 0 : dummy_soln, rcond ) ! Intent(out)
454 :
455 : !$acc update device( rcond )
456 :
457 : end if
458 :
459 :
460 405213 : if ( tridiag_solve_method == lapack ) then
461 :
462 : !$acc update host( lhs, rhs )
463 :
464 : ! Perform LU decomp and solve system (LAPACK)
465 : call lapack_tridiag_solve( solve_name, ndim, 1, ngrdcol, & ! Intent(in)
466 : lhs, rhs, & ! Intent(inout)
467 405213 : soln ) ! Intent(out)
468 :
469 : !$acc update device( soln )
470 :
471 0 : else if ( tridiag_solve_method == tridiag_lu ) then
472 :
473 : ! Solve the system with a tri-diagonal specific LU decomp
474 : call tridiag_lu_solve( ndim, ngrdcol, & ! Intent(in)
475 : lhs, rhs, & ! Intent(in)
476 0 : soln ) ! Intent(out)
477 : else
478 :
479 : ! The solve method should match one of the above
480 0 : if ( clubb_at_least_debug_level( 0 ) ) then
481 0 : write(fstderr,*) "Error in tridiag_solve_single_rhs_multiple_lhs: "
482 0 : write(fstderr,*) " no case for tridiag_solve_method = ", tridiag_solve_method
483 0 : err_code = clubb_fatal_error
484 : end if
485 :
486 : end if
487 :
488 405213 : return
489 :
490 : end subroutine tridiag_solve_single_rhs_multiple_lhs
491 :
492 44675280 : subroutine tridiag_solve_multiple_rhs_lhs( &
493 : solve_name, tridiag_solve_method, & ! Intent(in)
494 : ngrdcol, ndim, nrhs, & ! Intent(in)
495 44675280 : lhs, rhs, & ! Intent(inout)
496 44675280 : soln, rcond ) ! Intent(out)
497 :
498 : use lapack_wrap, only: &
499 : lapack_tridiag_solve, & ! Procedure(s)
500 : lapack_tridiag_solvex
501 :
502 : use tridiag_lu_solvers, only: &
503 : tridiag_lu_solve ! Procedure(s)
504 :
505 : implicit none
506 :
507 : ! ----------------------- Input Variables -----------------------
508 : character(len=*), intent(in) :: &
509 : solve_name
510 :
511 : integer, intent(in) :: &
512 : tridiag_solve_method
513 :
514 : integer, intent(in) :: &
515 : ngrdcol, & ! Number of grid columns
516 : ndim, & ! The order of the LHS Matrix, i.e. the # of linear equations
517 : nrhs ! Number of RHS's to back substitute for
518 :
519 : real( kind = core_rknd ), dimension(3,ngrdcol,ndim), intent(inout) :: &
520 : lhs ! Left hand side
521 :
522 : real( kind = core_rknd ), dimension(ngrdcol,ndim,nrhs), intent(inout) :: &
523 : rhs ! Right hand side(s)
524 :
525 : ! ----------------------- Output Variables -----------------------
526 : real( kind = core_rknd ), dimension(ngrdcol,ndim,nrhs), intent(out) :: &
527 : soln
528 :
529 : ! ----------------------- Optional Out -----------------------
530 :
531 : ! The estimate of the reciprocal condition number of matrix
532 : ! after equilibration (if done).
533 : real( kind = core_rknd ), optional, dimension(ngrdcol), intent(out) :: &
534 : rcond
535 :
536 : ! ----------------------- Local Variables -----------------------
537 : real( kind = core_rknd ), dimension(3,ngrdcol,ndim) :: &
538 89350560 : lhs_copy ! Copy of left hand side
539 :
540 : real( kind = core_rknd ), dimension(ngrdcol,ndim,nrhs) :: &
541 89350560 : rhs_copy ! Copy of right hand side
542 :
543 : real( kind = core_rknd ), dimension(ngrdcol,ndim,nrhs) :: &
544 89350560 : dummy_soln
545 :
546 : integer :: i
547 :
548 : ! ----------------------- Begin Code -----------------------
549 :
550 44675280 : if ( present(rcond) ) then
551 :
552 : !$acc update host( lhs, rhs )
553 :
554 : ! Lapack overwrites lhs and rhs, so we'll give it copies of them.
555 0 : lhs_copy = lhs
556 0 : rhs_copy = rhs
557 :
558 : ! Perform LU decomp and solve system (LAPACK with diagnostics)
559 : call lapack_tridiag_solvex( solve_name, ndim, nrhs, ngrdcol, & ! Intent(in)
560 : lhs_copy, rhs_copy, & ! Intent(inout)
561 0 : dummy_soln, rcond ) ! Intent(out)
562 :
563 : !$acc update device( rcond )
564 :
565 : end if
566 :
567 :
568 44675280 : if ( tridiag_solve_method == lapack ) then
569 :
570 : !$acc update host( lhs, rhs )
571 :
572 : ! Perform LU decomp and solve system (LAPACK)
573 : call lapack_tridiag_solve( solve_name, ndim, nrhs, ngrdcol, & ! Intent(in)
574 : lhs, rhs, & ! Intent(inout)
575 44675280 : soln ) ! Intent(out)
576 :
577 : !$acc update device( soln )
578 :
579 0 : else if ( tridiag_solve_method == tridiag_lu ) then
580 :
581 : ! Solve the system with a tri-diagonal specific LU decomp
582 : call tridiag_lu_solve( ndim, nrhs, ngrdcol, & ! Intent(in)
583 : lhs, rhs, & ! Intent(in)
584 0 : soln ) ! Intent(out)
585 : else
586 :
587 : ! The solve method should match one of the above
588 0 : if ( clubb_at_least_debug_level( 0 ) ) then
589 0 : write(fstderr,*) "Error in tridiag_solve_multiple_rhs_lhs: "
590 0 : write(fstderr,*) " no case for tridiag_solve_method = ", tridiag_solve_method
591 0 : err_code = clubb_fatal_error
592 : end if
593 :
594 : end if
595 :
596 44675280 : return
597 :
598 : end subroutine tridiag_solve_multiple_rhs_lhs
599 :
600 :
601 : end module matrix_solver_wrapper
|