Octopus
|
This module is intended to contain "only mathematical" functions and procedures. More...
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... | |
|
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.
|
private |
Definition at line 399 of file solvers.F90.
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)
[in] | np | number of points |
[in,out] | x | the initial guess and the result |
[in] | b | the right side |
op | the matrix A as operator | |
dotu | the dot product (must be \( x^T*y \), not daggered) | |
nrm2 | the 2-norm of the vector x | |
prec | the preconditioner | |
[in,out] | iter | (in) the maximum number of iterations, (out) used iterations |
[out] | residue | the residue = abs(Ax-b) |
[in] | threshold | convergence threshold |
[in] | showprogress | should there be a progress bar |
[out] | converged | has the algorithm converged |
[in] | use_initial_guess | do we have a guess or not |
Definition at line 486 of file solvers.F90.
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:
[in] | np | number of points |
[in,out] | x | initial guess and result |
[in] | b | right side |
op | the matrix A as operator: y <- A*x | |
opt | the transposed matrix A as operator: y <- A^T*x | |
dotu | the dot product (must be \( x^T*y \), not daggered) | |
nrm2 | the 2-norm of the vector x | |
prec | preconditioner | |
prect | transposed preconditioner | |
[in,out] | iter | (in) the maximum number of iterations, (out) used iterations |
[out] | residue | the residue = abs(Ax-b) |
[in] | threshold | convergence threshold |
[in] | showprogress | should there be a progress bar |
[out] | converged | has the algorithm converged |
Definition at line 710 of file solvers.F90.
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.
|
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.
|
private |
Definition at line 1628 of file solvers.F90.
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)
[in] | np | number of points |
[in,out] | x | the initial guess and the result |
[in] | b | the right side |
op | the matrix A as operator | |
dotu | the dot product (must be \( x^T*y \), not daggered) | |
nrm2 | the 2-norm of the vector x | |
prec | the preconditioner | |
[in,out] | iter | (in) the maximum number of iterations, (out) used iterations |
[out] | residue | the residue = abs(Ax-b) |
[in] | threshold | convergence threshold |
[in] | showprogress | should there be a progress bar |
[out] | converged | has the algorithm converged |
[in] | use_initial_guess | do we have a guess or not |
Definition at line 1715 of file solvers.F90.
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:
[in] | np | number of points |
[in,out] | x | initial guess and result |
[in] | b | right side |
op | the matrix A as operator: y <- A*x | |
opt | the transposed matrix A as operator: y <- A^T*x | |
dotu | the dot product (must be \( x^T*y \), not daggered) | |
nrm2 | the 2-norm of the vector x | |
prec | preconditioner | |
prect | transposed preconditioner | |
[in,out] | iter | (in) the maximum number of iterations, (out) used iterations |
[out] | residue | the residue = abs(Ax-b) |
[in] | threshold | convergence threshold |
[in] | showprogress | should there be a progress bar |
[out] | converged | has the algorithm converged |
Definition at line 1939 of file solvers.F90.
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.