my_gmres_schur.m

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

function [x,flag,relres,iter,resvec] = my_gmres_schur(L,U,B,C,b,restart,tol,maxit)

% to solve the schure complement system C*inv(L*U)*B' x = y;
% Author: Liang Li, UESTC
% Date  : 2008-01-03

[m,n] = size(B);
% 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
outer = maxit;
inner = restart;
x = zeros(m,1);
x0 = x;

% 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_umsym(L,U,B,C,x); % Zero-th residual
temp_vec = B'*x; temp_vec = L \ temp_vec; 
temp_vec = U \ temp_vec; temp_vec = C*temp_vec;
r = b - temp_vec;
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;
   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_umsym(L,U,B,C,V(:,iinner)); %(C*inv(L*U)*B')
       u = B'*V(:,iinner); u = L \ u; u = U \ u; u = C*u;
      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_umsym(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_umsym(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