1function test22(nmat)
2%TEST22 test pos.def and indef. matrices
3% Example:
4% test22(nmat)
5%
6% if nmat <= 0, just test problematic matrices
7% See also cholmod_test
8
9% Copyright 2007, Timothy A. Davis, http://www.suitesparse.com
10
11fprintf ('=================================================================\n');
12fprintf ('test22: test pos.def and indef. matrices\n') ;
13
14index = ssget ;
15
16[ignore f] = sort (index.nrows) ;
17
18if (nargin > 0)
19    problematic = (nmat <= 0) ;
20    if (problematic)
21	% Matrices for which MATLAB and CHOLMOD differ, for which the error
22	% is high, or other issues arose during debugging
23	fprintf ('testing matrices for which MATLAB and CHOLMOD differ\n') ;
24	f = [
25	    186 % HB/jgl011
26	    109 % HB/curtis54
27	    793 % Qaplib/lp_nug05
28	    607 % LPnetlib/lp_brandy
29	    707 % LPnetlib/lp_bore3d
30	    231 % HB/plskz362
31	    794 % Qaplib/lp_nug06
32	    673 % LPnetlib/lp_scorpion
33	   1156 % Sandia/oscil_dcop_45 (also 1157:1168)
34	    795 % Qaplib/lp_nug07
35	    796 % Qaplib/lp_nug08
36	    260 % HB/well1033
37	    261 % HB/well1850
38	    230 % HB/plsk1919
39	    649 % LPnetlib/lp_pds_02
40	    660 % LPnetlib/lp_qap12
41	    609 % LPnetlib/lp_cre_a
42	    619 % LPnetlib/lp_dfl001
43	    661 % LPnetlib/lp_qap15
44	    650 % LPnetlib/lp_pds_06
45	    379 % Cote/vibrobox
46	    638 % LPnetlib/lp_ken_11
47	    799 % Qaplib/lp_nug20
48	    ]' ;
49    else
50	nmat = max (0,nmat) ;
51	nmat = min (nmat, length (f)) ;
52	f = f (1:nmat) ;
53    end
54end
55
56skip = [
57    811, ...		% Simon/appu, which is a random matrix
58    937:939, ...	% large ND/ problems
59    1157:1168 ...	% duplicates
60    799			% rather large: Qaplib/lp_nug20
61    ] ;
62
63tlimit = 0.1 ;
64fprintf ('test22: chol and chol2 are repeated so each take >= %g sec\n',tlimit);
65
66klimit = 1 ;
67
68% warmup, for more accurate timing
69[R,p] = chol (sparse (1)) ;						    %#ok
70[R,p] = chol2 (sparse (1)) ;						    %#ok
71clear R p
72
73for i = f
74
75    if (any (i == skip))
76	continue ;
77    end
78
79    Problem = ssget (i) ;
80    A = Problem.A ;
81    [m n] = size (A) ;
82    fprintf ('\n================== %4d: Problem: %s  m: %d n: %d nnz: %d', ...
83	i, Problem.name, m, n, nnz (A)) ;
84    clear Problem ;
85
86    try	% create a symmetric version of the matrix
87	if (m == n)
88	    if (nnz (A-A') > 0)
89		A = A+A' ;
90	    end
91	else
92	    A = A*A' ;
93	end
94    catch
95	fprintf ('skip\n') ;
96	continue
97    end
98
99    fprintf (' %d\n', nnz (A)) ;
100
101    p = amd2 (A) ;
102    A = A (p,p) ;
103    anorm = norm (A,1) ;
104
105    % Run each code for at least 'tlimit' seconds
106
107    % MATLAB
108    k = 0 ;
109    t1 = 0 ;
110    while (t1 < tlimit & k < klimit)					    %#ok
111	tic ;
112	[R1,p1] = chol (A) ;
113	t = toc ;
114	t1 = t1 + t ;
115	k = k + 1 ;
116    end
117    t1 = t1 / k ;
118
119    % CHOLMOD
120    k = 0 ;
121    t2 = 0 ;
122    while (t2 < tlimit & k < klimit)					    %#ok
123	tic ;
124	[R2,p2] = chol2 (A) ;
125	t = toc ;
126	t2 = t2 + t ;
127	k = k + 1 ;
128    end
129    t2 = t2 / k ;
130
131    if (klimit == 1)
132	rmin = full (min (abs (diag (R2)))) ;
133	rmax = full (max (abs (diag (R2)))) ;
134	if (p2 ~= 0 | isnan (rmin) | isnan (rmax) | rmax == 0)		    %#ok
135	    rcond = 0 ;
136	else
137	    rcond = rmin / rmax ;
138	end
139	fprintf ('rcond: %30.20e\n', rcond) ;
140    end
141
142    if (p1 == 1)
143	% MATLAB does not follow its own definitions.  If p is 1, then R is
144	% supposed to be 0-by-n, not 0-by-0.  CHOLMOD fixes this bug.
145	% Here, A is m-by-m
146	R1 = sparse (0,m) ;
147    end
148
149    kerr = 0 ;
150    if (p1 ~= p2)
151	% MATLAB and CHOLMOD don't agree.  See if both are correct,
152	% because differences in roundoff errors can make one go
153	% a little farther than the other.
154
155	% if p1 is zero, it means MATLAB was fully successful
156	k1 = p1 ;
157	if (k1 == 0)
158	    k1 = n ;
159	end
160
161	% if p2 is zero, it means CHOLMOD was fully successful
162	k2 = p2 ;
163	if (k2 == 0)
164	    k2 = n ;
165	end
166
167	if (k1 > k2)
168	    % MATLAB went further than CHOLMOD. This is OK if MATLAB found
169	    % a small entry where CHOLMOD stopped.
170	    k = k2 ;
171	    kerr = R1 (k,k) ;
172	    % now reduce R1 in size, to compare with R2
173	    R1 = R1 (1:k2-1,:) ;
174	else
175	    % CHOLMOD went further than MATLAB. This is OK if CHOLMOD found
176	    % a small entry where MATLAB stopped.
177	    k = k1 ;
178	    kerr = R2 (k,k) ;
179	    % now reduce R2 in size, to compare with R1
180	    R2 = R2 (1:k1-1,:) ;
181	end
182    end
183
184    err = norm (R1-R2,1) / max (anorm,1) ;
185    fprintf ('p: %6d %6d MATLAB: %10.4f CHOLMOD: %10.4f speedup %6.2f err:', ...
186	p1, p2, t1, t2, t1/t2) ;
187
188    if (err == 0)
189	fprintf ('      0') ;
190    else
191	fprintf (' %6.0e', err) ;
192    end
193
194    if (kerr == 0)
195	fprintf ('      0\n') ;
196    else
197	fprintf (' %6.0e\n', kerr) ;
198    end
199
200%    if (err > 1e-6)
201%	error ('!') ;
202%    end
203
204% pause
205    clear R1 R2 p1 p2 p A
206
207end
208
209fprintf ('test22: all tests passed\n') ;
210