>> help Jsym_band_Lanczos J-symmetric band Lanczos method ----------------------------------------------------------------------------- Notes: 1) The band Lanczos method simplifies when A is J-symmetric, i.e., J A = A^T J, and the block of left starting vectors is the J-multiple of the block of right starting vectors, i.e., L = J R. Here, A^T (= A.') denotes the transpose of A and J = J^T is assumed to be a nonsingular symmetric matrix. In this case, the left and right Lanczos vectors are J-multiples of each other: w_k = J v_k for all k. The J-symmetric 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 the matrix-vector products with J, the remaining computational costs and the storage requirements of the function "Jsym_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-symmetric, where J is a nonsingular symmetric matrix. 3) The "orthogonalizations" performed in the algorithm are with respect to the bilinear form (w,v) = w^T J v (= w.' * J * v), where w^T (= w.') denotes the transpose of w; note that in the complex case, this bilinear form is not an inner product and in the real case, this bilinear form is an inner product only if the symmetric 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 = Jsym_band_Lanczos(A,J,R,nmax) result = Jsym_band_Lanczos(A,J,R,nmax,sflag) result = Jsym_band_Lanczos(A,J,R,nmax,sflag,tol) result = Jsym_band_Lanczos(A,J,R,nmax,sflag,tol,n,result) where J is a symmetric matrix and A is a J-symmetric 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 = Jsym_band_Lanczos(@(x) afun(x,...),J,R,nmax) result = Jsym_band_Lanczos(A,@(x) Jfun(x,...),R,nmax,sflag,tol) result = Jsym_band_Lanczos(@(x) afun(x,...),Jfun(x,...),R, ... nmax,sflag,tol,n,result) ----------------------------------------------------------------------------- Required inputs: A = a J-symmetric matrix or an (anonymous) function that computes matrix-vector products with A J = a nonsingular symmetric 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-symmetric matrix, where J is a nonsingular symmetric 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 symmetric. It is not checked if A is J-symmetric. ----------------------------------------------------------------------------- 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 "Jsym_band_Lanczos" has generated n Lanczos vectors and that the current call to "Jsym_band_Lanczos" resumes the iteration at step n+1; n = 0 means that the band J-symmetric Lanczos method is started from scratch result = the output structure of the previous call to "Jsym_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 "Jsym_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 = Jsym_band_Lanczos(A,J,R,nmax0,sflag,tol) n = result.n result = Jsym_band_Lanczos(A,J,R,nmax,sflag,tol,n,result) The first call to the function "Jsym_band_Lanczos" runs the J-symmetric band Lanczos method from scratch and generates n Lanczos vectors. The second call to the function of "Jsym_band_Lanczos" resumes the iteration at step n+1. Example 2 (J-symmetric band Lanczos method, run one step at a time): result = Jsym_band_Lanczos(A,J,R,1,sflag,tol,0,[]) for n = 1 : nmax - 1, result = Jsym_band_Lanczos(A,J,R,n+1,[],[],result.n,result) end This will run the J-symmetric 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 -----------------------------------------------------------------------------