pdsxxx.m

来自「％ Atomizer Main Directory, Version .802 」· M 代码 · 共 105 行

function y = pdsxxx( mode, mlsqr, nlsqr, x, Aname, rw )
%
% pdsxxx.m is required by pdsco.m (when it calls lsqr.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 'pdsmat'.
%
% rw contains parameters [explicit LSproblem LSmethod LSdamp]
% from pdsco.m to say which least-squares subproblem is being solved.
%
% global AAA DDD1 DDD2 provides A and 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 AAA DDD1 DDD2 now used for efficiency.
% 11 Feb 2000: Added diagonal preconditioning for LSQR, LSproblem = 1.
%-----------------------------------------------------------------------

global AAA DDD1 DDD2;

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 = DDD2 .* x; end
       t = feval( Aname, 2, m, n, x );   % Ask 'aprod' to form t = A'x.
       y = [ (DDD1.*t); (LSdamp*x) ];
    else
       t = DDD1.*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 = DDD2 .* y; end
    end

elseif LSproblem == 2
    % The operator M is just [ DA' ].
    m = nlsqr;
    n = mlsqr;
    if mode == 1
       t = feval( Aname, 2, m, n, x );   % Ask 'aprod' to form t = A'x.
       y = DDD1.*t;
    else
       t = DDD1.*x;
       y = feval( Aname, 1, m, n, t );   % Ask 'aprod' to form y = At.
    end

elseif LSproblem == 11
    % The operator M is [ AD;  LSdamp*I ].
    m   = mlsqr - nlsqr;
    n   = nlsqr;

    if mode == 1
       t = DDD1.*x;
       y = feval( Aname, 1, m, n, t );   % y = At.
       y = [ y; (LSdamp*x) ];
    else
       t = x(1:m);
       y = feval( Aname, 2, m, n, t );   % y = A't.
       y = DDD1.*y + LSdamp*x(m+1:mlsqr);
    end

elseif LSproblem == 12
    % The operator M is just AD.
    m   = mlsqr;
    n   = nlsqr;

    if mode == 1
       t = DDD1.*x;
       y = feval( Aname, 1, m, n, t );   % y = At.
    else
       t = x(1:m);
       y = feval( Aname, 2, m, n, t );   % y = A't.
       y = DDD1.*y;
    end

elseif LSproblem == 21
    % The operator M is [ A; DDD1 ], where DDD1 = LSdamp*Dinv.
    m   = mlsqr - nlsqr;
    n   = nlsqr;

    if mode == 1
       y = feval( Aname, 1, m, n, x );   % y = Ax.
       y = [ y; (DDD1.*x) ];
    else
       t = x(1:m);
       y = feval( Aname, 2, m, n, t );   % y = A't.
       y = y + DDD1.*x(m+1:mlsqr);
    end
end

%-----------------------------------------------------------------------
% End of pdsxxx.m
%-----------------------------------------------------------------------