my_gmres_schur_block.asv

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

function [x,y,flag,relres,iter,resvec] = my_gmres_schur_block(A,f,g,restart,tol,maxit,L,U,B,C,x0,y0,alf,tol_inner,ptype)
%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: B*inv(L*U)*B'+C, 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); % m 是(2,2)块的大小,n is the size of block (1,1)
% Check for all zero right hand side vector => all zero solution
b = [f;g];
n2b = norm(b);                      % Norm of rhs vector, b
if (n2b == 0)                       % if    rhs vector is all zeros
   x = zeros(n,1);                  % then  solution is all zeros
   y = zeros(m,1);
   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
x = x0;y = y0;
% 
% 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;y];                          % 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
r1 = f - A*x - B'*y;
r2 = g - C*x;
% r = b - schur_Ax(L,U,B,C,x);                     % Zero-th residual
normr = norm([r1;r2]);                   % 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(n+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(n+m,inner);            % W = V*inv(R)
   iinner = 0;                    % "inner" iteration counter
   vh1 = r1; vh2 = r2;
   % precondition
   if ptype == 2 % 对角预条件
       vh1 = L \ vh1;
       vh1 = U \ vh1;
       vh2 = my_gmres_schur(L,U,B,C,vh2,10,tol_inner,500);
   elseif ptype == 3 % 三角预条件,alpha 在 block (2,2)
       vh2 = my_gmres_schur(L,U,B,C,vh2/alf,10,tol_inner,500);
       vh1 = vh1 - B'*vh2;
       vh1 = L \ vh1;
       vh1 = U \ vh1;
   elseif ptype == 4 % 三角预条件,alpha 在 block (1,2)
       vh2 = my_gmres_schur(L,U,B,C,vh2/alf,10,tol_inner,500);
       vh1 = vh1 - alf*B'*vh2;
       vh1 = L \ vh1;
       vh1 = U \ vh1;
   end
   h(1) = norm([vh1;vh2);
   V(:,1) = vh / h(1);
   QT(1,1) = 1;
   phibar = h(1);
   
   for iinner = 1 : inner  % inner iteration
       u = schur_Ax(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(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(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