Choral
krylov_mod.f90
Go to the documentation of this file.
1 
24 
25 module krylov_mod
26 
28  use real_type
29  use basic_tools
30  use abstract_interfaces, only: rntorn
31  use r1d_mod, only: reorder, mult
32  use cg_mod
33  use gmres_mod
34  use csr_mod
35  use precond_mod
36 
37  implicit none
38  private
39 
40  ! %----------------------------------------%
41  ! | |
42  ! | PUBLIC DATA |
43  ! | |
44  ! %----------------------------------------%
45  !! TYPE
46  public :: krylov
47 
48  !! GENERIC
49  public :: clear
50  public :: solve !! TESTED
51  public :: print
52 
53 
54  ! %----------------------------------------%
55  ! | |
56  ! | DERIVED TYPE |
57  ! | |
58  ! %----------------------------------------%
63  type krylov
64 
65  !! parameters provided by the constructor
66  !!
68  integer :: type=kry_cg
70  real(RP) :: tol = 1e-8_rp
72  integer :: itmax=1000
74  integer :: restart=15
76  integer :: verb=0
77 
78  !! Parameter defined by the constructor
79  !!
81  character(10) :: name='CG'
82 
83  !! Parameter returned after solving
84  !!
86  real(RP) :: res
88  integer :: iter
90  integer :: aeval=0
92  logical :: ierr = .false.
93 
94  !! Data used when solving
95  !!
97  real(RP), dimension(:), allocatable :: xr
99  real(RP), dimension(:), allocatable :: br
100 
101  contains
102 
104  final :: krylov_clear
105 
106  end type krylov
107 
108  ! %----------------------------------------%
109  ! | |
110  ! | GENERIC INTERFACES |
111  ! | |
112  ! %----------------------------------------%
114  interface print
115  module procedure krylov_print
116  end interface print
117 
119  interface clear
120  module procedure krylov_clear
121  end interface clear
122 
123  !! constructor
124  interface krylov
125  module procedure krylov_create
126  end interface krylov
127 
129  interface solve
130  !!
131  !! For real arrays
132  !!
133  module procedure krylov_solve_raw ! TESTED
134  module procedure krylov_solve_csr ! TESTED
135  module procedure krylov_solve_pc ! TESTED
136  module procedure krylov_solve_prec ! TESTED
137 
138  end interface solve
139 
140 
141  ! %----------------------------------------%
142  ! | |
143  ! | GLOBAL PARAMETERSS |
144  ! | |
145  ! %----------------------------------------%
146  real(DP) :: time_start, time_end
147 
148 
149 contains
150 
151 
154  subroutine krylov_clear(k)
155  type(krylov), intent(inout) :: k
156 
157  call freemem(k%xr)
158  call freemem(k%br)
159 
160  end subroutine krylov_clear
161 
162 
163 
181  function krylov_create(type, tol, itMax, restart, verb) &
182  & result(k)
183  type(krylov) :: k
184  integer , intent(in), optional :: type
185  real(RP), intent(in), optional :: tol
186  integer , intent(in), optional :: itMax, restart, verb
187 
188  if (present(type)) then
189  k%type = type
190  select case(k%type)
191  case(kry_cg)
192  k%name = 'CG'
193  case(kry_gmres)
194  k%name = 'GMRES'
195  case default
196  call quit("krylov_mod: krylov_create: unknown type")
197  end select
198  end if
199 
200  if(present(itmax )) k%itMax = itmax
201  if(present(restart)) k%restart = restart
202  if(present(tol )) k%tol = tol
203  if(present(verb )) k%verb = verb
204 
205  end function krylov_create
206 
208  subroutine krylov_print (k)
209  type(krylov), intent(in) :: k
210 
211  write(*,*)"krylov_mod : krylov_print"
212  write(*,*)" name = "// trim(k%name)
213  write(*,*)" tolerance =", k%tol
214  write(*,*)" itMax =", k%itMax
215  write(*,*)" verbosity =", k%verb
216  write(*,*)" restart =", k%restart
217 
218  end subroutine krylov_print
219 
220 
221 
224  subroutine solve_start(k)
225  type(krylov), intent(inout) :: k
226 
227  !! initialise clock
228  time_start = clock()
229 
230  !! initialise krylov output parameters
231  k%ierr = .false.
232  k%res = 1._rp
233  k%iter = 0
234  k%AEval = 0
235 
236  end subroutine solve_start
237 
240  subroutine solve_end(k)
241  type(krylov), intent(inout) :: k
242 
243  if (k%iter<0) k%ierr = .true.
244 
245  if (k%verb>1) then
246  write(*,*) " Iterations =", k%iter
247  write(*,*) " Performed x-->A*x =", k%AEval
248  write(*,*) " Residual =", k%res
249  time_end = clock()
250  write(*,*) ' CPU =',&
251  & real(time_end - time_start, sp)
252  end if
253 
254  end subroutine solve_end
255 
256 
259  subroutine krylov_solve_raw(x, k, b, A)
260  real(RP), dimension(:), intent(inout) :: x
261  type(krylov) , intent(inout) :: k
262 
263  real(RP), dimension(:), intent(in) :: b
264  procedure(rntorn) :: A
265 
266  if (k%verb>0) write(*,*)"krylov_mod : solve_raw = ",&
267  & trim(k%name)
268 
269  call solve_start(k)
270 
271  select case(k%type)
272  case(kry_cg)
273 
274  call cg(x, k%iter, k%res, &
275  & b, a, k%tol, k%itmax, k%verb)
276  k%AEval = k%iter + 1
277 
278  case(kry_gmres)
279  call gmres(x, k%iter, k%AEval, k%res, &
280  & b, a, k%tol, k%itmax, k%restart, k%verb)
281 
282  case default
283  call quit("krylov_mod: krylov_solve_raw:&
284  & incorrect solver type")
285 
286  end select
287 
288  call solve_end(k)
289 
290  end subroutine krylov_solve_raw
291 
292 
297  subroutine krylov_solve_csr(x, k, b, mat)
299  real(RP), dimension(:), intent(inout) :: x
300  type(krylov) , intent(inout) :: k
301 
302  real(RP), dimension(:), intent(in) :: b
303  type(csr) , intent(in) :: mat
304 
305  call krylov_solve_raw(x, k, b, a)
306 
307  contains
308 
309  subroutine a(y, x)
310  real(RP), dimension(:), intent(out) :: y
311  real(RP), dimension(:), intent(in) :: x
312 
313  call matvecprod(y, mat, x)
314 
315  end subroutine a
316 
317  end subroutine krylov_solve_csr
318 
321  subroutine krylov_solve_pc(x, k, b, A, pc)
322  real(RP), dimension(:), intent(inout) :: x
323  type(krylov) , intent(inout) :: k
324  real(RP), dimension(:), intent(in) :: b
325  procedure(rntorn) :: A, pc
326 
327  if (k%verb>0) write(*,*)"krylov_mod : solve_pc = ",&
328  & trim(k%name)
329 
330  call solve_start(k)
331 
332  select case(k%type)
333  case(kry_cg)
334  call pcg(x, k%iter, k%res, &
335  & b, a, pc, k%tol, k%itmax, k%verb)
336  k%AEval = k%iter + 1
337 
338  case default
339  call quit("krylov_mod: krylov_solve_pc:&
340  & incorrect solver type")
341 
342  end select
343 
344  call solve_end(k)
345 
346  end subroutine krylov_solve_pc
347 
350  subroutine krylov_solve_prec(x, k, b, mat, pc)
352  real(RP), dimension(:), intent(inout) :: x
353  type(krylov) , intent(inout) :: k
354  real(RP), dimension(:), intent(in) :: b
355  type(csr) , intent(in) :: mat
356  type(precond) , intent(inout) :: pc
357 
358  logical :: bool
359 
360  bool = .true.
361  select case(pc%type)
362 
363  case(pc_0)
364  call krylov_solve_raw(x, k, b, a)
365 
366  case(pc_jacobi)
367  call krylov_solve_pc(x, k, b, a, prec_jac)
368 
369  case(pc_icc0)
370 
371  !! checks allocation and size for the auxiliary
372  !! data xr and br
373  if (.NOT.allocated(k%xr)) call allocmem(k%xr, mat%nl)
374  if (size(k%xr,1)/=mat%nl) call allocmem(k%xr, mat%nl)
375  if (.NOT.allocated(k%br)) call allocmem(k%br, mat%nl)
376  if (size(k%br,1)/=mat%nl) call allocmem(k%br, mat%nl)
377 
378  !! reorder rhs and initial guess
379  call reorder(k%br, pc%perm, b)
380  call reorder(k%xr, pc%perm, x)
381 
382  !! solve
383  call krylov_solve_pc(k%xr, k, k%br, a, prec_icc0)
384 
385  !! reorder the computed solution
386  call reorder(x, pc%permInv, k%xr)
387 
388  case default
389  call quit("krylov_mod: krylov_solve_prec: &
390  & incorrect preconditioning type")
391 
392  end select
393 
394  contains
395 
396  subroutine a(y, x)
397  real(RP), dimension(:), intent(out) :: y
398  real(RP), dimension(:), intent(in) :: x
399 
400  call matvecprod(y, mat, x)
401 
402  end subroutine a
403 
404  subroutine prec_jac(y, x)
405  real(RP), dimension(:), intent(out) :: y
406  real(RP), dimension(:), intent(in) :: x
407 
408  call mult(y, pc%invD , x)
409 
410  end subroutine prec_jac
411 
412  subroutine prec_icc0(y, x)
413  real(RP), dimension(:), intent(out) :: y
414  real(RP), dimension(:), intent(in) :: x
415 
416  call invllt(y, pc%LLT, x)
417 
418  end subroutine prec_icc0
419 
420  end subroutine krylov_solve_prec
421 
422 end module krylov_mod
423 
subroutine krylov_clear(k)
Destructor.
Definition: krylov_mod.f90:155
DEFINITION OF ABSTRACT_INTERFACES FOR THE LIBRARY CHORAL
subroutine, public pcg(x, iter, res, b, A, pc, tol, itmax, verb)
Preconditioned conjugate gradient.
Definition: cg_mod.f90:155
BASIC TOOLS
Definition: basic_tools.F90:8
deallocate memory for real(RP) arrays
Definition: real_type.F90:163
double precision function, public clock()
Returns the internal clock time
Definition: basic_tools.F90:94
subroutine solve_start(k)
Reset krylov parameters before solving.
Definition: krylov_mod.f90:225
The type precond defines preconditioning for linear systems.
Definition: precond_mod.F90:56
integer, parameter kry_gmres
GmRes linear solver.
subroutine, public cg(x, iter, res, b, A, tol, itmax, verb)
Conjugate gradient (no preconditioning)
Definition: cg_mod.f90:60
real(dp) time_start
Definition: krylov_mod.f90:146
subroutine krylov_solve_csr(x, k, b, mat)
SOLVE : KRYLOV no preconditioning.
Definition: krylov_mod.f90:298
GMRES LINEAR SOLVER
Definition: gmres_mod.f90:9
The type krylov defines the settings of a linear solver.
Definition: krylov_mod.f90:63
x solution of LLT x = b, L lower diagonal LT = transpose(L)
Definition: csr_mod.F90:172
subroutine, public quit(message)
Stop program execution, display an error messahe.
integer, parameter, public sp
simple precision for real numbers
Definition: real_type.F90:60
z = x*y
Definition: R1d_mod.F90:87
subroutine, public gmres(x, iter, nbPrd, res, b, A, tol, itmax, rst, verb)
GMRES (no preconditioning)
Definition: gmres_mod.f90:49
integer, parameter br
real(dp) time_end
Definition: krylov_mod.f90:146
REAL NUMBERS PRECISION IN CHORAL: selects simple/double/quad
Definition: real_type.F90:33
OPENMP OPERATIONS ON 1-DIMENSIONAL REAL ARRAYS
Definition: R1d_mod.F90:10
CG LINEAR SOLVER = Conjugate Gradient
Definition: cg_mod.f90:14
real(rp) function e(x, v1, v2)
CHORAL CONSTANTS
Derived type for sparse matrices in CSR format.
Definition: csr_mod.F90:62
integer, parameter pc_jacobi
Jacobi diagonal preconditioning.
subroutine krylov_solve_pc(x, k, b, A, pc)
SOLVE : KRYLOV with preconditioning.
Definition: krylov_mod.f90:322
allocate memory for real(RP) arrays
Definition: real_type.F90:158
subroutine krylov_solve_raw(x, k, b, A)
SOLVE : KRYLOV no preconditioning.
Definition: krylov_mod.f90:260
subroutine krylov_solve_prec(x, k, b, mat, pc)
KRYLOV with preconditioning defined with a prec type.
Definition: krylov_mod.f90:351
integer, parameter pc_icc0
Incomplete Cholesky, order 0.
type(krylov) function krylov_create(type, tol, itMax, restart, verb)
Constructor for the type krylov
Definition: krylov_mod.f90:183
subroutine k(y, x)
subroutine solve_end(k)
ends up solving
Definition: krylov_mod.f90:241
subroutine prec_jac(y, x)
Definition: krylov_mod.f90:405
DERIVED TYPE precond: preconditioning for linear systems
Definition: precond_mod.F90:25
subroutine prec_icc0(y, x)
Definition: krylov_mod.f90:413
DERIVED TYPE csr for sparse matrices
Definition: csr_mod.F90:13
DERIVED TYPE krylov: for the resolution of linear systems
Definition: krylov_mod.f90:25
integer, parameter pc_0
void preconditioning
array entry permutation
Definition: R1d_mod.F90:72
subroutine krylov_print(k)
Print a short description.
Definition: krylov_mod.f90:209
subroutine a(y, x)
Definition: krylov_mod.f90:310
integer, parameter kry_cg
CG linear solver.
print a short description
Definition: krylov_mod.f90:114