krylov_mod.f90

Go to the documentation of this file.

module krylov_mod

  use choral_constants
  use real_type
  use basic_tools
  use abstract_interfaces, only: rntorn
  use r1d_mod, only: reorder, mult
  use cg_mod
  use gmres_mod
  use csr_mod
  use precond_mod

  implicit none
  private

  !       %----------------------------------------%
  !       |                                        |
  !       |          PUBLIC DATA                   |
  !       |                                        |
  !       %----------------------------------------%
  !! TYPE
  public :: krylov

  !! GENERIC
  public :: clear
  public :: solve     !! TESTED 
  public :: print


  !       %----------------------------------------%
  !       |                                        |
  !       |          DERIVED TYPE                  |
  !       |                                        |
  !       %----------------------------------------%
  type krylov

     !! parameters provided by the constructor
     !!
     integer  :: type=kry_cg
     real(RP) :: tol = 1e-8_rp
     integer  :: itmax=1000
     integer  :: restart=15
     integer  :: verb=0

     !! Parameter defined by the constructor
     !!
     character(10) :: name='CG'

     !! Parameter returned after solving
     !!
     real(RP) :: res
     integer  :: iter
     integer  :: aeval=0
     logical  :: ierr = .false.

     !! Data used when solving
     !!
     real(RP), dimension(:), allocatable :: xr
     real(RP), dimension(:), allocatable :: br

   contains

     final :: krylov_clear

  end type krylov

  !       %----------------------------------------%
  !       |                                        |
  !       |          GENERIC INTERFACES            |
  !       |                                        |
  !       %----------------------------------------%
  interface print
     module procedure krylov_print
  end interface print

  interface clear
     module procedure krylov_clear
  end interface clear

  !! constructor
  interface krylov
     module procedure krylov_create
  end interface krylov

  interface solve
     !!
     !! For real arrays
     !!
     module procedure krylov_solve_raw         ! TESTED
     module procedure krylov_solve_csr         ! TESTED 
     module procedure krylov_solve_pc          ! TESTED
     module procedure krylov_solve_prec        ! TESTED

  end interface solve

  
  !       %----------------------------------------%
  !       |                                        |
  !       |          GLOBAL PARAMETERSS            |
  !       |                                        |
  !       %----------------------------------------%
  real(DP) :: time_start, time_end


contains


  subroutine krylov_clear(k)
    type(krylov), intent(inout) :: k
    
    call freemem(k%xr)
    call freemem(k%br)

  end subroutine krylov_clear



  function krylov_create(type, tol, itMax, restart, verb) &
       & result(k)
    type(krylov)                   :: k
    integer , intent(in), optional :: type
    real(RP), intent(in), optional :: tol
    integer , intent(in), optional :: itMax, restart, verb

    if (present(type)) then
       k%type = type
       select case(k%type)
       case(kry_cg)
          k%name = 'CG'
       case(kry_gmres)
          k%name = 'GMRES'
       case default
          call quit("krylov_mod: krylov_create: unknown type")
       end select
    end if

    if(present(itmax  )) k%itMax   = itmax 
    if(present(restart)) k%restart = restart
    if(present(tol    )) k%tol     = tol
    if(present(verb   )) k%verb    = verb

  end function krylov_create

  subroutine krylov_print (k)
    type(krylov), intent(in) :: k

    write(*,*)"krylov_mod      : krylov_print"
    write(*,*)"  name                           = "// trim(k%name)
    write(*,*)"  tolerance                      =", k%tol    
    write(*,*)"  itMax                          =", k%itMax  
    write(*,*)"  verbosity                      =", k%verb
    write(*,*)"  restart                        =", k%restart

  end subroutine krylov_print



  subroutine solve_start(k)
    type(krylov), intent(inout) :: k

    !! initialise clock
    time_start = clock()

    !! initialise krylov output parameters
    k%ierr  = .false.
    k%res   = 1._rp
    k%iter  = 0
    k%AEval = 0

  end subroutine solve_start

  subroutine solve_end(k)
    type(krylov), intent(inout) :: k

    if (k%iter<0)  k%ierr = .true.

    if (k%verb>1) then
       write(*,*) "  Iterations                     =", k%iter
       write(*,*) "  Performed x-->A*x              =", k%AEval
       write(*,*) "  Residual                       =", k%res
       time_end = clock()
       write(*,*) '  CPU                            =',&
            & real(time_end - time_start, sp)
    end if

  end subroutine solve_end


  subroutine krylov_solve_raw(x, k, b, A)
    real(RP), dimension(:), intent(inout) :: x
    type(krylov)          , intent(inout) :: k

    real(RP), dimension(:), intent(in)    :: b
    procedure(rntorn)                     :: A

    if (k%verb>0) write(*,*)"krylov_mod      : solve_raw      = ",&
         & trim(k%name)

    call solve_start(k)

    select case(k%type)
    case(kry_cg)

       call cg(x, k%iter, k%res, &
            &  b, a, k%tol, k%itmax, k%verb)
       k%AEval = k%iter + 1

    case(kry_gmres)
       call gmres(x, k%iter, k%AEval, k%res, &
            &      b, a, k%tol, k%itmax, k%restart, k%verb)

    case default
       call quit("krylov_mod: krylov_solve_raw:&
            & incorrect solver type")

    end select

    call solve_end(k)

  end subroutine krylov_solve_raw


  subroutine krylov_solve_csr(x, k, b, mat)

    real(RP), dimension(:), intent(inout) :: x
    type(krylov)          , intent(inout) :: k

    real(RP), dimension(:), intent(in)    :: b
    type(csr)             , intent(in)    :: mat

    call krylov_solve_raw(x, k, b, a)

  contains

    subroutine a(y, x)
      real(RP), dimension(:), intent(out) :: y
      real(RP), dimension(:), intent(in)  :: x

      call matvecprod(y, mat, x)

    end subroutine a

  end subroutine krylov_solve_csr

  subroutine krylov_solve_pc(x, k, b, A, pc)
    real(RP), dimension(:), intent(inout) :: x
    type(krylov)          , intent(inout) :: k
    real(RP), dimension(:), intent(in)    :: b
    procedure(rntorn)                     :: A, pc

    if (k%verb>0) write(*,*)"krylov_mod      : solve_pc       = ",&
         & trim(k%name)

    call solve_start(k)

    select case(k%type)
    case(kry_cg)
       call pcg(x, k%iter, k%res, &
            &  b, a, pc, k%tol, k%itmax, k%verb)
       k%AEval = k%iter + 1

    case default
       call quit("krylov_mod: krylov_solve_pc:&
            & incorrect solver type")

    end select

    call solve_end(k)

  end subroutine krylov_solve_pc

  subroutine krylov_solve_prec(x, k, b, mat, pc)

    real(RP), dimension(:), intent(inout) :: x
    type(krylov)          , intent(inout) :: k
    real(RP), dimension(:), intent(in)    :: b
    type(csr)             , intent(in)    :: mat
    type(precond)         , intent(inout) :: pc

    logical :: bool

    bool = .true.
    select case(pc%type)

    case(pc_0)
       call krylov_solve_raw(x, k, b, a)

    case(pc_jacobi)
       call krylov_solve_pc(x, k, b, a, prec_jac)

    case(pc_icc0)

       !! checks allocation and size for the auxiliary
       !! data xr and br
       if (.NOT.allocated(k%xr)) call allocmem(k%xr, mat%nl)
       if (size(k%xr,1)/=mat%nl) call allocmem(k%xr, mat%nl)
       if (.NOT.allocated(k%br)) call allocmem(k%br, mat%nl)
       if (size(k%br,1)/=mat%nl) call allocmem(k%br, mat%nl)

       !! reorder rhs and initial guess
       call reorder(k%br, pc%perm, b)
       call reorder(k%xr, pc%perm, x)

       !! solve
       call krylov_solve_pc(k%xr, k, k%br, a, prec_icc0)

       !! reorder the computed solution
       call reorder(x, pc%permInv, k%xr)

    case default
       call quit("krylov_mod: krylov_solve_prec: &
            & incorrect preconditioning type")

    end select

  contains

    subroutine a(y, x)
      real(RP), dimension(:), intent(out) :: y
      real(RP), dimension(:), intent(in)  :: x

      call matvecprod(y, mat, x)

    end subroutine a

    subroutine prec_jac(y, x)
      real(RP), dimension(:), intent(out) :: y
      real(RP), dimension(:), intent(in)  :: x

      call mult(y, pc%invD , x)

    end subroutine prec_jac

    subroutine prec_icc0(y, x)
      real(RP), dimension(:), intent(out) :: y
      real(RP), dimension(:), intent(in)  :: x

      call invllt(y, pc%LLT, x)

    end subroutine prec_icc0

  end subroutine krylov_solve_prec

end module krylov_mod