function [X, out] = rgmres_bl_sh(A, B, shifts, opt)
% [X, out] = rgmres_bl_sh(A, B, shifts, opts)
%
% Simple implementation of the block shifted recycled GMRES algorithm which builds
% a basis for K_m(P*A,R(:,1)) + K_m(P*A,R(:,2)) + ...
% + K_m(P*A,R(:,end))
% where U has columns spanning a recycled subspace, C=A*U has orthonormal columns
% P = eye(n) - C*C, and R(:.i) = B(:,i) - A*X0(:,i) is orthogonal to range(C)
% for X0 = opt.X0 some initial approximate solutions
% and solves the GMRES minimization for all shifted systems
% (A + shifts(1)*I)X(:,1) = B(:,1)
% (A + shifts(2)*I)X(:,2) = B(:,2)
% .
% .
% .
% (A + shifts(end)*I)X(:,end) = B(:,end)
% using the block Arnoldi relation arising from the building of the basis.
%
% This algorithm is described in...
%
% INPUTS:
% A Base coefficient matrix or a matvec procedure
% B Columns of B are the right-hand sides for each shifted system
% shifts Shifts for each system, i.e., we are solving
% (A + shifts(i))x = B(:,i). Note, if there is a "base system",
% i.e., Ax=B(:,i) then there should be a shift shifts(i) = 0.
% !!!!!!!IMPORTANT!!!!!!
% In other words, there is no zero shift assumed in this code.
% opt Optional inputs
% opt.tol Optional convergence tolerance. Default: 1e-8
% opt.deptol Optional tolerance for detection of dependent vectors.
% Default: 1e-8
% opt.nrestarts Optional maximum restart cycles. Default: 10
% opt.m Optional maximum iterations per cycle. Default: 50
% opt.X0 Optional initial approxations. Size(opt.X0) =
% [n,length(shifts)]. Default: zeros(n,length(shifts))
% opt.U A recycled subspace from a previous iteration. Default: []
% opt.k The maximum dimension of recycled subspace passed from cycle to cycle
% default: 20
% opt.isOutNMV Do we output the total number of matvecs nmv?
% opt.isOutResvecFro Do we output the vector of block residual Frobenius
% norms?
% opt.isOutResvec2norm Do we output the individual 2-norm residual
% vectors for each right-hand side.
% opt.rhsCollinear Boolean indicating if we are starting with collinear
% initial residuals
% opt.rhsCoeffs If the initial residuals are collinear, we need to
% know how they are related. This is a vector of
% scalars describing each residual as a scalar multiple
% if the first residual
% opt.isGMRESCustomBackSub Do we use Matlab's in-built back
% substitution or the custom one we built
% to conform to Trilinos limitations in
% regards to back substitution. Default: false
% opt.isOutU Do we output U as out.U?
% opt.isOutNMV Do we output the total number of matvecs nmv?
% opt.recycleMethod At the end of the cycle, we need to discard part
% of the search space and pass the rest to the next
% cycle. What method is used. Default: ritz
% Possibilities are
% 'ritz': Compute Ritz vectors and take ones
% associated to smallest Ritz values.
% 'harmritz' Compute some harmonic ritz vectors
% similar to the ritz vector procedure.
% opt.stopNorm Which norm do we use to stop? We can evaluate either
% the block Frobenius norm or the max over all
% column-wise 2-norms. There are two possibilities:
% 'fro' or '2norm'. Default: 'fro'.
% opt.isComplex Does the matrix have complex entries? Default: false
% OUTPUTS:
% X Columns of X are the approximations for each shifted system
% out Optional outputs
%
% Kirk Soodhalter
% 22-Dec-2014 07:22:34
% kirk.math.soodhalter.com
n = size(B,1);
s = length(shifts);
isAProc = isa(A,'function_handle');
if ~exist('opt','var')
opt = struct(); %empty structure
end
if isfield(opt,'tol')
tol = opt.tol;
else
tol = 1e-8;
end
if isfield(opt,'deptol')
deptol = opt.deptol;
else
deptol = 1e-8;
end
if isfield(opt,'nrestarts')
maxcycles = opt.nrestarts;
else
maxcycles = 10;
end
if isfield(opt,'X0')
X0 = opt.X0;
else
X0 = zeros(n,s);
end
if isfield(opt,'m')
m = opt.m;
else
m = 50;
end
if isfield(opt,'isStab')
isStab = opt.isStab;
else
isStab = false;
end
if isfield(opt,'rhsCollinear')
rhsCollinear = opt.rhsCollinear;
else
rhsCollinear = false;
end
if isfield(opt,'stopNorm')
stopNorm = opt.stopNorm;
else
stopNorm = 'fro';
end
if isfield(opt,'isGMRESCustomBackSub')
isGMRESCustomBackSub = opt.isGMRESCustomBackSub;
else
isGMRESCustomBackSub = false;
end
if isfield(opt,'k')
k = opt.k;
else
k = 10;
end
if isfield(opt,'recycleMethod');
recycleMethod = opt.recycleMethod;
validRecycleMethods = {'ritz','harmritz'};
if ~ismember(recycleMethod,validRecycleMethods)
error('%s is not a valid value for opt.recycleMethod', recycleMethod);
end
else
recycleMethod = 'ritz';
end
if isfield(opt,'U') && ~isempty(opt.U)
U = opt.U;
else
U = [];
end
if isfield(opt,'isComplex')
isComplex = opt.isComplex;
else
isComplex = false;
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Logical array with true/false values %
% for various possible opt outputs %
% so we can test if we need to initialize %
% opt %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
isOpt = true(0);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Check to see if we %
% should output the updated U %
% and the option is turned off %
% for all subsequent recursive calls %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if isfield(opt,'isOutNMV')
isOutNMV = opt.isOutNMV;
else
isOutNMV = true;
end
isOpt = [isOpt isOutNMV];
if isfield(opt,'isOutResvecFro')
isOutResvecFro = opt.isOutResvecFro;
else
isOutResvecFro = true;
%resvecs = zeros(maxcycles*m,s);
end
isOpt = [isOpt isOutResvecFro];
if isfield(opt,'isOutResvec2norm')
isOutResvec2norm = opt.isOutResvec2norm;
else
isOutResvec2norm = true;
%resvecs = zeros(maxcycles*m,s);
end
isOpt = [isOpt isOutResvec2norm];
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Check to see if we %
% should output the updated U %
% and the option is turned off %
% for all subsequent recursive calls %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if isfield(opt,'isOutU')
isOutU = opt.isOutU;
isOpt = [isOpt isOutU];
opt.isOutU = false;
else
isOutU = false;
end
if sum(isOpt) > 0
out = struct();
%put initializations here
else
out = [];
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Construct anon. function %
% to apply shift operator %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%
D = diag(shifts);
if isAProc
Amv_sh = @(x) A(x) + x*D;
else
Amv_sh = @(x) A*x + x*D;
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Construct anon. function %
% to apply base operator %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if isAProc
Amv = @(x) A(x);
else
Amv = @(x) A*x;
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Calculate initial residual. %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
R = B - Amv_sh(X0);
Bnorm = norm(R);
%%%%%%%%%%%%%%%%%
% Main Function %
%%%%%%%%%%%%%%%%%
[X,resvecFro,resvec2norm,nmv,U] = srgmres(Amv,X0,R,shifts,tol*Bnorm,deptol,m,k,maxcycles,U,stopNorm,rhsCollinear,...
isStab,recycleMethod,isGMRESCustomBackSub,isComplex);
if isOutNMV
out.nmv = nmv;
end
if isOutU
%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Outermost function call %
% Output U for future use %
%%%%%%%%%%%%%%%%%%%%%%%%%%%
out.U = U;
end
if isOutResvecFro
out.resvecFro = resvecFro;
end
if isOutResvec2norm
out.resvec2norm = resvec2norm;
end
if isOutNMV
out.nmv = nmv;
end
end
function [X,resvecFro,resvec2norm,nmv,U] = srgmres(Amv,X,R,shifts,tol,deptol,m,k,maxcycles,U,stopNorm,rhsCollinear,...
isStab,recycleMethod,isGMRESCustomBackSub,isComplex)
[n,s] = size(X);
%%%%%%%%%%%%%%%%%%%%%%%%
% % initialize outputs %
%%%%%%%%%%%%%%%%%%%%%%%%
resvecFro = zeros(m*maxcycles, 1);
resvec2norm = zeros(m*maxcycles, s);
nmv = 0;
cycle = 0;
T = zeros(size(X));
C = zeros(size(X,1),k+1);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Calculate initial residual norms %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
resvecFro(1) = norm(R,'fro');
fprintf('||R|| = %e\t\tnmv = %d\n',resvecFro(1),1);
for i = 1:s
resvec2norm(1,i) = norm(R(:,i));
end
resvecPointer = 2;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% If we have a previous recycled %
% subspace, project residuals, if %
% not then perform a cycle of %
% GMRES to get a subspace U %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if isempty(U)
if rhsCollinear
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Render collinear residuals uncollinear %
% using one-step GMRES projection %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
beta = zeros(1,s);
for i = 1:s
beta(i) = norm(R(:,i));
end
V(:,1) = R(:,1)/beta(1);
w = Amv(V(:,1));
nmv = nmv + 1;
H(1,1) = V(:,1)'*w;
w = w - H(1,1)*V(:,1);
H(2,1) = norm(w);
V(:,2) = w/H(2,1);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% GMRES for each shifted system %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for i = 1:s
Hsig = (H + shifts(i)*[1;0]);
y = Hsig\(beta(i)*[1;0]);
T(:,i) = T(:,i) + V(:,1)*y;
R(:,i) = R(:,i) - V*(Hsig*y);
end
resvecFro(resvecPointer) = norm(R,'fro');
for i = 1:s
resvec2norm(resvecPointer,i) = norm(R(:,i));
end
fprintf('||R|| = %e (after one-step projection)\n',resvecFro(resvecPointer));
resvecPointer = resvecPointer + 1;
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Cycle of block shifted GMRES %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[T, R, W, H, iter, resvecFro_cycle, resvec2norm_cycle] = sblock_gmres_cycle(Amv, T, R, m, tol, deptol, shifts, isStab, isGMRESCustomBackSub, stopNorm);
nmv = nmv + iter;
resvecFro(resvecPointer:resvecPointer+iter-1) = resvecFro_cycle;
resvec2norm(resvecPointer:resvecPointer+iter-1,:) = resvec2norm_cycle;
resvecPointer = resvecPointer + iter;
cycle = cycle + 1;
fprintf('||R||/||R0|| = %e\t\tnmv = %d\n',resvecFro(resvecPointer-1)/resvecFro(1),nmv);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Build a recycled subspace %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if k >= iter*s %i.e. keep everything
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% form new U (the subspace to recycle) %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
U(:,1:iter*s) = W(:,1:iter*s);
keff = iter*s;
else
%%%%%%%%%%%%%%%%%%%%%%%%%%%
% downselect to k vectors %
%%%%%%%%%%%%%%%%%%%%%%%%%%%
switch recycleMethod
case 'ritz'
[P,keff] = getRitzVecsKryl(iter*s,k,H(1:(iter+1)*s,1:iter*s),isComplex);
case 'harmritz'
[P,keff] = getHarmVecsKryl(iter*s,k,H(1:(iter+1)*s,1:iter*s),isComplex);
end
U(:,1:keff) = W(:,1:iter*s) * P;
end
if iter < m
X = X + T;
U = U(:,1:keff);
resvec2norm = resvec2norm(1:resvecPointer-1,:);
resvecFro = resvecFro(1:resvecPointer-1);
return;
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Continuing to another cycle; Form orthonormalized C and adjust U accordingly so that C = A*U %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if k >= iter*s %i.e. keep everything
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% form new U (the subspace to recycle) %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[QQ,F] = qr(H,0);
C(:,1:keff) = W(:,1:(iter+1)*s) * QQ;
U(:,1:keff) = U(:,1:keff) / F;
else
%%%%%%%%%%%%%%%%%%%%%%%%%%%
% downselect to k vectors %
%%%%%%%%%%%%%%%%%%%%%%%%%%%
switch recycleMethod
case {'ritz','harmritz'}
[QQ,F] = qr(H*P,0);
C(:,1:keff) = W(:,1:(iter+1)*s) * QQ;
U(:,1:keff) = U(:,1:keff) / F;
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% After cycle of block shifted GMRES, %
% residuals are no longer all orthogonal to C, %
% so we need to do oblique projection of non-base residuals %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
CTU = C(:,1:keff)'*U(:,1:keff);
CTR = C(:,1:keff)'*R;
for i = 1:length(shifts)
rProj = ( eye(size(U,2)) + shifts(i)*CTU ) \ CTR(:,i);
T(:,i) = T(:,i) + U(:,1:keff)*rProj;
R(:,i) = R(:,i) - C(:,1:keff)*rProj - U(:,1:keff)*(shifts(i)*rProj);
end
else
% Set size of effective k
keff = size(U,2);
%grow U to correct size for allocation purposes
if keff < k + 1
U(:,keff+1:k+1) = zeros(n,k-keff+1);
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Get image of U under action of A %
% Get orth. basis for C and overwrite %
% Scale U so that AU=C %
% Note: this maps using base matrix A %
% might be better to use %
% shifted system (**) %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
C(:,1:keff) = Amv(U(:,1:keff));
[C(:,1:keff),RR] = qr(C(:,1:keff),0);
U(:,1:keff) = U(:,1:keff) / RR;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Perform oblique projection of each residual %
% so it is orthogonal to C %
% Projection orthogonal for shift=0 due to (**) %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
CTU = C(:,1:keff)'*U(:,1:keff);
CTR = C(:,1:keff)'*R;
for i = 1:length(shifts)
rProj = ( eye(keff) + shifts(i)*CTU ) \ CTR(:,i);
T(:,i) = T(:,i) + U(:,1:keff)*rProj;
R(:,i) = R(:,i) - C(:,1:keff)*rProj - U(:,1:keff)*(shifts(i)*rProj);
end
end
%%%%%%%%%%%%%%%%%%%%%%%
% Main Body of Solver %
%%%%%%%%%%%%%%%%%%%%%%%
switch stopNorm
case 'fro'
currNorm = resvecFro(resvecPointer-1);
case'2norm'
currNorm = max(resvec2norm(resvecPointer-1,:));
end
while(currNorm > tol && cycle <= maxcycles)
[T, R, W, H, B, iter, resvecFro_cycle, resvec2norm_cycle] = sblock_rgmres_cycle(Amv, T, R, shifts, m, keff, U, C, tol, deptol, isStab, stopNorm);
resvecFro(resvecPointer:resvecPointer+iter-1) = resvecFro_cycle;
resvec2norm(resvecPointer:resvecPointer+iter-1,:) = resvec2norm_cycle;
resvecPointer = resvecPointer + iter;
switch stopNorm
case 'fro'
currNorm = resvecFro(resvecPointer-1);
case'2norm'
currNorm = max(resvec2norm(resvecPointer-1,:));
end
nmv = nmv + iter;
cycle = cycle + 1;
fprintf('||R||/||R0|| = %e\t\tnmv = %d\n',resvecFro(resvecPointer-1)/resvecFro(1),nmv);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Rescale U and store the inverses of the norms %
% of its columns in the diagonals of D. %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[U(:,1:keff),D] = scaleCols(U(:,1:keff),keff);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Build a recycled subspace %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if k >= iter*s + keff
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Keep all Krylov vectors, update keff %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
U(:,keff+1:keff+iter*s) = W(:,1:iter*s);
keff_old = keff;
keff = keff_old + iter*s;
else
%%%%%%%%%%%%%%%%%%%%%%%%%%%
% downselect to k vectors %
%%%%%%%%%%%%%%%%%%%%%%%%%%%
switch recycleMethod
case 'ritz'
%%%%%%%%%%%%%%%%
% Form large H %
%%%%%%%%%%%%%%%%
H2 = modArnoldiHessenberg(D,B,H,keff,iter*s,s);
keff_old = keff;
[P,keff] = getRitzVecsAug(iter*s,k,keff_old,H2,isComplex);
case 'harmritz'
%%%%%%%%%%%%%%%%
% Form large H %
%%%%%%%%%%%%%%%%
H2 = modArnoldiHessenberg(D,B,H,keff,iter*s,s);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Calculate Harmonic Ritz vectors %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
keff_old = keff;
[P,keff] = getHarmVecsAug(iter*s,k,keff_old,s,H2,W,U(:,1:keff_old),C(:,1:keff_old),isComplex);
end
%%%%%%%%%%%%%%
% Form new U %
%%%%%%%%%%%%%%
U = U(:,1:keff_old) * P(1:keff_old,:) + W(:,1:iter*s) * P(keff_old+1:keff_old+iter*s,:);
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% If iter < m, early convergence of GMRES %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if iter < m
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Assign all (unassigned) outgoing values and return %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
X = X + T;
U = U(:,1:keff);
resvec2norm = resvec2norm(1:resvecPointer-1,:);
resvecFro = resvecFro(1:resvecPointer-1);
return;
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Continuing to another cycle; Form orthonormalized C and adjust U accordingly so that C = A*U %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if k >= iter*s + keff %i.e. keep everything
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% form new U (the subspace to recycle) %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
H2 = modArnoldiHessenberg(D,B,H,keff_old,iter*s,s);
[Q,F] = qr(H2,0);
C(:,1:keff) = C(:,1:keff_old) * Q(1:keff_old,:) + W * Q(keff_old+1:keff_old+(iter+1)*s,:);
U(:,1:keff) = U(:,1:keff) / F;
else
%%%%%%%%%%%%%%%%%%%%%%%%%%%
% downselect to k vectors %
%%%%%%%%%%%%%%%%%%%%%%%%%%%
switch recycleMethod
case {'harmritz', 'ritz'}
%%%%%%%%%%%%%%%%%%%%%%%%%
% H2, U already updated %
%%%%%%%%%%%%%%%%%%%%%%%%%
[Q,F] = qr(H2*P,0);
C(:,1:keff) = C(:,1:keff_old) * Q(1:keff_old,:) + W * Q(keff_old+1:keff_old+(iter+1)*s,:);%[C V] * Q;
U(:,1:keff) = U(:,1:keff) / F;
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% After cycle of block shifted rGMRES, %
% residuals are no longer all orthogonal to C, %
% so we need to do oblique projection of non-base residuals %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
CTU = C(:,1:keff)'*U(:,1:keff);
CTR = C(:,1:keff)'*R;
for i = 1:length(shifts)
rProj = ( eye(size(U(:,1:keff),2)) + shifts(i)*CTU ) \ CTR(:,i);
T(:,i) = T(:,i) + U(:,1:keff)*rProj;
R(:,i) = R(:,i) - C(:,1:keff)*rProj - U(:,1:keff)*(shifts(i)*rProj);
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% We have performed max number of cycles, no conv. %
% Assign all (unassigned) outgoing values and return %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
X = X + T;
U = U(:,1:keff);
resvec2norm = resvec2norm(1:resvecPointer-1,:);
resvecFro = resvecFro(1:resvecPointer-1);
end
function [X, R, W, H, B, iter, resvecFro, resvec2norm] = sblock_rgmres_cycle(Amv, X, R, shifts, m, k, U, C, tol, deptol, isStab, stopNorm)
[n,s] = size(R);
resvecFro = zeros(m,1);
resvec2norm = zeros(m,s);
H = zeros((m+1)*s,m*s);
W = zeros(n,(m+1)*s);
B = zeros(k,m*s);
CTU = C(:,1:k)'*U(:,1:k);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% If there are dependent vectors created %
% we hold them for orthogonalization purposes %
% even though they are replaced with random vectors %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
VdepOrth = zeros(n,s);
VdepOrthPointer = 1;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Get first block Arnoldi vector %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[W(:,1:s),S0] = qr(R,0);
%%%%%%%%%%%%%%%%%%%%%%%%%
% Initialize RHS for %
% least squares problem %
%%%%%%%%%%%%%%%%%%%%%%%%%
RHS = zeros((m+1)*s,s);
%keyboard
iter = 0;
%%%%%%%%%%%%%%%%%%%%
% Iterate by block %
%%%%%%%%%%%%%%%%%%%%
for pp = 1:s:m*s
iter = iter + 1;
V = Amv(W(:,pp:pp+s-1));
W(:,pp+s:pp+2*s-1) = V;
%%%%%%%%%%%%%%%%%%%%%
% Orthog against C %
% Store coeffs in B %
%%%%%%%%%%%%%%%%%%%%%
B(:,pp:pp+s-1) = C(:,1:k)'*W(:,pp+s:pp+2*s-1);
W(:,pp+s:pp+2*s-1) = W(:,pp+s:pp+2*s-1) - C(:,1:k)*B(:,pp:pp+s-1);
%%%%%%%%%%%%%%%%%
% Block Arnoldi %
% Process %
%%%%%%%%%%%%%%%%%
for j=1:s:pp
H(j:j+s-1,pp:pp+s-1) = W(:,j:j+s-1)'*W(:,pp+s:pp+2*s-1);
W(:,pp+s:pp+2*s-1) = W(:,pp+s:pp+2*s-1) - W(:,j:j+s-1)*H(j:j+s-1,pp:pp+s-1);
end
[W(:,pp+s:pp+2*s-1),H(pp+s:pp+2*s-1,pp:pp+s-1)] = qr(W(:,pp+s:pp+2*s-1),0);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% If we have vectors which were replaced due to %
% dependence of Arnoldi vectors %
% we still hold them for orthogonalization purposes %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if VdepOrthPointer > 1
HH = VdepOrth(:,1:VdepOrthPointer-1)'*W(:,pp+s:pp+2*s-1);
W(:,pp+s:pp+2*s-1) = W(:,pp+s:pp+2*s-1) - VdepOrth(:,1:VdepOrthPointer-1)*HH;
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Check if any dependent vectors to replace %
% If so, do the replacement %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
qrDiag = diag(H(pp+s:pp+2*s-1,pp:pp+s-1));
idx = find(abs(qrDiag) < deptol, 1);
if ~isempty(idx)
for j=1:s
W(:,iter*s+j) = V(:,j);
for p=1:iter*s+j-1
H(p,(iter-1)*s+j) = W(:,p)'*W(:,iter*s+j);
W(:,iter*s+j) = W(:,iter*s+j) - H(p,(iter-1)*s+j)*W(:,p);
end
H(iter*s+j,(iter-1)*s+j) = norm(W(:,iter*s+j));
if H(iter*s+j,(iter-1)*s+j) < tol
H(iter*s+j,(iter-1)*s+j) = 0;
VdepOrth(:,VdepOrthPointer) = W(:,iter*s+j);
W(:,iter*s+j) = rand(n,1);
for p=1:iter*s+j-1
HH = W(:,p)'*W(:,iter*s+j);
W(:,iter*s+j) = W(:,iter*s+j) - HH*W(:,p);
end
W(:,iter*s+j) = W(:,iter*s+j)/norm(W(:,iter*s+j));
else
W(:,iter*s+j) = W(:,iter*s+j)/H(iter*s+j,(iter-1)*s+j);
end
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Optional reorthogonalization for stabilization %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if isStab
for j=1:s:pp
HH = W(:,j:j+s-1)'*W(:,pp+s:pp+2*s-1);
W(:,pp+s:pp+2*s-1) = W(:,pp+s:pp+2*s-1) - W(:,j:j+s-1)*HH;
end
[W(:,pp+s:pp+2*s-1),~] = qr(W(:,pp+s:pp+2*s-1),0);
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Using technique from de Sturler and Kilmer %
% we construct a block matrix which is derived from %
% the minimization of the shifted systems %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
WTU = W(:,1:(iter+1)*s)'*U(:,1:k);
Uhat = U(:,1:k) - C(:,1:k)*CTU - W(:,1:(iter+1)*s)*WTU;
Vhat = orth(Uhat);
VTU = Vhat'*U(:,1:k);
kr = size(Vhat,2);
RES = zeros((iter+1)*s+k+kr,s);
RHS = eye((iter+1)*s+k+kr,s);
RHS(1:s,:) = S0;
RHS(end-kr+1:end,:) = Vhat'*R;
Y = zeros(iter*s+k,s);
HH = @(sigma) H(1:(iter+1)*s,1:iter*s) + sigma*eye((iter+1)*s,iter*s);
G = @(sigma) [HH(sigma) sigma*WTU; B(:,1:iter*s) eye(k)+sigma*CTU; zeros(kr,iter*s) sigma*VTU];
for p=1:s
if shifts(p) ~= 0
Gsig = G(shifts(p));
Y(:,p) = Gsig\RHS(:,p);
RES(:,p) = RHS(:,p) - Gsig*Y(:,p);
else
rhstemp = zeros((iter+1)*s,1);
rhstemp(1:s) = S0(:,p);
ytemp = H(1:(iter+1)*s,1:iter*s)\rhstemp;
Y(1:iter*s,p) = ytemp;
Y(iter*s+1:end,p) = - B(:,1:iter*s)*ytemp;
restemp = rhstemp - H(1:(iter+1)*s,1:iter*s)*ytemp;
RES(1:(iter+1)*s,p) = restemp;
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%
% Compute residual norm %
% Stop loop if converged %
%%%%%%%%%%%%%%%%%%%%%%%%%%
resvecFro(iter) = norm(RES,'fro'); %COMPUTE FROBENIUS NORM
for p = 1:s
resvec2norm(iter,p) = norm(RES(:,p));
end
if strcmp(stopNorm,'fro')
if resvecFro(iter) < tol
break;
end
elseif strcmp(stopNorm,'2norm')
if max(resvec2norm(iter,:)) < tol
break;
end
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Compute approximation at end of cycle %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
X = X + W(:,1:iter*s)*Y(1:iter*s,:);
X = X + U(:,1:k)*Y(iter*s+1:end,:);
R = W(:,1:(iter+1)*s)*RES(1:(iter+1)*s,:);
R = R + C(:,1:k)*RES((iter+1)*s+1:(iter+1)*s+k,:);
R = R + Vhat*RES((iter+1)*s+k+1:end,:);
resvecFro = resvecFro(1:iter);
resvec2norm = resvec2norm(1:iter,:);
end
function [X, R, W, H, iter, resvecFro, resvec2norm] = sblock_gmres_cycle(Amv, X, R, m, tol, deptol, shifts, isStab, isGMRESCustomBackSub, stopNorm)
[n,s] = size(R);
resvecFro = zeros(m,1);
resvec2norm = zeros(m,s);
H = zeros((m+1)*s,m*s);
W = zeros(n,(m+1)*s);
House = cell(s,1);
G = cell(s,1);
for p= 1:s
G{p} = zeros((m+1)*s,m*s);
end
alpha = zeros(s,1);
vrefl = cell(s,1);
nvs = zeros(s,1);
w = zeros(s,1);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% If there are dependent vectors created %
% we hold them for orthogonalization purposes %
% even though they are replaced with random vectors %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
VdepOrth = zeros(n,s);
VdepOrthPointer = 1;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Get first block Arnoldi vector %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[W(:,1:s),S0] = qr(R,0);
%%%%%%%%%%%%%%%%%%%%%%%%%
% Initialize RHS for %
% least squares problem %
%%%%%%%%%%%%%%%%%%%%%%%%%
RHS = zeros((m+1)*s,s);
RHS(1:s,1:s) = S0;
%keyboard
iter = 0;
%%%%%%%%%%%%%%%%%%%%
% Iterate by block %
%%%%%%%%%%%%%%%%%%%%
for k = 1:s:m*s
iter = iter + 1;
V = Amv(W(:,k:k+s-1));
W(:,k+s:k+2*s-1) = V;
%%%%%%%%%%%%%%%%%
% Block Arnoldi %
% Process %
%%%%%%%%%%%%%%%%%
for j=1:s:k
H(j:j+s-1,k:k+s-1) = W(:,j:j+s-1)'*W(:,k+s:k+2*s-1);
for p = 1:s
G{p}(j:j+s-1,k:k+s-1) = H(j:j+s-1,k:k+s-1);
end
W(:,k+s:k+2*s-1) = W(:,k+s:k+2*s-1) - W(:,j:j+s-1)*H(j:j+s-1,k:k+s-1);
end
[W(:,k+s:k+2*s-1),H(k+s:k+2*s-1,k:k+s-1)] = qr(W(:,k+s:k+2*s-1),0);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% If we have vectors which were replaced due to %
% dependence of Arnoldi vectors %
% we still hold them for orthogonalization purposes %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if VdepOrthPointer > 1
HH = VdepOrth(:,1:VdepOrthPointer-1)'*W(:,k+s:k+2*s-1);
W(:,k+s:k+2*s-1) = W(:,k+s:k+2*s-1) - VdepOrth(:,1:VdepOrthPointer-1)*HH;
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Check if any dependent vectors to replace %
% If so, do the replacement %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
qrDiag = diag(H(k+s:k+2*s-1,k:k+s-1));
idx = find(abs(qrDiag) < deptol, 1);
if ~isempty(idx)
for j=1:s
W(:,iter*s+j) = V(:,j);
for p=1:iter*s+j-1
H(p,(iter-1)*s+j) = W(:,p)'*W(:,iter*s+j);
W(:,iter*s+j) = W(:,iter*s+j) - H(p,(iter-1)*s+j)*W(:,p);
end
H(iter*s+j,(iter-1)*s+j) = norm(W(:,iter*s+j));
if H(iter*s+j,(iter-1)*s+j) < tol
H(iter*s+j,(iter-1)*s+j) = 0;
VdepOrth(:,VdepOrthPointer) = W(:,iter*s+j);
W(:,iter*s+j) = rand(n,1);
for p=1:iter*s+j-1
HH = W(:,p)'*W(:,iter*s+j);
W(:,iter*s+j) = W(:,iter*s+j) - HH*W(:,p);
end
W(:,iter*s+j) = W(:,iter*s+j)/norm(W(:,iter*s+j));
else
W(:,iter*s+j) = W(:,iter*s+j)/H(iter*s+j,(iter-1)*s+j);
end
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Optional reorthogonalization for stabilization %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if isStab
for j=1:s:k
HH = W(:,j:j+s-1)'*W(:,k+s:k+2*s-1);
W(:,k+s:k+2*s-1) = W(:,k+s:k+2*s-1) - W(:,j:j+s-1)*HH;
end
[W(:,k+s:k+2*s-1),~] = qr(W(:,k+s:k+2*s-1),0);
end
lastdiag = iter*s;
for p = 1:s
G{p}(k+s:k+2*s-1,k:k+s-1) = H(k+s:k+2*s-1,k:k+s-1);
G{p}(lastdiag-s+1:lastdiag,lastdiag-s+1:lastdiag) = ...
G{p}(lastdiag-s+1:lastdiag,lastdiag-s+1:lastdiag) + shifts(p)*eye(s);
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Compute newest column of R from QR %
% factorization of H %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Perform previous reflections %
% on new column %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for j = 1:iter-1
Growrange = (j-1)*s+1:(j+1)*s;
Gcolrange = k:k+s-1;
for p = 1:s
G{p}(Growrange,Gcolrange) = House{p}{j}*G{p}(Growrange,Gcolrange);
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Calculate new refletions %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for p = 1:s
House{p}{iter} = eye(2*s,2*s);
end
for j=1:s
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Calculate current reflection %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
curcol = k+j-1;
Growrange = curcol:curcol+s;
Gcolrange = curcol+1:curcol+(s-j);
for p = 1:s
alpha(p) = -exp(1i*angle(G{p}(curcol,curcol)))*norm(G{p}(curcol:curcol+s,curcol));
vrefl{p} = G{p}(curcol:curcol+s,curcol);
vrefl{p}(1) = vrefl{p}(1) - alpha(p);
nvs(p) = norm(vrefl{p})^2;
w(p) = (G{p}(curcol:curcol+s,curcol)'*vrefl{p})/(vrefl{p}'*G{p}(curcol:curcol+s,curcol));
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Apply new reflector to: %
% 1. The current column (implicitely) %
% 2. To subsequent columns of H in current block %
% 3. To House{iter} to store product of reflections for this column %
% 4. RHS, the right-hand side of the least squares problem %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for p = 1:s
G{p}(curcol:curcol,curcol) = alpha(p);%1
G{p}(curcol+1:curcol+s,curcol) = zeros(s,1);%1
G{p}(Growrange,Gcolrange) =...
G{p}(Growrange,Gcolrange) -...
((1+w(p))/nvs(p))*vrefl{p}*(vrefl{p}'*G{p}(Growrange,Gcolrange));%2
House{p}{iter}(j:j+s,:) = House{p}{iter}(j:j+s,:) - ((1+w(p))/nvs(p))*vrefl{p}*(vrefl{p}'*House{p}{iter}(j:j+s,:));%3
RHS(Growrange,p) = RHS(Growrange,p) - ((1+w(p))/nvs(p))*vrefl{p}*(vrefl{p}'*RHS(Growrange,p));%4
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%
% Compute residual norm %
% Stop loop if converged %
%%%%%%%%%%%%%%%%%%%%%%%%%%
resvecFro(iter) = norm(RHS(k+s:k+2*s-1,:),'fro'); %COMPUTE FROBENIUS NORM
for p = 1:s
resvec2norm(iter,p) = norm(RHS(k+s:k+2*s-1,p));
end
if strcmp(stopNorm,'fro')
if resvecFro(iter) < tol
break;
end
elseif strcmp(stopNorm,'2norm')
if max(resvec2norm(iter,:)) < tol
break;
end
end
end
%%%%%%%%%%%%%%%%%%%%%%%
% Solve Least Squares %
% Calculate Solution %
%%%%%%%%%%%%%%%%%%%%%%%
Y = zeros(iter*s,s);
for p = 1:s
if isGMRESCustomBackSub
Y(:,p) = UTSolve_block(G{p}, RHS(:,p), iter*s, 1);
else
Y(:,p) = G{p}(1:(iter+1)*s,1:iter*s)\RHS(1:(iter+1)*s,p);
end
end
X = X + W(:,1:k+s-1)*Y;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Compute residual by applying plane rotations/reflections %
% in reverse order to modified-lease-squares residual %
% This allows us to recover H*y without doing the matvec %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
RES = zeros((iter+1)*s,s);
for p = 1:s
RES(iter*s+1:(iter+1)*s,p) = -RHS(iter*s+1:(iter+1)*s,p);
%%%%%%%%%%%%%%%%%%%%%
% Apply reflections %
%%%%%%%%%%%%%%%%%%%%%
for j = iter:-1:1
RESrowrange = (j-1)*s+1:(j+1)*s;
RES(RESrowrange,p) = House{p}{j}'*RES(RESrowrange,p);
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%
% Now we have H*Y-E1*S0 %
% Need to get rid of -S0 %
% in first p x p block %
%%%%%%%%%%%%%%%%%%%%%%%%%%
RES(1:s,:) = RES(1:s,:) + S0;
R = R - W(:,1:(iter+1)*s)*RES;
resvecFro = resvecFro(1:iter);
resvec2norm = resvec2norm(1:iter,:);
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% *****************Basis Downselection for Recycling****************** %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [PP,keff] = getHarmVecsKryl(m,k,H, isComplex)
% Build matrix for eigenvalue problem
harmRitzMat = H(1:m,:)' \ speye(m);
harmRitzMat(1:m,1:m-1) = 0;
harmRitzMat = H(1:m,:) + H(m+1,m)^2 * harmRitzMat;
% Compute k or k+1 smallest harmonic Ritz pairs
% We get reciprocals of the computed eigenvalues
% so that getRecycleVecs() will sort them correctly
% to get vectors associated to largest values.
[harmVecs, harmVals] = eig(harmRitzMat);
harmVals = 1./diag(harmVals);
% Downselect to k or k+1 Harmonic Ritz vectors
[PP,keff] = getRecycleVecs(harmVecs,harmVals,k,m,0,isComplex);
end
function [PP,keff] = getRitzVecsKryl(m,k,H, isComplex)
% Compute k smallest harmonic Ritz pairs
[ritzVecs, ritzVals] = eig(H(1:m,1:m));
ritzVals = diag(ritzVals);
% Downselect to k or k+1 Ritz vectors
[PP,keff] = getRecycleVecs(ritzVecs,ritzVals,k,m,0,isComplex);
end
function [PP,keff] = getRitzVecsAug(m,k,keff,G, isComplex)
% Compute k smallest Ritz pairs
[ritzVecs, ritzVals] = eig(G(1:m+keff,1:m+keff));
ritzVals = diag(ritzVals);
% Downselect to k or k+1 Ritz vectors
[PP,keff] = getRecycleVecs(ritzVecs,ritzVals,k,m,keff,isComplex);
end
function [PP,keff] = getHarmVecsAug(m,k,keff,s,G,V,U,C, isComplex)
B = G' * G;
% A = | C'*U 0 |
% | V_{m+1}'*U I |
A = zeros(m+keff+s,m+keff);
A(1:keff,1:keff) = C(:,1:keff)' * U(:,1:keff);
A(keff+1:m+keff+s,1:keff) = V(:,1:m+s)' * U(:,1:keff);
A(keff+1:m+keff,keff+1:m+keff) = eye(m);
A = G' * A;
% Compute k smallest harmonic Ritz pairs
[harmVecs, harmVals] = eig(A,B);
harmVals = diag(harmVals);
% Downselect to k or k+1 Harmonic Ritz vectors
[PP,keff] = getRecycleVecs(harmVecs,harmVals,k,m,keff,isComplex);
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% *******************Matrix Modification Functions******************** %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [PP,keff] = getRecycleVecs(vecs,vals,k,m,keff,isComplex)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% This function expects that if vals %
% is ordered with ascending values, %
% the vecs corresponding to the %
% smallest vals are the ones we want.%
% Thus when building a recycling %
% procedure, design accordingly %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Construct magnitude of each value
w = abs(vals);
[~,iperm] = sort(w,1,'ascend');
% k smallest values
% Optionally store k+1 vectors to capture complex conjugate pair
PP = zeros(m+keff,k+1);
idx = 1;
while(idx <= k)
if (isreal(vals(iperm(idx)))) || isComplex
PP(:,idx) = vecs(:,iperm(idx));
idx = idx + 1;
else
PP(:,idx) = real(vecs(:,iperm(idx)));
PP(:,idx+1) = imag(vecs(:,iperm(idx)));
idx = idx + 2;
end
end
% Return number of vectors selected
keff = idx-1;
if keff == k
PP = PP(:,1:k);
end
end
function G = modArnoldiHessenberg(D,B,H,k,p,s)
G = zeros(p+k+s,p+k);
G(1:k,1:k) = D;
G(1:k,k+1:p+k) = B(1:k,1:p);
G(k+1:p+k+s,k+1:p+k) = H(1:p+s,1:p);
end
function [U,D] = scaleCols(U,k)
d = zeros(1,k);
for i = 1:k
d(i) = norm(U(:,i));
U(:,i) = U(:,i) / d(i);
end
D = diag(1 ./ d);
end
function Y = UTSolve_block(R, RHS, n, p)
% Y = UTSolve_block(R, RHS, n)
%
% Solve via backward substitution the upper triangular system
% R*Y = RHS(1:n,:) which solves a least-squares problem
%
% DATA DICTIONARY
% INPUTS:
% R n+p x n upper triangular system
% RHS n+p x p block right hand side
% n the size of the linear system we solve to get the least squares
% solution
% p Number of right-hand sides.
% OUTPUTS:
% Y The block solution to the least squares problem.
Y = zeros(n,p);
Y(n,:) = RHS(n,:) / R(n,n);
for i = n-1:-1:1,
Y(i,:) = RHS(i,:);
for j = i+1:n,
Y(i,:) = Y(i,:) - R(i,j) * Y(j,:);
end
Y(i,:) = Y(i,:) / R(i,i);
end
end