1function umfpack_test (nmat)
2%UMFPACK_TEST for testing umfpack (requires ssget)
3%
4% Example:
5%   umfpack_test
6%   umfpack_test (100)  % runs the first 100 matrices
7% See also umfpack
8
9% Copyright 1995-2007 by Timothy A. Davis.
10
11index = ssget ;
12
13f = find (index.nrows == index.ncols) ;
14[ignore, i] = sort (index.nrows (f)) ;
15f = f (i) ;
16
17if (nargin < 1)
18    nmat = length (f) ;
19else
20    nmat = min (nmat, length (f)) ;
21end
22nmat = max (nmat, 1) ;
23f = f (1:nmat) ;
24
25Control = umfpack ;
26Control.prl = 0 ;
27
28clf
29
30h = waitbar (0, 'UMFPACK test') ;
31
32for k = 1:nmat
33
34    i = f (k) ;
35    waitbar (k/nmat, h, 'UMFPACK test') ;
36
37%    try
38
39        fprintf ('\nmatrix: %s %s %d\n', ...
40            index.Group{i}, index.Name{i}, index.nrows(i)) ;
41
42        Prob = ssget (i) ;
43        A = Prob.A ;
44        n = size (A,1) ;
45
46        b = rand (1,n) ;
47        c = b' ;
48
49	%-----------------------------------------------------------------------
50	% symbolic factorization
51	%-----------------------------------------------------------------------
52
53	[P1, Q1, Fr, Ch, Info] = umfpack (A, 'symbolic') ;
54	subplot (2,2,1)
55	spy (A)
56	title ('A')
57
58	subplot (2,2,2)
59	treeplot (Fr (1:end-1,2)') ;
60	title ('supercolumn etree')
61
62	%-----------------------------------------------------------------------
63	% P(R\A)Q = LU
64	%-----------------------------------------------------------------------
65
66	[L,U,P,Q,R,Info] = umfpack (A, struct ('details',1)) ;
67	err = lu_normest (P*(R\A)*Q, L, U) ;
68	fprintf ('norm est PR\\AQ-LU: %g relative: %g\n', ...
69	    err, err / norm (A,1)) ;
70
71	subplot (2,2,3)
72	spy (P*A*Q)
73	title ('PAQ') ;
74
75	cs = Info.number_of_column_singletons ; % (57) ;
76	rs = Info.number_of_row_singletons ; % (58) ;
77
78	subplot (2,2,4)
79	hold off
80        try
81            spy (L+U)
82        catch
83        end
84	hold on
85	if (cs > 0)
86	    plot ([0 cs n n 0] + .5, [0 cs cs 0 0]+.5, 'c') ;
87	end
88	if (rs > 0)
89	    plot ([0 rs rs 0 0] + cs +.5, [cs cs+rs n n cs]+.5, 'r') ;
90	end
91
92	title ('LU factors')
93	drawnow
94
95	%-----------------------------------------------------------------------
96	% PAQ = LU
97	%-----------------------------------------------------------------------
98
99	[L,U,P,Q] = umfpack (A) ;
100	err = lu_normest (P*A*Q, L, U) ;
101	fprintf ('norm est PAQ-LU:   %g relative: %g\n', ...
102	    err, err / norm (A,1)) ;
103
104	%-----------------------------------------------------------------------
105	% solve
106	%-----------------------------------------------------------------------
107
108	x1 = b/A ;
109	y1 = A\c ;
110	m1 = norm (b-x1*A) ;
111	m2 = norm (A*y1-c) ;
112
113	% factor the transpose
114	Control.irstep = 2 ;
115	[x, info] = umfpack (A', '\', c, Control) ;
116	lunz0 = info.nnz_in_L_plus_U ;
117	r = norm (A'*x-c) ;
118
119	fprintf (':: %8.2e  matlab: %8.2e %8.2e\n',  r, m1, m2) ;
120
121	% factor the original matrix and solve xA=b
122	for ir = 0:4
123            Control.irstep = ir ;
124	    [x, info] = umfpack (b, '/', A, Control) ;
125	    r = norm (b-x*A) ;
126	    if (ir == 0)
127		lunz1 = info.nnz_in_L_plus_U ;
128	    end
129	    fprintf ('%d: %8.2e %d\n', ir, r, info.iterative_refinement_steps) ;
130	end
131
132	% factor the original matrix and solve Ax=b
133	for ir = 0:4
134            Control.irstep = ir ;
135	    [x, info] = umfpack (A, '\', c, Control) ;
136	    r = norm (A*x-c) ;
137	    fprintf ('%d: %8.2e %d\n', ir, r, info.iterative_refinement_steps) ;
138	end
139
140	fprintf (...
141	    'lunz trans %12d    no trans: %12d  trans/notrans: %10.4f\n', ...
142	    lunz0, lunz1, lunz0 / lunz1) ;
143
144	%-----------------------------------------------------------------------
145	% get the determinant
146	%-----------------------------------------------------------------------
147
148	det1 = det (A) ;
149	det2 = umfpack (A, 'det') ;
150	[det3 dexp3] = umfpack (A, 'det') ;
151	err = abs (det1-det2) ;
152	err3 = abs (det1 - (det3 * 10^dexp3)) ;
153	denom = det1 ;
154	if (denom == 0)
155	    denom = 1 ;
156	end
157	err = err / denom ;
158	err3 = err3 / denom ;
159	fprintf ('det:  %20.12e + (%20.12e)i MATLAB\n', ...
160	    real(det1), imag(det1)) ;
161	fprintf ('det:  %20.12e + (%20.12e)i umfpack\n', ...
162	    real(det2), imag(det2)) ;
163	fprintf ('det: (%20.12e + (%20.12e)i) * 10^(%g) umfpack\n', ...
164	    real(det3), imag(det3), dexp3) ;
165	fprintf ('diff %g %g\n', err, err3) ;
166
167%    catch
168%        % out-of-memory is OK, other errors are not
169%        disp (lasterr) ;
170%        if (isempty (strfind (lasterr, 'Out of memory')))
171%            error (lasterr) ;                                               %#ok
172%        else
173%            fprintf ('test terminated early, but otherwise OK\n') ;
174%        end
175%    end
176
177end
178
179close (h) ;     % close the waitbar
180