function[R,pvec] = hprecon(H,upperbandw);
%HPRECON Sparse Cholesky factor of H-preconditioner
%
%   [R,PVEC] = HPRECON(H,UPPERBANDW) computes the 
%   sparse Cholesky factor (transpose of a (usually) banded 
%   preconditioner of square matrix H. So R'R  
%   approximates H. If 0 < UPPERBANDW <  n then the upper bandwidth of 
%   R is UPPERBANDW. If UPPERBANDW >= n then the structure of R 
%   corresponds to a sparse Cholesky factorization of H
%   using the symmmd ordering (the ordering is returned in PVEC).

%   Copyright (c) 1990-98 by The MathWorks, Inc.
%   $Revision: 1.2 $  $Date: 1998/03/20 01:59:08 $

if nargin < 1, 
   error('hprecon requires at least 1 input parameter.'); 
end
if isempty(H), 
   error('Input parameter H cannot be empty in function hprecon.'); 
end
n = length(H); 
pvec = (1:n);
epsi = .0001*ones(n,1);
info = 1;
if nargin <2, 
   upperbandw = 0; 
end

if upperbandw >= n-1 % Try complete approximation to H
   pvec = symmmd(H);
   ddiag = diag(H);
   mind = min(ddiag);
   lambda = 0;
   if mind < 0, 
      lambda = -mind + .001; 
   end
   while info > 0
      H = H + lambda*speye(n);
      [R,info] = chol(H(pvec,pvec));
      lambda = lambda + 10;
   end
elseif (upperbandw > 0) & ( upperbandw < n-1) % Banded approximation to H
   % Banded approximation
   lambda = 0;
   ddiag = diag(H);
   mind = min(ddiag);
   if mind < 0, 
      lambda = -mind + .001; 
   end
   H = tril(triu(H,-upperbandw),upperbandw);
   while info > 0
      H = H + lambda*speye(n);
      [R,info] = chol(H);
      lambda = 4*lambda;
      if lambda <= .001, 
         lambda = 1; 
      end
   end
elseif upperbandw == 0 % diagonal approximation for H
   dnrms = sqrt(sum(H.*H))';
   d = max(sqrt(dnrms),epsi);
   R = sparse(1:n,1:n,full(d));
   pvec = (1:n);
else
    error('upperbandw must be >= 0.')
end

