[f,out] = funm_kryl (A,b,param)For the two implementations of this restarted Krylov method we refer to [1] and [5].
Download funm_kryl (about 20 kB, version 20090409)To use funm_kryl extract the zip-archive on your computer and change to the directory funm_kryl in Matlab. Run the example files to see how it works.
param = param_init ()The fields of param can now be modified as required. Note that not every parameter configuration makes sense, e.g., you cannot use param.H_full==false if param.function is not a structure of type parfrac. Therefore you should rerun
param = param_init (param)which will check whether your settings are valid. If necessary, param_init returns a modified valid param structure.
Field of param | Type | Values | Examples |
function | function handle | a handle to a Matlab function which takes a square matrix argument H and returns f(H) | param.function = @sinm; param.function = @myfunm; |
---|---|---|---|
struct | a structure representing a partial fraction, see parfrac() for details | param.function = parfrac('expBA',16); | |
V_full | boolean | false: (default) do not store the full Krylov basis | --- |
true: keep the full Krylov basis in memory and return it in out.V_full | --- | ||
H_full | boolean | false: use recursive scheme to compute updates of Arnoldi approximations (only available if param.function is a structure as returned by parfrac) | --- |
true: (default) evaluate matrix function for full Hessenberg matrix H_full which grows in size with each restart cycle, out.H_full will be returned | --- | ||
restart_length | integer >= 1 | if the dimension of the Krylov space exceeds param.restart_length the Krylov approximation will be updated and a new Krylov space is generated | param.restart_length = 20; |
max_restarts | integer >= 1 | number of cycles of the restarting algorithm (if the number is 1 no restart is done) | param.max_restarts = 10; |
hermitian | boolean | false: (default) A is not Hermitian (the Arnoldi process is used) | --- |
true: A is Hermitian (the Lanczos process is used) | --- | ||
reorth_number | 0 (default), 1 or 2 | number reorthonormalizations in the Gram-Schmidt orthonormalization (has no effect in the Hermitian case) | --- |
thick | [] | (default) No thick-restarts are performed. | --- |
function handle | Thick-restarts are performed. This requires a function that reorderes the Schur form U*T*U' of a Hessenberg matrix H with eigenvalues D such that wanted eigenvalues occur in the upper left block of the reordered Schur matrix. The wanted eigenvalues are returned in out.thick_replaced and the unwanted eigenvalues are returned in out.thick_interpol. | param.thick = @thick; | |
harmonic_target | inf (default) or number | Target for harmonic Ritz extraction after each cycle. The value inf corresponds to the standard Ritz extraction. | --- |
exact | vector | exact solution of f(A)b; if provided, out.err is a vector of the absolute error in the 2-norm of the Krylov approximation for each restart cycle | param.exact = expm(A)*b; param.exact = []; |
bound | boolean | compute lower and upper error bounds (in non-Hermitan case these may be useful as error indicators) | --- |
stopping_accuracy | positive number | the restart iteration stops if the upper error bound or the absolute error of the Krylov approximation falls below this treshold (default accuracy is 1e-12) | --- |
min_decay | positive number | the restart iteration stops if the linear convergence rate of the upper error bound or the absolute error of the Krylov approximation is greater than this number (default linear decay rate is 0.95) | --- |
waitbar | boolean | show waitbar | --- |
[f,out] = funm_kryl (A,b,param) report(param.out) % show simple reportThe restarted Krylov approxiation is returned in f.
Field | Type | Values |
out.update | vector | norm of the updates of the Krylov approximations for each restart cycle |
out.err | vector | absolute error of the Krylov approximations for each restart cycle (only available if param.exact is non-empty) |
out.time | vector | computation time (i.e., Arnoldi/Lanczos decomposition + update of the approximation) for each restart cycle (sum(out.time) is the overall computation time) |
out.H_full | matrix | block-bidiagonal upper Hessenberg matrix, only returned if param.H_full == 1 |
out.V_full | matrix | the basis vectors of the full Krylov space, blockwise orthonormal, only returned if param.V_full == 1 |
out.bound_1, out.bound_2 | vector | lower and upper error estimators, only if param.bound == 1 (requires eigenvalue calculation of a matrix of size (param.restart_length times param.restart_length) in each restart cycle). For A Hermitian these bounds are based on Gauss-Radau and Gauss-Lobatto quadrature, respectively, and would be correct for A's smallest and largest eigenvalue to be known. Since we use the smallest and largest Ritz values instead, these bounds improve for later cycles. These bounds are based on an extension of an idea given in [8] and [7]. |
out.stop_condition | string | string that tells which stop condition was satisfied at termination of the algorithm |
out.thick_replaced, out.thick_interpol | cell array | if thick-restarts are used these arrays contain the wanted resp. unwanted eigenvalues which occured during the restart method. |
r = parfrac (type,degree,param1) param.parfrac = r x = 3; eval_parfrac(r,x) % test evaluation at point xThe structure r is a representation of a partial fraction having degree poles. The following types are implemented, further details may be found in [9]:
Type | Available degrees | Meaning | param1 |
expCF | 2,4,...,16(?) | quasi-best rational approximation to the exponential function on (-inf,0] via Caratheodory-Fejer | n.a. |
expWP | 1,2,... | rational approximation to the exponential function on (-inf,0] with poles on the Weidemann parabolic contour | n.a. |
expTW | 2,4,... | rational approximation to the exponential function on (-inf,0] with poles on the Talbot-Weidemann cotangent contour | n.a. |
expWH | 1,2,... | rational approximation to the exponential function on (-inf,0] with poles on the Weideman hyperbolic contour | n.a. |
expBA | 2,3,...,18 | best rational approximation to the exponential function on (-inf,0] (Cody-Meinardus-Varga approximation), table lookup with coefficients given in [3], [4] | n.a. |
cfud | 1,2,...,8 | best rational approximation to function defined by param1 on the unit disk (CF-method) | handle to function analytic on unit disk, e.g., param1=@exp, param1=@cos |
sqrtZ | --- | sqrtZ has been removed. Use sqrtinvZ to approximate x/sqrt(x). | --- |
sqrtinvZ | 1,2,...,50 | Zolotarovs relative best rational approximation to 1/sqrt(x) on the interval [1,param1] | right-end of approximation interval, e.g., param1=1e3 |
signZ | 1,2,...,50 | Zolotarovs best rational approximation to sign(x) on the set [-param1,param1]\[-1,1] | end-point of approximation interval, e.g., param1=1e3 |
inv | 1 | the function 1/z (solution of linear systems of equations) | n.a. |