module LightKrylov_IterativeSolvers !! This module provides some of the most important computational routines provided by !! `LightKrylov`. These include: !! !! - `eigs` : Compute the leading eigenpairs of a square linear operator \( A \). !! - `eighs` : Compute the leading eigenpairs of a symmetric positive definite !! operator \( A \). !! - `svds` : Compute the leading singular triplets of a linear operator \( A \). !! - `gmres` : Solve the linear system \( Ax = b \) using the *generalized minimum residual method*. !! - `cg` : Solve the linear system \( Ax = b \) where \( A \) is symmetric !! positive definite using the *Conjugate Gradient* method. !! !! It also provides abstract interfaces to pass user-defined solvers and preconditioners !! to `LightKrylov`. Note that these features are still experimental for now. !-------------------------------------------- !----- Fortran Standard Library ----- !-------------------------------------------- use iso_fortran_env, only: output_unit use stdlib_sorting, only: sort_index use stdlib_optval, only: optval use stdlib_io_npy, only: save_npy use stdlib_stats, only: median !------------------------------- !----- LightKrylov ----- !------------------------------- use LightKrylov_Constants use LightKrylov_Utils use LightKrylov_Logger, only: log_information, log_message, log_debug, & & type_error, stop_error, check_info use LightKrylov_Timing, only: timer => global_lightkrylov_timer, time_lightkrylov use LightKrylov_AbstractVectors use LightKrylov_AbstractLinops use LightKrylov_BaseKrylov implicit none(type, external) private character(len=*), parameter :: this_module = 'LK_Solvers' character(len=*), parameter :: this_module_long = 'LightKrylov_IterativeSolvers' character(len=*), parameter :: eigs_output = 'eigs_output.txt' public :: abstract_linear_solver_rsp public :: abstract_linear_solver_rdp public :: abstract_linear_solver_csp public :: abstract_linear_solver_cdp public :: save_eigenspectrum public :: eigs public :: eighs public :: svds public :: gmres public :: gmres_rsp public :: gmres_rdp public :: gmres_csp public :: gmres_cdp public :: fgmres public :: fgmres_rsp public :: fgmres_rdp public :: fgmres_csp public :: fgmres_cdp public :: cg public :: cg_rsp public :: cg_rdp public :: cg_csp public :: cg_cdp public :: write_results_rsp public :: write_results_rdp public :: write_results_csp public :: write_results_cdp !------------------------------------------------------ !----- ABSTRACT PRECONDITIONER DEFINITION ----- !------------------------------------------------------ type, abstract, public :: abstract_precond_rsp contains private procedure(abstract_apply_rsp), pass(self), public, deferred :: apply end type abstract_precond_rsp abstract interface subroutine abstract_apply_rsp(self, vec, iter, current_residual, target_residual) !! Abstract interface to apply a preconditioner in `LightKrylov`. import abstract_precond_rsp, abstract_vector_rsp import sp implicit none(type, external) class(abstract_precond_rsp), intent(inout) :: self !! Preconditioner. class(abstract_vector_rsp), intent(inout) :: vec !! Input/Output vector. integer, optional, intent(in) :: iter !! Current iteration number. real(sp), optional, intent(in) :: current_residual real(sp), optional, intent(in) :: target_residual end subroutine abstract_apply_rsp end interface type, abstract, public :: abstract_precond_rdp contains private procedure(abstract_apply_rdp), pass(self), public, deferred :: apply end type abstract_precond_rdp abstract interface subroutine abstract_apply_rdp(self, vec, iter, current_residual, target_residual) !! Abstract interface to apply a preconditioner in `LightKrylov`. import abstract_precond_rdp, abstract_vector_rdp import dp implicit none(type, external) class(abstract_precond_rdp), intent(inout) :: self !! Preconditioner. class(abstract_vector_rdp), intent(inout) :: vec !! Input/Output vector. integer, optional, intent(in) :: iter !! Current iteration number. real(dp), optional, intent(in) :: current_residual real(dp), optional, intent(in) :: target_residual end subroutine abstract_apply_rdp end interface type, abstract, public :: abstract_precond_csp contains private procedure(abstract_apply_csp), pass(self), public, deferred :: apply end type abstract_precond_csp abstract interface subroutine abstract_apply_csp(self, vec, iter, current_residual, target_residual) !! Abstract interface to apply a preconditioner in `LightKrylov`. import abstract_precond_csp, abstract_vector_csp import sp implicit none(type, external) class(abstract_precond_csp), intent(inout) :: self !! Preconditioner. class(abstract_vector_csp), intent(inout) :: vec !! Input/Output vector. integer, optional, intent(in) :: iter !! Current iteration number. real(sp), optional, intent(in) :: current_residual real(sp), optional, intent(in) :: target_residual end subroutine abstract_apply_csp end interface type, abstract, public :: abstract_precond_cdp contains private procedure(abstract_apply_cdp), pass(self), public, deferred :: apply end type abstract_precond_cdp abstract interface subroutine abstract_apply_cdp(self, vec, iter, current_residual, target_residual) !! Abstract interface to apply a preconditioner in `LightKrylov`. import abstract_precond_cdp, abstract_vector_cdp import dp implicit none(type, external) class(abstract_precond_cdp), intent(inout) :: self !! Preconditioner. class(abstract_vector_cdp), intent(inout) :: vec !! Input/Output vector. integer, optional, intent(in) :: iter !! Current iteration number. real(dp), optional, intent(in) :: current_residual real(dp), optional, intent(in) :: target_residual end subroutine abstract_apply_cdp end interface !-------------------------------------------------------- !----- GENERIC INTERFACE FOR LINEAR SOLVERS ----- !-------------------------------------------------------- abstract interface subroutine abstract_linear_solver_rsp(A, b, x, info, rtol, atol, preconditioner, options, transpose, meta) !! Abstract interface to use a user-defined linear solver in `LightKrylov`. import abstract_linop_rsp, abstract_vector_rsp, abstract_opts, abstract_metadata, abstract_precond_rsp, sp implicit none(type, external) class(abstract_linop_rsp), intent(inout) :: A !! Linear operator to invert. class(abstract_vector_rsp), intent(in) :: b !! Right-hand side vector. class(abstract_vector_rsp), intent(inout) :: x !! Solution vector. integer, intent(out) :: info !! Information flag. In case of successful exit, the flag should return the number of iterations required for convergence. real(sp), optional, intent(in) :: rtol !! Relative solver tolerance real(sp), optional, intent(in) :: atol !! Absolute solver tolerance class(abstract_precond_rsp), optional, intent(inout) :: preconditioner !! Preconditioner. class(abstract_opts), optional, intent(in) :: options !! Options passed to the linear solver. logical, optional, intent(in) :: transpose !! Determine whether \(\mathbf{A}\) (`.false.`) or \(\mathbf{A}^T\) (`.true.`) is being used. class(abstract_metadata), optional, intent(out) :: meta !! Metadata. end subroutine abstract_linear_solver_rsp subroutine abstract_linear_solver_rdp(A, b, x, info, rtol, atol, preconditioner, options, transpose, meta) !! Abstract interface to use a user-defined linear solver in `LightKrylov`. import abstract_linop_rdp, abstract_vector_rdp, abstract_opts, abstract_metadata, abstract_precond_rdp, dp implicit none(type, external) class(abstract_linop_rdp), intent(inout) :: A !! Linear operator to invert. class(abstract_vector_rdp), intent(in) :: b !! Right-hand side vector. class(abstract_vector_rdp), intent(inout) :: x !! Solution vector. integer, intent(out) :: info !! Information flag. In case of successful exit, the flag should return the number of iterations required for convergence. real(dp), optional, intent(in) :: rtol !! Relative solver tolerance real(dp), optional, intent(in) :: atol !! Absolute solver tolerance class(abstract_precond_rdp), optional, intent(inout) :: preconditioner !! Preconditioner. class(abstract_opts), optional, intent(in) :: options !! Options passed to the linear solver. logical, optional, intent(in) :: transpose !! Determine whether \(\mathbf{A}\) (`.false.`) or \(\mathbf{A}^T\) (`.true.`) is being used. class(abstract_metadata), optional, intent(out) :: meta !! Metadata. end subroutine abstract_linear_solver_rdp subroutine abstract_linear_solver_csp(A, b, x, info, rtol, atol, preconditioner, options, transpose, meta) !! Abstract interface to use a user-defined linear solver in `LightKrylov`. import abstract_linop_csp, abstract_vector_csp, abstract_opts, abstract_metadata, abstract_precond_csp, sp implicit none(type, external) class(abstract_linop_csp), intent(inout) :: A !! Linear operator to invert. class(abstract_vector_csp), intent(in) :: b !! Right-hand side vector. class(abstract_vector_csp), intent(inout) :: x !! Solution vector. integer, intent(out) :: info !! Information flag. In case of successful exit, the flag should return the number of iterations required for convergence. real(sp), optional, intent(in) :: rtol !! Relative solver tolerance real(sp), optional, intent(in) :: atol !! Absolute solver tolerance class(abstract_precond_csp), optional, intent(inout) :: preconditioner !! Preconditioner. class(abstract_opts), optional, intent(in) :: options !! Options passed to the linear solver. logical, optional, intent(in) :: transpose !! Determine whether \(\mathbf{A}\) (`.false.`) or \(\mathbf{A}^T\) (`.true.`) is being used. class(abstract_metadata), optional, intent(out) :: meta !! Metadata. end subroutine abstract_linear_solver_csp subroutine abstract_linear_solver_cdp(A, b, x, info, rtol, atol, preconditioner, options, transpose, meta) !! Abstract interface to use a user-defined linear solver in `LightKrylov`. import abstract_linop_cdp, abstract_vector_cdp, abstract_opts, abstract_metadata, abstract_precond_cdp, dp implicit none(type, external) class(abstract_linop_cdp), intent(inout) :: A !! Linear operator to invert. class(abstract_vector_cdp), intent(in) :: b !! Right-hand side vector. class(abstract_vector_cdp), intent(inout) :: x !! Solution vector. integer, intent(out) :: info !! Information flag. In case of successful exit, the flag should return the number of iterations required for convergence. real(dp), optional, intent(in) :: rtol !! Relative solver tolerance real(dp), optional, intent(in) :: atol !! Absolute solver tolerance class(abstract_precond_cdp), optional, intent(inout) :: preconditioner !! Preconditioner. class(abstract_opts), optional, intent(in) :: options !! Options passed to the linear solver. logical, optional, intent(in) :: transpose !! Determine whether \(\mathbf{A}\) (`.false.`) or \(\mathbf{A}^T\) (`.true.`) is being used. class(abstract_metadata), optional, intent(out) :: meta !! Metadata. end subroutine abstract_linear_solver_cdp end interface !------------------------------------------------------- !----- ----- !----- GENERALIZED MINIMUM RESIDUAL METHOD ----- !----- ----- !------------------------------------------------------- !----- Options and Metadata ----- type, extends(abstract_opts), public :: gmres_sp_opts !! GMRES options. integer :: kdim = 30 !! Dimension of the Krylov subspace (default: 30). integer :: maxiter = 10 !! Maximum number of `gmres` restarts (default: 10). logical :: if_print_metadata = .false. !! Print iteration metadata on exit (default: .false.). logical :: sanity_check = .true. !! Performs extra matrix-vector product for sanity check. end type gmres_sp_opts type, extends(abstract_metadata), public :: gmres_sp_metadata !! GMRES metadata. integer :: n_iter = 0 !! Total iteration counter. integer :: n_inner = 0 !! Number of inner iterations. integer :: n_outer = 0 !! Number of outer iterations (i.e. restarts) real(sp), dimension(:), allocatable :: res !! Residual history. logical :: converged = .false. !! Convergence flag. integer :: info = 0 !! Copy of the information flag for completeness. contains procedure, pass(self), public :: print => print_gmres_sp procedure, pass(self), public :: reset => reset_gmres_sp end type gmres_sp_metadata interface module subroutine print_gmres_sp(self, reset_counters, verbose) implicit none(type, external) class(gmres_sp_metadata), intent(inout) :: self logical, optional, intent(in) :: reset_counters !! Reset all counters to zero after printing? logical, optional, intent(in) :: verbose !! Print the full residual history? end subroutine print_gmres_sp module subroutine reset_gmres_sp(self) implicit none(type, external) class(gmres_sp_metadata), intent(inout) :: self end subroutine reset_gmres_sp end interface type, extends(abstract_opts), public :: gmres_dp_opts !! GMRES options. integer :: kdim = 30 !! Dimension of the Krylov subspace (default: 30). integer :: maxiter = 10 !! Maximum number of `gmres` restarts (default: 10). logical :: if_print_metadata = .false. !! Print iteration metadata on exit (default: .false.). logical :: sanity_check = .true. !! Performs extra matrix-vector product for sanity check. end type gmres_dp_opts type, extends(abstract_metadata), public :: gmres_dp_metadata !! GMRES metadata. integer :: n_iter = 0 !! Total iteration counter. integer :: n_inner = 0 !! Number of inner iterations. integer :: n_outer = 0 !! Number of outer iterations (i.e. restarts) real(dp), dimension(:), allocatable :: res !! Residual history. logical :: converged = .false. !! Convergence flag. integer :: info = 0 !! Copy of the information flag for completeness. contains procedure, pass(self), public :: print => print_gmres_dp procedure, pass(self), public :: reset => reset_gmres_dp end type gmres_dp_metadata interface module subroutine print_gmres_dp(self, reset_counters, verbose) implicit none(type, external) class(gmres_dp_metadata), intent(inout) :: self logical, optional, intent(in) :: reset_counters !! Reset all counters to zero after printing? logical, optional, intent(in) :: verbose !! Print the full residual history? end subroutine print_gmres_dp module subroutine reset_gmres_dp(self) implicit none(type, external) class(gmres_dp_metadata), intent(inout) :: self end subroutine reset_gmres_dp end interface !----- Interfaces for the GMRES solvers ----- interface gmres !! ### Description !! !! Solve a square linear system of equations !! !! \[ !! Ax = b !! \] !! !! using the *Generalized Minimum RESidual* (GMRES) method. !! !! **References** !! !! - Saad Y. and Schultz M. H. "GMRES: A generalized minimal residual algorithm for !! solving nonsymmetric linear systems." SIAM Journal on Scientific and Statistical !! Computing, 7(3), 1986. !! !! ### Syntax !! !! ```fortran !! call gmres(A, b, x, info [, rtol] [, atol] [, preconditioner] [, options] [, transpose]) !! ``` !! !! ### Arguments !! !! - `A` : Linear operator derived from one of the `abstract_linop` types provided by the !! `AbstractLinops` module. It is an `intent(inout)` argument. !! !! - `b` : Right-hand side vector derived from one the `abstract_vector` types provided !! by the `AbstractVectors` module. It needs to have the same type and kind as `A`. !! It is an `intent(in)` argument. !! !! - `x` : On entry, initial guess for the solution. On exit, the solution computed by !! gmres. It is a vector derived from one the `abstract_vector` types provided by the !! `AbstractVectors` module. It needs to have the same type and kind as `A`. It is !! an `intent(inout)` argument. !! !! - `info` : `integer` information flag. !! !! - `rtol` (optional) : `real` relative tolerance for the solver. !! !! - `atol` (optional) : `real` absolute tolerance for the solver. !! !! - `preconditioner` (optional) : Right preconditioner used to solve the system. It needs !! to be consistent with the `abstract_preconditioner` interface. !! It is an optional `intent(in)` argument. !! !! - `options` (optional) : Container for the gmres options given by the `gmres_opts` type. !! It is an optional `intent(in)` argument. !! !! - `transpose` (optional): `logical` flag controlling whether \( Ax = b\) or \( A^H x = b \) is being solved. !! !! - `meta` (optional) : Container for the gmres metadata. It needs to be of type `gmres_metadata`. ! --- Interface for GMRES with Abstract Linops and Abstract Vectors --- module subroutine gmres_rsp(A, b, x, info, rtol, atol, preconditioner, options, transpose, meta) implicit none(type, external) class(abstract_linop_rsp), intent(inout) :: A !! Linear operator to be inverted. class(abstract_vector_rsp), intent(in) :: b !! Right-hand side vector. class(abstract_vector_rsp), intent(inout) :: x !! Solution vector. integer, intent(out) :: info !! Information flag. real(sp), optional, intent(in) :: rtol !! Relative solver tolerance real(sp), optional, intent(in) :: atol !! Absolute solver tolerance class(abstract_precond_rsp), optional, intent(inout) :: preconditioner !! Preconditioner (optional). class(abstract_opts), optional, intent(in) :: options !! GMRES options. logical, optional, intent(in) :: transpose !! Whether \(\mathbf{A}\) or \(\mathbf{A}^H\) is being used. class(abstract_metadata), optional, intent(out) :: meta !! Metadata. end subroutine gmres_rsp module subroutine dense_gmres_rsp(A, b, x, info, rtol, atol, preconditioner, options, transpose, meta) implicit none(type, external) real(sp), intent(in) :: A(:, :) !! Linear operator to be inverted. real(sp), intent(in) :: b(:) !! Right-hand side vector. real(sp), intent(inout) :: x(:) !! Solution vector. integer, intent(out) :: info !! Information flag. real(sp), optional, intent(in) :: rtol !! Relative solver tolerance real(sp), optional, intent(in) :: atol !! Absolute solver tolerance class(abstract_precond_rsp), optional, intent(inout) :: preconditioner !! Preconditioner (optional). class(abstract_opts), optional, intent(in) :: options !! GMRES options. logical, optional, intent(in) :: transpose !! Whether \(\mathbf{A}\) or \(\mathbf{A}^H\) is being used. class(abstract_metadata), optional, intent(out) :: meta !! Metadata. end subroutine dense_gmres_rsp module subroutine gmres_rdp(A, b, x, info, rtol, atol, preconditioner, options, transpose, meta) implicit none(type, external) class(abstract_linop_rdp), intent(inout) :: A !! Linear operator to be inverted. class(abstract_vector_rdp), intent(in) :: b !! Right-hand side vector. class(abstract_vector_rdp), intent(inout) :: x !! Solution vector. integer, intent(out) :: info !! Information flag. real(dp), optional, intent(in) :: rtol !! Relative solver tolerance real(dp), optional, intent(in) :: atol !! Absolute solver tolerance class(abstract_precond_rdp), optional, intent(inout) :: preconditioner !! Preconditioner (optional). class(abstract_opts), optional, intent(in) :: options !! GMRES options. logical, optional, intent(in) :: transpose !! Whether \(\mathbf{A}\) or \(\mathbf{A}^H\) is being used. class(abstract_metadata), optional, intent(out) :: meta !! Metadata. end subroutine gmres_rdp module subroutine dense_gmres_rdp(A, b, x, info, rtol, atol, preconditioner, options, transpose, meta) implicit none(type, external) real(dp), intent(in) :: A(:, :) !! Linear operator to be inverted. real(dp), intent(in) :: b(:) !! Right-hand side vector. real(dp), intent(inout) :: x(:) !! Solution vector. integer, intent(out) :: info !! Information flag. real(dp), optional, intent(in) :: rtol !! Relative solver tolerance real(dp), optional, intent(in) :: atol !! Absolute solver tolerance class(abstract_precond_rdp), optional, intent(inout) :: preconditioner !! Preconditioner (optional). class(abstract_opts), optional, intent(in) :: options !! GMRES options. logical, optional, intent(in) :: transpose !! Whether \(\mathbf{A}\) or \(\mathbf{A}^H\) is being used. class(abstract_metadata), optional, intent(out) :: meta !! Metadata. end subroutine dense_gmres_rdp module subroutine gmres_csp(A, b, x, info, rtol, atol, preconditioner, options, transpose, meta) implicit none(type, external) class(abstract_linop_csp), intent(inout) :: A !! Linear operator to be inverted. class(abstract_vector_csp), intent(in) :: b !! Right-hand side vector. class(abstract_vector_csp), intent(inout) :: x !! Solution vector. integer, intent(out) :: info !! Information flag. real(sp), optional, intent(in) :: rtol !! Relative solver tolerance real(sp), optional, intent(in) :: atol !! Absolute solver tolerance class(abstract_precond_csp), optional, intent(inout) :: preconditioner !! Preconditioner (optional). class(abstract_opts), optional, intent(in) :: options !! GMRES options. logical, optional, intent(in) :: transpose !! Whether \(\mathbf{A}\) or \(\mathbf{A}^H\) is being used. class(abstract_metadata), optional, intent(out) :: meta !! Metadata. end subroutine gmres_csp module subroutine dense_gmres_csp(A, b, x, info, rtol, atol, preconditioner, options, transpose, meta) implicit none(type, external) complex(sp), intent(in) :: A(:, :) !! Linear operator to be inverted. complex(sp), intent(in) :: b(:) !! Right-hand side vector. complex(sp), intent(inout) :: x(:) !! Solution vector. integer, intent(out) :: info !! Information flag. real(sp), optional, intent(in) :: rtol !! Relative solver tolerance real(sp), optional, intent(in) :: atol !! Absolute solver tolerance class(abstract_precond_csp), optional, intent(inout) :: preconditioner !! Preconditioner (optional). class(abstract_opts), optional, intent(in) :: options !! GMRES options. logical, optional, intent(in) :: transpose !! Whether \(\mathbf{A}\) or \(\mathbf{A}^H\) is being used. class(abstract_metadata), optional, intent(out) :: meta !! Metadata. end subroutine dense_gmres_csp module subroutine gmres_cdp(A, b, x, info, rtol, atol, preconditioner, options, transpose, meta) implicit none(type, external) class(abstract_linop_cdp), intent(inout) :: A !! Linear operator to be inverted. class(abstract_vector_cdp), intent(in) :: b !! Right-hand side vector. class(abstract_vector_cdp), intent(inout) :: x !! Solution vector. integer, intent(out) :: info !! Information flag. real(dp), optional, intent(in) :: rtol !! Relative solver tolerance real(dp), optional, intent(in) :: atol !! Absolute solver tolerance class(abstract_precond_cdp), optional, intent(inout) :: preconditioner !! Preconditioner (optional). class(abstract_opts), optional, intent(in) :: options !! GMRES options. logical, optional, intent(in) :: transpose !! Whether \(\mathbf{A}\) or \(\mathbf{A}^H\) is being used. class(abstract_metadata), optional, intent(out) :: meta !! Metadata. end subroutine gmres_cdp module subroutine dense_gmres_cdp(A, b, x, info, rtol, atol, preconditioner, options, transpose, meta) implicit none(type, external) complex(dp), intent(in) :: A(:, :) !! Linear operator to be inverted. complex(dp), intent(in) :: b(:) !! Right-hand side vector. complex(dp), intent(inout) :: x(:) !! Solution vector. integer, intent(out) :: info !! Information flag. real(dp), optional, intent(in) :: rtol !! Relative solver tolerance real(dp), optional, intent(in) :: atol !! Absolute solver tolerance class(abstract_precond_cdp), optional, intent(inout) :: preconditioner !! Preconditioner (optional). class(abstract_opts), optional, intent(in) :: options !! GMRES options. logical, optional, intent(in) :: transpose !! Whether \(\mathbf{A}\) or \(\mathbf{A}^H\) is being used. class(abstract_metadata), optional, intent(out) :: meta !! Metadata. end subroutine dense_gmres_cdp end interface !---------------------------------------------------------------- !----- ----- !----- FLEXIBLE GENERALIZED MINIMUM RESIDUAL METHOD ----- !----- ----- !---------------------------------------------------------------- !----- Options and Metadata ----- type, extends(abstract_opts), public :: fgmres_sp_opts !! FGMRES options. integer :: kdim = 30 !! Dimension of the Krylov subspace (default: 30). integer :: maxiter = 10 !! Maximum number of `fgmres` restarts (default: 10). logical :: if_print_metadata = .false. !! Print iteration metadata on exit (default: .false.). logical :: sanity_check = .true. !! Performs extra matrix-vector product for sanity check. end type fgmres_sp_opts type, extends(abstract_metadata), public :: fgmres_sp_metadata !! FGMRES metadata. integer :: n_iter = 0 !! Total iteration counter. integer :: n_inner = 0 !! Number of inner iterations. integer :: n_outer = 0 !! Number of outer iterations (i.e. restarts) real(sp), dimension(:), allocatable :: res !! Residual history. logical :: converged = .false. !! Convergence flag. integer :: info = 0 !! Copy of the information flag for completeness. contains procedure, pass(self), public :: print => print_fgmres_sp procedure, pass(self), public :: reset => reset_fgmres_sp end type fgmres_sp_metadata interface module subroutine print_fgmres_sp(self, reset_counters, verbose) implicit none(type, external) class(fgmres_sp_metadata), intent(inout) :: self logical, optional, intent(in) :: reset_counters !! Reset all counters to zero after printing? logical, optional, intent(in) :: verbose !! Print the full residual history? end subroutine print_fgmres_sp module subroutine reset_fgmres_sp(self) implicit none(type, external) class(fgmres_sp_metadata), intent(inout) :: self end subroutine reset_fgmres_sp end interface type, extends(abstract_opts), public :: fgmres_dp_opts !! FGMRES options. integer :: kdim = 30 !! Dimension of the Krylov subspace (default: 30). integer :: maxiter = 10 !! Maximum number of `fgmres` restarts (default: 10). logical :: if_print_metadata = .false. !! Print iteration metadata on exit (default: .false.). logical :: sanity_check = .true. !! Performs extra matrix-vector product for sanity check. end type fgmres_dp_opts type, extends(abstract_metadata), public :: fgmres_dp_metadata !! FGMRES metadata. integer :: n_iter = 0 !! Total iteration counter. integer :: n_inner = 0 !! Number of inner iterations. integer :: n_outer = 0 !! Number of outer iterations (i.e. restarts) real(dp), dimension(:), allocatable :: res !! Residual history. logical :: converged = .false. !! Convergence flag. integer :: info = 0 !! Copy of the information flag for completeness. contains procedure, pass(self), public :: print => print_fgmres_dp procedure, pass(self), public :: reset => reset_fgmres_dp end type fgmres_dp_metadata interface module subroutine print_fgmres_dp(self, reset_counters, verbose) implicit none(type, external) class(fgmres_dp_metadata), intent(inout) :: self logical, optional, intent(in) :: reset_counters !! Reset all counters to zero after printing? logical, optional, intent(in) :: verbose !! Print the full residual history? end subroutine print_fgmres_dp module subroutine reset_fgmres_dp(self) implicit none(type, external) class(fgmres_dp_metadata), intent(inout) :: self end subroutine reset_fgmres_dp end interface !----- Interfaces for the FGMRES solvers ----- interface fgmres !! ### Description !! !! Solve a square linear system of equations !! !! \[ !! Ax = b !! \] !! !! using the *Flexible Generalized Minimum RESidual* (FGMRES) method. !! !! **References** !! !! - Saad Y. and Schultz M. H. "GMRES: A generalized minimal residual algorithm for !! solving nonsymmetric linear systems." SIAM Journal on Scientific and Statistical !! Computing, 7(3), 1986. !! !! ### Syntax !! !! ```fortran !! call fgmres(A, b, x, info [, rtol] [, atol] [, preconditioner] [, options] [, transpose]) !! ``` !! !! ### Arguments !! !! - `A` : Linear operator derived from one of the `abstract_linop` types provided by the !! `AbstractLinops` module. It is an `intent(inout)` argument. !! !! - `b` : Right-hand side vector derived from one the `abstract_vector` types provided !! by the `AbstractVectors` module. It needs to have the same type and kind as `A`. !! It is an `intent(in)` argument. !! !! - `x` : On entry, initial guess for the solution. On exit, the solution computed by !! gmres. It is a vector derived from one the `abstract_vector` types provided by the !! `AbstractVectors` module. It needs to have the same type and kind as `A`. It is !! an `intent(inout)` argument. !! !! - `info` : `integer` information flag. !! !! - `rtol` (optional) : `real` relative tolerance for the solver. !! !! - `atol` (optional) : `real` absolute tolerance for the solver. !! !! - `preconditioner` (optional) : Right preconditioner used to solve the system. It needs to be !! consistent with the `abstract_preconditioner` interface. !! It is an optional `intent(in)` argument. !! !! - `options` (optional) : Container for the gmres options given by the `gmres_opts` type. !! It is an optional `intent(in)` argument. !! !! - `transpose` (optional): `logical` flag controlling whether \( Ax = b\) or !! \( A^H x = b \) is being solved. module subroutine fgmres_rsp(A, b, x, info, rtol, atol, preconditioner, options, transpose, meta) implicit none(type, external) class(abstract_linop_rsp), intent(inout) :: A !! Linear operator to be inverted. class(abstract_vector_rsp), intent(in) :: b !! Right-hand side vector. class(abstract_vector_rsp), intent(inout) :: x !! Solution vector. integer, intent(out) :: info !! Information flag. real(sp), optional, intent(in) :: rtol !! Relative solver tolerance real(sp), optional, intent(in) :: atol !! Absolute solver tolerance class(abstract_precond_rsp), optional, intent(inout) :: preconditioner !! Preconditioner (optional). class(abstract_opts), optional, intent(in) :: options !! GMRES options. logical, optional, intent(in) :: transpose !! Whether \(\mathbf{A}\) or \(\mathbf{A}^H\) is being used. class(abstract_metadata), optional, intent(out) :: meta !! Metadata. end subroutine fgmres_rsp module subroutine dense_fgmres_rsp(A, b, x, info, rtol, atol, preconditioner, options, transpose, meta) implicit none(type, external) real(sp), intent(in) :: A(:, :) !! Linear operator to be inverted. real(sp), intent(in) :: b(:) !! Right-hand side vector. real(sp), intent(inout) :: x(:) !! Solution vector. integer, intent(out) :: info !! Information flag. real(sp), optional, intent(in) :: rtol !! Relative solver tolerance real(sp), optional, intent(in) :: atol !! Absolute solver tolerance class(abstract_precond_rsp), optional, intent(inout) :: preconditioner !! Preconditioner (optional). class(abstract_opts), optional, intent(in) :: options !! GMRES options. logical, optional, intent(in) :: transpose !! Whether \(\mathbf{A}\) or \(\mathbf{A}^H\) is being used. class(abstract_metadata), optional, intent(out) :: meta !! Metadata. end subroutine dense_fgmres_rsp module subroutine fgmres_rdp(A, b, x, info, rtol, atol, preconditioner, options, transpose, meta) implicit none(type, external) class(abstract_linop_rdp), intent(inout) :: A !! Linear operator to be inverted. class(abstract_vector_rdp), intent(in) :: b !! Right-hand side vector. class(abstract_vector_rdp), intent(inout) :: x !! Solution vector. integer, intent(out) :: info !! Information flag. real(dp), optional, intent(in) :: rtol !! Relative solver tolerance real(dp), optional, intent(in) :: atol !! Absolute solver tolerance class(abstract_precond_rdp), optional, intent(inout) :: preconditioner !! Preconditioner (optional). class(abstract_opts), optional, intent(in) :: options !! GMRES options. logical, optional, intent(in) :: transpose !! Whether \(\mathbf{A}\) or \(\mathbf{A}^H\) is being used. class(abstract_metadata), optional, intent(out) :: meta !! Metadata. end subroutine fgmres_rdp module subroutine dense_fgmres_rdp(A, b, x, info, rtol, atol, preconditioner, options, transpose, meta) implicit none(type, external) real(dp), intent(in) :: A(:, :) !! Linear operator to be inverted. real(dp), intent(in) :: b(:) !! Right-hand side vector. real(dp), intent(inout) :: x(:) !! Solution vector. integer, intent(out) :: info !! Information flag. real(dp), optional, intent(in) :: rtol !! Relative solver tolerance real(dp), optional, intent(in) :: atol !! Absolute solver tolerance class(abstract_precond_rdp), optional, intent(inout) :: preconditioner !! Preconditioner (optional). class(abstract_opts), optional, intent(in) :: options !! GMRES options. logical, optional, intent(in) :: transpose !! Whether \(\mathbf{A}\) or \(\mathbf{A}^H\) is being used. class(abstract_metadata), optional, intent(out) :: meta !! Metadata. end subroutine dense_fgmres_rdp module subroutine fgmres_csp(A, b, x, info, rtol, atol, preconditioner, options, transpose, meta) implicit none(type, external) class(abstract_linop_csp), intent(inout) :: A !! Linear operator to be inverted. class(abstract_vector_csp), intent(in) :: b !! Right-hand side vector. class(abstract_vector_csp), intent(inout) :: x !! Solution vector. integer, intent(out) :: info !! Information flag. real(sp), optional, intent(in) :: rtol !! Relative solver tolerance real(sp), optional, intent(in) :: atol !! Absolute solver tolerance class(abstract_precond_csp), optional, intent(inout) :: preconditioner !! Preconditioner (optional). class(abstract_opts), optional, intent(in) :: options !! GMRES options. logical, optional, intent(in) :: transpose !! Whether \(\mathbf{A}\) or \(\mathbf{A}^H\) is being used. class(abstract_metadata), optional, intent(out) :: meta !! Metadata. end subroutine fgmres_csp module subroutine dense_fgmres_csp(A, b, x, info, rtol, atol, preconditioner, options, transpose, meta) implicit none(type, external) complex(sp), intent(in) :: A(:, :) !! Linear operator to be inverted. complex(sp), intent(in) :: b(:) !! Right-hand side vector. complex(sp), intent(inout) :: x(:) !! Solution vector. integer, intent(out) :: info !! Information flag. real(sp), optional, intent(in) :: rtol !! Relative solver tolerance real(sp), optional, intent(in) :: atol !! Absolute solver tolerance class(abstract_precond_csp), optional, intent(inout) :: preconditioner !! Preconditioner (optional). class(abstract_opts), optional, intent(in) :: options !! GMRES options. logical, optional, intent(in) :: transpose !! Whether \(\mathbf{A}\) or \(\mathbf{A}^H\) is being used. class(abstract_metadata), optional, intent(out) :: meta !! Metadata. end subroutine dense_fgmres_csp module subroutine fgmres_cdp(A, b, x, info, rtol, atol, preconditioner, options, transpose, meta) implicit none(type, external) class(abstract_linop_cdp), intent(inout) :: A !! Linear operator to be inverted. class(abstract_vector_cdp), intent(in) :: b !! Right-hand side vector. class(abstract_vector_cdp), intent(inout) :: x !! Solution vector. integer, intent(out) :: info !! Information flag. real(dp), optional, intent(in) :: rtol !! Relative solver tolerance real(dp), optional, intent(in) :: atol !! Absolute solver tolerance class(abstract_precond_cdp), optional, intent(inout) :: preconditioner !! Preconditioner (optional). class(abstract_opts), optional, intent(in) :: options !! GMRES options. logical, optional, intent(in) :: transpose !! Whether \(\mathbf{A}\) or \(\mathbf{A}^H\) is being used. class(abstract_metadata), optional, intent(out) :: meta !! Metadata. end subroutine fgmres_cdp module subroutine dense_fgmres_cdp(A, b, x, info, rtol, atol, preconditioner, options, transpose, meta) implicit none(type, external) complex(dp), intent(in) :: A(:, :) !! Linear operator to be inverted. complex(dp), intent(in) :: b(:) !! Right-hand side vector. complex(dp), intent(inout) :: x(:) !! Solution vector. integer, intent(out) :: info !! Information flag. real(dp), optional, intent(in) :: rtol !! Relative solver tolerance real(dp), optional, intent(in) :: atol !! Absolute solver tolerance class(abstract_precond_cdp), optional, intent(inout) :: preconditioner !! Preconditioner (optional). class(abstract_opts), optional, intent(in) :: options !! GMRES options. logical, optional, intent(in) :: transpose !! Whether \(\mathbf{A}\) or \(\mathbf{A}^H\) is being used. class(abstract_metadata), optional, intent(out) :: meta !! Metadata. end subroutine dense_fgmres_cdp end interface !--------------------------------------------- !----- ----- !----- CONJUGATE GRADIENT METHOD ----- !----- ----- !--------------------------------------------- !----- Options and Metadata ----- type, extends(abstract_opts), public :: cg_sp_opts !! Conjugate gradient options. integer :: maxiter = 100 !! Maximum number of `cg` iterations (default: 100). logical :: if_print_metadata = .false. !! Print interation metadata on exit (default = .false.) end type cg_sp_opts type, extends(abstract_metadata), public :: cg_sp_metadata !! Conjugate gradient metadata. integer :: n_iter = 0 !! Iteration counter real(sp), dimension(:), allocatable :: res !! Residual history logical :: converged = .false. !! Convergence flag integer :: info = 0 !! Copy of the information flag for completeness contains procedure, pass(self), public :: print => print_cg_sp procedure, pass(self), public :: reset => reset_cg_sp end type cg_sp_metadata interface module subroutine print_cg_sp(self, reset_counters, verbose) implicit none(type, external) class(cg_sp_metadata), intent(inout) :: self logical, optional, intent(in) :: reset_counters !! Reset all counters to zero after printing? logical, optional, intent(in) :: verbose !! Print the residual full residual history? end subroutine print_cg_sp module subroutine reset_cg_sp(self) implicit none(type, external) class(cg_sp_metadata), intent(inout) :: self end subroutine reset_cg_sp end interface type, extends(abstract_opts), public :: cg_dp_opts !! Conjugate gradient options. integer :: maxiter = 100 !! Maximum number of `cg` iterations (default: 100). logical :: if_print_metadata = .false. !! Print interation metadata on exit (default = .false.) end type cg_dp_opts type, extends(abstract_metadata), public :: cg_dp_metadata !! Conjugate gradient metadata. integer :: n_iter = 0 !! Iteration counter real(dp), dimension(:), allocatable :: res !! Residual history logical :: converged = .false. !! Convergence flag integer :: info = 0 !! Copy of the information flag for completeness contains procedure, pass(self), public :: print => print_cg_dp procedure, pass(self), public :: reset => reset_cg_dp end type cg_dp_metadata interface module subroutine print_cg_dp(self, reset_counters, verbose) implicit none(type, external) class(cg_dp_metadata), intent(inout) :: self logical, optional, intent(in) :: reset_counters !! Reset all counters to zero after printing? logical, optional, intent(in) :: verbose !! Print the residual full residual history? end subroutine print_cg_dp module subroutine reset_cg_dp(self) implicit none(type, external) class(cg_dp_metadata), intent(inout) :: self end subroutine reset_cg_dp end interface !----- Interfaces for the Conjugate Gradient solvers ----- interface cg !! ### Description !! !! Given a symmetric (positive definite) matrix \( A \), solves the linear system !! !! \[ !! Ax = b !! \] !! !! using the *Conjugate Gradient* method. !! !! **References** !! !! - Hestenes, M. R., and Stiefel, E. (1952). "Methods of Conjugate Gradients for Solving !! Linear Systems," Journal of Research of the National Bureau of Standards, !! 49(6), 409–436. !! !! ### Syntax !! !! ```fortran !! call cg(A, b, x, info [, rtol] [, atol] [, preconditioner] [, options]) !! ``` !! !! ### Arguments !! !! - `A` : Linear operator derived from one of the `abstract_sym_linop` or `abstract_hermitian_linop` !! types provided by the `AbstractLinops` module. It is an `intent(inout)` argument. !! !! - `b` : Right-hand side vector derived from one the `abstract_vector` types provided !! by the `AbstractVectors` module. It needs to have the same type and kind as `A`. !! It is an `intent(in)` argument. !! !! - `x` : On entry, initial guess for the solution. On exit, the solution computed by !! cg. It is a vector derived from one the `abstract_vector` types provided by the !! `AbstractVectors` module. It needs to have the same type and kind as `A`. It is !! an `intent(inout)` argument. !! !! - `info` : `integer` information flag. !! !! - `rtol` (optional) : `real` relative tolerance for the solver. !! !! - `atol` (optional) : `real` absolute tolerance for the solver. !! !! - `preconditioner` (optional) : Right preconditioner used to solve the system. It needs to be !! consistent with the `abstract_preconditioner` interface. !! It is an optional `intent(in)` argument. !! !! - `options` (optional) : Container for the gmres options given by the `cg_opts` type. !! It is an optional `intent(in)` argument. module subroutine cg_rsp(A, b, x, info, rtol, atol, preconditioner, options, meta) implicit none(type, external) class(abstract_sym_linop_rsp), intent(inout) :: A !! Linear operator to be inverted. class(abstract_vector_rsp), intent(in) :: b !! Right-hand side vector. class(abstract_vector_rsp), intent(inout) :: x !! Solution vector. integer, intent(out) :: info !! Information flag. real(sp), optional, intent(in) :: rtol !! Relative solver tolerance real(sp), optional, intent(in) :: atol !! Absolute solver tolerance class(abstract_precond_rsp), optional, intent(inout) :: preconditioner !! Preconditioner (not yet supported). type(cg_sp_opts), optional, intent(in) :: options !! Options for the conjugate gradient solver. class(abstract_metadata), optional, intent(out) :: meta !! Metadata. end subroutine cg_rsp module subroutine cg_rdp(A, b, x, info, rtol, atol, preconditioner, options, meta) implicit none(type, external) class(abstract_sym_linop_rdp), intent(inout) :: A !! Linear operator to be inverted. class(abstract_vector_rdp), intent(in) :: b !! Right-hand side vector. class(abstract_vector_rdp), intent(inout) :: x !! Solution vector. integer, intent(out) :: info !! Information flag. real(dp), optional, intent(in) :: rtol !! Relative solver tolerance real(dp), optional, intent(in) :: atol !! Absolute solver tolerance class(abstract_precond_rdp), optional, intent(inout) :: preconditioner !! Preconditioner (not yet supported). type(cg_dp_opts), optional, intent(in) :: options !! Options for the conjugate gradient solver. class(abstract_metadata), optional, intent(out) :: meta !! Metadata. end subroutine cg_rdp module subroutine cg_csp(A, b, x, info, rtol, atol, preconditioner, options, meta) implicit none(type, external) class(abstract_hermitian_linop_csp), intent(inout) :: A !! Linear operator to be inverted. class(abstract_vector_csp), intent(in) :: b !! Right-hand side vector. class(abstract_vector_csp), intent(inout) :: x !! Solution vector. integer, intent(out) :: info !! Information flag. real(sp), optional, intent(in) :: rtol !! Relative solver tolerance real(sp), optional, intent(in) :: atol !! Absolute solver tolerance class(abstract_precond_csp), optional, intent(inout) :: preconditioner !! Preconditioner (not yet supported). type(cg_sp_opts), optional, intent(in) :: options !! Options for the conjugate gradient solver. class(abstract_metadata), optional, intent(out) :: meta !! Metadata. end subroutine cg_csp module subroutine cg_cdp(A, b, x, info, rtol, atol, preconditioner, options, meta) implicit none(type, external) class(abstract_hermitian_linop_cdp), intent(inout) :: A !! Linear operator to be inverted. class(abstract_vector_cdp), intent(in) :: b !! Right-hand side vector. class(abstract_vector_cdp), intent(inout) :: x !! Solution vector. integer, intent(out) :: info !! Information flag. real(dp), optional, intent(in) :: rtol !! Relative solver tolerance real(dp), optional, intent(in) :: atol !! Absolute solver tolerance class(abstract_precond_cdp), optional, intent(inout) :: preconditioner !! Preconditioner (not yet supported). type(cg_dp_opts), optional, intent(in) :: options !! Options for the conjugate gradient solver. class(abstract_metadata), optional, intent(out) :: meta !! Metadata. end subroutine cg_cdp end interface !------------------------------------------ !----- ----- !----- SINGULAR VALUE SOLVERS ----- !----- ----- !------------------------------------------ interface svds !! ### Description !! !! Computes the leading singular triplets of an arbitrary linear operator \(A\) !! using the Lanczos iterative process. Given a linear operator \(A\), it finds !! the leading singular values and singular vectors such that: !! !! \[ !! \begin{aligned} !! Av & = \sigma u \\ !! A^H u & = \sigma v. !! \end{aligned} !! \] !! !! The subspaces \(U\) and \(V\) are constructed via Lanczos factorization, resulting in !! a bidiagonal matrix \(B\). The singular values of \(A\) are approximated by those of !! \(B\) and the singular vectors are computed accordingly. !! !! **References** !! !! - Golub, G. H., & Kahan, W. (1965). "Calculating the Singular Values and !! Pseudo-Inverse of a Matrix." !! - Baglama, J., & Reichel, L. (2005). "Augmented implicitly restarted Lanczos !! bidiagonalization methods." !! - R. M. Larsen. "Lanczos bidiagonalization with partial reorthogonalization." !! Technical Report, 1998. !! !! ### Syntax !! !! ```fortran !! call svds(A, U, S, V, residuals, info [, kdim] [,tolerance]) !! ``` !! !! ### Arguments !! !! - `A` : Linear operator derived from `abstract_sym_linop_rsp`, `abstract_sym_linop_rdp`, !! `abstract_hermitian_linop_csp` or `abstract_hermitian_linop_cdp` whose leading !! eigenpairs need to be computed. It is an `intent(inout)` argument. !! !! - `U` : Array of `abstract_vectors` with the same type and kind as `A`. On exit, it !! contains the left singular vectors of `A`. Note that the dimension of `U` fixes !! the number of eigenpairs computed. It is an `intent(out)` argument. !! !! - `S` : Rank-1 array of `real` numbers. On exit, it contains the leading !! singular values of `A`. It is an `intent(out)` argument. !! !! - `V` : Array of `abstract_vectors` with the same type and kind as `A`. On exit, it !! contains the left singular vectors of `A`. Note that the dimension of `U` fixes !! the number of eigenpairs computed. It is an `intent(out)` argument. !! !! - `residuals` : Rank-1 array of `real` numbers. On exit, it contains the residuals !! associated with each singular triplet. It is an `intent(out)` argument. !! !! - `info` : `integer` Information flag. !! !! - `kdim` (*optional*) : `integer`, maximum dimension of the Krylov subspace used to !! approximate the leading singular triplets. It is an optional !! `intent(in)` argument. By default, `kdim = 4*size(X)`. !! !! - `tolerance` (*optional*) : `real` tolerance below which a triplet is considered as being converged. It is an !! optional `intent(in)` argument. By default, tolerance = rtol_sp` or `tolerance = rtol_dp`. !! @note !! This implementation does not currently include an automatic restarting procedure !! such as `krylov_schur` for `eigs`. This is work in progress. !! @endnote module subroutine svds_rsp(A, U, S, V, residuals, info, u0, kdim, tolerance, write_intermediate) implicit none(type, external) class(abstract_linop_rsp), intent(inout) :: A !! Linear operator whose leading singular triplets need to be computed. class(abstract_vector_rsp), intent(out) :: U(:) !! Leading left singular vectors. real(sp), allocatable, intent(out) :: S(:) !! Leading singular values. class(abstract_vector_rsp), intent(out) :: V(:) !! Leading right singular vectors. real(sp), allocatable, intent(out) :: residuals(:) !! Residuals associated to each Ritz eigenpair. integer, intent(out) :: info !! Information flag. class(abstract_vector_rsp), optional, intent(in) :: u0 integer, optional, intent(in) :: kdim !! Desired number of eigenpairs. real(sp), optional, intent(in) :: tolerance !! Tolerance. logical, optional, intent(in) :: write_intermediate !! Write intermediate eigenvalues to file during iteration? end subroutine svds_rsp module subroutine svds_rdp(A, U, S, V, residuals, info, u0, kdim, tolerance, write_intermediate) implicit none(type, external) class(abstract_linop_rdp), intent(inout) :: A !! Linear operator whose leading singular triplets need to be computed. class(abstract_vector_rdp), intent(out) :: U(:) !! Leading left singular vectors. real(dp), allocatable, intent(out) :: S(:) !! Leading singular values. class(abstract_vector_rdp), intent(out) :: V(:) !! Leading right singular vectors. real(dp), allocatable, intent(out) :: residuals(:) !! Residuals associated to each Ritz eigenpair. integer, intent(out) :: info !! Information flag. class(abstract_vector_rdp), optional, intent(in) :: u0 integer, optional, intent(in) :: kdim !! Desired number of eigenpairs. real(dp), optional, intent(in) :: tolerance !! Tolerance. logical, optional, intent(in) :: write_intermediate !! Write intermediate eigenvalues to file during iteration? end subroutine svds_rdp module subroutine svds_csp(A, U, S, V, residuals, info, u0, kdim, tolerance, write_intermediate) implicit none(type, external) class(abstract_linop_csp), intent(inout) :: A !! Linear operator whose leading singular triplets need to be computed. class(abstract_vector_csp), intent(out) :: U(:) !! Leading left singular vectors. real(sp), allocatable, intent(out) :: S(:) !! Leading singular values. class(abstract_vector_csp), intent(out) :: V(:) !! Leading right singular vectors. real(sp), allocatable, intent(out) :: residuals(:) !! Residuals associated to each Ritz eigenpair. integer, intent(out) :: info !! Information flag. class(abstract_vector_csp), optional, intent(in) :: u0 integer, optional, intent(in) :: kdim !! Desired number of eigenpairs. real(sp), optional, intent(in) :: tolerance !! Tolerance. logical, optional, intent(in) :: write_intermediate !! Write intermediate eigenvalues to file during iteration? end subroutine svds_csp module subroutine svds_cdp(A, U, S, V, residuals, info, u0, kdim, tolerance, write_intermediate) implicit none(type, external) class(abstract_linop_cdp), intent(inout) :: A !! Linear operator whose leading singular triplets need to be computed. class(abstract_vector_cdp), intent(out) :: U(:) !! Leading left singular vectors. real(dp), allocatable, intent(out) :: S(:) !! Leading singular values. class(abstract_vector_cdp), intent(out) :: V(:) !! Leading right singular vectors. real(dp), allocatable, intent(out) :: residuals(:) !! Residuals associated to each Ritz eigenpair. integer, intent(out) :: info !! Information flag. class(abstract_vector_cdp), optional, intent(in) :: u0 integer, optional, intent(in) :: kdim !! Desired number of eigenpairs. real(dp), optional, intent(in) :: tolerance !! Tolerance. logical, optional, intent(in) :: write_intermediate !! Write intermediate eigenvalues to file during iteration? end subroutine svds_cdp end interface !------------------------------------------------ !----- ----- !----- HERMITIAN EIGENVALUE SOLVERS ----- !----- ----- !------------------------------------------------ interface eighs !! ### Description !! !! Computes the leading eigenpairs of a symmetric operator \(A\) using the !! Lanczos iterative process. Given a square linear operator \(A\), it finds !! the leading eigvalues and eigvectors such that: !! !! \[ !! Ax = \lambda x !! \] !! !! The subspace \(X\) is constructed via Lanczos factorization, resulting in a symmetric !! tridiagonal matrix \(T\). The eigenvalues of \(A\) are approximated by those of \(T\) !! and the eigenvectors are computed accordingly. !! !! **References** !! !! - Lanczos, C. (1950). "An Iteration Method for the Solution of the Eigenvalue Problem !! of Linear Differential and Integral Operators". United States Governm. Press Office. !! !! ### Syntax !! !! ```fortran !! call eighs(A, X, eigvals, residuals, info [, kdim] [,tolerance]) !! ``` !! !! ### Arguments !! !! - `A` : Linear operator derived from `abstract_sym_linop_rsp`, `abstract_sym_linop_rdp`, !! `abstract_hermitian_linop_csp` or `abstract_hermitian_linop_cdp` whose leading !! eigenpairs need to be computed. It is an `intent(inout)` argument. !! !! - `X` : Array of `abstract_vectors` with the same type and kind as `A`. On exit, it !! contains the leading eigenvectors of `A`. Note that the dimension of `X` fixes !! the number of eigenpairs computed. It is an `intent(out)` argument. !! !! - `eigvals` : Rank-1 array of `real` numbers. On exit, it contains the leading !! eigenvalues of `A`. It is an `intent(out)` argument. !! !! - `residuals` : Rank-1 array of `real` numbers. On exit, it contains the residuals !! associated with each eigenpairs. It is an `intent(out)` argument. !! !! - `info` : `integer` Information flag. !! !! - `kdim` (*optional*) : `integer`, maximum dimension of the Krylov subspace used to !! approximate the leading eigenpairs. It is an optional `intent(in)` !! argument. By default, `kdim = 4*size(X)`. !! !! - `tolerance` (*optional*) : `real` tolerance below which an eigenpair is considered as !! being converged. It is an optional `intent(in)` argument. By default, !! `tolerance = rtol_sp` or `tolerance = rtol_dp`. !! @note !! This implementation does not currently include an automatic restarting procedure !! such as `krylov_schur` for `eigs`. This is work in progress. !! @endnote module subroutine eighs_rsp(A, X, eigvals, residuals, info, x0, kdim, tolerance, write_intermediate) implicit none(type, external) class(abstract_sym_linop_rsp), intent(inout) :: A !! Linear operator whose leading eigenpairs need to be computed. class(abstract_vector_rsp), intent(out) :: X(:) !! Leading eigevectors of \( \mathbf{A} \). real(sp), allocatable, intent(out) :: eigvals(:) !! Leading eigenvalues of \( \mathbf{A} \). real(sp), allocatable, intent(out) :: residuals(:) !! Residuals associated to each Ritz eigenpairs. integer, intent(out) :: info !! Information flag. class(abstract_vector_rsp), optional, intent(in) :: x0 !! Optional starting vector to generate the Krylov subspace. integer, optional, intent(in) :: kdim !! Desired number of eigenpairs. real(sp), optional, intent(in) :: tolerance !! Tolerance logical, optional, intent(in) :: write_intermediate !! Write intermediate eigenvalues to file during iteration? end subroutine eighs_rsp module subroutine eighs_rdp(A, X, eigvals, residuals, info, x0, kdim, tolerance, write_intermediate) implicit none(type, external) class(abstract_sym_linop_rdp), intent(inout) :: A !! Linear operator whose leading eigenpairs need to be computed. class(abstract_vector_rdp), intent(out) :: X(:) !! Leading eigevectors of \( \mathbf{A} \). real(dp), allocatable, intent(out) :: eigvals(:) !! Leading eigenvalues of \( \mathbf{A} \). real(dp), allocatable, intent(out) :: residuals(:) !! Residuals associated to each Ritz eigenpairs. integer, intent(out) :: info !! Information flag. class(abstract_vector_rdp), optional, intent(in) :: x0 !! Optional starting vector to generate the Krylov subspace. integer, optional, intent(in) :: kdim !! Desired number of eigenpairs. real(dp), optional, intent(in) :: tolerance !! Tolerance logical, optional, intent(in) :: write_intermediate !! Write intermediate eigenvalues to file during iteration? end subroutine eighs_rdp module subroutine eighs_csp(A, X, eigvals, residuals, info, x0, kdim, tolerance, write_intermediate) implicit none(type, external) class(abstract_hermitian_linop_csp), intent(inout) :: A !! Linear operator whose leading eigenpairs need to be computed. class(abstract_vector_csp), intent(out) :: X(:) !! Leading eigevectors of \( \mathbf{A} \). real(sp), allocatable, intent(out) :: eigvals(:) !! Leading eigenvalues of \( \mathbf{A} \). real(sp), allocatable, intent(out) :: residuals(:) !! Residuals associated to each Ritz eigenpairs. integer, intent(out) :: info !! Information flag. class(abstract_vector_csp), optional, intent(in) :: x0 !! Optional starting vector to generate the Krylov subspace. integer, optional, intent(in) :: kdim !! Desired number of eigenpairs. real(sp), optional, intent(in) :: tolerance !! Tolerance logical, optional, intent(in) :: write_intermediate !! Write intermediate eigenvalues to file during iteration? end subroutine eighs_csp module subroutine eighs_cdp(A, X, eigvals, residuals, info, x0, kdim, tolerance, write_intermediate) implicit none(type, external) class(abstract_hermitian_linop_cdp), intent(inout) :: A !! Linear operator whose leading eigenpairs need to be computed. class(abstract_vector_cdp), intent(out) :: X(:) !! Leading eigevectors of \( \mathbf{A} \). real(dp), allocatable, intent(out) :: eigvals(:) !! Leading eigenvalues of \( \mathbf{A} \). real(dp), allocatable, intent(out) :: residuals(:) !! Residuals associated to each Ritz eigenpairs. integer, intent(out) :: info !! Information flag. class(abstract_vector_cdp), optional, intent(in) :: x0 !! Optional starting vector to generate the Krylov subspace. integer, optional, intent(in) :: kdim !! Desired number of eigenpairs. real(dp), optional, intent(in) :: tolerance !! Tolerance logical, optional, intent(in) :: write_intermediate !! Write intermediate eigenvalues to file during iteration? end subroutine eighs_cdp end interface !------------------------------------- !----- Utility functions ----- !------------------------------------- interface save_eigenspectrum !! ### Description !! !! Utility function to save the eigenspectrum computed from the Arnoldi factorization. !! It outpost a .npy file. !! !! ### Syntax !! !! ```fortran !! call save_eigenspectrum(eigvals, residuals, fname) !! ``` !! !! ### Arguments !! !! - `eigvals` : `complex` rank-1 array containing the eigenvalues. It is an !! `intent(in)` argument. !! !! - `residuals` : `real` rank-1 array containing the residuals associated to each !! eigenvalues. It is an `intent(in)` argument. !! !! `fname` : Name of the file to save the eigenspectrum. It is an `intent(in)` argument. module subroutine save_eigenspectrum_rsp(lambda, residuals, fname) !! Saves the eigenspectrum and corresponding residuals to disk use the `npy` binary format. implicit none(type, external) real(sp), intent(in) :: lambda(:) !! Eigenalues. real(sp), intent(in) :: residuals(:) !! Residual of the corresponding Ritz eigenpairs. character(len=*), intent(in) :: fname !! Name of the output file. end subroutine save_eigenspectrum_rsp module subroutine save_eigenspectrum_rdp(lambda, residuals, fname) !! Saves the eigenspectrum and corresponding residuals to disk use the `npy` binary format. implicit none(type, external) real(dp), intent(in) :: lambda(:) !! Eigenalues. real(dp), intent(in) :: residuals(:) !! Residual of the corresponding Ritz eigenpairs. character(len=*), intent(in) :: fname !! Name of the output file. end subroutine save_eigenspectrum_rdp module subroutine save_eigenspectrum_csp(lambda, residuals, fname) !! Saves the eigenspectrum and corresponding residuals to disk use the `npy` binary format. implicit none(type, external) complex(sp), intent(in) :: lambda(:) !! Eigenalues. real(sp), intent(in) :: residuals(:) !! Residual of the corresponding Ritz eigenpairs. character(len=*), intent(in) :: fname !! Name of the output file. end subroutine save_eigenspectrum_csp module subroutine save_eigenspectrum_cdp(lambda, residuals, fname) !! Saves the eigenspectrum and corresponding residuals to disk use the `npy` binary format. implicit none(type, external) complex(dp), intent(in) :: lambda(:) !! Eigenalues. real(dp), intent(in) :: residuals(:) !! Residual of the corresponding Ritz eigenpairs. character(len=*), intent(in) :: fname !! Name of the output file. end subroutine save_eigenspectrum_cdp end interface interface eigs !! ### Description !! !! Computes the leading eigenpairs of a square linear operator \(A\) !! using the Arnoldi iterative process. Given a square linear operator \(A\), it finds !! the leading eigenvalues and eigenvectors such that: !! !! \[ !! Ax = \lambda x !! \] !! !! or !! !! \[ !! A^H x = \lambda x. !! \] !! !! The subspace \(X\) is constructed via Arnoldi factorization, resulting in an upper !! Hessenberg matrix \(H\). The eigenvalues of \(A\) are approximated by those of \(H\) !! and the eigenvectors are computed accordingly. !! !! **References** !! !! - Arnoldi, W. E. (1951). "The Principle of Minimized Iterations in the Solution of !! the Matrix Eigenvalue Problem." Quarterly of Applied Mathematics, 9(1), 17–29. !! !! ### Syntax !! !! ```fortran !! call eigs(A, X, eigvals, residuals, info [, kdim] [, select] [,tolerance] [, transpose]) !! ``` !! !! ### Arguments !! !! - `A` : Linear operator derived from `abstract_sym_linop_rsp`, `abstract_sym_linop_rdp`, !! `abstract_hermitian_linop_csp` or `abstract_hermitian_linop_cdp` whose leading !! eigenpairs need to be computed. It is an `intent(inout)` argument. !! !! - `X` : Array of `abstract_vectors` with the same type and kind as `A`. On exit, it !! contains the leading eigenvectors of `A`. Note that the dimension of `X` fixes !! the number of eigenpairs computed. It is an `intent(out)` argument. !! !! - `eigvals` : Rank-1 array of `real` numbers. On exit, it contains the leading !! eigenvalues of `A`. It is an `intent(out)` argument. !! !! - `residuals` : Rank-1 array of `real` numbers. On exit, it contains the residuals !! associated with each eigenpairs. It is an `intent(out)` argument. !! !! - `info` : `integer` Information flag. !! !! - `kdim` (*optional*) : `integer`, maximum dimension of the Krylov subspace used to !! approximate the leading eigenpairs. It is an optional `intent(in)` !! argument. By default, `kdim = 4*size(X)`. !! !! - `select` (*optional*) : Function to select which eigenvalues to compute. !! !! - `tolerance` (*optional*) : `real` tolerance below which an eigenpair is considered as !! being converged. It is an optional `intent(in)` argument. !! By default, `tolerance = rtol_sp` or `tolerance = rtol_dp`. !! !! - `transpose` (*optional*) : `logical` flag determining whether the eigenvalues of \(A\) !! or \(A^H\) need to be computed. It is an optional `intent(in)` !! argument. module procedure eigs_rsp module procedure eigs_rdp module procedure eigs_csp module procedure eigs_cdp end interface contains !------------------------------------- !----- UTILITY FUNCTIONS ----- !------------------------------------- subroutine write_results_rsp(filename, vals, res, tol) !! Prints the intermediate results of iterative eigenvalue/singular value decompositions implicit none(type, external) character(len=*), intent(in) :: filename !! Output filename. This file will be overwritten real(sp), intent(in) :: vals(:) !! Intermediate values real(sp), intent(inout) :: res(:) !! Residuals real(sp), intent(in) :: tol !! Convergence tolerance ! internals integer :: i, k, idx integer, allocatable :: indices(:) character(len=*), parameter :: fmt = '(I6,2(2X,E16.9),2X,L4)' k = size(vals) if (io_rank()) then ! only IO rank writes allocate(indices(k)) call sort_index(res, indices) ! res is returned in sorted order open (1234, file=filename, status='replace', action='write') write (1234, '(A6,2(A18),A6)') 'Iter', 'value', 'residual', 'conv' do i = 1, k idx = indices(i) write (1234, fmt) k, vals(idx), res(i), res(i) < tol end do close (1234) end if end subroutine write_results_rsp subroutine write_results_rdp(filename, vals, res, tol) !! Prints the intermediate results of iterative eigenvalue/singular value decompositions implicit none(type, external) character(len=*), intent(in) :: filename !! Output filename. This file will be overwritten real(dp), intent(in) :: vals(:) !! Intermediate values real(dp), intent(inout) :: res(:) !! Residuals real(dp), intent(in) :: tol !! Convergence tolerance ! internals integer :: i, k, idx integer, allocatable :: indices(:) character(len=*), parameter :: fmt = '(I6,2(2X,E16.9),2X,L4)' k = size(vals) if (io_rank()) then ! only IO rank writes allocate(indices(k)) call sort_index(res, indices) ! res is returned in sorted order open (1234, file=filename, status='replace', action='write') write (1234, '(A6,2(A18),A6)') 'Iter', 'value', 'residual', 'conv' do i = 1, k idx = indices(i) write (1234, fmt) k, vals(idx), res(i), res(i) < tol end do close (1234) end if end subroutine write_results_rdp subroutine write_results_csp(filename, vals, res, tol) !! Prints the intermediate results of iterative eigenvalue/singular value decompositions implicit none(type, external) character(len=*), intent(in) :: filename !! Output filename. This file will be overwritten complex(sp), intent(in) :: vals(:) !! Intermediate values real(sp), intent(inout) :: res(:) !! Residuals real(sp), intent(in) :: tol !! Convergence tolerance ! internals integer :: i, k, idx integer, allocatable :: indices(:) real(sp) :: modulus character(len=*), parameter :: fmt = '(I6,4(2X,E16.9),2X,L4)' k = size(vals) if (io_rank()) then ! only IO rank writes allocate(indices(k)) call sort_index(res, indices) ! res is returned in sorted order open (1234, file=filename, status='replace', action='write') write (1234, '(A6,4(A18),A6)') 'Iter', 'Re', 'Im', 'modulus', 'residual', 'conv' do i = 1, k idx = indices(i) modulus = sqrt(vals(idx)%re**2 + vals(idx)%im**2) write (1234, fmt) k, vals(idx)%re, vals(idx)%im, modulus, res(i), res(i) < tol end do close (1234) end if end subroutine write_results_csp subroutine write_results_cdp(filename, vals, res, tol) !! Prints the intermediate results of iterative eigenvalue/singular value decompositions implicit none(type, external) character(len=*), intent(in) :: filename !! Output filename. This file will be overwritten complex(dp), intent(in) :: vals(:) !! Intermediate values real(dp), intent(inout) :: res(:) !! Residuals real(dp), intent(in) :: tol !! Convergence tolerance ! internals integer :: i, k, idx integer, allocatable :: indices(:) real(dp) :: modulus character(len=*), parameter :: fmt = '(I6,4(2X,E16.9),2X,L4)' k = size(vals) if (io_rank()) then ! only IO rank writes allocate(indices(k)) call sort_index(res, indices) ! res is returned in sorted order open (1234, file=filename, status='replace', action='write') write (1234, '(A6,4(A18),A6)') 'Iter', 'Re', 'Im', 'modulus', 'residual', 'conv' do i = 1, k idx = indices(i) modulus = sqrt(vals(idx)%re**2 + vals(idx)%im**2) write (1234, fmt) k, vals(idx)%re, vals(idx)%im, modulus, res(i), res(i) < tol end do close (1234) end if end subroutine write_results_cdp elemental pure function compute_residual_rsp(beta, x) result(residual) !! Computes the residual associated with a Ritz eigenpair. real(sp), intent(in) :: beta !! Norm of the residual Krylov vector. real(sp), intent(in) :: x !! Last entry of the low-dimensional Ritz eigenvector. real(sp) :: residual !! Residual associated to the corresponding Ritz eigenpair. residual = abs(beta*x) end function compute_residual_rsp elemental pure function compute_residual_rdp(beta, x) result(residual) !! Computes the residual associated with a Ritz eigenpair. real(dp), intent(in) :: beta !! Norm of the residual Krylov vector. real(dp), intent(in) :: x !! Last entry of the low-dimensional Ritz eigenvector. real(dp) :: residual !! Residual associated to the corresponding Ritz eigenpair. residual = abs(beta*x) end function compute_residual_rdp elemental pure function compute_residual_csp(beta, x) result(residual) !! Computes the residual associated with a Ritz eigenpair. complex(sp), intent(in) :: beta !! Norm of the residual Krylov vector. complex(sp), intent(in) :: x !! Last entry of the low-dimensional Ritz eigenvector. real(sp) :: residual !! Residual associated to the corresponding Ritz eigenpair. residual = abs(beta*x) end function compute_residual_csp elemental pure function compute_residual_cdp(beta, x) result(residual) !! Computes the residual associated with a Ritz eigenpair. complex(dp), intent(in) :: beta !! Norm of the residual Krylov vector. complex(dp), intent(in) :: x !! Last entry of the low-dimensional Ritz eigenvector. real(dp) :: residual !! Residual associated to the corresponding Ritz eigenpair. residual = abs(beta*x) end function compute_residual_cdp module procedure save_eigenspectrum_rsp ! Internal variables. real(sp) :: array(size(lambda), 2) array(:, 1) = lambda ; array(:, 2) = residuals ! Save the eigenspectrum to disk. call save_npy(fname, array) end procedure save_eigenspectrum_rsp module procedure save_eigenspectrum_rdp ! Internal variables. real(dp) :: array(size(lambda), 2) array(:, 1) = lambda ; array(:, 2) = residuals ! Save the eigenspectrum to disk. call save_npy(fname, array) end procedure save_eigenspectrum_rdp module procedure save_eigenspectrum_csp ! Internal variables. real(sp) :: array(size(lambda), 3) array(:, 1) = lambda%re ; array(:, 2) = lambda%im ; array(:, 3) = residuals ! Save the eigenspectrum to disk. call save_npy(fname, array) end procedure save_eigenspectrum_csp module procedure save_eigenspectrum_cdp ! Internal variables. real(dp) :: array(size(lambda), 3) array(:, 1) = lambda%re ; array(:, 2) = lambda%im ; array(:, 3) = residuals ! Save the eigenspectrum to disk. call save_npy(fname, array) end procedure save_eigenspectrum_cdp !--------------------------------------------------- !----- GENERAL EIGENVALUE COMPUTATIONS ----- !--------------------------------------------------- subroutine eigs_rsp(A, X, eigvals, residuals, info, x0, kdim, tolerance, transpose, write_intermediate) implicit none(type, external) class(abstract_linop_rsp), intent(inout) :: A !! Linear operator whose leading eigenpairs need to be computed. class(abstract_vector_rsp), intent(out) :: X(:) !! Leading eigenvectors of \(\mathbf{A}\). complex(sp), allocatable, intent(out) :: eigvals(:) !! Leading eigenvalues of \(\mathbf{A}\). real(sp), allocatable, intent(out) :: residuals(:) !! Residuals associated to each Ritz eigenpair. integer, intent(out) :: info !! Information flag. class(abstract_vector_rsp), optional, intent(in) :: x0 !! Optional starting vector for generating the Krylov subspace. integer, optional, intent(in) :: kdim !! Maximum dimension of the Krylov subspace (optional). real(sp), optional, intent(in) :: tolerance !! Tolerance. logical, optional, intent(in) :: transpose !! Determine whether \(\mathbf{A}\) or \(\mathbf{A}^H\) is being used. logical, optional, intent(in) :: write_intermediate !! Write intermediate eigenvalues to file during iteration? !-------------------------------------- !----- Internal variables ----- !-------------------------------------- ! Krylov subspace and Krylov subspace dimension. class(abstract_vector_rsp), allocatable :: Xwrk(:) integer :: kdim_, kstart ! Hessenberg matrix. real(sp), allocatable :: H(:, :) ! Working arrays for the eigenvectors and eigenvalues. real(sp), allocatable :: eigvecs_wrk(:, :) complex(sp), allocatable :: eigvals_wrk(:) real(sp), allocatable :: residuals_wrk(:) ! Miscellaneous. character(len=*), parameter :: this_procedure = 'eigs_rsp' integer :: nev, conv integer :: i, j, k, niter, krst real(sp) :: tol, x0_norm real(sp) :: beta real(sp) :: alpha logical :: outpost character(len=256) :: msg if (time_lightkrylov()) call timer%start(this_procedure) ! Deals with optional parameters. nev = size(X) kdim_ = optval(kdim, 4*nev) tol = optval(tolerance, rtol_sp) outpost = optval(write_intermediate, .true.) ! Allocate eigenvalues. allocate(eigvals(nev)) ; eigvals = 0.0_sp ! Allocate working variables. allocate(Xwrk(kdim_+1), source=X(1)) ; call zero_basis(Xwrk) if (present(x0)) then call copy(Xwrk(1), x0) x0_norm = x0%norm(); call Xwrk(1)%scal(one_rsp/x0_norm) else call Xwrk(1)%rand(.true.) endif allocate(H(kdim_+1, kdim_)) ; H = 0.0_sp allocate(eigvecs_wrk(kdim_, kdim_)) ; eigvecs_wrk = 0.0_sp allocate(eigvals_wrk(kdim_)) ; eigvals_wrk = 0.0_sp allocate(residuals_wrk(kdim_)) ; residuals_wrk = 0.0_sp ! Ritz eigenpairs computation. H = 0.0_sp kstart = 1 ; conv = 0 ; niter = 0 ; krst = 1 krylovschur: do while (conv < nev) arnoldi_factorization: do k = kstart, kdim_ ! Arnoldi step. call arnoldi(A, Xwrk, H, info, kstart=k, kend=k, transpose=transpose) call check_info(info, 'arnoldi', this_module, this_procedure) ! Spectral decomposition of the k x k Hessenberg matrix. eigvals_wrk = 0.0_sp ; eigvecs_wrk = 0.0_sp if (time_lightkrylov()) call timer%start('eig') call eig(H(:k, :k), eigvecs_wrk(:k, :k), eigvals_wrk(:k)) if (time_lightkrylov()) call timer%stop('eig') ! Compute residuals. beta = H(k+1, k) do i = 1, k if (eigvals_wrk(i)%im > 0) then alpha = abs(cmplx(eigvecs_wrk(k, i), eigvecs_wrk(k, i+1), kind=sp)) else if (eigvals_wrk(i)%im < 0) then alpha = abs(cmplx(eigvecs_wrk(k, i-1), eigvecs_wrk(k, i), kind=sp)) else alpha = abs(eigvecs_wrk(k, i)) endif residuals_wrk(i) = compute_residual_rsp(beta, alpha) enddo ! Check convergence. niter = niter + 1 conv = count(residuals_wrk(:k) < tol) write(msg,'(I0,A,I0,A,I0,A)') conv, '/', nev, ' eigenvalues converged after ', niter, & & ' steps of the Arnoldi process.' call log_information(msg, this_module, this_procedure) if (outpost) call write_results_csp(eigs_output, eigvals_wrk(:k), residuals_wrk(:k), tol) if (conv >= nev) exit arnoldi_factorization enddo arnoldi_factorization write(msg,'(I0,A,I0,A,I0,A)') conv, '/', nev, ' eigenvalues converged after ', krst, & & ' Krylov-Schur restarts of the Arnoldi process.' call log_information(msg, this_module, this_procedure) ! Krylov-Schur restarting procedure. krst = krst + 1 call krylov_schur(kstart, Xwrk, H, median_selector) ; kstart = kstart + 1 end do krylovschur !-------------------------------- !----- POST-PROCESS ----- !-------------------------------- block integer :: indices(kdim_) real(sp) :: abs_eigvals(kdim_) ! Re-compute eigenvalues and eigenvectors. k = min(k, kdim_); if (time_lightkrylov()) call timer%start('eig') call eig(H(:k, :k), eigvecs_wrk(:k, :k), eigvals_wrk(:k)) if (time_lightkrylov()) call timer%stop('eig') ! Sort eigenvalues. abs_eigvals = abs(eigvals_wrk) ; call sort_index(abs_eigvals, indices, reverse=.true.) eigvals_wrk = eigvals_wrk(indices) ; eigvecs_wrk = eigvecs_wrk(:, indices) residuals_wrk = residuals_wrk(indices) ! Store converged eigenvalues. eigvals = eigvals_wrk(:nev) ; residuals = residuals_wrk(:nev) end block ! Construct eigenvectors. do i = 1, nev call X(i)%zero() do j = 1, k call X(i)%axpby(eigvecs_wrk(j, i), Xwrk(j), one_rsp) enddo enddo info = niter if (time_lightkrylov()) call timer%stop(this_procedure) contains function median_selector(lambda) result(selected) implicit none(type, external) complex(sp), intent(in) :: lambda(:) logical, allocatable :: selected(:) selected = abs(lambda) > median(abs(lambda)) end function median_selector end subroutine eigs_rsp subroutine eigs_rdp(A, X, eigvals, residuals, info, x0, kdim, tolerance, transpose, write_intermediate) implicit none(type, external) class(abstract_linop_rdp), intent(inout) :: A !! Linear operator whose leading eigenpairs need to be computed. class(abstract_vector_rdp), intent(out) :: X(:) !! Leading eigenvectors of \(\mathbf{A}\). complex(dp), allocatable, intent(out) :: eigvals(:) !! Leading eigenvalues of \(\mathbf{A}\). real(dp), allocatable, intent(out) :: residuals(:) !! Residuals associated to each Ritz eigenpair. integer, intent(out) :: info !! Information flag. class(abstract_vector_rdp), optional, intent(in) :: x0 !! Optional starting vector for generating the Krylov subspace. integer, optional, intent(in) :: kdim !! Maximum dimension of the Krylov subspace (optional). real(dp), optional, intent(in) :: tolerance !! Tolerance. logical, optional, intent(in) :: transpose !! Determine whether \(\mathbf{A}\) or \(\mathbf{A}^H\) is being used. logical, optional, intent(in) :: write_intermediate !! Write intermediate eigenvalues to file during iteration? !-------------------------------------- !----- Internal variables ----- !-------------------------------------- ! Krylov subspace and Krylov subspace dimension. class(abstract_vector_rdp), allocatable :: Xwrk(:) integer :: kdim_, kstart ! Hessenberg matrix. real(dp), allocatable :: H(:, :) ! Working arrays for the eigenvectors and eigenvalues. real(dp), allocatable :: eigvecs_wrk(:, :) complex(dp), allocatable :: eigvals_wrk(:) real(dp), allocatable :: residuals_wrk(:) ! Miscellaneous. character(len=*), parameter :: this_procedure = 'eigs_rdp' integer :: nev, conv integer :: i, j, k, niter, krst real(dp) :: tol, x0_norm real(dp) :: beta real(dp) :: alpha logical :: outpost character(len=256) :: msg if (time_lightkrylov()) call timer%start(this_procedure) ! Deals with optional parameters. nev = size(X) kdim_ = optval(kdim, 4*nev) tol = optval(tolerance, rtol_dp) outpost = optval(write_intermediate, .true.) ! Allocate eigenvalues. allocate(eigvals(nev)) ; eigvals = 0.0_dp ! Allocate working variables. allocate(Xwrk(kdim_+1), source=X(1)) ; call zero_basis(Xwrk) if (present(x0)) then call copy(Xwrk(1), x0) x0_norm = x0%norm(); call Xwrk(1)%scal(one_rdp/x0_norm) else call Xwrk(1)%rand(.true.) endif allocate(H(kdim_+1, kdim_)) ; H = 0.0_dp allocate(eigvecs_wrk(kdim_, kdim_)) ; eigvecs_wrk = 0.0_dp allocate(eigvals_wrk(kdim_)) ; eigvals_wrk = 0.0_dp allocate(residuals_wrk(kdim_)) ; residuals_wrk = 0.0_dp ! Ritz eigenpairs computation. H = 0.0_dp kstart = 1 ; conv = 0 ; niter = 0 ; krst = 1 krylovschur: do while (conv < nev) arnoldi_factorization: do k = kstart, kdim_ ! Arnoldi step. call arnoldi(A, Xwrk, H, info, kstart=k, kend=k, transpose=transpose) call check_info(info, 'arnoldi', this_module, this_procedure) ! Spectral decomposition of the k x k Hessenberg matrix. eigvals_wrk = 0.0_dp ; eigvecs_wrk = 0.0_dp if (time_lightkrylov()) call timer%start('eig') call eig(H(:k, :k), eigvecs_wrk(:k, :k), eigvals_wrk(:k)) if (time_lightkrylov()) call timer%stop('eig') ! Compute residuals. beta = H(k+1, k) do i = 1, k if (eigvals_wrk(i)%im > 0) then alpha = abs(cmplx(eigvecs_wrk(k, i), eigvecs_wrk(k, i+1), kind=dp)) else if (eigvals_wrk(i)%im < 0) then alpha = abs(cmplx(eigvecs_wrk(k, i-1), eigvecs_wrk(k, i), kind=dp)) else alpha = abs(eigvecs_wrk(k, i)) endif residuals_wrk(i) = compute_residual_rdp(beta, alpha) enddo ! Check convergence. niter = niter + 1 conv = count(residuals_wrk(:k) < tol) write(msg,'(I0,A,I0,A,I0,A)') conv, '/', nev, ' eigenvalues converged after ', niter, & & ' steps of the Arnoldi process.' call log_information(msg, this_module, this_procedure) if (outpost) call write_results_cdp(eigs_output, eigvals_wrk(:k), residuals_wrk(:k), tol) if (conv >= nev) exit arnoldi_factorization enddo arnoldi_factorization write(msg,'(I0,A,I0,A,I0,A)') conv, '/', nev, ' eigenvalues converged after ', krst, & & ' Krylov-Schur restarts of the Arnoldi process.' call log_information(msg, this_module, this_procedure) ! Krylov-Schur restarting procedure. krst = krst + 1 call krylov_schur(kstart, Xwrk, H, median_selector) ; kstart = kstart + 1 end do krylovschur !-------------------------------- !----- POST-PROCESS ----- !-------------------------------- block integer :: indices(kdim_) real(dp) :: abs_eigvals(kdim_) ! Re-compute eigenvalues and eigenvectors. k = min(k, kdim_); if (time_lightkrylov()) call timer%start('eig') call eig(H(:k, :k), eigvecs_wrk(:k, :k), eigvals_wrk(:k)) if (time_lightkrylov()) call timer%stop('eig') ! Sort eigenvalues. abs_eigvals = abs(eigvals_wrk) ; call sort_index(abs_eigvals, indices, reverse=.true.) eigvals_wrk = eigvals_wrk(indices) ; eigvecs_wrk = eigvecs_wrk(:, indices) residuals_wrk = residuals_wrk(indices) ! Store converged eigenvalues. eigvals = eigvals_wrk(:nev) ; residuals = residuals_wrk(:nev) end block ! Construct eigenvectors. do i = 1, nev call X(i)%zero() do j = 1, k call X(i)%axpby(eigvecs_wrk(j, i), Xwrk(j), one_rdp) enddo enddo info = niter if (time_lightkrylov()) call timer%stop(this_procedure) contains function median_selector(lambda) result(selected) implicit none(type, external) complex(dp), intent(in) :: lambda(:) logical, allocatable :: selected(:) selected = abs(lambda) > median(abs(lambda)) end function median_selector end subroutine eigs_rdp subroutine eigs_csp(A, X, eigvals, residuals, info, x0, kdim, tolerance, transpose, write_intermediate) implicit none(type, external) class(abstract_linop_csp), intent(inout) :: A !! Linear operator whose leading eigenpairs need to be computed. class(abstract_vector_csp), intent(out) :: X(:) !! Leading eigenvectors of \(\mathbf{A}\). complex(sp), allocatable, intent(out) :: eigvals(:) !! Leading eigenvalues of \(\mathbf{A}\). real(sp), allocatable, intent(out) :: residuals(:) !! Residuals associated to each Ritz eigenpair. integer, intent(out) :: info !! Information flag. class(abstract_vector_csp), optional, intent(in) :: x0 !! Optional starting vector for generating the Krylov subspace. integer, optional, intent(in) :: kdim !! Maximum dimension of the Krylov subspace (optional). real(sp), optional, intent(in) :: tolerance !! Tolerance. logical, optional, intent(in) :: transpose !! Determine whether \(\mathbf{A}\) or \(\mathbf{A}^H\) is being used. logical, optional, intent(in) :: write_intermediate !! Write intermediate eigenvalues to file during iteration? !-------------------------------------- !----- Internal variables ----- !-------------------------------------- ! Krylov subspace and Krylov subspace dimension. class(abstract_vector_csp), allocatable :: Xwrk(:) integer :: kdim_, kstart ! Hessenberg matrix. complex(sp), allocatable :: H(:, :) ! Working arrays for the eigenvectors and eigenvalues. complex(sp), allocatable :: eigvecs_wrk(:, :) complex(sp), allocatable :: eigvals_wrk(:) real(sp), allocatable :: residuals_wrk(:) ! Miscellaneous. character(len=*), parameter :: this_procedure = 'eigs_csp' integer :: nev, conv integer :: i, j, k, niter, krst real(sp) :: tol, x0_norm complex(sp) :: beta logical :: outpost character(len=256) :: msg if (time_lightkrylov()) call timer%start(this_procedure) ! Deals with optional parameters. nev = size(X) kdim_ = optval(kdim, 4*nev) tol = optval(tolerance, rtol_sp) outpost = optval(write_intermediate, .true.) ! Allocate eigenvalues. allocate(eigvals(nev)) ; eigvals = 0.0_sp ! Allocate working variables. allocate(Xwrk(kdim_+1), source=X(1)) ; call zero_basis(Xwrk) if (present(x0)) then call copy(Xwrk(1), x0) x0_norm = x0%norm(); call Xwrk(1)%scal(one_csp/x0_norm) else call Xwrk(1)%rand(.true.) endif allocate(H(kdim_+1, kdim_)) ; H = 0.0_sp allocate(eigvecs_wrk(kdim_, kdim_)) ; eigvecs_wrk = 0.0_sp allocate(eigvals_wrk(kdim_)) ; eigvals_wrk = 0.0_sp allocate(residuals_wrk(kdim_)) ; residuals_wrk = 0.0_sp ! Ritz eigenpairs computation. H = 0.0_sp kstart = 1 ; conv = 0 ; niter = 0 ; krst = 1 krylovschur: do while (conv < nev) arnoldi_factorization: do k = kstart, kdim_ ! Arnoldi step. call arnoldi(A, Xwrk, H, info, kstart=k, kend=k, transpose=transpose) call check_info(info, 'arnoldi', this_module, this_procedure) ! Spectral decomposition of the k x k Hessenberg matrix. eigvals_wrk = 0.0_sp ; eigvecs_wrk = 0.0_sp if (time_lightkrylov()) call timer%start('eig') call eig(H(:k, :k), eigvecs_wrk(:k, :k), eigvals_wrk(:k)) if (time_lightkrylov()) call timer%stop('eig') ! Compute residuals. beta = H(k+1, k) residuals_wrk(:k) = compute_residual_csp(beta, eigvecs_wrk(k,:k)) ! Check convergence. niter = niter + 1 conv = count(residuals_wrk(:k) < tol) write(msg,'(I0,A,I0,A,I0,A)') conv, '/', nev, ' eigenvalues converged after ', niter, & & ' steps of the Arnoldi process.' call log_information(msg, this_module, this_procedure) if (outpost) call write_results_csp(eigs_output, eigvals_wrk(:k), residuals_wrk(:k), tol) if (conv >= nev) exit arnoldi_factorization enddo arnoldi_factorization write(msg,'(I0,A,I0,A,I0,A)') conv, '/', nev, ' eigenvalues converged after ', krst, & & ' Krylov-Schur restarts of the Arnoldi process.' call log_information(msg, this_module, this_procedure) ! Krylov-Schur restarting procedure. krst = krst + 1 call krylov_schur(kstart, Xwrk, H, median_selector) ; kstart = kstart + 1 end do krylovschur !-------------------------------- !----- POST-PROCESS ----- !-------------------------------- block integer :: indices(kdim_) real(sp) :: abs_eigvals(kdim_) ! Re-compute eigenvalues and eigenvectors. k = min(k, kdim_); if (time_lightkrylov()) call timer%start('eig') call eig(H(:k, :k), eigvecs_wrk(:k, :k), eigvals_wrk(:k)) if (time_lightkrylov()) call timer%stop('eig') ! Sort eigenvalues. abs_eigvals = abs(eigvals_wrk) ; call sort_index(abs_eigvals, indices, reverse=.true.) eigvals_wrk = eigvals_wrk(indices) ; eigvecs_wrk = eigvecs_wrk(:, indices) residuals_wrk = residuals_wrk(indices) ! Store converged eigenvalues. eigvals = eigvals_wrk(:nev) ; residuals = residuals_wrk(:nev) end block ! Construct eigenvectors. do i = 1, nev call X(i)%zero() do j = 1, k call X(i)%axpby(eigvecs_wrk(j, i), Xwrk(j), one_csp) enddo enddo info = niter if (time_lightkrylov()) call timer%stop(this_procedure) contains function median_selector(lambda) result(selected) implicit none(type, external) complex(sp), intent(in) :: lambda(:) logical, allocatable :: selected(:) selected = abs(lambda) > median(abs(lambda)) end function median_selector end subroutine eigs_csp subroutine eigs_cdp(A, X, eigvals, residuals, info, x0, kdim, tolerance, transpose, write_intermediate) implicit none(type, external) class(abstract_linop_cdp), intent(inout) :: A !! Linear operator whose leading eigenpairs need to be computed. class(abstract_vector_cdp), intent(out) :: X(:) !! Leading eigenvectors of \(\mathbf{A}\). complex(dp), allocatable, intent(out) :: eigvals(:) !! Leading eigenvalues of \(\mathbf{A}\). real(dp), allocatable, intent(out) :: residuals(:) !! Residuals associated to each Ritz eigenpair. integer, intent(out) :: info !! Information flag. class(abstract_vector_cdp), optional, intent(in) :: x0 !! Optional starting vector for generating the Krylov subspace. integer, optional, intent(in) :: kdim !! Maximum dimension of the Krylov subspace (optional). real(dp), optional, intent(in) :: tolerance !! Tolerance. logical, optional, intent(in) :: transpose !! Determine whether \(\mathbf{A}\) or \(\mathbf{A}^H\) is being used. logical, optional, intent(in) :: write_intermediate !! Write intermediate eigenvalues to file during iteration? !-------------------------------------- !----- Internal variables ----- !-------------------------------------- ! Krylov subspace and Krylov subspace dimension. class(abstract_vector_cdp), allocatable :: Xwrk(:) integer :: kdim_, kstart ! Hessenberg matrix. complex(dp), allocatable :: H(:, :) ! Working arrays for the eigenvectors and eigenvalues. complex(dp), allocatable :: eigvecs_wrk(:, :) complex(dp), allocatable :: eigvals_wrk(:) real(dp), allocatable :: residuals_wrk(:) ! Miscellaneous. character(len=*), parameter :: this_procedure = 'eigs_cdp' integer :: nev, conv integer :: i, j, k, niter, krst real(dp) :: tol, x0_norm complex(dp) :: beta logical :: outpost character(len=256) :: msg if (time_lightkrylov()) call timer%start(this_procedure) ! Deals with optional parameters. nev = size(X) kdim_ = optval(kdim, 4*nev) tol = optval(tolerance, rtol_dp) outpost = optval(write_intermediate, .true.) ! Allocate eigenvalues. allocate(eigvals(nev)) ; eigvals = 0.0_dp ! Allocate working variables. allocate(Xwrk(kdim_+1), source=X(1)) ; call zero_basis(Xwrk) if (present(x0)) then call copy(Xwrk(1), x0) x0_norm = x0%norm(); call Xwrk(1)%scal(one_cdp/x0_norm) else call Xwrk(1)%rand(.true.) endif allocate(H(kdim_+1, kdim_)) ; H = 0.0_dp allocate(eigvecs_wrk(kdim_, kdim_)) ; eigvecs_wrk = 0.0_dp allocate(eigvals_wrk(kdim_)) ; eigvals_wrk = 0.0_dp allocate(residuals_wrk(kdim_)) ; residuals_wrk = 0.0_dp ! Ritz eigenpairs computation. H = 0.0_dp kstart = 1 ; conv = 0 ; niter = 0 ; krst = 1 krylovschur: do while (conv < nev) arnoldi_factorization: do k = kstart, kdim_ ! Arnoldi step. call arnoldi(A, Xwrk, H, info, kstart=k, kend=k, transpose=transpose) call check_info(info, 'arnoldi', this_module, this_procedure) ! Spectral decomposition of the k x k Hessenberg matrix. eigvals_wrk = 0.0_dp ; eigvecs_wrk = 0.0_dp if (time_lightkrylov()) call timer%start('eig') call eig(H(:k, :k), eigvecs_wrk(:k, :k), eigvals_wrk(:k)) if (time_lightkrylov()) call timer%stop('eig') ! Compute residuals. beta = H(k+1, k) residuals_wrk(:k) = compute_residual_cdp(beta, eigvecs_wrk(k,:k)) ! Check convergence. niter = niter + 1 conv = count(residuals_wrk(:k) < tol) write(msg,'(I0,A,I0,A,I0,A)') conv, '/', nev, ' eigenvalues converged after ', niter, & & ' steps of the Arnoldi process.' call log_information(msg, this_module, this_procedure) if (outpost) call write_results_cdp(eigs_output, eigvals_wrk(:k), residuals_wrk(:k), tol) if (conv >= nev) exit arnoldi_factorization enddo arnoldi_factorization write(msg,'(I0,A,I0,A,I0,A)') conv, '/', nev, ' eigenvalues converged after ', krst, & & ' Krylov-Schur restarts of the Arnoldi process.' call log_information(msg, this_module, this_procedure) ! Krylov-Schur restarting procedure. krst = krst + 1 call krylov_schur(kstart, Xwrk, H, median_selector) ; kstart = kstart + 1 end do krylovschur !-------------------------------- !----- POST-PROCESS ----- !-------------------------------- block integer :: indices(kdim_) real(dp) :: abs_eigvals(kdim_) ! Re-compute eigenvalues and eigenvectors. k = min(k, kdim_); if (time_lightkrylov()) call timer%start('eig') call eig(H(:k, :k), eigvecs_wrk(:k, :k), eigvals_wrk(:k)) if (time_lightkrylov()) call timer%stop('eig') ! Sort eigenvalues. abs_eigvals = abs(eigvals_wrk) ; call sort_index(abs_eigvals, indices, reverse=.true.) eigvals_wrk = eigvals_wrk(indices) ; eigvecs_wrk = eigvecs_wrk(:, indices) residuals_wrk = residuals_wrk(indices) ! Store converged eigenvalues. eigvals = eigvals_wrk(:nev) ; residuals = residuals_wrk(:nev) end block ! Construct eigenvectors. do i = 1, nev call X(i)%zero() do j = 1, k call X(i)%axpby(eigvecs_wrk(j, i), Xwrk(j), one_cdp) enddo enddo info = niter if (time_lightkrylov()) call timer%stop(this_procedure) contains function median_selector(lambda) result(selected) implicit none(type, external) complex(dp), intent(in) :: lambda(:) logical, allocatable :: selected(:) selected = abs(lambda) > median(abs(lambda)) end function median_selector end subroutine eigs_cdp end module LightKrylov_IterativeSolvers