>> help JHerm_band_Lanczos J-Hermitian band Lanczos method ----------------------------------------------------------------------------- Notes: 1) The band Lanczos method simplifies when A is J-Hermitian, i.e., J A = A^H J, and the block of left starting vectors is the complex conjugate of the J-multiple of the block of the right starting vectors, i.e., L = conj(J R). Here, A^H (= A') denotes the complex conjugate transpose of A and J = J^H is assumed to be a nonsingular Hermitian matrix. In this case, the left Lanczos vectors are complex conjugates of the J-multiples of the right Lanczos vectors: w_k = conj(J v_k) for all k. The J-Hermitian band Lanczos method exploits this fact and generates only the right Lanczos vectors. To this end, one matrix-vector product y = J x with J needs to be computed in each Lanczos iteration. Excluding these matrix-vector products with J, the remaining computational costs and the storage requirements of the function "JHerm_band_Lanczos" are only roughly half of those of the function "band_Lanczos". 2) The matrices A, J, and R are allowed to be complex in general. The matrix A is assumed to be J-Hermitian, where J is a nonsingular Hermitian matrix. 3) The "orthogonalizations" performed in the algorithm are with respect to the bilinear form (w,v) = w^H J v (= w' * J * v); note that this bilinear form is an inner product only if the Hermitian matrix J is positive definite. 4) This algorithm explicitly enforces only the minimum possible amount of "orthogonalization". Moreover, only the entries below and on the diagonal of the Lanczos matrix T are computed directly via inner products and norms of candidate vectors and previously deflated vectors. 5) This algorithm has no built-in look-ahead procedure to continue in the case of a breakdown or near-breakdown, and instead, the algorithm simply stops in this case. ----------------------------------------------------------------------------- Usage: result = JHerm_band_Lanczos(A,J,R,nmax) result = JHerm_band_Lanczos(A,J,R,nmax,sflag) result = JHerm_band_Lanczos(A,J,R,nmax,sflag,tol) result = JHerm_band_Lanczos(A,J,R,nmax,sflag,tol,n,result) where J is a Hermitian matrix and A is a J-Hermitian matrix. The matrices A and/or J in these call can be replaced by functions "afun" and/or "Jfun" such that y = afun(x,...) computes the matrix-vector product y = A x (= A * x) with A and/or y = Jfun(x,...) computes the matrix-vector product y = J x (= J * x) with J. Some examples of such usage: result = JHerm_band_Lanczos(@(x) afun(x,...),J,R,nmax) result = JHerm_band_Lanczos(A,@(x) Jfun(x,...),R,nmax,sflag,tol) result = JHerm_band_Lanczos(@(x) afun(x,...),Jfun(x,...),R, ... nmax,sflag,tol,n,result) ----------------------------------------------------------------------------- Required inputs: A = a J-Hermitian matrix or an (anonymous) function that computes matrix-vector products with A J = a nonsingular Hermitian matrix or an (anonymous) function that computes matrix-vector products with J R = matrix the m columns of which are the starting vectors nmax = maximum number of Lanczos vectors to be generated It is assumed that A is a J-Hermitian matrix, where J is a nonsingular Hermitian matrix, and that A, J, and R have the same number of rows. It is checked if A, J, and R have the same number of rows when the inputs A and J are matrices, but not when A or J is a function. When A is a matrix, it is checked if A is square. When J is a matrix, it is checked if J is Hermitian. It is not checked if A is J-Hermitian. ----------------------------------------------------------------------------- Optional inputs: sflag = an integer that controls the storage of the Lanczos vectors: sflag = 0 means that only the vectors needed for the "orthogonalizations" are stored; a nonzero value of sflag means that all vectors are stored tol = a structure that contains tolerances and parameters for the deflation procedure and the breakdown check n = a nonnegative integer; n > 0 means that a previous call to the function "JHerm_band_Lanczos" has generated n Lanczos vectors and that the current call to "JHerm_band_Lanczos" resumes the iteration at step n+1; n = 0 means that the band J-Hermitian Lanczos method is started from scratch result = the output structure of the previous call to "JHerm_band_Lanczos" if n > 0; if n = 0, this input is ignored If sflag is not provided as an input, sflag = 0 is used. If "tol" is provided as an input, it needs to contain tol.defl_tol = unscaled deflation tolerance tol.brk_tol = tolerance for breakdown check and values of the following three flags: tol.defl_flag = 0 (use unscaled deflation tolerance) 1 (use scaled deflation tolerance) 2 (use scaled deflation tolerance only for the starting block R) 3 (use scaled deflation tolerance except for the starting block R) tol.normA_flag = 1 (an estimate for the norm of A is generated within the algorithm) 2 (an estimate for the norm of A is provided as tol.normA) tol.normJ_flag = 1 (an estimate for the norm of J is generated within the algorithm) 2 (an estimate for the norm of J is provided as tol.normJ) If tol.normA_flag = 2, then an estimate for the norm of A needs to be provided as tol.normA If tol.normJ_flag = 2, then an estimate for the norm of J needs to be provided as tol.normJ If "tol" is not provided as an input, the following default values are used: tol.defl_tol = sqrt(eps) (eps = machine precision) tol.brk_tol = 100 * eps; tol.defl_flag = 1 tol.normA_flag = 1 If n is not provided as an input, n = 0 is used. If n > 0, the input "result" needs to be provided. It is assumed, but not checked, that "result" is the output structure of a previous call to the function "JHerm_band_Lanczos" (applied to the same matrices A, J, and R). In this case, the inputs "sflag" and "tol" are ignored and tol = result.tol and sflag = result.sflag are used instead. ----------------------------------------------------------------------------- On return, "result" is a structure the fields of which include the following quantities: result.n = number of Lanczos vectors that were generated result.V = matrix V the columns of which are the stored Lanczos vectors; if sflag = 0, V has at most m+1 columns and these columns are the Lanczos vectors that are needed for orthogonalization; if sflag is nonzero, V contains the first n Lanczos vectors result.Vh_defl = matrix Vh_defl the columns of which are the candidate vectors for the next mc Lanczos vectors and the m-mc deflated vectors result.D = vector D the entries of which are the products (V(:,k)' * J * V(:,k), k = 1, 2, ..., n result.T = the n x n matrix T that represents the oblique projection of the matrix A onto the subspace spanned by the columns of V, and orthogonally to the subspace spanned by the columns of J * V; A, J, V, and T are connected via the relation V' * J * A * V = (V' * J * V) * T result.rho = the matrix rho that contains the coefficients used to turn the starting vectors (in R) into the first Lanczos vectors; R, J, V, and rho are connected via the relation V' * J * R = (V' * J * V) * rho ----------------------------------------------------------------------------- This routine can be run in incremental fashion. Example 1: result = JHerm_band_Lanczos(A,J,R,nmax0,sflag,tol) n = result.n result = JHerm_band_Lanczos(A,J,R,nmax,sflag,tol,n,result) The first call to the function "JHerm_band_Lanczos" runs the J-Hermitian band Lanczos method from scratch and generates n Lanczos vectors. The second call to the function of "JHerm_band_Lanczos" resumes the iteration at step n+1. Example 2 (J-Hermitian band Lanczos method, run one step at a time): result = JHerm_band_Lanczos(A,J,R,1,sflag,tol,0,[]) for n = 1 : nmax - 1, result = JHerm_band_Lanczos(A,J,R,n+1,[],[],result.n,result) end This will run the J-Hermitian band Lanczos method for nmax steps. ----------------------------------------------------------------------------- BANDITS: a Matlab Package of Band Krylov Subspace Iterations Copyright (c) 2018-2019 Roland W. Freund See LICENSE.txt for license -----------------------------------------------------------------------------