pdsco.m

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

r       = rlin  -  LSdamp2 * y;
t       = tlin  +  (gamma2 * x  +  grad);
v       = mu    -  Xz;

% Initialize other things.

infinity  = 1.0e+20;
itn       = 0;
converged = 0;
atol      = atol1;
atol2     = max( atol2, atolmin );

%  Iteration log.

f       = [r; t; v];
normf   = norm(f);            % Merit function for linesearch.
normr   = norm(r);
normt   = norm(t);
normv   = norm(v);
normx   = norm(x);
normz   = norm(z);
objtrue = obj;
stepx   = 0;
stepz   = 0;
nf      = 0;
itncg   = 0;
nfail   = 0;
center  = max(Xz) / min(Xz);

head1   = '\n\nItn   mu stepx stepz  Pinf  Dinf';
head2   = ' mu-xz   Objective    nf center';
if direct,
    head3 = '';
else
    head3 =['  atol ' solver];
end
fprintf( [ head1 head2 head3 ] )

fprintf('\n%3g                 ', itn         )
fprintf('%6.1f%6.1f' , log10(normr), log10(normt))
fprintf('%6.1f%15.7e', log10(normv), objtrue     )
fprintf('   %7.1f'   , center                    )


if kminor
    fprintf('\n\nStart of first minor itn...\n')
    keyboard
end

%-----------------------------------------------------------------------
% Main loop.
%-----------------------------------------------------------------------

while ~converged
    itn = itn + 1;
%  x1  = x;  y1 = y;  z1 = z;  % Save for debugging at end of iteration.
%  t1  = t;  r1 = r;  v1 = v;

%-----------------------------------------------------------------------
%  Define a damped Newton iteration for solving
%     ( r, t, v ) = 0,
%  keeping  x, z > 0.  We eliminate dz
%  to obtain the system
%
%     [-H   A'] [ dx ] = [ w ],    d2I = delta^2 I,   w = t - v./x;
%     [ A  d2I] [ dy ] = [ r ]
%
%  which is equivalent to
%
%     [-dI  DA'] [ s  ] = [  Dw   ],    dI = delta I,
%     [ AD  dI ] [ dy ] = [r/delta]     D  = H^{-1/2),  dx = delta D s,
%
%  and also equivalent to the least-squares problem
%
%     min || [ DA']dy  -  [  Dw   ] ||.                                (*)
%         || [ dI ]       [r/delta] ||
%
% 17 Mar 1998: We solve the latter as the damped least-squares problem
%
%     min || [ DA']dybar  -  [D wbar] ||,  wbar = w - (A'r)/delta^2,   (**)
%         || [ dI ]          [  0   ] ||     dy = dybar + r/delta^2,
%
%    to allow lsqr to work with the smaller operator DA'.
%
% 30 Mar 1998: LSproblem = 1 or 2 selects (*) or (**) respectively.
%              (*)  seems safer if delta is small, but
%              (**) is more efficient (and safe) if delta = 1, say.
%
% 31 Mar 1998: LSproblem = 11 or 12 selects alternative LS problem
%
%     min || [ AD ]s  -  [    r    ] ||,     dx = Ds,                 (***)
%         || [ dI ]      [-delta Dw] ||      dy = (r - Adx) / delta^2
%
% and associated damped LS problem:
%
%     min || [ AD ]sbar  -  [ rbar ] ||,      s = sbar - Dw,         (****)
%         || [ dI ]         [  0   ] ||    rbar = r + AD^2 w.
%
% 06 Apr 1998: LSproblem = 21 selects equivalent LS problem
%
%     min || [    A     ]dx  -  [    r    ] ||,                     (*****)
%         || [delta Dinv]       [-delta Dw] ||     dy = (r - Adx) / delta^2
%-----------------------------------------------------------------------
    Q       = hess  +  gamma2;
    xinv    = 1 ./ x;
    H       = Q  +  xinv.*z;
    w       = t  -  xinv.*v;
    Hinv    = 1 ./ H;
    D       = sqrt(Hinv);
    DDD1    = D;

    rw      = [explicit  LSproblem  LSmethod  LSdamp m n 0];
              % Passed to LSQR, SYMMLQ.

    if LSproblem == 1
       %-----------------------------------------------------------------
       %  Solve (*) for dy.
       %-----------------------------------------------------------------
       Dw      = D.*w;

       if direct
          DD   = sparse( 1:n, 1:n, D, n, n );
          AD   = AAA*DD;

          if useChol
             d2I  = LSdamp2 * em;
             d2I  = sparse( 1:m, 1:m, d2I, m, m );
             ADDA = AD * AD'  +  d2I;

             if itn==1,  P = symmmd(ADDA);  end  % Do ordering only once.

             [R,indef] = chol(ADDA(P,P));
             if indef
                fprintf('\n\n   Chol says AD^2A'' is not positive definite')
                break
             end

%           R    = full(R);  %%% Matlab 4.0 must have a bug ---
%                            %%% Needed for R \ (R' \ rhs) below.

             rhs  = r  +  AAA * (Hinv.*w);
%           dy   = ADDA \ rhs;
             dy   = R \ (R' \ rhs(P));   dy(P) = dy;
          else % useQR
             dI   = LSdamp * em;
             rhs  = [ Dw; r/LSdamp ];
             dy   = [ AD'; diag(dI) ] \ rhs;
          end

          Ady  = AAA'*dy;
          dx   = Hinv .* (Ady - w);
          Adx  = AAA*dx;
       else
          precon  = explicit;
       %  precon  = false;    % TEST
          if precon           % Construct diagonal preconditioner for LSQR
             rw(7)= precon;
             DD   = sparse( 1:n, 1:n, D, n, n );
             AD   = AAA*DD;
             AD2  = AD.^2;
             wD   = sum( (AD2') )';   % Sum of squares of each row of AD
             wD   = sqrt( wD + LSdamp2 );
             DDD2 = 1 ./ wD;
          end
          rhs    = [ Dw; r/LSdamp ];
          damp   = 0;
          [ dy, istop, itncg ] = ...
              lsqr( nb, m, 'pdsxxx', Aname, rw, rhs, damp, ...
                    atol, btol, conlim, itnlim, show );

          if precon, dy = DDD2 .* dy; end
          Ady = feval( Aname, 2, m, n, dy );
          dx  = Hinv .* (Ady - w);
          Adx = feval( Aname, 1, m, n, dx );
       end

    elseif LSproblem == 2
       %-----------------------------------------------------------------
       % Solve (**) for dybar (damped least-squares problem).
       %-----------------------------------------------------------------
       yshift = (1/LSdamp2)*r;

       if direct
          wshift = AAA'*yshift;
          wbar   = w - wshift;
          rhs    = D.*wbar;
          DD     = sparse( 1:n, 1:n, D, n, n );
          dI     = LSdamp * em;
          dybar  = [ DD*(AAA'); diag(dI) ] \ [rhs; zeros(m,1)];

          dy     = dybar + yshift;
          Ady    = AAA'*dy;
          dx     = Hinv .* (Ady - w);
          Adx    = AAA*dx;
       else
          wshift = feval( Aname, 2, m, n, yshift );
          wbar   = w - wshift;
          rhs    = D.*wbar;
          damp   = LSdamp;
          [ dybar, istop, itncg ] = ...
             lsqr( n, m, 'pdsxxx', Aname, rw, rhs, damp, ...
                   atol, btol, conlim, itnlim, show );

          dy     = dybar + yshift;
          Ady = feval( Aname, 2, m, n, dy );
          dx  = Hinv .* (Ady - w);
          Adx = feval( Aname, 1, m, n, dx );
       end

    elseif LSproblem == 11
       %-----------------------------------------------------------------
       % Solve (***) for s (then dx = Ds).
       %-----------------------------------------------------------------
       Dw      = D.*w;
       rhs     = [ r; (-LSdamp * Dw) ];

       if direct
          DD     = sparse( 1:n, 1:n, D, n, n );
          dI     = LSdamp * en;
          s      = [ (AAA*DD); diag(dI) ] \ rhs;

          dx     = D.*s;
          Adx    = AAA*dx;
          dy     = (r - Adx) / LSdamp2;
          Ady    = AAA'*dy;
       else
          damp   = 0;
          [ s, istop, itncg ] = ...
             lsqr( nb, n, 'pdsxxx', Aname, rw, rhs, damp, ...
                   atol, btol, conlim, itnlim, show );

          dx     = D.*s;
          Adx    = feval( Aname, 1, m, n, dx );
          dy     = (r - Adx) / LSdamp2;
          Ady    = feval( Aname, 2, m, n, dy );
       end

    elseif LSproblem == 12
       %-----------------------------------------------------------------
       % Solve (****) for sbar.
       %-----------------------------------------------------------------
       sshift = D.*w;
       wbar   = D.*sshift;

       if direct
          rshift = AAA*wbar;
          rhs    = r + rshift;
          DD     = sparse( 1:n, 1:n, D, n, n );
          dI     = LSdamp * en;
          sbar   = [ AAA*DD; diag(dI) ] \ [rhs; zeros(n,1)];

          dx     = D.*(sbar - sshift);
          Adx    = AAA*dx;
          dy     = (r - Adx) / LSdamp2;
          Ady    = AAA'*dy;
       else
          rshift = feval( Aname, 1, m, n, wbar );
          rhs    = r + rshift;
          damp   = LSdamp;
          [ sbar, istop, itncg ] = ...
             lsqr( m, n, 'pdsxxx', Aname, rw, rhs, damp, ...
                   atol, btol, conlim, itnlim, show );

          dx     = D.*(sbar - sshift);
          Adx    = feval( Aname, 1, m, n, dx );
          dy     = (r - Adx) / LSdamp2;
          Ady    = feval( Aname, 2, m, n, dy );
       end

    elseif LSproblem == 21
       %-----------------------------------------------------------------
       % Solve (*****) for dx.
       %-----------------------------------------------------------------
       Dw      = D.*w;
       rhs     = [ r; (-LSdamp * Dw) ];

       if direct
          dx     = [ AAA; diag(LSdamp./D) ] \ rhs;
          Adx    = AAA*dx;
          dy     = (r - Adx) / LSdamp2;
          Ady    = AAA'*dy;
       else
          damp   = 0;
          DDD1   = LSdamp ./ D;
          [ dx, istop, itncg ] = ...
             lsqr( nb, n, 'pdsxxx', Aname, rw, rhs, damp, ...
                   atol, btol, conlim, itnlim, show );

          Adx    = feval( Aname, 1, m, n, dx );
          dy     = (r - Adx) / LSdamp2;
          Ady    = feval( Aname, 2, m, n, dy );
       end

    elseif LSproblem == 31
       %-----------------------------------------------------------------
       % Solve 3x3 system symmetrized by Z^{1/2}.
       %-----------------------------------------------------------------
       Zroot = sqrt(z);
       K3    = [ diag(x)      diag(Zroot)   zeros(n,m)
                 diag(Zroot)  diag(- Q)         AAA'
                 zeros(m,n)       AAA     diag(LSdamp2*em) ];
       rhs   = [(v ./ Zroot); t; r];
       dwxy  = K3 \ rhs;
       dz    = Zroot .* dwxy(1:n);
       dx    = dwxy(  n+1:n+n  );
       dy    = dwxy(n+n+1:n+n+m);
       Adx   = AAA*dx;
       Ady   = AAA'*dy;