pdsco.m

来自「％ Atomizer Main Directory, Version .802 」· M 代码 · 共 1,246 行 · 第 1/4 页

     'The least-squares solution is good enough, given atol     '
     'The estimate of cond(Abar) has exceeded conlim            '
     'Ax - b is small enough for this machine                   '
     'The least-squares solution is good enough for this machine'
     'Cond(Abar) seems to be too large for this machine         '
     'The iteration limit has been reached                      '];

%wantvar= nargout >= 6;
%if wantvar, var = zeros(n,1); end
  
itn    = 0;		istop  = 0;		nstop  = 0;
ctol   = 0;		if conlim > 0, ctol = 1/conlim; end;
anorm  = 0;		acond  = 0;
dampsq = damp^2;	ddnorm = 0;		res2   = 0;
xnorm  = 0;		xxnorm = 0;		z      = 0;
cs2    = -1;		sn2    = 0;

% Set up the first vectors u and v for the bidiagonalization.
% These satisfy  beta*u = b,  alfa*v = A'u.

u      = b(1:m);	x    = zeros(n,1);
alfa   = 0;		beta = norm( u );
if beta > 0
   u = (1/beta) * u;	v = feval( aprodname, 2, m, n, u, iw, rw );
   alfa = norm( v );
end
if alfa > 0
   v = (1/alfa) * v;    w = v;
end

arnorm = alfa * beta;   if arnorm == 0, disp(msg(1,:)); return, end

rhobar = alfa;		phibar = beta;		bnorm  = beta;
rnorm  = beta;
head1  = '   Itn      x(1)      Function';
head2  = ' Compatible   LS      Norm A   Cond A';

if show
   disp(' ')
   disp([head1 head2])
   test1  = 1;		test2  = alfa / beta;
   str1   = sprintf( '%6g %12.5e %10.3e',   itn, x(1), rnorm );
   str2   = sprintf( '  %8.1e %8.1e',       test1, test2 );
   disp([str1 str2])
end

%------------------------------------------------------------------
%     Main iteration loop.
%------------------------------------------------------------------
while itn < itnlim
      itn = itn + 1;
%     Perform the next step of the bidiagonalization to obtain the
%     next  beta, u, alfa, v.  These satisfy the relations
%                beta*u  =  a*v   -  alfa*u,
%                alfa*v  =  A'*u  -  beta*v.

      u    = feval( aprodname, 1, m, n, v, iw, rw )  -  alfa*u;
      beta = norm( u );
      if beta > 0
         u     = (1/beta) * u;
         anorm = norm([anorm alfa beta damp]);
         v     = feval( aprodname, 2, m, n, u, iw, rw )  -  beta*v;
         alfa  = norm( v );
         if alfa > 0,  v = (1/alfa) * v; end
      end

%     Use a plane rotation to eliminate the damping parameter.
%     This alters the diagonal (rhobar) of the lower-bidiagonal matrix.

      rhobar1 = norm([rhobar damp]);
      cs1     = rhobar / rhobar1;
      sn1     = damp   / rhobar1;
      psi     = sn1 * phibar;
      phibar  = cs1 * phibar;

%     Use a plane rotation to eliminate the subdiagonal element (beta)
%     of the lower-bidiagonal matrix, giving an upper-bidiagonal matrix.

      rho     = norm([rhobar1 beta]);
      cs      =   rhobar1/ rho;
      sn      =   beta   / rho;
      theta   =   sn * alfa;
      rhobar  = - cs * alfa;
      phi     =   cs * phibar;
      phibar  =   sn * phibar;
      tau     =   sn * phi;

%     Update x and w.

      t1      =   phi  /rho;
      t2      = - theta/rho;
      dk      =   (1/rho)*w;

      x       = x      +  t1*w;
      w       = v      +  t2*w;
      ddnorm  = ddnorm +  norm(dk)^2;
%     if wantvar, var = var  +  dk.*dk; end

%     Use a plane rotation on the right to eliminate the
%     super-diagonal element (theta) of the upper-bidiagonal matrix.
%     Then use the result to estimate  norm(x).

      delta   =   sn2 * rho;
      gambar  = - cs2 * rho;
      rhs     =   phi  -  delta * z;
      zbar    =   rhs / gambar;
      xnorm   =   sqrt(xxnorm + zbar^2);
      gamma   =   norm([gambar theta]);
      cs2     =   gambar / gamma;
      sn2     =   theta  / gamma;
      z       =   rhs    / gamma;
      xxnorm  =   xxnorm  +  z^2;

%     Test for convergence.
%     First, estimate the condition of the matrix  Abar,
%     and the norms of  rbar  and  Abar'rbar.

      acond   =   anorm * sqrt( ddnorm );
      res1    =   phibar^2;
      res2    =   res2  +  psi^2;
      rnorm   =   sqrt( res1 + res2 );
      arnorm  =   alfa * abs( tau );

%     Now use these norms to estimate certain other quantities,
%     some of which will be small near a solution.

      test1   =   rnorm / bnorm;
      test2   =   arnorm/( anorm * rnorm );
      test3   =       1 / acond;
      t1      =   test1 / (1    +  anorm * xnorm / bnorm);
      rtol    =   btol  +  atol *  anorm * xnorm / bnorm;

%     The following tests guard against extremely small values of
%     atol, btol  or  ctol.  (The user may have set any or all of
%     the parameters  atol, btol, conlim  to 0.)
%     The effect is equivalent to the normal tests using
%     atol = eps,  btol = eps,  conlim = 1/eps.

      if itn >= itnlim,   istop = 7; end
      if 1 + test3  <= 1, istop = 6; end
      if 1 + test2  <= 1, istop = 5; end
      if 1 + t1     <= 1, istop = 4; end

%     Allow for tolerances set by the user.

      if  test3 <= ctol,  istop = 3; end
      if  test2 <= atol,  istop = 2; end
      if  test1 <= rtol,  istop = 1; end

%-----------------------------------------------------------------------
%     SPECIAL TEST THAT DEPENDS ON pdsco.m.
%     Aname in pdsco  is  iw in lsqr.
%     dy              is  x
%     Other stuff     is in info.
%     We allow for diagonal preconditioning in pdsDDD2.
%-----------------------------------------------------------------------
      if istop > 0
	 r3new     = arnorm;
	 r3ratio   = r3new / info.r3norm;
         atolold   = atol;
	 atolnew   = atol;
	 
	 if atol > info.atolmin
	    if     r3ratio <= 0.1     % dy seems good
	       % Relax
	    elseif r3ratio <= 0.5     % Accept dy but make next one more accurate.
	       atolnew = atolnew * 0.1;
	    else                      % Recompute dy more accurately
	       fprintf('\n                                ')
	       fprintf('                                ')
	       fprintf(' %5.1f%7g%7.3f', log10(atolold), itn, r3ratio)
	       atol    = atol * 0.1;
	       atolnew = atol;
	       istop   = 0;
	    end
	 end

	 outfo.atolold = atolold;
	 outfo.atolnew = atolnew;
      	 outfo.r3ratio = r3ratio;
      end
      
%-----------------------------------------------------------------------
%     See if it is time to print something.
%-----------------------------------------------------------------------
      prnt = 0;
      if n     <= 40       , prnt = 1; end
      if itn   <= 10       , prnt = 1; end
      if itn   >= itnlim-10, prnt = 1; end
      if rem(itn,10) == 0  , prnt = 1; end
      if test3 <=  2*ctol  , prnt = 1; end
      if test2 <= 10*atol  , prnt = 1; end
      if test1 <= 10*rtol  , prnt = 1; end
      if istop ~=  0       , prnt = 1; end

      if prnt == 1
         if show
            str1 = sprintf( '%6g %12.5e %10.3e',   itn, x(1), rnorm );
            str2 = sprintf( '  %8.1e %8.1e',       test1, test2 );
            str3 = sprintf( ' %8.1e %8.1e',        anorm, acond );
            disp([str1 str2 str3])
         end
      end
      if istop > 0, break, end
end

%     End of iteration loop.
%     Print the stopping condition.

if show
   disp(' ')
   disp('LSQR finished')
   disp(msg(istop+1,:))
   disp(' ')
   str1 = sprintf( 'istop  =%8g   itn    =%8g',      istop, itn    );
   str2 = sprintf( 'anorm  =%8.1e   acond  =%8.1e',  anorm, acond  );
   str3 = sprintf( 'rnorm  =%8.1e   arnorm =%8.1e',  rnorm, arnorm );
   str4 = sprintf( 'bnorm  =%8.1e   xnorm  =%8.1e',  bnorm, xnorm  );
   disp([str1 '   ' str2])
   disp([str3 '   ' str4])
   disp(' ')
end

%-----------------------------------------------------------------------
% End private function pdsxxxlsqr
%-----------------------------------------------------------------------


function y = pdsxxxlsqrmat( mode, mlsqr, nlsqr, x, Aname, rw )
%
% pdsxxxlsqrmat is required by pdsco.m (when it calls pdsxxxlsqr.m).
% It forms Mx or M'x for some operator M that depends on LSproblem below.
%
% mlsqr, nlsqr  are the dimensions of the LS problem that lsqr is solving.
%
% Aname is the name of the user's (Ax,A'y) routine
% or the default routine 'pdsxxxmat'.
%
% rw contains parameters [explicit LSproblem LSmethod LSdamp]
% from pdsco.m to say which least-squares subproblem is being solved.
%
% global pdsDDD1 pdsDDD2 provides various diagonal matrices
% for each value of LSproblem.

%-----------------------------------------------------------------------
% 17 Mar 1998: First version to go with pdsco.m and lsqr.m.
% 01 Apr 1998: global pdsDDD1 pdsDDD2 now used for efficiency.
% 11 Feb 2000: Added diagonal preconditioning for LSQR, LSproblem = 1.
% 14 Dec 2000: Added diagonal preconditioning for LSQR, LSproblem = 12.
% 30 Jan 2001: Added diagonal preconditioning for LSQR, LSproblem = 21.
% 12 Feb 2001: Included in pdsco.m as private function.
%              Specialized to allow only LSproblem = 1.
%-----------------------------------------------------------------------

global pdsDDD1 pdsDDD2;

explicit  = rw(1);
LSproblem = rw(2);
LSmethod  = rw(3);
LSdamp    = rw(4);
precon    = rw(7);

if LSproblem == 1
    % The operator M is [ DA';  LSdamp*I ].
    m = nlsqr;
    n = mlsqr - m;
    if mode == 1
       if precon, x = pdsDDD2 .* x; end
       t = feval( Aname, 2, m, n, x );   % Ask 'aprod' to form t = A'x.
       y = [ (pdsDDD1.*t); (LSdamp*x) ];
    else
       t = pdsDDD1.*x(1:n);
       y = feval( Aname, 1, m, n, t );   % Ask 'aprod' to form y = At.
       y = y  +  LSdamp * x(n+1:mlsqr);
       if precon, y = pdsDDD2 .* y; end
    end

else
    disp('Error in pdsxxx: Only LSproblem = 1 is allowed')
    break
end

%-----------------------------------------------------------------------
% End private function pdsxxxlsqrmat
%-----------------------------------------------------------------------


function y = pdsxxxmat( mode, m, n, x )
%        y = pdsxxxmat( mode, m, n, x )
%    computes y = Ax (mode=1) or A'x (mode=2)
%    for some matrix A, for use with pdsco.m.

%-----------------------------------------------------------------------
% 04 Apr 1998: Default A*x and A'*y function for pdsco.m.
%              It assumes A is the global matrix pdsAAA created by pdsco.m
%              from the user's input parameter A.
%-----------------------------------------------------------------------

global pdsAAA;

if mode == 1
   y = pdsAAA*x;
else
   y = pdsAAA'*x;
end

%-----------------------------------------------------------------------
% End private function pdsxxxmat
%-----------------------------------------------------------------------