OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
linear_solver_mod.F
Go to the documentation of this file.
1Copyright> OpenRadioss
2Copyright> Copyright (C) 1986-2025 Altair Engineering Inc.
3Copyright>
4Copyright> This program is free software: you can redistribute it and/or modify
5Copyright> it under the terms of the GNU Affero General Public License as published by
6Copyright> the Free Software Foundation, either version 3 of the License, or
7Copyright> (at your option) any later version.
8Copyright>
9Copyright> This program is distributed in the hope that it will be useful,
10Copyright> but WITHOUT ANY WARRANTY; without even the implied warranty of
11Copyright> MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12Copyright> GNU Affero General Public License for more details.
13Copyright>
14Copyright> You should have received a copy of the GNU Affero General Public License
15Copyright> along with this program. If not, see <https://www.gnu.org/licenses/>.
16Copyright>
17Copyright>
18Copyright> Commercial Alternative: Altair Radioss Software
19Copyright>
20Copyright> As an alternative to this open-source version, Altair also offers Altair Radioss
21Copyright> software under a commercial license. Contact Altair to discuss further if the
22Copyright> commercial version may interest you: https://www.altair.com/radioss/.
23!||====================================================================
24!|| linear_solver_mod ../engine/share/modules/linear_solver_mod.F
25!||--- called by ------------------------------------------------------
26!|| diffusion_mod ../engine/share/modules/diffusion_mod.F
27!|| init_diffusion ../engine/share/modules/diffusion_mod.F
28!|| radioss2 ../engine/source/engine/radioss2.F
29!||--- calls -----------------------------------------------------
30!|| get_global_dim ../engine/share/modules/linear_solver_mod.F
31!||--- uses -----------------------------------------------------
32!|| matrix_mod ../common_source/linearalgebra/matrix_mod.F
33!|| vector_mod ../common_source/linearalgebra/vector_mod.F
34!||====================================================================
35 MODULE linear_solver_mod
36 USE matrix_mod
37 USE vector_mod
38 implicit none
39#ifdef MUMPS5
40#include "dmumps_struc.h"
41#endif
42#include "my_real.inc"
43
44! .___________________________. !
45! | | !
46! | ********************** | !
47! | ** Type definitions ** | !
48! | ********************** | !
49! |___________________________| !
50
51! ********************* !
52! Generic Linear Solver !
53! ********************* !
54
55 TYPE t_linear_solver
56 INTEGER, PRIVATE :: global_dimension
57 CONTAINS
58 PROCEDURE, pass :: init_solver
59 PROCEDURE, pass :: set_matrix
60 PROCEDURE, pass :: set_rhs
61 PROCEDURE, pass :: solve
62 PROCEDURE, pass :: terminate
63 PROCEDURE, NOPASS :: error
64 PROCEDURE, pass :: get_global_dim
65 END TYPE t_linear_solver
66
67! ******************** !
68! MUMPS Linear Solver !
69! ******************** !
70#ifdef MUMPS5
71 TYPE, EXTENDS(t_linear_solver) :: t_mumps_solver
72 TYPE(dmumps_struc), PRIVATE :: mumps_par
73 LOGICAL :: job_1_done = .false.
74 CONTAINS
75 PROCEDURE, pass :: init_solver_mumps
76 PROCEDURE, pass :: set_matrix_mumps
77 PROCEDURE, pass :: set_rhs_mumps
78 PROCEDURE, pass :: solve_mumps
79 PROCEDURE, pass :: terminate_mumps
80 END TYPE t_mumps_solver
81#endif
82
83! ************************* !
84! Conjugate gradient solver !
85! ************************* !
86 TYPE, EXTENDS(t_linear_solver) :: t_cg_solver
87 TYPE(t_cfs_matrix), POINTER :: mat
88 TYPE(t_vector), POINTER :: rhs
89 TYPE(t_vector) :: sol_vec, r, rnew, temp, p
90 my_real, DIMENSION(:), ALLOCATABLE :: diag
91 INTEGER :: nrhs
92 CONTAINS
93 PROCEDURE, pass :: init_solver_cg
94 PROCEDURE, pass :: set_matrix_cg
95 PROCEDURE, pass :: set_rhs_cg
96 PROCEDURE, pass :: solve_cg
97 PROCEDURE, pass :: terminate_cg
98 END TYPE t_cg_solver
99
100 CONTAINS
101
102! ._____________________________. !
103! | | !
104! | ************************ | !
105! | ** Interface routines ** | !
106! | ************************ | !
107! |_____________________________| !
108
109! **************** !
110! Generic routines !
111! **************** !
112
113! Error
114! -----
115!||====================================================================
116!|| error ../engine/share/modules/linear_solver_mod.F
117!||====================================================================
118 SUBROUTINE error()
119 print*, "ERROR"
120 END SUBROUTINE error
121
122!||====================================================================
123!|| get_global_dim ../engine/share/modules/linear_solver_mod.F
124!||--- called by ------------------------------------------------------
125!|| get_solution ../engine/share/modules/diffusion_mod.F
126!|| linear_solver_mod ../engine/share/modules/linear_solver_mod.F
127!|| solve_diffusion ../engine/share/modules/diffusion_mod.F
128!||--- uses -----------------------------------------------------
129!|| spmd_comm_world_mod ../engine/source/mpi/spmd_comm_world.F90
130!||====================================================================
131 FUNCTION get_global_dim(this)
132C-----------------------------------------------
133C I m p l i c i t T y p e s
134C-----------------------------------------------
135 USE spmd_comm_world_mod, ONLY : spmd_comm_world
136#include "implicit_f.inc"
137 class(t_linear_solver), INTENT(IN) :: this
138 INTEGER :: get_global_dim
139 get_global_dim = this%GLOBAL_DIMENSION
140 END FUNCTION get_global_dim
141
142! Solver initialization
143! ---------------------
144!||====================================================================
145!|| init_solver ../engine/share/modules/linear_solver_mod.F
146!||--- uses -----------------------------------------------------
147!|| spmd_comm_world_mod ../engine/source/mpi/spmd_comm_world.F90
148!||====================================================================
149 SUBROUTINE init_solver(this, MAT_DIM)
150C-----------------------------------------------
151C I m p l i c i t T y p e s
152C-----------------------------------------------
153 USE spmd_comm_world_mod, ONLY : spmd_comm_world
154#include "implicit_f.inc"
155 class(t_linear_solver), INTENT(INOUT) :: this
156 INTEGER, INTENT(IN) :: mat_dim
157 this%GLOBAL_DIMENSION = mat_dim
158 SELECT TYPE(this)
159#ifdef MUMPS5
160 TYPE IS (t_mumps_solver)
161 CALL this%INIT_SOLVER_MUMPS(mat_dim)
162#endif
163 TYPE IS (t_cg_solver)
164 CALL this%INIT_SOLVER_CG(mat_dim)
165 CLASS DEFAULT
166 CALL this%ERROR()
167 END SELECT
168 END SUBROUTINE init_solver
169
170! Set matrix
171! ----------
172!||====================================================================
173!|| set_matrix ../engine/share/modules/linear_solver_mod.F
174!||--- uses -----------------------------------------------------
175!|| matrix_mod ../common_source/linearalgebra/matrix_mod.F
176!|| spmd_comm_world_mod ../engine/source/mpi/spmd_comm_world.F90
177!||====================================================================
178 SUBROUTINE set_matrix(this, MAT)
179C-----------------------------------------------
180C M o d u l e s
181C-----------------------------------------------
182 USE matrix_mod
183C-----------------------------------------------
184C I m p l i c i t T y p e s
185C-----------------------------------------------
186 USE spmd_comm_world_mod, ONLY : spmd_comm_world
187#include "implicit_f.inc"
188 class(t_linear_solver), INTENT(INOUT) :: this
189 TYPE(t_cfs_matrix), INTENT(INOUT) :: MAT
190 SELECT TYPE(this)
191#ifdef MUMPS5
192 TYPE IS (t_mumps_solver)
193 CALL this%SET_MATRIX_MUMPS(mat)
194#endif
195 TYPE IS (t_cg_solver)
196 CALL this%SET_MATRIX_CG(mat)
197 CLASS DEFAULT
198 CALL this%ERROR()
199 END SELECT
200 END SUBROUTINE set_matrix
201
202! Set right hand side
203! -------------------
204!||====================================================================
205!|| set_rhs ../engine/share/modules/linear_solver_mod.F
206!||--- uses -----------------------------------------------------
207!|| spmd_comm_world_mod ../engine/source/mpi/spmd_comm_world.F90
208!|| vector_mod ../common_source/linearalgebra/vector_mod.F
209!||====================================================================
210 SUBROUTINE set_rhs(this, NRHS, RHS)
211C-----------------------------------------------
212C M o d u l e s
213C-----------------------------------------------
214 USE vector_mod
215C-----------------------------------------------
216C I m p l i c i t T y p e s
217C-----------------------------------------------
218 USE spmd_comm_world_mod, ONLY : spmd_comm_world
219#include "implicit_f.inc"
220 class(t_linear_solver), INTENT(INOUT) :: this
221 INTEGER, INTENT(IN) :: NRHS
222 TYPE(t_vector), INTENT(INOUT) :: RHS
223 SELECT TYPE(this)
224#ifdef MUMPS5
225 TYPE IS (t_mumps_solver)
226 CALL this%SET_RHS_MUMPS(nrhs, rhs)
227#endif
228 TYPE IS (t_cg_solver)
229 CALL this%SET_RHS_CG(nrhs, rhs)
230 CLASS DEFAULT
231 CALL this%ERROR()
232 END SELECT
233 END SUBROUTINE set_rhs
234
235! Solve the linear system
236! -----------------------
237!||====================================================================
238!|| solve ../engine/share/modules/linear_solver_mod.F
239!||--- uses -----------------------------------------------------
240!|| spmd_comm_world_mod ../engine/source/mpi/spmd_comm_world.F90
241!||====================================================================
242 SUBROUTINE solve(this, SOL, DIM)
243C-----------------------------------------------
244C I m p l i c i t T y p e s
245C-----------------------------------------------
246 USE spmd_comm_world_mod, ONLY : spmd_comm_world
247#include "implicit_f.inc"
248 class(t_linear_solver), INTENT(INOUT) :: this
249 INTEGER, INTENT(IN) :: DIM
250 DOUBLE PRECISION, DIMENSION(DIM), INTENT(OUT) :: SOL
251 SELECT TYPE(this)
252#ifdef MUMPS5
253 TYPE IS (t_mumps_solver)
254 CALL this%SOLVE_MUMPS(sol, dim)
255#endif
256 TYPE IS (t_cg_solver)
257 CALL this%SOLVE_CG(sol, dim)
258 CLASS DEFAULT
259 CALL this%ERROR()
260 END SELECT
261 END SUBROUTINE solve
262
263! End the solver instance
264! -----------------------
265!||====================================================================
266!|| terminate ../engine/share/modules/linear_solver_mod.F
267!||--- uses -----------------------------------------------------
268!|| spmd_comm_world_mod ../engine/source/mpi/spmd_comm_world.F90
269!||====================================================================
270 SUBROUTINE terminate(this)
271C-----------------------------------------------
272C I m p l i c i t T y p e s
273C-----------------------------------------------
274 USE spmd_comm_world_mod, ONLY : spmd_comm_world
275#include "implicit_f.inc"
276 class(t_linear_solver), INTENT(INOUT) :: this
277 SELECT TYPE(this)
278#ifdef MUMPS5
279 TYPE IS (t_mumps_solver)
280 CALL this%TERMINATE_MUMPS()
281#endif
282 TYPE IS (t_cg_solver)
283 CALL this%TERMINATE_CG()
284 CLASS DEFAULT
285 CALL this%ERROR()
286 END SELECT
287 END SUBROUTINE terminate
288
289! ************** !
290! MUMPS routines !
291! ************** !
292
293! Solver initialization
294! ---------------------
295#ifdef MUMPS5
296!||====================================================================
297!|| init_solver_mumps ../engine/share/modules/linear_solver_mod.F
298!||--- calls -----------------------------------------------------
299!||--- uses -----------------------------------------------------
300!|| spmd_comm_world_mod ../engine/source/mpi/spmd_comm_world.F90
301!||====================================================================
302 SUBROUTINE init_solver_mumps(this, MAT_DIM)
303C-----------------------------------------------
304C I m p l i c i t T y p e s
305C-----------------------------------------------
306 USE spmd_comm_world_mod, ONLY : spmd_comm_world
307#include "implicit_f.inc"
308C-----------------------------------------------
309C M e s s a g e P a s s i n g
310C-----------------------------------------------
311#include "spmd.inc"
312C-----------------------------------------------
313 class(t_mumps_solver), INTENT(INOUT) :: this
314 INTEGER, INTENT(IN) :: MAT_DIM
315 this%MUMPS_PAR%PAR = 1
316#ifdef MPI
317 this%MUMPS_PAR%COMM = spmd_comm_world
318#else
319 this%MUMPS_PAR%COMM = -1
320#endif
321 this%mumps_par%job = -1
322 this%mumps_par%sym = 0
323 call dmumps(this%mumps_par)
324
325! matrice globale
326 this%MUMPS_PAR%ICNTL(5) = 0
327! matrice distribuee
328 this%MUMPS_PAR%ICNTL(18) = 3
329! taille de la matrice
330 this%MUMPS_PAR%N = mat_dim
331! matrice sym trique d finie positive
332 this%MUMPS_PAR%SYM = 1
333! distributed rhs
334 this%MUMPS_PAR%ICNTL(20) = 10
335! distributed solution
336 !mumps_par%icntl(21) = 1
337! un-distributed solution
338 this%MUMPS_PAR%ICNTL(21) = 0
339! info on solution
340 this%MUMPS_PAR%ICNTL(11) = 1
341! dump matrice
342! this%mumps_par%write_problem = "mat"
343! verbosity
344 this%MUMPS_PAR%ICNTL(4) = 0
345
346 END SUBROUTINE init_solver_mumps
347
348! Set matrix
349! ----------
350!||====================================================================
351!|| set_matrix_mumps ../engine/share/modules/linear_solver_mod.F
352!||--- uses -----------------------------------------------------
353!|| matrix_mod ../common_source/linearalgebra/matrix_mod.F
354!|| spmd_comm_world_mod ../engine/source/mpi/spmd_comm_world.F90
355!||====================================================================
356 SUBROUTINE set_matrix_mumps(this, MAT)
357C-----------------------------------------------
358C M o d u l e s
359C-----------------------------------------------
360 USE matrix_mod
361C-----------------------------------------------
362C I m p l i c i t T y p e s
363C-----------------------------------------------
364 USE spmd_comm_world_mod, ONLY : spmd_comm_world
365#include "implicit_f.inc"
366 class(t_mumps_solver), INTENT(INOUT) :: this
367 TYPE(t_cfs_matrix), INTENT(INOUT) :: MAT
368 CALL mat%MATRIX_ASSOCIATE(this%MUMPS_PAR%IRN_LOC, this%MUMPS_PAR%JCN_LOC, this%MUMPS_PAR%A_LOC)
369 this%MUMPS_PAR%NNZ_LOC = mat%GET_DIM()
370 END SUBROUTINE set_matrix_mumps
371
372! Set right hand side
373! -------------------
374!||====================================================================
375!|| set_rhs_mumps ../engine/share/modules/linear_solver_mod.F
376!||--- uses -----------------------------------------------------
377!|| spmd_comm_world_mod ../engine/source/mpi/spmd_comm_world.F90
378!|| vector_mod ../common_source/linearalgebra/vector_mod.F
379!||====================================================================
380 SUBROUTINE set_rhs_mumps(this, NRHS, RHS)
381C-----------------------------------------------
382C M o d u l e s
383C-----------------------------------------------
384 USE vector_mod
385C-----------------------------------------------
386C I m p l i c i t T y p e s
387C-----------------------------------------------
388 USE spmd_comm_world_mod, ONLY : spmd_comm_world
389#include "implicit_f.inc"
390 class(t_mumps_solver), INTENT(INOUT) :: this
391 INTEGER, INTENT(IN) :: NRHS
392 TYPE(t_vector), INTENT(INOUT) :: RHS
393C-----------------------------------------------
394C D u m m y a r g u m e n t s
395C-----------------------------------------------
396 INTEGER :: DIM
397
398 CALL rhs%ASSOCIATE(this%MUMPS_PAR%IRHS_LOC, this%MUMPS_PAR%RHS_LOC)
399 dim = rhs%GET_DIM() / nrhs
400 this%MUMPS_PAR%NRHS = nrhs
401 this%MUMPS_PAR%NLOC_RHS = dim
402 this%MUMPS_PAR%LRHS_LOC = dim
403 END SUBROUTINE set_rhs_mumps
404
405! Solve the linear system
406! -----------------------
407!||====================================================================
408!|| solve_mumps ../engine/share/modules/linear_solver_mod.F
409!||--- calls -----------------------------------------------------
410!||--- uses -----------------------------------------------------
411!|| spmd_comm_world_mod ../engine/source/mpi/spmd_comm_world.F90
412!||====================================================================
413 SUBROUTINE solve_mumps(this, SOL, DIM)
414C-----------------------------------------------
415C I m p l i c i t T y p e s
416C-----------------------------------------------
417 USE spmd_comm_world_mod, ONLY : spmd_comm_world
418#include "implicit_f.inc"
419#include "com01_c.inc"
420#include "task_c.inc"
421C-----------------------------------------------
422C M e s s a g e P a s s i n g
423C-----------------------------------------------
424#include "spmd.inc"
425C-----------------------------------------------
426 class(t_mumps_solver), INTENT(INOUT) :: this
427 INTEGER, INTENT(IN) :: DIM
428! my_real, dimension(dim), intent(out), target :: sol
429 DOUBLE PRECISION, DIMENSION(DIM), INTENT(OUT), TARGET :: SOL
430C-----------------------------------------------
431C L o c a l V a r i a b l e s
432C-----------------------------------------------
433 INTEGER :: IERR
434 IF (dim /= this%MUMPS_PAR%N * this%MUMPS_PAR%NRHS) THEN
435 print*, "*** DIMENSION MISMATCH IN SOLUTION VECTOR"
436 RETURN
437 ELSE
438 IF (.NOT. this%JOB_1_DONE) THEN
439! analysis
440 this%MUMPS_PAR%JOB = 1
441 CALL dmumps(this%MUMPS_PAR)
442 this%JOB_1_DONE = .true.
443 ENDIF
444! factorization
445 this%MUMPS_PAR%JOB = 2
446 CALL dmumps(this%MUMPS_PAR)
447
448! solve
449 this%MUMPS_PAR%RHS => sol
450 this%MUMPS_PAR%LRHS = this%MUMPS_PAR%N
451 this%MUMPS_PAR%JOB = 3
452 CALL dmumps(this%MUMPS_PAR)
453
454#ifdef MPI
455! Sent to all procs
456 IF (nspmd > 1) THEN
457 CALL mpi_bcast(sol, dim, real, 0, spmd_comm_world, ierr)
458 ENDIF
459#endif
460 ENDIF
461 END SUBROUTINE solve_mumps
462
463! End the solver instance
464! -----------------------
465!||====================================================================
466!|| terminate_mumps ../engine/share/modules/linear_solver_mod.F
467!||--- calls -----------------------------------------------------
468!||--- uses -----------------------------------------------------
469!|| spmd_comm_world_mod ../engine/source/mpi/spmd_comm_world.F90
470!||====================================================================
471 SUBROUTINE terminate_mumps(this)
472C-----------------------------------------------
473C I m p l i c i t T y p e s
474C-----------------------------------------------
475 USE spmd_comm_world_mod, ONLY : spmd_comm_world
476#include "implicit_f.inc"
477C-----------------------------------------------
478C M e s s a g e P a s s i n g
479C-----------------------------------------------
480#include "spmd.inc"
481C-----------------------------------------------
482 class(t_mumps_solver), INTENT(INOUT) :: this
483
484 this%MUMPS_PAR%JOB = -2
485 CALL dmumps(this%MUMPS_PAR)
486
487 END SUBROUTINE terminate_mumps
488#endif
489C END OF MUMPS5 SPECIFIC CODE
490
491! *********** !
492! CG routines !
493! *********** !
494! Solver initialization
495! ---------------------
496!||====================================================================
497!|| init_solver_cg ../engine/share/modules/linear_solver_mod.F
498!||--- uses -----------------------------------------------------
499!|| spmd_comm_world_mod ../engine/source/mpi/spmd_comm_world.F90
500!||====================================================================
501 SUBROUTINE init_solver_cg(this, MAT_DIM)
502C-----------------------------------------------
503C I m p l i c i t T y p e s
504C-----------------------------------------------
505 USE spmd_comm_world_mod, ONLY : spmd_comm_world
506#include "implicit_f.inc"
507 class(t_cg_solver), INTENT(INOUT) :: this
508 INTEGER, INTENT(IN) :: MAT_DIM
509
510 CALL this%SOL_VEC%CREATE(mat_dim)
511 CALL this%R%CREATE(mat_dim)
512 CALL this%RNEW%CREATE(mat_dim)
513 CALL this%TEMP%CREATE(mat_dim)
514 CALL this%P%CREATE(mat_dim)
515 ALLOCATE(this%DIAG(mat_dim))
516
517 END SUBROUTINE init_solver_cg
518
519! Set matrix
520! ----------
521!||====================================================================
522!|| set_matrix_cg ../engine/share/modules/linear_solver_mod.F
523!||--- uses -----------------------------------------------------
524!|| matrix_mod ../common_source/linearalgebra/matrix_mod.F
525!|| spmd_comm_world_mod ../engine/source/mpi/spmd_comm_world.F90
526!||====================================================================
527 SUBROUTINE set_matrix_cg(this, MAT)
528C-----------------------------------------------
529C M o d u l e s
530C-----------------------------------------------
531 USE matrix_mod
532C-----------------------------------------------
533C I m p l i c i t T y p e s
534C-----------------------------------------------
535 USE spmd_comm_world_mod, ONLY : spmd_comm_world
536#include "implicit_f.inc"
537 class(t_cg_solver), INTENT(INOUT) :: this
538 TYPE(t_cfs_matrix), INTENT(INOUT), TARGET :: MAT
539
540 this%MAT => mat
541
542 END SUBROUTINE set_matrix_cg
543
544! Set right hand side
545! -------------------
546!||====================================================================
547!|| set_rhs_cg ../engine/share/modules/linear_solver_mod.F
548!||--- uses -----------------------------------------------------
549!|| spmd_comm_world_mod ../engine/source/mpi/spmd_comm_world.F90
550!|| vector_mod ../common_source/linearalgebra/vector_mod.F
551!||====================================================================
552 SUBROUTINE set_rhs_cg(this, NRHS, RHS)
553C-----------------------------------------------
554C M o d u l e s
555C-----------------------------------------------
556 USE vector_mod
557C-----------------------------------------------
558C I m p l i c i t T y p e s
559C-----------------------------------------------
560 USE spmd_comm_world_mod, ONLY : spmd_comm_world
561#include "implicit_f.inc"
562 class(t_cg_solver), INTENT(INOUT) :: this
563 INTEGER, INTENT(IN) :: NRHS
564 TYPE(t_vector), INTENT(INOUT), TARGET :: RHS
565
566 this%RHS => rhs
567 this%NRHS = nrhs
568
569
570 END SUBROUTINE set_rhs_cg
571
572! Solve the linear system
573! -----------------------
574!||====================================================================
575!|| solve_cg ../engine/share/modules/linear_solver_mod.F
576!||--- calls -----------------------------------------------------
577!|| prod_vec ../common_source/linearalgebra/matrix_mod.F
578!||--- uses -----------------------------------------------------
579!|| spmd_comm_world_mod ../engine/source/mpi/spmd_comm_world.F90
580!|| vector_mod ../common_source/linearalgebra/vector_mod.F
581!||====================================================================
582 SUBROUTINE solve_cg(this, SOL, DIM)
583 USE vector_mod
584C-----------------------------------------------
585C I m p l i c i t T y p e s
586C-----------------------------------------------
587 USE spmd_comm_world_mod, ONLY : spmd_comm_world
588#include "implicit_f.inc"
589#include "com01_c.inc"
590#include "task_c.inc"
591C-----------------------------------------------
592C M e s s a g e P a s s i n g
593C-----------------------------------------------
594#include "spmd.inc"
595C-----------------------------------------------
596 class(t_cg_solver), INTENT(INOUT) :: this
597 INTEGER, INTENT(IN) :: DIM
598 DOUBLE PRECISION, DIMENSION(DIM), INTENT(OUT), TARGET :: SOL
599C-----------------------------------------------
600C L o c a l V a r i a b l e s
601C-----------------------------------------------
602 INTEGER :: ITER, SYSTEM_SIZE, IRHS, II, I, J, MAT_NNZ
603 my_real :: error, norm_init
604 my_real :: alpha, beta
605 ! a sortir sous forme de parametre
606 INTEGER :: MAX_ITER
607 my_real :: tol
608
609 tol = 1.d-8
610
611 system_size = dim / this%NRHS
612 max_iter = system_size
613
614! diaginal matrix made of inverse of square root of diagonal elements of the
615! system matrix
616 mat_nnz = this%MAT%GET_DIM()
617 DO ii = 1, mat_nnz
618 i = this%MAT%IROW(ii)
619 j = this%MAT%JCOL(ii)
620 IF (i == j) THEN
621 this%DIAG(i) = one / sqrt(this%MAT%VAL(ii))
622 ENDIF
623 ENDDO
624
625
626 DO irhs = 1, this%NRHS
627 this%SOL_VEC%VAL(1:system_size) = zero
628! initialisation du solver
629 CALL prod_vec(this%MAT, this%SOL_VEC, this%TEMP)
630 this%R%VAL(1:system_size) = this%RHS%VAL(system_size * (irhs - 1) + 1 : system_size * (irhs - 1) + system_size) -
631 . this%TEMP%VAL(1:system_size)
632 this%P%VAL(1:system_size) = this%R%VAL(1:system_size)
633 norm_init = this%R%NORM()
634 error = norm_init
635 iter = 0
636 DO WHILE (iter <= max_iter .AND. error > tol)
637 iter = iter + 1
638 CALL prod_vec(this%MAT, this%P, this%TEMP)
639 alpha = dot_product(this%R%VAL(1:system_size),this% R%VAL(1:system_size)) /
640 . dot_product(this%TEMP%VAL(1:system_size), this%P%VAL(1:system_size))
641 DO ii = 1, system_size
642 this%SOL_VEC%VAL(ii) = this%SOL_VEC%VAL(ii) + alpha * this%P%VAL(ii)
643 this%RNEW%VAL(ii) = this%R%VAL(ii) - alpha * this%TEMP%VAL(ii)
644 ENDDO
645 beta = dot_product(this%RNEW%VAL(1:system_size), this%RNEW%VAL(1:system_size)) /
646 . dot_product(this%R%VAL(1:system_size), this%R%VAL(1:system_size))
647 DO ii = 1, system_size
648 this%P%VAL(ii) = this%RNEW%VAL(ii) + beta * this%P%VAL(ii)
649 this%R%VAL(ii) = this%RNEW%VAL(ii)
650 ENDDO
651 error = this%R%NORM() / norm_init
652 ENDDO
653 sol(system_size * (irhs - 1) + 1:system_size * (irhs - 1) + system_size) =
654 . this%SOL_VEC%VAL(1:system_size)
655 ENDDO
656
657
658 END SUBROUTINE solve_cg
659
660! End the solver instance
661! -----------------------
662!||====================================================================
663!|| terminate_cg ../engine/share/modules/linear_solver_mod.F
664!||--- uses -----------------------------------------------------
665!|| spmd_comm_world_mod ../engine/source/mpi/spmd_comm_world.F90
666!||====================================================================
667 SUBROUTINE terminate_cg(this)
668C-----------------------------------------------
669C I m p l i c i t T y p e s
670C-----------------------------------------------
671 USE spmd_comm_world_mod, ONLY : spmd_comm_world
672#include "implicit_f.inc"
673 class(t_cg_solver), INTENT(INOUT) :: this
674 CALL this%SOL_VEC%DESTROY()
675 CALL this%R%DESTROY()
676 CALL this%TEMP%DESTROY()
677 CALL this%P%DESTROY()
678 CALL this%RNEW%DESTROY()
679 DEALLOCATE(this%DIAG)
680 END SUBROUTINE terminate_cg
681
682 END MODULE linear_solver_mod
my_real
#define my_real
Definition cppsort.cpp:32
alpha
#define alpha
Definition eval.h:35
mpi_bcast
subroutine mpi_bcast(buffer, cnt, datatype, root, comm, ierr)
Definition mpi.f:205
linear_solver_mod
Definition linear_solver_mod.F:35
linear_solver_mod::get_global_dim
integer function get_global_dim(this)
Definition linear_solver_mod.F:132
linear_solver_mod::set_matrix_cg
subroutine set_matrix_cg(this, mat)
Definition linear_solver_mod.F:528
linear_solver_mod::set_matrix
subroutine set_matrix(this, mat)
Definition linear_solver_mod.F:179
linear_solver_mod::set_rhs_cg
subroutine set_rhs_cg(this, nrhs, rhs)
Definition linear_solver_mod.F:553
linear_solver_mod::terminate
subroutine terminate(this)
Definition linear_solver_mod.F:271
linear_solver_mod::set_rhs
subroutine set_rhs(this, nrhs, rhs)
Definition linear_solver_mod.F:211
linear_solver_mod::error
subroutine error()
Definition linear_solver_mod.F:119
linear_solver_mod::init_solver
subroutine init_solver(this, mat_dim)
Definition linear_solver_mod.F:150
linear_solver_mod::solve
subroutine solve(this, sol, dim)
Definition linear_solver_mod.F:243
linear_solver_mod::terminate_cg
subroutine terminate_cg(this)
Definition linear_solver_mod.F:668
linear_solver_mod::solve_cg
subroutine solve_cg(this, sol, dim)
Definition linear_solver_mod.F:583
linear_solver_mod::init_solver_cg
subroutine init_solver_cg(this, mat_dim)
Definition linear_solver_mod.F:502
matrix_mod
Definition matrix_mod.F:34
vector_mod
Definition vector_mod.F:37
linear_solver_mod::t_cg_solver
Definition linear_solver_mod.F:86
linear_solver_mod::t_linear_solver
Definition linear_solver_mod.F:55
matrix_mod::t_cfs_matrix
Definition matrix_mod.F:42
vector_mod::t_vector
Definition vector_mod.F:45