my_gmres_schur_unsym.m

求解线性系统的Krylov方法的工具箱

function [x,flag,relres,iter,resvec] = my_gmres_schur_unsym(L,U,B,C,b,restart,tol,maxit,M1,D,M2,x0)
%GMRES  Generalized Minimum Residual Method. 
% function [x,flag,relres,iter,resvec] = my_gmres_schur(L,U,B,C,b,restart,tol,maxit,M1,D,M2,x0)
% the A*x is replaced by a schur operation
% on input: C*inv(L*U)*B', b---the linear system; restart---restart GMRES after restart
% times; tol---stopping criterion; maxit---maximum iteration times;
% M1,D,M2---preconditioners; x0---initial guess
% on output: x---the approximate solution of the linear system;
% flag:
%    0 GMRES converged to the desired tolerance TOL within MAXIT
%      iterations;
%    1 GMRES iterated MAXIT times but did not converge;
%    2 preconditioner M was ill-conditioned;
%    3 GMRES stagnated (two consecutive iterates were the same);
% relres---the relative residual NORM(b-A*x)/NORM(b), if flag is 0, relres
%          <= tol;
% iter---both the outer and inner iteration numbers at which x was
%        computed: 0 <= iter(1) <= maxit and 0 <= iter(2) <= restart;
% resvec---a vector of the residual norms at each inner iteration, including
%          NORM(b-A*x0).

% Author: Liang Li, UESTC
% Date  : 2008-01-03

if (nargin < 5)
   error('Not enough input arguments.');
end

[m,n] = size(C);
% Check for all zero right hand side vector => all zero solution
n2b = norm(b);                      % Norm of rhs vector, b
if (n2b == 0)                       % if    rhs vector is all zeros
   x = zeros(m,1);                  % then  solution is all zeros
   flag = 0;                        % a valid solution has been obtained
   relres = 0;                      % the relative residual is actually 0/0
   iter = [0 0];                    % no iterations need be performed
   resvec = 0;                      % resvec(1) = norm(b-A*x) = norm(0)      
   return                           % no need to proceed
end

% Assign default values to unspecified parameters
if (nargin < 6) | isempty(restart) | (restart == n)
    restarted = 0;
else
    restarted = 1;
end
if (nargin < 7) | isempty(tol)
   tol = 1e-6;
end
if (nargin < 8 | isempty(maxit))
    if restarted
        maxit = min(ceil(n/restart),10);
    else
        maxit = min(n,10);
    end
end

if restarted
    outer = maxit;
    inner = restart;
else
    outer = 1;
    inner = maxit;
end
if ((nargin >= 9) & ~isempty(M1))
   existM1 = 1;
else
   existM1 = 0;
end

if ((nargin >= 10) & ~isempty(D))
   existD = 1;
else
   existD = 0;
end

if ((nargin >= 11) & ~isempty(M2))
   existM2 = 1;
else
   existM2 = 0;
end

if ((nargin >= 12) & ~isempty(x0))
   if ~isequal(size(x0),[m,1])
      es = sprintf(['Initial guess must be a column vector of' ...
            ' length %d to match the problem size.'],n);
      error(es);
   end
else
   x0 = zeros(m,1);
end
x = x0;

% Set up for the method
flag = 1;
xmin = x;                          % Iterate which has minimal residual so far
imin = 0;                          % "Outer" iteration at which xmin was computed
jmin = 0;                          % "Inner" iteration at which xmin was computed
tolb = tol * n2b;                  % Relative tolerance
r = b - schur_Ax_unsym(L,U,B,C,x);       % Zero-th residual
normr = norm(r);                   % Norm of residual

if (normr <= tolb)                 % Initial guess is a good enough solution
   flag = 0;
   relres = normr / n2b;
   iter = [0 0];
   resvec = normr;
   return                         % The initial guess is good enough.
end


resvec = zeros(inner*outer+1,1); % Preallocate vector for norm of residuals
resvec(1) = normr;                 % resvec(1) = norm(b-A*x0)
normrmin = normr;                  % Norm of residual from xmin
rho = 1;
stag = 0;                          % stagnation of the method

% loop over "outer" iterations (unless convergence or failure)

for iouter = 1 : outer
   V = zeros(m,inner+1);          % Arnoldi vectors
   h = zeros(inner+1,1);          % upper Hessenberg st A*V = V*H ...
   QT = zeros(inner+1,inner+1);   % orthogonal factor st QT*H = R
   R = zeros(inner,inner);        % upper triangular factor st H = Q*R
   f = zeros(inner,1);            % y = R\f => x = x0 + V*y
   W = zeros(m,inner);            % W = V*inv(R)
   iinner = 0;                    % "inner" iteration counter
   vh = r;
   if existM1
       vh = M1 \ vh;
       if isinf(norm(vh,inf))
         flag = 2;
         break
      end
   end
   if existD
       vh = D \ vh;
       if isinf(norm(vh,inf))
         flag = 2;
         break
      end
   end
   if existM2
       vh = M2 \ vh;
      if isinf(norm(vh,inf))
         flag = 2;
         break
      end
   end
   h(1) = norm(vh);
   V(:,1) = vh / h(1);
   QT(1,1) = 1;
   phibar = h(1);
   
   for iinner = 1 : inner  % inner iteration
       u = schur_Ax_unsym(L,U,B,C,V(:,iinner));
      if existM1
          u = M1 \ u;
         if isinf(norm(u,inf))
            flag = 2;
            break
         end
      end
      if existD
          u = D \ u;
          if isinf(norm(u,inf))
            flag = 2;
            break
         end
      end
      if existM2
          u = M2 \ u;
         if isinf(norm(u,inf))
            flag = 2;
            break
         end
      end
      for k = 1 : iinner  % orthogonalization
         h(k) = V(:,k)' * u;
         u = u - h(k) * V(:,k);
      end
      h(iinner+1) = norm(u);
      V(:,iinner+1) = u / h(iinner+1);
      R(1:iinner,iinner) = QT(1:iinner,1:iinner) * h(1:iinner);
      rt = R(iinner,iinner);
      
      % find cos(theta) and sin(theta) of Givens rotation
      if h(iinner+1) == 0
         c = 1.0;                      % theta = 0
         s = 0.0;
      elseif abs(h(iinner+1)) > abs(rt)
         temp = rt / h(iinner+1);
         s = 1.0 / sqrt(1.0 + abs(temp)^2); % pi/4 < theta < 3pi/4
         c = - temp * s;
      else
         temp = h(iinner+1) / rt;
         c = 1.0 / sqrt(1.0 + abs(temp)^2); % -pi/4 <= theta < 0 < theta <= pi/4
         s = - temp * c;
      end
      R(iinner,iinner) = c' * rt - s' * h(iinner+1);

      q = QT(iinner,1:iinner);
      QT(iinner,1:iinner) = c' * q;
      QT(iinner+1,1:iinner) = s * q;
      QT(iinner,iinner+1) = -s';
      QT(iinner+1,iinner+1) = c;
      f(iinner) = c' * phibar;
      phibar = s * phibar;
      
      if iinner < inner
         if f(iinner) == 0                 % stagnation of the method
            stag = 1;
         end
         W(:,iinner) = (V(:,iinner) - W(:,1:iinner-1) * R(1:iinner-1,iinner))/ R(iinner,iinner);
         % Check for stagnation of the method
         if stag == 0
            stagtest = zeros(n,1);
            ind = (x ~= 0);
            if ~(iouter==1 & iinner==1)
               stagtest(ind) = W(ind,iinner) ./ x(ind);
               stagtest(~ind & W(:,iinner) ~= 0) = Inf;
               if abs(f(iinner))*norm(stagtest,inf) < eps
                  stag = 1;
               end
            end
         end
         x = x + f(iinner) * W(:,iinner);       % form the new inner iterate
      else % j == inner
         vrf = V(:,1:iinner)*(R(1:iinner,1:iinner)\f(1:iinner));
         % Check for stagnation of the method
         if stag == 0
            stagtest = zeros(n,1);
            ind = (x0 ~= 0);
            stagtest(ind) = vrf(ind) ./ x0(ind);
            stagtest(~ind & vrf ~= 0) = Inf;
            if norm(stagtest,inf) < eps
               stag = 1;
            end
         end
         x = x0 + vrf;                % form the new outer iterate
      end 
      r = b - schur_Ax_unsym(L,U,B,C,x);
      normr = norm(r);
      resvec((iouter-1)*inner+iinner+1) = normr;
      
      if normr <= tolb               % check for convergence
         if iinner < inner
            y = R(1:iinner,1:iinner) \ f(1:iinner);
            x = x0 + V(:,1:iinner) * y;     % more accurate computation of xj
%             r = b - A * x;
%             normr = norm(r);
            resvec((iouter-1)*inner+iinner+1) = normr;
         end
         if normr <= tolb             % check using more accurate xj
            flag = 0;
            iter = [iouter iinner];
            break
         end
      end
      
      if stag == 1
         flag = 3;
         break
      end
      
      if normr < normrmin            % update minimal norm quantities
         normrmin = normr;
         xmin = x;
         imin = iouter;
         jmin = iinner;
      end
   end                              % for j = 1 : inner
   
   if flag == 1
      x0 = x;                        % save for the next outer iteration
      r = b - schur_Ax_unsym(L,U,B,C,x0);
  else
      break
   end
   
end                                % for i = 1 : maxit

% returned solution is that with minimum residual
if flag == 0
   relres = normr / n2b;
else
   x = xmin;
   iter = [imin jmin];
   relres = normrmin / n2b;
end

% truncate the zeros from resvec
if flag <= 1 | flag == 3
   resvec = resvec(1:(iouter-1)*inner+iinner+1);
else
   if iinner == 0
      resvec = resvec(1:(iouter-1)*inner+1);
   else
      resvec = resvec(1:(iouter-1)*inner+iinner);
   end
end