my_minres.asv

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

   end
else % no preconditioner
   v = u;
end
beta1 = vold' * v;
if (beta1 <= 0)
   flag = 5;
   relres = normr / n2b;
   iter = 0;
   resvec = resvec(1);
   resveccg = resveccg(1);
   if nargout < 2
      itermsg('minres',tol,maxit,0,flag,iter,relres);
   end
   return
end
beta1 = sqrt(beta1);
snprod = beta1;
vv = v / beta1;
if isequal(atype,'matrix')
   v = A * vv;
else     
   v = iterapp(afun,atype,afcnstr,vv,varargin{:});
end
alpha = vv' * v;
v = v - (alpha/beta1) * vold;

% Local reorthogonalization
numer = vv' * v;
denom = vv' * vv;
v = v - (numer/denom) * vv;
volder = vold;
vold = v;

if existM1
   if isequal(m1type,'matrix')
      u = M1 \ vold;
   else
      u = iterapp(m1fun,m1type,m1fcnstr,vold,varargin{:});
   end
   if isinf(norm(v,inf))
      flag = 2;
      relres = normr / n2b;
      iter = 0;
      resvec = resvec(1);
      resveccg = resveccg(1);
      if nargout < 2
         itermsg('minres',tol,maxit,0,flag,iter,relres);
      end
      return
   end
else % no preconditioner
   u = vold;
end
if existM2
   if isequal(m2type,'matrix')
      v = M2 \ u;
   else
      v = iterapp(m2fun,m2type,m2fcnstr,u,varargin{:});
   end
      if isinf(norm(v,inf))
      flag = 2;
      relres = normr / n2b;
      iter = 0;
      resvec = resvec(1);
      resveccg = resveccg(1);
      if nargout < 2
         itermsg('minres',tol,maxit,0,flag,iter,relres);
      end
      return
   end
else % no preconditioner
   v = u;
end
betaold = beta1;
beta = vold' * v;
if (beta < 0)
   flag = 5;
   relres = normr / n2b;
   iter = 0;
   resvec = resvec(1);
   resveccg = resveccg(1);
   if nargout < 2
      itermsg('minres',tol,maxit,0,flag,iter,relres);
   end
   return
end
beta = sqrt(beta);
gammabar = alpha;
epsilon = 0;
deltabar = beta;
gamma = sqrt(gammabar^2 + beta^2);
mold = zeros(n,1);
m = vv / gamma;
cs = gammabar / gamma;
sn = beta / gamma;
x = x + snprod * cs * m;
snprod = snprod * sn;
% This recurrence produces CG iterates
xcg = x + snprod * (sn/cs) * m;

if existM
   if isequal(atype,'matrix')
      r = b - A * x;
   else
      r = b - iterapp(afun,atype,afcnstr,x,varargin{:});
   end
   normr = norm(r);
else
   normr = abs(snprod);
end
resvec(2,1) = normr;
if (cs == 0)
   % It's possible that this cs value is zero (CG iterate does not exist)
   normrcg = Inf;
else
   normrcg = normr / abs(cs);
end
resveccg(2,1) = normrcg;

stag = 0;                          % stagnation of the method

% loop over maxit iterations (unless convergence or failure)

for i = 1 : maxit
   
   vv = v / beta;
   if isequal(atype,'matrix')
      v = A * vv;
   else     
      v = iterapp(afun,atype,afcnstr,vv,varargin{:});
   end
   v = v - (beta / betaold) * volder;
   alpha = vv' * v;
   v = v - (alpha / beta) * vold;
   volder = vold;
   vold = v;
   if existM1
      if isequal(m1type,'matrix')
         u = M1 \ vold;
      else
         u = iterapp(m1fun,m1type,m1fcnstr,vold,varargin{:});
      end
      if isinf(norm(u,inf))
         flag = 2;
         break
      end
   else % no preconditioner
      u = vold;
   end
   if existM2
      if isequal(m2type,'matrix')
         v = M2 \ u;
      else
         v = iterapp(m2fun,m2type,m2fcnstr,u,varargin{:});
      end
      if isinf(norm(v,inf))
         flag = 2;
         break
      end
   else % no preconditioner
      v = u;
   end
   betaold = beta;
   beta = vold' * v;
   if (beta < 0)
      flag = 5;
      break
   end
   beta = sqrt(beta);
   delta = cs * deltabar + sn * alpha;
   molder = mold;
   mold = m;
   m = vv - delta * mold - epsilon * molder;
   gammabar = sn * deltabar - cs * alpha;
   epsilon = sn * beta;
   deltabar = - cs * beta;
   gamma = sqrt(gammabar^2 + beta^2);
   m = m / gamma;
   %	csold = cs;
   cs = gammabar / gamma;
   sn = beta / gamma;
   % Check for stagnation of the method
   stagtest = zeros(n,1);
   ind = (x ~= 0);
   stagtest(ind) = m(ind) ./ x(ind);
   stagtest(~ind & (m ~= 0)) = Inf;
   if (snprod*cs == 0) | (abs(snprod*cs)*norm(stagtest,inf) < eps)
      % increment the number of consecutive iterates which are the same
      stag = stag + 1;
   else
      % this iterate is not the same as the previous one
      stag = 0;
   end
   x = x + snprod * cs * m;
   snprod = snprod * sn;
   % This recurrence produces CG iterates
   xcg = x + snprod * (sn/cs) * m;
   
   if existM
      if isequal(atype,'matrix')
         r = b - A * x;
      else
         r = b - iterapp(afun,atype,afcnstr,x,varargin{:});
      end
      normr = norm(r);
   else
      normr = abs(snprod);
   end
   resvec(i+2,1) = normr;
   if (cs == 0)
      % It's possible that this cs value is zero (CG iterate does not exist)
      normrcg = Inf;
   else
      normrcg = normr / abs(cs);
   end
   resveccg(i+2,1) = normrcg;
   
   if (normr <= tolb) % check for convergence
      % double check residual norm is less than tolerance
      if isequal(atype,'matrix')
         r = b - A * x;
      else
         r = b - iterapp(afun,atype,afcnstr,x,varargin{:});
      end
      normr = norm(r);
      if (normr <= tolb)
         flag = 0;
         iter = i;
         break
      end
   end
   if (normrcg <= tolb)
      % Conjugate Gradients solution
      xcg = x + snprod * (sn/cs) * m;
      % double check residual norm is less than tolerance
      if isequal(atype,'matrix')
         r = b - A * xcg;                 % CG residual
      else
         r = b - iterapp(afun,atype,afcnstr,xcg,varargin{:});
      end
      normrcg = norm(r);
      if (normrcg <= tolb)
         x = xcg;
         flag = 0;
         iter = i;
         break
      end
   end
   if (stag >= 2)                  % 3 consecutive iterates are the same
      flag = 3;
      break
   end
   if (normr < normrmin)           % update minimal norm quantities
      normrmin = normr;
      xmin = x;
      imin = i;
   end   
end                                % for i = 1 : maxit

% returned solution is first with minimal residual
if (flag == 0)
   relres = normr / n2b;
else
   x = xmin;
   iter = imin;
   relres = normrmin / n2b;
end

% truncate the zeros from resvec
if ((flag <= 1) | (flag == 3))
   resvec = resvec(1:i+1);
   resveccg = resveccg(1:i+1);
else
   resvec = resvec(1:i);
   resveccg = resveccg(1:i);
end

% only display a message if the output flag is not used
if (nargout < 2)
   itermsg('minres',tol,maxit,i,flag,iter,relres);
end