function [A,xy,B,hbtype,rhstype] = hbo(file)
%HBO	Read and process a Harwell-Boeing sparse matrix file.
%	A = HBO('matfile') gets the sparse matrix from the specified .mat file.
%
%	[A,xy,B,hbtype,rhstype] = HBO('matfile') also gets:
%       xy: node coordinates (if present)
%       B:  right-hand sides, starting guesses, or exact solutions (if present)
%       hbtype: matrix type, which will be one of the following three-letter
%	        codes: CSA, PRA, PSA, PSE, PUA, RRA, RSA, RUA, RZA, UNK
%       rhstype: a description of B, for which see the user's guide.
%
%	In addition to reading the file, HBO assembles any unassembled matrices,
%	symmetrizes any symmetric matrices, and implicitizes any explicit zeros.
%
%	HBO .mat files are created from Harwell-Boeing data files by the
%	stand-alone Fortran process, HBO2MAT.  These .mat files contain:
%
%	    For assembled, type xxA, matrices:
%	        A       - the sparse matrix.
%	        hbtitle - The first 72 characters of the first "card".
%	        hbname  - The matrix name, same as file name without .mat.
%	        hbtype  - One of those three letter codes.
%	        hbfill  - If present, the value inserted in pattern matrices.
%	        hbzero  - If present, the value inserted for explicit zeros.
%	        rhstype - If present, the right hand side type.
%               xy      - If present, a matrix whose rows are node coordinates.
%               B       - If present, a matrix whose columns are right-hand
%                         sides, etc.
%
%	    For unassembled, xxE, matrices:
%	        hbtitle - The first 72 characters of the first "card".
%	        hbname  - The matrix name, same as file name without .mat.
%	        hbtype  - Only type PSE exists in current collection.
%	        varind and elptr - The indices specifying the locations
%	            of the constituent elements.
%
%  Say "help harwell" for more information about the collection.
%  See also HBOLIST, HBOFIND.

%	Cleve Moler, The MathWorks, 4/2/94.

load(file)

if ~exist('hbtype')
   hbtype = 'UNK';
   A = A;

% PSE  - Pattern symmetric unassembled
elseif strcmp(hbtype,'PSE')
   n1 = length(elptr);
   elptr(n1) = [];
   J = varind;
   I = zeros(size(J));
   I(elptr) = ones(size(elptr));
   I = cumsum(I);
   A = sparse(I,J,1);
   A = A'*A;

% RSA  - Real symmetric
elseif strcmp(hbtype,'RSA')
   A = A + A' - diag(diag(A));

% RZA  - Real skew symmetric
elseif strcmp(hbtype,'RZA')
   A = A - A';

% RUA  - Real unsymmetric
elseif strcmp(hbtype,'RUA')
   A = A;

% RRA  - Real rectangular
elseif strcmp(hbtype,'RRA')
   A = A;

% CSA  - Complex symmetric
elseif strcmp(hbtype,'CSA')
   A = A + A' - diag(diag(A));

% PSA  - Pattern symmetric
elseif strcmp(hbtype,'PSA')
   A = A + A' - diag(diag(A));

% PUA  - Pattern unsymmetric
elseif strcmp(hbtype,'PUA')
   A = A;

% PRA  - Pattern rectangular
elseif strcmp(hbtype,'PRA')
   A = A;

else
   error(['Harwell-Boeing type ' hbtype ' unexpected.'])
end

% Remove any explict zeros

if exist('hbzero')
   k = find(A == hbzero);
   A(k) = sparse(length(k),1);
end

% Get any right-hand side vectors or coordinates.

if exist('xyz')
    xy = xyz;
end;

if ~exist('xy')
    xy = [];
end;

if ~exist('B')
    B = [];
    rhstype = 'NON';
end;

if ~exist('rhstype')
    rhstype = 'UNK';
end;