1function [A, Z, title, key, mtype] = RBfix (filename)
2%RBFIX read a possibly corrupted matrix from a R/B file
3% (assembled format only).  Usage:
4%
5% [A Z title key mtype] = RBfix (filename)
6%
7% The Rutherford/Boeing format stores a sparse matrix in a file in compressed-
8% column form, using 3 arrays: Ap, Ai, and Ax.  The row indices of entries in
9% A(:,j) are in Ai(p1:p2) and the corresponding numerical values are Ax(p1:p2),
10% where p1 = Ap(j) and p2 = Ap(j+1)-1.  The row indices ought to be sorted, and
11% no duplicates should appear, but this function ignores that requirement.
12% Duplicate entries are summed if they exist, and A is returned with sorted
13% columns.  Symmetric matrices are stored with just their lower triangular
14% parts in the file.  Normally, it is an error if entries are present in the
15% upper triangular part of a matrix that is declared in the file to be
16% symmetric.  This function simply ignores those entries.
17%
18% If CHOLMOD is installed, this function is faster and uses less memory.
19%
20% Example:
21%
22%   load west0479
23%   RBwrite ('mywest', west0479, [ ], 'My west0479 file', 'west0479') ;
24%   [A Z title key mtype] = RBfix ('mywest') ;
25%   isequal (A, west0479)
26%   title, key, mtype
27%
28% See also mread, RBread, RBwrite, RBreade, sparse2.
29
30% Optionally uses the CHOLMOD sparse2 mexFunction.
31
32% Copyright 2009, Timothy A. Davis
33
34%-------------------------------------------------------------------------------
35% read in the raw contents of the Rutherford/Boeing file
36%-------------------------------------------------------------------------------
37
38[mtype Ap Ai Ax title key nrow] = RBraw (filename) ;
39mtype = lower (mtype) ;
40
41%-------------------------------------------------------------------------------
42% determine dimension, number of entries, and convert numerical entries
43%-------------------------------------------------------------------------------
44
45% number of columns
46ncol = length (Ap) - 1 ;
47
48% number of entries
49nz = length (Ai) ;
50
51% check column pointers
52if (any (Ap ~= sort (Ap)) | (Ap (1) ~= 1) | (Ap (ncol+1) - 1 ~= nz))	    %#ok
53    error ('invalid column pointers') ;
54end
55
56% check row indices
57if ((double (max (Ai)) > nrow) | double (min (Ai)) < 1)			    %#ok
58    error ('invalid row indices') ;
59end
60
61% Ax can be empty, for a p*a matrix
62if (~isempty (Ax))
63    if (mtype (1) == 'c')
64	% Ax is real, with real/imaginary parts interleaved
65	if (2 * nz ~= length (Ax))
66	    error ('invalid matrix') ;
67	end
68	Ax = Ax (1:2:end) + (1i * Ax (2:2:end)) ;
69    elseif (mtype (1) == 'i')
70	Ax = double (Ax) ;
71    end
72    % numerical values must be of the right size
73    if (nz ~= length (Ax))
74	error ('invalid matrix') ;
75    end
76end
77
78%-------------------------------------------------------------------------------
79% create the triplet form
80%-------------------------------------------------------------------------------
81
82% construct column indices
83Aj = zeros (nz,1) ;
84for j = 1:ncol
85    p1 = Ap (j) ;
86    p2 = Ap (j+1) - 1 ;
87    Aj (p1:p2) = j ;
88end
89
90%-------------------------------------------------------------------------------
91% create the sparse matrix form
92%-------------------------------------------------------------------------------
93
94if (exist ('sparse2') == 3)						    %#ok
95    % Use sparse2 in CHOLMOD.  It's faster, allows integer Ai and Aj, and
96    % returns the Z matrix as the 2nd output argument.
97    if (isempty (Ax))
98	Ax = 1 ;
99    end
100    % numerical matrix
101    [A Z] = sparse2 (Ai, Aj, Ax, nrow, ncol) ;
102else
103    % stick with MATLAB, without CHOLMOD.  This is slower and takes more memory.
104    Ai = double (Ai) ;
105    Aj = double (Aj) ;
106    if (isempty (Ax))
107	% pattern-only matrix
108	A = spones (sparse (Ai, Aj, 1, nrow, ncol)) ;
109	Z = sparse (nrow, ncol) ;
110    else
111	% numerical matrix
112	A = sparse (Ai, Aj, Ax, nrow, ncol) ;
113	% determine the pattern of explicit zero entries
114	S = spones (sparse (Ai, Aj, 1, nrow, ncol)) ;
115	Z = S - spones (A) ;
116    end
117end
118
119% check for entries in upper part
120if (any (mtype (2) == 'shz') & nnz (triu (A,1) > 0))			    %#ok
121    fprintf ('entries in upper triangular part of %s matrix ignored\n', mtype);
122end
123
124% add the upper triangular part
125if (mtype (2) == 's')
126    A = A + tril (A,-1).' ;
127    Z = Z + tril (Z,-1)' ;
128elseif (mtype (2) == 'h')
129    A = A + tril (A,-1)' ;
130    Z = Z + tril (Z,-1)' ;
131elseif (mtype (2) == 'z')
132    A = A - tril (A,-1).' ;
133    Z = Z + tril (Z,-1)' ;
134end
135
136