sparskit.F90

Go to the documentation of this file.

!! Copyright (C) 2005-2006 Heiko Appel
!!
!! This program is free software; you can redistribute it and/or modify
!! it under the terms of the GNU General Public License as published by
!! the Free Software Foundation; either version 2, or (at your option)
!! any later version.
!!
!! This program is distributed in the hope that it will be useful,
!! but WITHOUT ANY WARRANTY; without even the implied warranty of
!! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
!! GNU General Public License for more details.
!!
!! You should have received a copy of the GNU General Public License
!! along with this program; if not, write to the Free Software
!! Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
!! 02110-1301, USA.
!!
 
#include "global.h"
 
module sparskit_oct_m
  use debug_oct_m
  use global_oct_m
  use, intrinsic :: iso_fortran_env
  use messages_oct_m
  use namespace_oct_m
  use parser_oct_m
  use profiling_oct_m
 
  implicit none
 
  private
 
  integer, public, parameter ::  &
    SK_CG      =  1,             &  !< Conjugate Gradient Method
    sk_cgnr    =  2,             &  
    sk_bcg     =  3,             &  
    sk_dbcg    =  4,             &  
    sk_bcgstab =  5,             &  
    sk_tfqmr   =  6,             &  
    sk_fom     =  7,             &  
    sk_gmres   =  8,             &  
    sk_fgmres  =  9,             &  
    sk_dqgmres = 10,             &  
    sk_minval  = sk_cg,          &
    sk_maxval  = sk_dqgmres
 
  public ::                      &
    sparskit_solver_t
 
  public ::                      &
    sparskit_solver_init,        &
    dsparskit_solver_run,        &
    zsparskit_solver_run,        &
    sparskit_solver_end,         &
    sparskit_solver_copy
 
  type sparskit_solver_t
    private
    logical :: is_complex
    integer :: size
    integer :: solver_type
    integer :: krylov_size
    integer :: preconditioning
    integer :: maxiter
    integer :: used_iter
    integer :: iter_out
    real(real64)   :: residual_norm
    real(real64)   :: rel_tolerance
    real(real64)   :: abs_tolerance
 
    real(real64), allocatable :: sk_work(:), sk_b(:), sk_y(:)
 
    integer :: ipar(16)
    real(real64)   :: fpar(16)
    logical :: verbose
  end type sparskit_solver_t
 
  interface
    subroutine cg(n, rhs, sol, ipar, fpar, w)
      integer n, ipar(16)
      real*8 rhs(n), sol(n), fpar(16), w(n,*)
    end subroutine cg
 
    subroutine cgnr(n, rhs, sol, ipar, fpar, wk)
      integer n, ipar(16)
      real*8 rhs(n),sol(n),fpar(16),wk(n,*)
    end subroutine cgnr
 
    subroutine bcg(n, rhs, sol, ipar, fpar, w)
      integer n, ipar(16)
      real*8  fpar(16), rhs(n), sol(n), w(n,*)
    end subroutine bcg
 
    subroutine dbcg(n, rhs, sol, ipar, fpar, w)
      integer n, ipar(16)
      real*8 rhs(n), sol(n), fpar(16), w(n,*)
    end subroutine dbcg
 
    subroutine bcgstab(n, rhs, sol, ipar, fpar, w)
      integer n, ipar(16)
      real*8 rhs(n), sol(n), fpar(16), w(n,8)
    end subroutine bcgstab
 
    subroutine tfqmr(n, rhs, sol, ipar, fpar, w)
      integer n, ipar(16)
      real*8 rhs(n), sol(n), fpar(16), w(n,*)
    end subroutine tfqmr
 
    subroutine fom(n, rhs, sol, ipar, fpar, w)
      integer n, ipar(16)
      real*8 rhs(n), sol(n), fpar(16), w(*)
    end subroutine fom
 
    subroutine gmres(n, rhs, sol, ipar, fpar, w)
      integer n, ipar(16)
      real*8 rhs(n), sol(n), fpar(16), w(*)
    end subroutine gmres
 
    subroutine dqgmres(n, rhs, sol, ipar, fpar, w)
      integer n, ipar(16)
      real*8 rhs(n), sol(n), fpar(16), w(*)
    end subroutine dqgmres
 
    subroutine fgmres(n, rhs, sol, ipar, fpar, w)
      integer n, ipar(16)
      real*8 rhs(n), sol(n), fpar(16), w(*)
    end subroutine fgmres
  end interface
 
contains
 
  ! ---------------------------------------------------------
  subroutine sparskit_solver_init(namespace, n, sk, is_complex)
    type(namespace_t),       intent(in)  :: namespace
    integer,                 intent(in)  :: n
    type(sparskit_solver_t), intent(out) :: sk
    logical,                 intent(in)  :: is_complex
 
    integer :: workspace_size, m
 
    push_sub(sparskit_solver_init)
 
    sk%is_complex = is_complex
    ! there might be some incompatibilities to check of real/complex and available methods?
 
    !%Variable SPARSKITSolver
    !%Type integer
    !%Default sk_bcg
    !%Section Math::SPARSKIT
    !%Description
    !% Specifies what kind of linear solver will be used.
    !%Option sk_cg 1
    !% Conjugate Gradient Method
    !%Option sk_cgnr 2
    !% Conjugate Gradient Method (Normal Residual equation)
    !%Option sk_bcg 3
    !% Bi-Conjugate Gradient Method
    !%Option sk_dbcg 4
    !% BCG with partial pivoting
    !%Option sk_bcgstab 5
    !% BCG stabilized
    !%Option sk_tfqmr 6
    !% Transpose-Free Quasi-Minimum Residual method
    !%Option sk_fom 7
    !% Full Orthogonalization Method
    !%Option sk_gmres 8
    !% Generalized Minimum Residual method
    !%Option sk_fgmres 9
    !% Flexible version of Generalized Minimum Residual method
    !%Option sk_dqgmres 10
    !% Direct versions of the Quasi-Generalized Minimum Residual method
    !%End
    call parse_variable(namespace, 'SPARSKITSolver', sk_bcg, sk%solver_type)
    if (sk%solver_type < sk_minval.or.sk%solver_type > sk_maxval) then
      call messages_input_error(namespace, 'SPARSKITSolver')
    end if
 
    !%Variable SPARSKITKrylovSubspaceSize
    !%Type integer
    !%Default 15
    !%Section Math::SPARSKIT
    !%Description
    !% Some of the SPARSKIT solvers are Krylov subspace methods.
    !% This variable determines what size the solver will use
    !% for the subspace.
    !%End
    call parse_variable(namespace, 'SPARSKITKrylovSubspaceSize', 15, sk%krylov_size)
 
    ! preconditioner not implemented
    sk%preconditioning = 0
 
    !%Variable SPARSKITMaxIter
    !%Type integer
    !%Default 50000
    !%Section Math::SPARSKIT
    !%Description
    !% This variable controls the maximum number of iteration steps that
    !% will be performed by the (iterative) linear solver.
    !%End
    call parse_variable(namespace, 'SPARSKITMaxIter', 5000, sk%maxiter)
 
    !%Variable SPARSKITIterOut
    !%Type integer
    !%Default -1
    !%Section Math::SPARSKIT
    !%Description
    !% Determines how often status info of the solver is printed.
    !% If <= 0, will never be printed.
    !%End
    call parse_variable(namespace, 'SPARSKITIterOut', -1, sk%iter_out)
 
    !%Variable SPARSKITRelTolerance
    !%Type float
    !%Default 1e-8
    !%Section Math::SPARSKIT
    !%Description
    !% Some SPARSKIT solvers use a relative tolerance as a stopping criterion
    !% for the iterative solution process. This variable can be used to
    !% specify the tolerance.
    !%End
    call parse_variable(namespace, 'SPARSKITRelTolerance', 1e-8_real64, sk%rel_tolerance)
 
    !%Variable SPARSKITAbsTolerance
    !%Type float
    !%Default 1e-10
    !%Section Math::SPARSKIT
    !%Description
    !% Some SPARSKIT solvers use an absolute tolerance as a stopping criterion
    !% for the iterative solution process. This variable can be used to
    !% specify the tolerance.
    !%End
    call parse_variable(namespace, 'SPARSKITAbsTolerance', 1e-10_real64, sk%abs_tolerance)
 
    !%Variable SPARSKITVerboseSolver
    !%Type logical
    !%Default no
    !%Section Math::SPARSKIT
    !%Description
    !% When set to yes, the SPARSKIT solver will write more detailed output.
    !%End
    call parse_variable(namespace, 'SPARSKITVerboseSolver', .false., sk%verbose)
 
    ! size of the problem
    if (is_complex) then
      sk%size = 2*n
    else
      sk%size = n
    end if
 
    ! initialize workspace size
    workspace_size = 0
 
    ! Krylov subspace size
    m = sk%krylov_size
    if (mod(m, 2) /= 0) m = m + 1
 
    select case (sk%solver_type)
    case (sk_cg)
      message(1) = 'Info: SPARSKIT solver type: Conjugate Gradient Method'
      workspace_size = 5*sk%size
    case (sk_cgnr)
      message(1) = 'Info: SPARSKIT solver type: Conjugate Gradient Method (Normal Residual equation)'
      workspace_size = 5*sk%size
    case (sk_bcg)
      message(1) = 'Info: SPARSKIT solver type: Bi-Conjugate Gradient Method'
      workspace_size = 7*sk%size
    case (sk_dbcg)
      message(1) = 'Info: SPARSKIT solver type: BCG with partial pivoting'
      workspace_size = 11*sk%size
    case (sk_bcgstab)
      message(1) = 'Info: SPARSKIT solver type: BCG stabilized'
      workspace_size = 8*sk%size
    case (sk_tfqmr)
      message(1) = 'Info: SPARSKIT solver type: Transpose-Free Quasi-Minimum Residual method'
      workspace_size = 11*sk%size
    case (sk_fom)
      message(1) = 'Info: SPARSKIT solver type: Full Orthogonalization Method'
      workspace_size = (sk%size+3)*(m+2) + (m+1)*m/2
    case (sk_gmres)
      message(1) = 'Info: SPARSKIT solver type: Generalized Minimum Residual method'
      workspace_size = (sk%size+3)*(m+2) + (m+1)*m/2
    case (sk_fgmres)
      message(1) = 'Info: SPARSKIT solver type: Flexible version of Generalized Minimum Residual method'
      workspace_size =  2*sk%size*(m+1) + (m+1)*m/2 + 3*m + 2
    case (sk_dqgmres)
      message(1) = 'Info: SPARSKIT solver type: Direct versions of Quasi-Generalized Minimum Residual method'
      workspace_size = sk%size + (m+1) * (2*sk%size+4)
    case default
      write(message(1), '(a,i4,a)') "Input: '", sk%solver_type, &
        "' is not a valid SPARSKIT Solver"
      message(2) = '( SPARSKIT Solver =  cg | cgnr | bcg | dbcg | bcgstab | tfqmr | fom | gmres | fgmres | dqgmres )'
      call messages_fatal(2, namespace=namespace)
    end select
    call messages_info(1, namespace=namespace)
 
    ! Now we initialize the arrays for the reverse communication protocol
    sk%ipar = 0
    sk%fpar = 0
 
    ! A call to the solver with ipar(1) == 0 will initialize the iterative solver.
    sk%ipar(1) = 0
 
    ! Stopping criteria; use convergence test scheme 2
    sk%ipar(3) = 2
 
    sk%ipar(4) = workspace_size
    sk%ipar(5) = sk%krylov_size
 
    ! Maximum number of matrix-vector multiplies
    sk%ipar(6) = sk%maxiter
 
    ! Relative tolerance
    sk%fpar(1) = sk%rel_tolerance
 
    ! Absolute tolerance
    sk%fpar(2) = sk%abs_tolerance
 
    ! allocate and initialize work arrays
    safe_allocate(sk%sk_b(1:sk%size))
    safe_allocate(sk%sk_y(1:sk%size))
    safe_allocate(sk%sk_work(1:workspace_size))
    sk%sk_work = m_zero
    sk%sk_y    = m_zero
    sk%sk_b    = m_zero
 
    pop_sub(sparskit_solver_init)
  end subroutine sparskit_solver_init
 
  ! ---------------------------------------------------------
  subroutine sparskit_solver_end(sk)
    type(sparskit_solver_t), intent(inout) :: sk
 
    push_sub(sparskit_solver_end)
 
    safe_deallocate_a(sk%sk_b)
    safe_deallocate_a(sk%sk_y)
    safe_deallocate_a(sk%sk_work)
 
    pop_sub(sparskit_solver_end)
  end subroutine sparskit_solver_end
 
  ! ---------------------------------------------------------
  subroutine sparskit_solver_copy(sko, ski)
    type(sparskit_solver_t), intent(inout) :: sko
    type(sparskit_solver_t), intent(in)    :: ski
 
    push_sub(sparskit_solver_end)
 
    sko%is_complex      = ski%is_complex
    sko%size            = ski%size
    sko%solver_type     = ski%solver_type
    sko%krylov_size     = ski%krylov_size
    sko%preconditioning = ski%preconditioning
    sko%maxiter         = ski%maxiter
    sko%used_iter       = ski%used_iter
    sko%iter_out        = ski%iter_out
    sko%residual_norm   = ski%residual_norm
    sko%rel_tolerance   = ski%rel_tolerance
    sko%abs_tolerance   = ski%abs_tolerance
    safe_allocate_source_a(sko%sk_work, ski%sk_work)
    safe_allocate_source_a(sko%sk_b,    ski%sk_b)
    safe_allocate_source_a(sko%sk_y,    ski%sk_y)
    sko%ipar            = ski%ipar
    sko%fpar            = ski%fpar
    sko%verbose         = ski%verbose
 
    pop_sub(sparskit_solver_end)
  end subroutine sparskit_solver_copy
 
#include "undef.F90"
#include "real.F90"
#include "sparskit_inc.F90"
 
#include "undef.F90"
#include "complex.F90"
#include "sparskit_inc.F90"
 
! distdot function for dot products is defined in mesh_function_oct_m
 
end module sparskit_oct_m
 
!! Local Variables:
!! mode: f90
!! coding: utf-8
!! End: