function[RPCMTX,ppvec] = aprecon(DM,DG,A,upperbandw);
%APRECON Banded preconditioner function for least-squares problems.
%
%   [RPCMTX,PPVEC] = APRECON(DM,DG,A,UPPERBW) produces the sparse
%   nonsingular upper triangular matrix RPCMTX such that
%   RPCMTX'*RPCMTX is a preconditioner for the
%   matrix M = DM*(A'*A)*DM + DG, where DM is a positive
%   diagonal matrix, DG is a non-negative diagonal matrix,
%   and A is sparse rectangular matrix with more rows than columns.
%   PPVEC is the associated permutation (row) vector and UPPERBW
%   specifies the upperbandwidth of RPCMTX.

%   Copyright (c) 1990-98 by The MathWorks, Inc.
%   $Revision: 1.2 $  $Date: 1998/07/21 22:50:21 $

% Initialization
if nargin < 3, 
   error('aprecon requires at least 3 input parameters.'), 
end
[m,n] = size(A);
if m < n
   error('A matrix must have at least as many rows as columns.')
end
n = length(DM);
RPCMTX = speye(n); 
ppvec = (1:n);
if isempty(A)
   RPCMTX = sparse([]);
   ppvec=[];
   return
end
epsi = sqrt(eps);
if nargin < 4
   upperbandw = 0; 
elseif isempty(upperbandw)
   upperbandw = 0;
end

% Form matrix M
TM = A*DM;

% Determine factor of preconditioner.
if upperbandw == 0 % Diag. preconditioner based on column 2-norms
   M = TM'*TM + DG;
   % M cannot be a row vector so no problem with sum.
   dnrms = sqrt(sum(M.*M))';
   d = max(sqrt(dnrms),epsi);
   RPCMTX = sparse(1:n,1:n,full(d));
   ppvec = (1:n);
   % But what if RPCMTX is singular?
elseif upperbandw >= n-1 % Attempt sparse QR
   dgg = sqrt(diag(DG)); 
   dgg = full(dgg);
   DDG = sparse(1:n,1:n,dgg);
   TM = [TM;DDG];
   p = colmmd(TM);
   RPCMTX = qr(TM(:,p)); 
   RPCMTX = RPCMTX(1:n,1:n);
   ppvec = p;
   
   %    Modify for singularity?
   mdiag = min(abs(diag(RPCMTX)));
   lambda = 1;
   while mdiag < sqrt(eps);
      TM = [A*DM; DDG + lambda*speye(n)];
      lambda = 4*lambda;
      p = colmmd(TM);
      RPCMTX = qr(TM(:,p));
      RPCMTX = RPCMTX(1:n,1:n);
      ppvec = p;
      mdiag = min(abs(diag(RPCMTX)));
   end
   
elseif (upperbandw > 0) & (upperbandw < n-1) % Band approximation.
   M = TM'*TM + DG;
   p = (1:n);
   M = tril(triu(M(p,p),-upperbandw),upperbandw);
   RPCMTX = sparse(n,n);
   [RPCMTX,info] = chol(M);
   lambda = 1;
   while info > 0
      M = M + lambda*speye(n);
      RPCMTX = sparse(n,n);
      [RPCMTX,info] = chol(M);
      lambda = 4*lambda;
   end
   ppvec = p;
else
   error('upperbandw must be >= 0.');
end


