Octopus
solvers_oct_m Module Reference

This module is intended to contain "only mathematical" functions and procedures. More...

Detailed Description

This module is intended to contain "only mathematical" functions and procedures.

Data Types

interface  dconjugate_gradients
 
interface  ddotp_i
 
interface  dnrm_i
 
interface  doperator_i
 
interface  zconjugate_gradients
 
interface  zdotp_i
 
interface  znrm_i
 
interface  zoperator_i
 

Functions/Subroutines

subroutine zsym_conjugate_gradients (np, x, b, op, dotp, iter, residue, threshold)
 The two following subroutines, sym_conjugate_gradients, and bi_conjugate_gradients, must be called under a common interface: conjugate_gradients. It provides an approximate solution to the linear system problem Ax = b. Solving a symmetric linear system, which is either real or complex symmetric or Hermitian the best choice is sym_conjugate_gradients (does not need \( A^\dagger \)) Solving a real unsymmetric or a complex non-Hermitian system bi_conjugate_gradients has to be chosen (where one does need \( A^\dagger \)). More...
 
subroutine zbi_conjugate_gradients (np, x, b, op, opt, dotp, iter, residue, threshold)
 
subroutine, public zqmr_sym_gen_dotu (np, x, b, op, dotu, nrm2, prec, iter, residue, threshold, showprogress, converged, use_initial_guess)
 for complex symmetric matrices W Chen and B Poirier, J Comput Phys 219, 198-209 (2006) More...
 
subroutine, public zqmr_gen_dotu (np, x, b, op, opt, dotu, nrm2, prec, prect, iter, residue, threshold, showprogress, converged)
 for general complex matrices taken from 'An Implementation of the QMR Method based on Coupled Two-Term Recurrences' by R. W. Freund and N. M. Nachtigal (page 25) http: More...
 
complex(real64) function, dimension(size(b, 1), size(b, 2)), public zidrs (b, s, preconditioner, matrixvector, ddotprod, zdotprod, tolerance, maximum_iterations, variant, flag, relres, iterations, x0, U0, omega, resvec, H)
 This is the "Induced Dimension Reduction", IDR(s) (for s=4). IDR(s) is a robust and efficient short recurrence Krylov subspace method for solving large nonsymmetric systems of linear equations. It is described in [Peter Sonneveld and Martin B. van Gijzen, SIAM J. Sci. Comput. 31, 1035 (2008)]. More...
 
subroutine dsym_conjugate_gradients (np, x, b, op, dotp, iter, residue, threshold)
 The two following subroutines, sym_conjugate_gradients, and bi_conjugate_gradients, must be called under a common interface: conjugate_gradients. It provides an approximate solution to the linear system problem Ax = b. Solving a symmetric linear system, which is either real or complex symmetric or Hermitian the best choice is sym_conjugate_gradients (does not need \( A^\dagger \)) Solving a real unsymmetric or a complex non-Hermitian system bi_conjugate_gradients has to be chosen (where one does need \( A^\dagger \)). More...
 
subroutine dbi_conjugate_gradients (np, x, b, op, opt, dotp, iter, residue, threshold)
 
subroutine, public dqmr_sym_gen_dotu (np, x, b, op, dotu, nrm2, prec, iter, residue, threshold, showprogress, converged, use_initial_guess)
 for complex symmetric matrices W Chen and B Poirier, J Comput Phys 219, 198-209 (2006) More...
 
subroutine, public dqmr_gen_dotu (np, x, b, op, opt, dotu, nrm2, prec, prect, iter, residue, threshold, showprogress, converged)
 for general complex matrices taken from 'An Implementation of the QMR Method based on Coupled Two-Term Recurrences' by R. W. Freund and N. M. Nachtigal (page 25) http: More...
 
real(real64) function, dimension(size(b, 1), size(b, 2)), public didrs (b, s, preconditioner, matrixvector, ddotprod, zdotprod, tolerance, maximum_iterations, variant, flag, relres, iterations, x0, U0, omega, resvec, H)
 This is the "Induced Dimension Reduction", IDR(s) (for s=4). IDR(s) is a robust and efficient short recurrence Krylov subspace method for solving large nonsymmetric systems of linear equations. It is described in [Peter Sonneveld and Martin B. van Gijzen, SIAM J. Sci. Comput. 31, 1035 (2008)]. More...
 

Function/Subroutine Documentation

◆ zsym_conjugate_gradients()

subroutine solvers_oct_m::zsym_conjugate_gradients ( integer, intent(in)  np,
complex(real64), dimension(:), intent(inout), contiguous  x,
complex(real64), dimension(:), intent(in)  b,
  op,
procedure(zdotp_i dotp,
integer, intent(inout)  iter,
real(real64), intent(out), optional  residue,
real(real64), intent(in), optional  threshold 
)
private

The two following subroutines, sym_conjugate_gradients, and bi_conjugate_gradients, must be called under a common interface: conjugate_gradients. It provides an approximate solution to the linear system problem Ax = b. Solving a symmetric linear system, which is either real or complex symmetric or Hermitian the best choice is sym_conjugate_gradients (does not need \( A^\dagger \)) Solving a real unsymmetric or a complex non-Hermitian system bi_conjugate_gradients has to be chosen (where one does need \( A^\dagger \)).

Note: the complex-valued versions (both CG and BiCG) work only with symmetric operators but they may be non-Hermitian. This is a property of the CG algorithm. A comment on this can be found in chapter 4.1 of ftp:

subroutine conjugate_gradients(np, x, b, op, [opt,] dotp, iter [, residue] [, threshold] ) integer, intent(in) :: np => The dimension of the problem. real(real64), intent(inout) :: x => On input, an estimate to the solution. => On output, the approximate solution. real(real64), intent(in) :: b => The inhomogeneous term of the equation. interface subroutine op(x, y) real(real64), intent(in) :: x(:) real(real64), intent(out) :: y(:) end subroutine op end interface => This should be a procedure that computes \( Ax = y \). interface subroutine opt(x, y) real(real64), intent(in) :: x(:) real(real64), intent(out) :: y(:) end subroutine opt end interface => If present, this should be a procedure that computes \( A^\dagger x = y \). Only useful for non-Hermitian operators. interface complex(real64) function dotp(x, y) complex(real64), intent(inout) :: x(:) complex(real64), intent(in) :: y(:) end function dotp => Calculates the \( <x | y> \). Depending on the matrix A one should choose: complex symmetric: \( <x | y> = x^T * y \) hermitian: \( <x | y> = x^\dagger * y \) general: \( <x | y> = x^\dagger * y \) integer, intent(inout) :: iter => On input, the maximum number of iterations that the procedure is allowed to take. On output, the iterations it actually did. real(real64), intent(out) :: residue => If present, it measures the final error: residue = \( <Ax - b | Ax - b> \). real(real64), intent(in) :: threshold => If present, it sets the required accuracy threshold for the algorithm to stop. If not present, this is set to 1.0e-6. [The algorithm stops when \( <Ax - b | Ax - b> <= threshold \), or iter iterations are reached.] end subroutine conjugate_gradients

(*) NOTE: The algorithm assumes that the vectors are given in an orthonormal basis.

Definition at line 330 of file solvers.F90.

◆ zbi_conjugate_gradients()

subroutine solvers_oct_m::zbi_conjugate_gradients ( integer, intent(in)  np,
complex(real64), dimension(:), intent(inout)  x,
complex(real64), dimension(:), intent(in)  b,
  op,
  opt,
procedure(zdotp_i dotp,
integer, intent(inout)  iter,
complex(real64), intent(out), optional  residue,
real(real64), intent(in), optional  threshold 
)
private

Definition at line 399 of file solvers.F90.

◆ zqmr_sym_gen_dotu()

subroutine, public solvers_oct_m::zqmr_sym_gen_dotu ( integer, intent(in)  np,
complex(real64), dimension(:), intent(inout), contiguous  x,
complex(real64), dimension(:), intent(in)  b,
procedure(zoperator_i op,
procedure(zdotp_i dotu,
procedure(znrm_i nrm2,
procedure(zoperator_i prec,
integer, intent(inout)  iter,
real(real64), intent(out), optional  residue,
real(real64), intent(in), optional  threshold,
logical, intent(in), optional  showprogress,
logical, intent(out), optional  converged,
logical, intent(in), optional  use_initial_guess 
)

for complex symmetric matrices W Chen and B Poirier, J Comput Phys 219, 198-209 (2006)

Parameters
[in]npnumber of points
[in,out]xthe initial guess and the result
[in]bthe right side
opthe matrix A as operator
dotuthe dot product (must be \( x^T*y \), not daggered)
nrm2the 2-norm of the vector x
precthe preconditioner
[in,out]iter(in) the maximum number of iterations, (out) used iterations
[out]residuethe residue = abs(Ax-b)
[in]thresholdconvergence threshold
[in]showprogressshould there be a progress bar
[out]convergedhas the algorithm converged
[in]use_initial_guessdo we have a guess or not

Definition at line 486 of file solvers.F90.

◆ zqmr_gen_dotu()

subroutine, public solvers_oct_m::zqmr_gen_dotu ( integer, intent(in)  np,
complex(real64), dimension(:), intent(inout)  x,
complex(real64), dimension(:), intent(in)  b,
procedure(zoperator_i op,
procedure(zoperator_i opt,
procedure(zdotp_i dotu,
procedure(znrm_i nrm2,
procedure(zoperator_i prec,
procedure(zoperator_i prect,
integer, intent(inout)  iter,
real(real64), intent(out), optional  residue,
real(real64), intent(in), optional  threshold,
logical, intent(in), optional  showprogress,
logical, intent(out), optional  converged 
)

for general complex matrices taken from 'An Implementation of the QMR Method based on Coupled Two-Term Recurrences' by R. W. Freund and N. M. Nachtigal (page 25) http:

Parameters
[in]npnumber of points
[in,out]xinitial guess and result
[in]bright side
opthe matrix A as operator: y <- A*x
optthe transposed matrix A as operator: y <- A^T*x
dotuthe dot product (must be \( x^T*y \), not daggered)
nrm2the 2-norm of the vector x
precpreconditioner
precttransposed preconditioner
[in,out]iter(in) the maximum number of iterations, (out) used iterations
[out]residuethe residue = abs(Ax-b)
[in]thresholdconvergence threshold
[in]showprogressshould there be a progress bar
[out]convergedhas the algorithm converged

Definition at line 710 of file solvers.F90.

◆ zidrs()

complex(real64) function, dimension (size(b,1), size(b,2)), public solvers_oct_m::zidrs ( complex(real64), dimension(:, :), intent(in)  b,
integer, intent(in)  s,
  preconditioner,
  matrixvector,
  ddotprod,
  zdotprod,
real(real64), intent(in), optional  tolerance,
integer, intent(in), optional  maximum_iterations,
character(len=8), intent(in), optional  variant,
integer, intent(out), optional  flag,
real(real64), intent(out), optional  relres,
integer, intent(out), optional  iterations,
complex(real64), dimension(:, :), intent(in), optional  x0,
complex(real64), dimension(:, :, :), intent(in), optional  U0,
complex(real64), dimension(:), intent(in), optional  omega,
real(real64), dimension(:), intent(out), optional  resvec,
complex(real64), dimension(:, :), intent(out), optional  H 
)

This is the "Induced Dimension Reduction", IDR(s) (for s=4). IDR(s) is a robust and efficient short recurrence Krylov subspace method for solving large nonsymmetric systems of linear equations. It is described in [Peter Sonneveld and Martin B. van Gijzen, SIAM J. Sci. Comput. 31, 1035 (2008)].

We have adapted the code released by M. B. van Gizjen [http: That code is licenced under the MIT licence, which allows to re-release this modified version under the GPL.

The original copyright notice, and licence, follows:

Copyright (c) July 2015, Martin van Gijzen

Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the "Software"), to deal in the Software without restriction, including without limitation the rights to use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of the Software, and to permit persons to whom the Software is furnished to do so, subject to the following conditions:

The above copyright notice and this permission notice shall be included in all copies or substantial portions of the Software.

THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.

Definition at line 924 of file solvers.F90.

◆ dsym_conjugate_gradients()

subroutine solvers_oct_m::dsym_conjugate_gradients ( integer, intent(in)  np,
real(real64), dimension(:), intent(inout), contiguous  x,
real(real64), dimension(:), intent(in)  b,
  op,
procedure(ddotp_i dotp,
integer, intent(inout)  iter,
real(real64), intent(out), optional  residue,
real(real64), intent(in), optional  threshold 
)
private

The two following subroutines, sym_conjugate_gradients, and bi_conjugate_gradients, must be called under a common interface: conjugate_gradients. It provides an approximate solution to the linear system problem Ax = b. Solving a symmetric linear system, which is either real or complex symmetric or Hermitian the best choice is sym_conjugate_gradients (does not need \( A^\dagger \)) Solving a real unsymmetric or a complex non-Hermitian system bi_conjugate_gradients has to be chosen (where one does need \( A^\dagger \)).

Note: the complex-valued versions (both CG and BiCG) work only with symmetric operators but they may be non-Hermitian. This is a property of the CG algorithm. A comment on this can be found in chapter 4.1 of ftp:

subroutine conjugate_gradients(np, x, b, op, [opt,] dotp, iter [, residue] [, threshold] ) integer, intent(in) :: np => The dimension of the problem. real(real64), intent(inout) :: x => On input, an estimate to the solution. => On output, the approximate solution. real(real64), intent(in) :: b => The inhomogeneous term of the equation. interface subroutine op(x, y) real(real64), intent(in) :: x(:) real(real64), intent(out) :: y(:) end subroutine op end interface => This should be a procedure that computes \( Ax = y \). interface subroutine opt(x, y) real(real64), intent(in) :: x(:) real(real64), intent(out) :: y(:) end subroutine opt end interface => If present, this should be a procedure that computes \( A^\dagger x = y \). Only useful for non-Hermitian operators. interface real(real64) function dotp(x, y) real(real64), intent(inout) :: x(:) real(real64), intent(in) :: y(:) end function dotp => Calculates the \( <x | y> \). Depending on the matrix A one should choose: complex symmetric: \( <x | y> = x^T * y \) hermitian: \( <x | y> = x^\dagger * y \) general: \( <x | y> = x^\dagger * y \) integer, intent(inout) :: iter => On input, the maximum number of iterations that the procedure is allowed to take. On output, the iterations it actually did. real(real64), intent(out) :: residue => If present, it measures the final error: residue = \( <Ax - b | Ax - b> \). real(real64), intent(in) :: threshold => If present, it sets the required accuracy threshold for the algorithm to stop. If not present, this is set to 1.0e-6. [The algorithm stops when \( <Ax - b | Ax - b> <= threshold \), or iter iterations are reached.] end subroutine conjugate_gradients

(*) NOTE: The algorithm assumes that the vectors are given in an orthonormal basis.

Definition at line 1559 of file solvers.F90.

◆ dbi_conjugate_gradients()

subroutine solvers_oct_m::dbi_conjugate_gradients ( integer, intent(in)  np,
real(real64), dimension(:), intent(inout)  x,
real(real64), dimension(:), intent(in)  b,
  op,
  opt,
procedure(ddotp_i dotp,
integer, intent(inout)  iter,
real(real64), intent(out), optional  residue,
real(real64), intent(in), optional  threshold 
)
private

Definition at line 1628 of file solvers.F90.

◆ dqmr_sym_gen_dotu()

subroutine, public solvers_oct_m::dqmr_sym_gen_dotu ( integer, intent(in)  np,
real(real64), dimension(:), intent(inout), contiguous  x,
real(real64), dimension(:), intent(in)  b,
procedure(doperator_i op,
procedure(ddotp_i dotu,
procedure(dnrm_i nrm2,
procedure(doperator_i prec,
integer, intent(inout)  iter,
real(real64), intent(out), optional  residue,
real(real64), intent(in), optional  threshold,
logical, intent(in), optional  showprogress,
logical, intent(out), optional  converged,
logical, intent(in), optional  use_initial_guess 
)

for complex symmetric matrices W Chen and B Poirier, J Comput Phys 219, 198-209 (2006)

Parameters
[in]npnumber of points
[in,out]xthe initial guess and the result
[in]bthe right side
opthe matrix A as operator
dotuthe dot product (must be \( x^T*y \), not daggered)
nrm2the 2-norm of the vector x
precthe preconditioner
[in,out]iter(in) the maximum number of iterations, (out) used iterations
[out]residuethe residue = abs(Ax-b)
[in]thresholdconvergence threshold
[in]showprogressshould there be a progress bar
[out]convergedhas the algorithm converged
[in]use_initial_guessdo we have a guess or not

Definition at line 1715 of file solvers.F90.

◆ dqmr_gen_dotu()

subroutine, public solvers_oct_m::dqmr_gen_dotu ( integer, intent(in)  np,
real(real64), dimension(:), intent(inout)  x,
real(real64), dimension(:), intent(in)  b,
procedure(doperator_i op,
procedure(doperator_i opt,
procedure(ddotp_i dotu,
procedure(dnrm_i nrm2,
procedure(doperator_i prec,
procedure(doperator_i prect,
integer, intent(inout)  iter,
real(real64), intent(out), optional  residue,
real(real64), intent(in), optional  threshold,
logical, intent(in), optional  showprogress,
logical, intent(out), optional  converged 
)

for general complex matrices taken from 'An Implementation of the QMR Method based on Coupled Two-Term Recurrences' by R. W. Freund and N. M. Nachtigal (page 25) http:

Parameters
[in]npnumber of points
[in,out]xinitial guess and result
[in]bright side
opthe matrix A as operator: y <- A*x
optthe transposed matrix A as operator: y <- A^T*x
dotuthe dot product (must be \( x^T*y \), not daggered)
nrm2the 2-norm of the vector x
precpreconditioner
precttransposed preconditioner
[in,out]iter(in) the maximum number of iterations, (out) used iterations
[out]residuethe residue = abs(Ax-b)
[in]thresholdconvergence threshold
[in]showprogressshould there be a progress bar
[out]convergedhas the algorithm converged

Definition at line 1939 of file solvers.F90.

◆ didrs()

real(real64) function, dimension (size(b,1), size(b,2)), public solvers_oct_m::didrs ( real(real64), dimension(:, :), intent(in)  b,
integer, intent(in)  s,
  preconditioner,
  matrixvector,
  ddotprod,
  zdotprod,
real(real64), intent(in), optional  tolerance,
integer, intent(in), optional  maximum_iterations,
character(len=8), intent(in), optional  variant,
integer, intent(out), optional  flag,
real(real64), intent(out), optional  relres,
integer, intent(out), optional  iterations,
real(real64), dimension(:, :), intent(in), optional  x0,
real(real64), dimension(:, :, :), intent(in), optional  U0,
real(real64), dimension(:), intent(in), optional  omega,
real(real64), dimension(:), intent(out), optional  resvec,
real(real64), dimension(:, :), intent(out), optional  H 
)

This is the "Induced Dimension Reduction", IDR(s) (for s=4). IDR(s) is a robust and efficient short recurrence Krylov subspace method for solving large nonsymmetric systems of linear equations. It is described in [Peter Sonneveld and Martin B. van Gijzen, SIAM J. Sci. Comput. 31, 1035 (2008)].

We have adapted the code released by M. B. van Gizjen [http: That code is licenced under the MIT licence, which allows to re-release this modified version under the GPL.

The original copyright notice, and licence, follows:

Copyright (c) July 2015, Martin van Gijzen

Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the "Software"), to deal in the Software without restriction, including without limitation the rights to use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of the Software, and to permit persons to whom the Software is furnished to do so, subject to the following conditions:

The above copyright notice and this permission notice shall be included in all copies or substantial portions of the Software.

THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.

Definition at line 2153 of file solvers.F90.