|
| |
Numerical Library
- CG
conjugate gradient
- stable and robust
- symmetric positive defnite
- PCG
preconditioned conjugate gradient
ICCG = incomplete Choleski decomposition
MICCG = modified incomplete Choleski decomposition (MODULE
miccg3DdV)
- CGNE
conjugate gradient on normal equations
- minimization of error
AAHx=AHb
Craig's method
- CGNR
AHAx=AHb
conjugate gradient on residual
- minimization of residual, equivalent to LSQR
(+ row protection method as preconditioning = with use of block Jacobi or
block SSOR preconditioning)
CGNR and CGNE converge slowly for nonhermitian A (therefore: GMRES, restarted
GMRES(m))
- CGS
conjugate gradient squared
- BiCG
biconjugate gradient
- BiCGSTAB (Vorst)
biconjugate gradient stabilized
- based on unsymmetric Lanczos
- GMRES
generalised minimal residual
- guaranteed to converge if A is definite
- GMRES(m)
restarted GMRES
- FMRES
CG with variable preconditioning Flexible GMRES
- GMRESR
- LSQR
- GMERR
- QMR (Freund)
quasi-minimal residual
- QMR uses look-ahead Lanczos process to generate basis vectors
- more stable and robust than GMRES, BiCGSTAB
- half the work and storage for hermitian/symmetric matrices
- + Sonneveld's trick
- QMR is mathematically equivalent to MINRES in the special case of A
Hermititan.
- MINRES and SYMMLQ require a symmetric positive definite preconditioner
wheresas QMR only a symmetric one.
- QMR will converge fast if the eigenvalues are clustered aound 1.
- QMRSTAB
- MINRES
- only symmetric positive definite preconditioner
- without preconditioning symmetric QMR and MINRES are mathematically
equivalent
- SYMMLQ
- only symmetric positive definite preconditioner
- less sensitive to ill-conditioned A than MINRES but often takes more steps
- Both SYMMLQ and MINRES use the Hermitian Lanczos recursion to generate an
orthonormal basis for the Krylov subspaces (3-term recurrence).
- SYMMBK
- slightly less expensive variant of SYMMLQ
- TFQMR (Freund)
transpose-free quasi-minimal residual
no need for AT, only A is need
- ORTHODIR
- guaranteed to converge if A is definite
- shift-inverse
- ORTHOMIN
- similar to GMRES
- similar to generalized CG residual method
- ORTHORES
- symmetric positive definite
guaranteed to converge: CG, PCG, ICCG, MICCG, ICCG/MICCG
- Choleski decomposition leads to trianguar decomposition
- symmetric indefinite
QMR, MINRES, SYMMLQ
- SYMMLQ less sensitive to ill-conditioned A than MINRES but often takes more
steps
- MINRES and SYMMLQ don't allow an indefinite preconditioner, but QMR does.
(CGNE, CGNR)
- SYMMLQ and MINRES are standard algorithms for symmetric indefinite systems
but they require a symmetric positive preconditioner M=M1M2
where M2=M1T. This is a problem for highly
indefinite systems as you need a preconditioner that is also highly indefinite
to better approximate A.
- Solution: use symmetric QMR. This is an adaption to the general
nonsymmetric QMR.
- M-Matrix
(aij<0 for i unequal j, A is nonsingular, A-1>0)
PCG (VORST method ICCG-MICCG)
- H-Matrix
(diagonally dominant (somewhat))
- hermitian positive definite
- hermitian indefinite
QMR, TFQMR, GMRES, BiCGSTAB
- complex symmetric
QMR, TFQMR, GMRES, BiCGSTAB
- definite
(i.e. the symmetric part (A+AT)/2 or its negative is
symmetric positive definite)
guaranteed to converge: ORTHODIR, GMRES
- unsymmetric
QMR, TFQMR, GMRES, GMRES(m), BiCGSTAB, CGR
- Jacobi
D-1
- SSOR
(L+D)D-1(D+LT)
- mixed preconditioning is more efficient because of Eisenstat's trick
- The CGNE/SSOR and CGNR/SSOR will not brake down if A is nonsingular since
then the matrices ATA or AAT are symmetric positive
definite.
- extend SSOR preconditioner to indefinite matrices
- CGNE/ILQ (Craig's method)
- CGNR/ILQ
- LU
lower-upper factorization
- for indefinite matrices the LU decomposition could lead to numerical
instabilities
- ILU
incomplete LU preconditioning with positional dropping ILU(k)
ILU for indefinite matrices:
- LU sometimes inaccurate, i.e. not close to A, due to small pivots
- L or U sometimes not diagonally dominant
- ILUT
incomplete LU preconditioning with numerical dropping ILUT(k)
- MILU
modified incomplete LU preconditioning with positional dropping MILU(k)
- incomplete QR (Saad)
- polynomial preconditioning
It is difficult to find a good polynomial in the case of indefinite matrices.
- use symmetric positive definite preconditioner to cluster eigenvalues
around 1 and -1
use symmetric indefinite preconditioner to cluster eigenvalues around 1
- Hermitian matrices have special spectral properties.
- Roughly speaking, Krylov subspace methods are most effective for
coefficient matrices A whose spectrum, except for possibly a few isolated
eigenvalues, is contained in a half-plane which excludes the origin of the
complex plane.
All Krylov.subspace iterations converge very slowly for highly indefinite
coefficient matrices.
Krylov space methods normally do not require any a priori spectral information
and there are methods in this class especially designed for indefinite and
nonsymmetric problems.
- We should point out, however, that the incomplete Cholesky factorisation
is not guaranteed to exist for an arbitrary positive definite symmetric
matrix. (should be somehow diagonal dominant)
- The convergence of conjugate gradients depends on the distribution of the
eigenvalues (spectrum) of the matrix.
- The main difficulty in indefinite preconditioning can be explained as
follows. The Krylov subspace methods converge fast when the eigenvalues are
clustered, say aroung 1. This means, that the preconditioner, often viewed as
an approximation to the inverse of the given matrix, must transfer the
eigenvalues to 1. It may easily happen that, due to approximation errors that
are intentionally made in order to keep the process efficient, eigenvalues are
transferred to values close to 0. In that case the convergence of Krylov
subspace methods will be slow, and it may happen that the unpreconditioned
iteration process converges faster than the preconditioned iteration.
Complex-valued problems can be treated as the problems discussed in this
paper, in particular Hermitian problems can be treated as real symmetric ones
by replacing AT by AH , and by the appropriate change of
innerproduct.
- The Lanczos method reduces to only one recursion for the case of A
Hermitian or complex symmetric.
The MODULE
fermi_functions contains several subroutines and functions for the evaluation of Fermi-Dirac
integrals
and the exponential function.
MODULE EXPONENT_FUNCTION
function exponentV(argV,typeC,factorlesshuge) result(valueV)
!-------------------------------------------------------------------------------
! Function allows one to call predefined functions of exp(argV) for any argument
! from -INF to +INF and returns an accurate value. Routine serves purpose to
! prevent overflow and roundoff-errors by using appropriate reordering or
! Talyor expansion before calling exp-function. Max. error is of the order
! of 10^-32 on Suns and Decs.
!
! depending on typeC, returns vector-valued function of exp(argV)
! input ..
! argV = V, vector, argument of exponential
! maybe scalar but result is vector of at least length 1
! typeC .. character that specifies function to be evaluated
! factorlesshuge .. optional argument relevant only for typeC='exp'
! ensures exp(argV) < HUGE * factorlesshuge
! where HUGE is largest representable number (10^308)
! example: factorlesshuge=1.d-20 causes upper bound of result
! to be HUGE * 1.d-20 = 10^(308-20)
! output ..
! valueV , vector of same length as argV, whereby
!
! typeC = 'exp' ... valueV = exp(V)
! typeC = 'fermi' ... valueV = 1/(1+exp(V))
! typeC = 'bernexp' ... valueV = V exp(V)/(exp(V)-1)
! typeC = 'bern' ... valueV = V /(exp(V)-1)
! typeC = 'lnexp' ... valueV = ln(1+exp(V))
!
-------------------------------------------------------------------------------------------------------------
!------------------------------------------------------------------------------
FUNCTION fermi_trap_levelV(argV,gV,ch) result(valueV)
!----------------------------------------------------------------
!
! calculates 1/(1+gV*exp(argV)) for ch='fc'
! -exp(argV)/(1+gV*exp(argV))^2 for ch='df'
!
!----------------------------------------------------------------
MODULE FERMI_FUNCTIONS
!--------------------------------------------------------------------------------
!
! 1. FUNCTION fermi3_2
!
! 2. FUNCTION fermi1_2
!
! 3. FUNCTION fermi_min_1_2
!
! 4. FUNCTION fermi_min_3_2
!
!--------------------------------------------------------------------------------
!---------------------------------------------------------------------------------
FUNCTION fermi1_2(x,ch,err) result(erg)
!---------------------------------------------------------------------------------
!
! calculates fermiintegral of order j=1/2
!
! /
j
! 2 |
z
! --------- | ------------------- dz
! sqrt(pi) | exp(z-x) + 1
! /
!
! in this case VAN HALEN/PULFREY is taken
! err is set to 1.0d10
!
! VAN HALEN attains a precision of only 1.0d-4 to 1.0d-5 between 2.0 and 4.0
! therefore modified version of WONG et alias
!---------------------------------------------------------------------------------
!
! input:
!
! x ... real(8)
!
! ch ... character(len=1)
! = '1' : precision 5.0d-3 (AYMERICH-HUMET)
! = '2' : precision 5.0d-6 (VAN HALEN, PULFREY)
!
! err ... real(8) : precision , if ch='3' is chosen
!
!---------------------------------------------------------------------------------
!---------------------------------------------------------------------------------
FUNCTION fermi_min_1_2(x,ch,err) result(erg)
!---------------------------------------------------------------------------------
!
! calculates fermiintegral of order j=-1/2
!
! /
j
! 1 |
z
! ---------- | ------------------ dz
! sqrt(pi) | exp(z-x) + 1
! /
!
!---------------------------------------------------------------------------------
!
! input:
!
! x ... real(8)
!
! ch ... character(len=1)
!
! = '2','1' : precision 6.1d-6 (VAN HALEN, PULFREY)
!
! err ... real(8) : precision , if ch='3' is chosen
!
! VAN HALEN : basis functions from his paper
! coefficients calculated by sth
!---------------------------------------------------------------------------------
The Newton method in subroutine
newton is used to calculate the roots of a functional.
!------------------------------------------------------------------------------------------------------------------
SUBROUTINE
newton(xV,ctrl,iter,tolf,tolmin,stpmx,ALF,newton_correction,calculate_function,calculate_gradient)
!------------------------------------------------------------------------------------------------------------------
! based on Numerical Recipes
! modified by sth 23.10.1998
!
! calculates root of N-dim. nonlinear equation by a globally convergent Newton's
method
! On entry, xV is the initial guess
! ctrl = .true. , if routine has converged to a local minimum. Restart with
different initial value
! = .false. on normal return
! iter (output) number of iterations
! (input) maximum number of iterations
!
! tolf ... convergence criterion on function values
! tolmin ... criterion for deciding whether spurious convergence to a minimum
has occurred
! tolx ... convergence criterion on dx
! stpmx ... scaled maximum step length allowed in line searches
!
!
! User must provide following SUBROUTINES:
!
! SUBROUTINE newton_correction(dxV,xV,fV)
! given xV,fV (input) this subroutine calculates one newton correction step dxV
:
! fdjac*dxV = -fV
!
! SUBROUTINE calculate_function(xV,fV)
! calculates value of function whose root should be calculated
!
! SUBROUTINE calculate_gradient(xV,fV,gradientV)
! calculates gradient fV*fdjac for linesearch
! note: fmin=0.5*fV*fV ==> fV=0 <=> 0 = grad fmin = fV * fdjacV ^=
matmult(fV,fdjac)
!
!------------------------------------------------------------------------------------------------------------------
!=========================================================
subroutine
cg (n, xV, iter, resid, info, brhs_V , matvec , psolve)
!-----------------------------------------------------------------------------------
! Univ. of Tennessee and Oak Ridge National Laboratory
! October 1, 1993
! Details of this algorithm are described in "Templates for the
! Solution of Linear Systems: Building Blocks for Iterative
! Methods", Barrett, Berry, Chan, Demmel, Donato, Dongarra,
! Eijkhout, Pozo, Romine, and van der Vorst, SIAM Publications,
! 1993. (ftp netlib2.cs.utk.edu; cd linalg; get templates.ps).
!
! converted to f90 and modified by pv, 3-29-98 , modified by sth, 10-30-98
!
! CG solves the linear system PmatrixM * xV = brhs_V using the
! Conjugate Gradient iterative method with preconditioning.
!
! Convergence test: ( norm( brhs_V - PmatrixM * xV ) / norm( brhs_V ) ) < TOL.
! For other measures, see the above reference.
! --Done in CGREVCOM.
! input/output ....
! n on entry, the dimension of the matrix. unchanged on exit.
! xv(n) on input, the initial guess. this is commonly set to the zero
! vector. on exit, if info = 0, the iterated approximate solution.
! set by cgrevcom.
! iter on input, the maximum iterations to be performed. on output, actual
! number of iterations performed. set by cgrevcom.
! resid on input, the allowable convergence measure for
! norm( PchargeV - PmatrixM*xV ) / norm( PchargeV ). on output,
! the final value of this measure.set by cgrevcom.
!
! info (output) . set by cgrevcom()
!
! brhs_V right handside of equation to be solved
!
!
! matvec (subroutine)
! The user must provide a subroutine to perform the matrix-vector product
! yV := alpha*PmatrixM*xV + beta*yV
! where alpha and beta are scalars, xV and yV are vectors, and PmatrixM
! is the Poisson matrix. Vector xV must remain unchanged. The solution
! is over-written on vector yV. The call is:
!
! call matvec( alpha, xV, beta, yV )
!
!
! psolve ( subroutine)
! The user must provide a subroutine to perform the preconditioner
! solve routine for the linear system
!
! M*xV = bV,
!
! where xV and bV are given vectors, and M a matrix. Vector bV must
! remain unchanged. The solution is over-written on vector xV.
! The call is:
!
! call psolve( xV, bV)
!
! ============================================================
!================================================================
SUBROUTINE
BICGSTAB (N, xV, ITER, RESID, INFO, brhs_V , matvec , psolve)
!-------------------------------------------------------------------------------
!
! -- Iterative template routine --
! Univ. of Tennessee and Oak Ridge National Laboratory
! October 1, 1993
! Details of this algorithm are described in "Templates for the
! Solution of Linear Systems: Building Blocks for Iterative
! Methods", Barrett, Berry, Chan, Demmel, Donato, Dongarra,
! Eijkhout, Pozo, Romine, and van der Vorst, SIAM Publications,
! 1993. (ftp netlib2.cs.utk.edu; cd linalg; get templates.ps).
!
! converted to f90 and modified by pv, 4-29-98
!
! BICGSTAB solves the linear system JmatrixM*xV = RecombV using the
! BiConjugate Gradient Stabilized iterative method with
! preconditioning.
!
! Convergence test: ( norm( RecombV - JmatrixM*xV ) / norm( RecombV ) ) < TOL.
! For other measures, see the above reference.
!
! Arguments
! =========
!
! N (input) INTEGER.
! On entry, the dimension of the matrix.
! Unchanged on exit.
! xV (input/output) DOUBLE PRECISION array, dimension N.
! On input, the initial guess. This is commonly set to
! the zero vector.
! On exit, if INFO = 0, the iterated approximate solution.
! ITER (input/output) INTEGER
! On input, the maximum iterations to be performed.
! On output, actual number of iterations performed.
! RESID (input/output) DOUBLE PRECISION
! On input, the allowable convergence measure for
! norm(RecombV - JmatrixM*xV ) / norm( RecombV ).
! On output, the final value of this measure.
! preconC chosen preconditioner, can be one of
! 'jacbi' .. Jacobi preconditioner
! 'ident' .. no preconditioner applied
! 'dilu' .. dilu preconditioner
! INFO (output) INTEGER
! = 0: Successful exit. Iterated approximate solution returned
! > 0: Convergence to tolerance not achieved. This will be
! set to the number of iterations performed.
! < 0: Illegal input parameter, or breakdown occurred
! during iteration.
! Illegal parameter:
! -1: matrix dimension N < 0
! -2: can no longer occur in f90 version (pv)
! -3: Maximum number of iterations ITER <= 0.
! BREAKDOWN: If parameters RHO or OMEGA become smaller
! than some tolerance, the program will terminate.
! Here we check against tolerance BREAKTOL.
! -10: RHO < BREAKTOL: RHO and RTLD have become
! orthogonal.
! -11: OMEGA < BREAKTOL: S and T have become
! orthogonal relative to T'*T.
! BREAKTOL is set in function GETBREAK.
!
! input via module currentmatrix (passed to subroutines of bicgstab) ....
! JdiagV(n), JsubdiagV(n),JsupdiagV(n) ..current matrix in form of
! main, sub- and super-diagonal
! RecombV(n) ... right hand side of current equation
!
! workM (workspace) real array, dimension (n,7)
! Workspace for residual, direction vector, etc.
! Note that vectors R and S shared the same workspace.
!
! routines called ..
!
! matvec (subroutine)
! The user must provide a subroutine to perform the
! matrix-vector product
!
! yV := alpha*JmatrixM*xV + beta*yV,
!
! where alpha and beta are scalars, xV and yV are vectors,
! and JmatrixM is the current matrix. Vector xV must remain unchanged.
! The solution is over-written on vector yV.
!
! The call is:
!
! CALL matvec( ALPHA, xV, BETA, YV )
!
!
! PSOLVE (subroutine)
! The user must provide a subroutine to perform the
! preconditioner solve routine for the linear system
!
! M*xV = bV,
!
! where xV and bV are vectors, and M a matrix. Vector bV must
! remain unchanged.
! The solution is over-written on vector bV.
!
! The call is:
! CALL PSOLVE( xV, bV )
!
! ==============================================================
subroutine
eigensolver_real_arpack(NDIM,EIGENVECTORL,WHICH_EIGVALC,MAX_EIGVAL,&
TOLERANCE,MAX_ITERATIONS,SIGMA,RESID_LIN_EQ,IPROBLEM,NC, &
eigenvalueV,eigenvectorM,lin_equation_solver_arnoldi,matvec)
!-----------------------------------------------------------------------------
! diagonalize symmetric Hamiltonian matrix with Arnoldi iteration
! compute max_eigval lowest eigenvalues
! input ...
! ndim .. actual size of Hamilton = number of grid points
! eigenvectorL .. T calculate eigenvalues+eigenvectors, F only former
! which_eigvalC(2) .. which eigenvalues are to be computed, can be
! 'LA','SA','LM','SM','BE' ... 1.letter:
! large/small/bothends,
! 2.letter:algebraic/in magnitude
! max_eigval .. compute that many eigenvalues <=ndim,
! starting with lowest
! tolerance ... convergence criterion, may be set to 0.d0
! (-> epsilon(0d0))
! max_iterations .. maximum number of iterations
! sigma ........ in shift-invert mode,
! eigenvalues near a point sigma_si are sought
!
! output ...
! eigenvalueV(1:max_eigval) ... eigenvalues in eV
! eigenvectorM(1:ndim, 1:max_eigval) ...eigenvectors psi,
! nc ... number of eigenvalues that converged
!
!-----------------------------------------------------------------------------
!
! USER must provide following routines as input:
!
! 1. SUBROUTINE lin_equation_solver_arnoldi(ndim,yV,xV,resid_lin_eq,ierr,&
! sigma)
! calculates yV = inv[A-sigma*I]*xV
! resid_lin_eq ... tolerance
! ierr .. output, if ierr/=0 routine failed
! sigma .. real number, see above
!
! 2. SUBROUTINE matvec(ndim,xV,yV)
!
! calculates yV = Matrix*xV
!
!------------------------------------------------
!
! uses arpack and lapack routines (contained in arpack_sg.f)
!
!-----------------------------------------------------------------------------
subroutine
eigensolver_realsym(ndim,num_ev,eigenvaluesV,eigenfunctionsM, &
init,itere_max,itsub_max,eps_itere,eps_itsub,info_itere,info_itsub, &
iwrite,matrix_vector,precond,inner_product)
!------------------------------------------------------------------------------
!
! eigensolver for real symmetric matrices (symmetric with respect to the
! inner product)
! calculates the num_ev LOWEST eigenvalues
! matrix is only accessed through the three user-defined subroutines
!
!------------------------------------------------------------------------------
!
! ndim (input) ... dimension of the matrix
! num_ev (input) ... number of eigenvalues to be calculated
!
! eigenvaluesV (input/output) ... real,dimension(num_ev)
! contains eigenvalues on return
! for init/=0 input values are taken as initial guess
! eigenfunctionsM (input/output) ... real,dimension(ndim,num_ev)
! contains eigenfunctions on return
!
! for init/=0 input values are taken as initial guess
! init (input) ... /=0 --> the program takes input values of eigenvaluesV and
! eigenfunctionsM as initial guess
! in this case eigenfunctions must be normalized
! =0 --> initial guess is generated by random numbers
!
! itere_max (input/output) ... maximum number of outer iterations
! contains actual number of iterations on return
! eps_itere (input/output) ... precision for outer loop
! contains actual precision on return
! info_itere (output) ... =0 ,if precision is reached within itere_max
! =1 ,oterwise
!
! itsub_max (input/output) ... maximum number of inner iterations
! contains actual number of iterations on return
! eps_itsub (input/output) ... precision for inner loop
! contains actual precision on return
! info_itsub (output) ... =0 ,if precision is reached within itsub_max
! =1 ,otherwise
!
! iwrite (input) ... =0 --> no information printed
! /=0 --> information printed
!
!-----------------------------------------------------------------------------
!
! subroutine matrix_vector(ndim,x_inV,x_outV)
! integer,intent(in) :: ndim
! real(8),dimension(ndim),intent(in) :: x_inV
! real(8),dimension(ndim),intent(out) :: x_outV
!
! calculates matrix-vector product : x_outV = Matrix * x_inV
!
!---------------------------------------------------------------
!
! subroutine precond(ndim,x_inV,x_outV)
! integer,intent(in) :: ndim
! real(8),dimension(ndim),intent(in) :: x_inV
! real(8),dimension(ndim),intent(out) :: x_outV
!
! preconditions the vector x_inV --> x_outV
!
!--------------------------------------------------------------
!
! subroutine inner_product(ndim,x_in1V,x_in2V,inn_prod)
! integer :: ndim
! real(8),dimension(ndim),intent(in) :: x_in1V,x_in2V
! real(8),intent(out) :: inn_prod
!
! calculates inner product: inn_prod = <x_in1V|x_in2V>
!
!------------------------------------------------------------------------------
SUBROUTINE
arnoldi_complex(ndim,eigenvectorL,which_eigvalC,max_eigval, &
tolerance,max_iterations,sigma,resid_lin_eq,iproblem, &
eigenvalueV,eigenvectorM,lin_equation_solver_arnoldi, &
matvec)
!-------------------------------------------------------------------------
!
! diagonalize complex matrix with Arnoldi iteration
! compute max_eigval lowest eigenvalues
!
! note: the ARPACK routines used can handle also nonhermitian matrices
! therefore eigenvalueV and eigenvectorV are complex
! If your matrix is hermitian, the same routines are recommended
! to be used, according to the ARPACK manual.
! Although the eigenvalues must be real in this case, there may be
! small imaginary parts due to round-off errors
!
! input ..
!
! ndim ... actual size of matrix
! eigenvectorL ... T calculate eigenvalues+eigenvectors, F former only
! which_eigvalC(2).. which eigenvalues are to be computed, can be
! 'LR','SR','LM','SM','LI','SI'
! 1.letter: large/small
! 2.letter: in magnitude/real/imaginary
! should be 'LM' in inverse shift mode
! max_eigval ... compute that many eigenvalues <= ndim, starting with 1
! must be <= ndim/3
! tolerance ...convergence criterion, may be set to 0.d0(-> epsilon(0d0))
! max_iterations .. maximum number of iterations
! sigma ... in shift-invert mode,
! eigenvalues near a point sigma are sought
! resid_lin_eq ... in shift-invert mode
! tolerance for solving the linear equation
! should be slightly smaller than tolerance above
! iproblem ... =1 normal mode
! =3 inverse shift mode
!
! output ..
!
! eigenvalueV(1:max_eigval) .. eigenvalues of matrix (real)
! sorted in descending order
! eigenvectorM(1:ndim,1:max_eigval) .. eigenvectors, if eigenvectorL = T
!
!---------------------------------------------
!
! USER must provide following routines as input:
!
! 1. SUBROUTINE lin_equation_solver_arnoldi(ndim,yV,xV,resid_lin_eq,ierr,&
! sigma)
! calculates yV = inv[A-sigma*I]*xV
! resid_lin_eq ... tolerance
! ierr .. output, if ierr/=0 routine failed
! sigma .. complex number, see above
!
! 2. SUBROUTINE matvec_arnoldi_complex(ndim,xV,yV)
!
! calculates yV = Matrix*xV
!
!------------------------------------------------
!
! uses arpack routines (contained in arpack_kp.f)
! and lapack routines
! included with -xlic_lib=sunperf on sun
!--------------------------------------------------------------------------------------------------------------------------
!-------------------------------------------------------------------------------------------------------------------
subroutine
eigensolver_hermitian(ndim,num_ev,eigenvaluesV,eigenfunctionsM, &
init,itere_max,itsub_max,eps_itere,eps_itsub,info_itere,info_itsub, &
iwrite,matrix_vector,precond,inner_product)
!--------------------------------------------------------------------------------------------------------------------
!
! eigensolver for complex,hermitian matrices (hermitian with respect to the
! inner product)
! calculates the num_ev LOWEST eigenvalues
! matrix is only accessed through the three user-defined subroutines
!
!------------------------------------------------------------------------------
!
! ndim (input) ... dimension of the matrix
! num_ev (input) ... number of eigenvalues to be calculated
!
! eigenvaluesV (input/output) ... real,dimension(num_ev)
! contains eigenvalues on return
! for init/=0 input values are taken as initial guess
! eigenfunctionsM (input/output) ... complex,dimension(ndim,num_ev)
! contains eigenfunctions on return
!
! for init/=0 input values are taken as initial guess
! init (input) ... /=0 --> the program takes input values of eigenvaluesV and
! eigenfunctionsM as initial guess
! in this case eigenfunctions must be normalized
! =0 --> initial guess is generated by random numbers
!
! itere_max (input/output) ... maximum number of outer iterations
! contains actual number of iterations on return
! eps_itere (input/output) ... precision for outer loop
! contains actual precision on return
! info_itere (output) ... =0 ,if precision is reached within itere_max
! =1 ,oterwise
!
! itsub_max (input/output) ... maximum number of inner iterations
! contains actual number of iterations on return
! eps_itsub (input/output) ... precision for inner loop
! contains actual precision on return
! info_itsub (output) ... =0 ,if precision is reached within itsub_max
! =1 ,otherwise
!
! iwrite (input) ... =0 --> no information printed
! /=0 --> information printed
!
!-----------------------------------------------------------------------------
!
! subroutine matrix_vector(ndim,x_inV,x_outV)
! integer,intent(in) :: ndim
! complex(8),dimension(ndim),intent(in) :: x_inV
! complex(8),dimension(ndim),intent(out) :: x_outV
!
! calculates matrix-vector product : x_outV = Matrix * x_inV
!
!---------------------------------------------------------------
!
! subroutine precond(ndim,x_inV,x_outV)
! integer,intent(in) :: ndim
! complex(8),dimension(ndim),intent(in) :: x_inV
! complex(8),dimension(ndim),intent(out) :: x_outV
!
! preconditions the vector x_inV --> x_outV
!
!--------------------------------------------------------------
!
! subroutine inner_product(ndim,x_in1V,x_in2V,inn_prod)
! integer :: ndim
! complex(8),dimension(ndim),intent(in) :: x_in1V,x_in2V
! complex(8),intent(out) :: inn_prod
!
! calculates inner product: inn_prod = <x_in1V|x_in2V>
!
!------------------------------------------------------------------------------
|
