1function test0 (nmat,f)
2%TEST0 test most CHOLMOD functions
3% Example:
4%   test0(nmat)
5% See also cholmod_test
6
7% Copyright 2007, Timothy A. Davis, http://www.suitesparse.com
8
9fprintf ('=================================================================\n');
10fprintf ('test0: test most CHOLMOD functions\n') ;
11
12% This test requires ssget, the MATLAB interface to the UF sparse matrix
13% collection.  You can obtain ssget from http://www.suitesparse.com.
14
15try
16    index = ssget ;
17catch
18    error ('Test aborted.  UF sparse matrix collection not available.\n') ;
19end
20
21if (nargin < 2)
22    f = find (index.posdef) ;
23    [ignore i] = sort (index.nrows (f)) ;
24    f = f (i) ;
25end
26
27rand ('state', 0) ;
28randn ('state', 0) ;
29
30doplots = 0 ;
31
32if (doplots)
33    clf
34end
35
36% skip = [937:939 1202:1211] ;
37skip = 937:939 ;
38if (nargin > 0)
39    nmat = max (0,nmat) ;
40    nmat = min (nmat, length (f)) ;
41    f = f (1:nmat) ;
42end
43
44% f= 229
45
46fprintf ('test matrices sorted by dimension:\n') ;
47for i = f
48    if (any (i == skip))
49	continue
50    end
51    fprintf ('%4d: %-20s %-20s %12d %d\n', i,  ...
52	index.Group {i}, index.Name {i}, index.nrows (i), index.posdef (i)) ;
53end
54
55% pause
56
57for i = f
58
59    if (any (i == skip))
60	continue
61    end
62
63    % try
64
65	Problem = ssget (i) ;
66	A = Problem.A ;
67	fprintf ('\n================== Problem: %d: %s  n: %d nnz: %d\n', ...
68	    i, Problem.name, size (A,1), nnz (A)) ;
69	fprintf ('title: %s\n', Problem.title) ;
70	clear Problem
71	n = size (A,1) ;
72
73        % use AMD from SuiteSparse
74        tic
75        p = amd2 (A) ;
76        t0 = toc ;
77        fprintf ('time: amd     %10.4f\n', t0) ;
78
79	S = A (p,p) ;
80
81	if (doplots)
82	    subplot (3,2,1) ;   spy (A) ;	title ('A original') ;
83	    subplot (3,2,2) ;   spy (S) ;	title ('A permuted') ;
84	    drawnow ;
85	end
86
87	% ensure chol, chol2, and lchol are loaded, for more accurate timing
88	R = chol2 (sparse (1)) ;	    %#ok
89	R = chol (sparse (1)) ;		    %#ok
90	R = lchol (sparse (1)) ;	    %#ok
91	R = ldlchol (sparse (1)) ;	    %#ok
92	R = ldlupdate (sparse (1), sparse (1)) ;	    %#ok
93	c = symbfact (sparse (1)) ;	    %#ok
94
95	tic ;
96	L = lchol (S) ;
97	t3 = toc ;
98	if (doplots)
99	    subplot (3,2,5) ;   spy (L) ;	title ('L=lchol') ;
100	    drawnow ;
101	end
102	fprintf ('CHOLMOD time: L=lchol  %10.4f  nnz(L): %d\n', t3, nnz (L)) ;
103	lnorm = norm (L, 1) ;
104
105	err = ldl_normest (S, L) / lnorm ;
106	if (err > 1e-6)
107	    error ('!') ;
108	end
109	clear L
110
111	tic ;
112	R = chol2 (S) ;
113	t2 = toc ;
114	if (doplots)
115	    subplot (3,2,3) ;   spy (R) ;	title ('R=chol2') ;
116	    drawnow ;
117	end
118	fprintf ('CHOLMOD time: R=chol2  %10.4f  nnz(R): %d\n', t2, nnz (R)) ;
119
120	err = ldl_normest (S, R') / lnorm ;
121	if (err > 1e-6)
122	    error ('!') ;
123	end
124	clear R
125
126	tic ;
127	R = chol (S) ;
128	t1 = toc ;
129	fprintf ('MATLAB time:  R=chol   %10.4f  nnz(R): %d\n', t1, nnz (R)) ;
130	if (doplots)
131	    subplot (3,2,4) ;   spy (R) ;	title ('chol') ;
132	    drawnow ;
133	end
134
135	err = ldl_normest (S, R') / lnorm ;
136	if (err > 1e-6)
137	    error ('!') ;
138	end
139	clear R
140
141	tic ;
142	[count,h,parent,post,R] = symbfact (S) ;
143	t7 = toc ;
144	fprintf ('MATLAB [..,R]=symbfact %10.4f  nnz(R): %d\n', t7, nnz (R)) ;
145
146	fprintf ('\nCHOLMOD speedup vs MATLAB chol:         R: %8.2f L: %8.2f\n\n', ...
147	    t1/t2, t1/t3) ;
148
149	fprintf ('\nCHOLMOD numeric lchol vs MATLAB symbfact:  %8.2f\n', t7/t3) ;
150
151	clear R S
152
153	% use AMD or METIS, doing the ordering in CHOLMOD
154	tic
155	[L,p,q] = lchol (A) ;
156	t4 = toc ;
157	fprintf ('CHOLMOD time: [L,,q]=lchol   %10.4f  nnz(L): %d\n', ...
158	    t4, nnz (L)) ;
159	if (doplots)
160	    subplot (3,2,6) ;   spy (L) ;	title ('[L,p,q]=lchol') ;
161	    drawnow ;
162	end
163
164	err = ldl_normest (A (q,q), L) / lnorm ;
165	if (err > 1e-6)
166	    error ('!') ;
167	end
168	clear L
169
170	% try an LDL' factorization, LD has LDL' factorization of S = A(q,q)
171	tic
172	[LD,p,q] = ldlchol (A) ;
173	t5 = toc ;
174	fprintf ('CHOLMOD time: [L,,q]=ldlchol %10.4f  nnz(L): %d\n', ...
175	    t5, nnz (LD)) ;
176	[L,D] = ldlsplit (LD) ;
177	S = A (q,q) ;
178
179	err = ldl_normest (S, L, D) / lnorm ;
180	if (err > 1e-6)
181	    error ('!') ;
182	end
183	clear L D A
184
185	% update the LDL' factorization (rank 1 to 8).  Pick a C that has
186	% the same pattern as a random set of columns of L, so no fill-in
187	% occurs.  Then add one arbitrary entry, to add some fill-in to L.
188	k = 1 + floor (rand (1) * 8) ;
189	cols = randperm (n) ;
190	cols = cols (1:k) ;
191	C = sprandn (LD (:,cols)) ;
192	row = 1 + floor (rand (1) * n) ;
193	C (row,1) = 1 ;
194
195	if (~isreal (C) | ~isreal (LD))					    %#ok
196	    fprintf ('skip update/downdate of complex matrix ...\n') ;
197	    continue ;
198	end
199
200	tic
201	LD2 = ldlupdate (LD, C) ;
202	t = toc ;
203	fprintf ('\nCHOLMOD time: rank-%d ldlupdate    %10.4f  nnz(L) %d', ...
204	    k, t, nnz (LD2)) ;
205
206	if (nnz (LD2) > nnz (LD))
207	    fprintf ('  with fill-in\n') ;
208	else
209	    fprintf ('  no fill-in\n') ;
210	end
211
212	% check the factorization, LD2 has LDL' factorization of S+C*C'
213	[L,D] = ldlsplit (LD2) ;
214	err = ldl_normest (S + C*C', L, D) / lnorm ;
215	if (err > 1e-6)
216	    error ('!') ;
217	end
218	clear L D
219
220	% downate the LDL' factorization, with just part of C
221	% no change to the pattern occurs.
222	k = max (1, floor (k/2)) ;
223	C1 = C (:, 1:k) ;
224	C2 = C (:, k+1:end) ;		%#ok
225	tic
226	LD3 = ldlupdate (LD2, C1, '-') ;
227	t = toc ;
228	clear LD2
229	fprintf ('CHOLMOD time: rank-%d ldldowndate  %10.4f  nnz(L) %d', ...
230	    k, t, nnz (LD3)) ;
231	fprintf ('  no fill-in\n') ;
232
233	% check the factorization, LD3 has LDL' factorization of A(q,q)+C2*C2'
234	[L,D] = ldlsplit (LD3) ;
235	S2 = S + C*C' - C1*C1' ;
236	err = ldl_normest (S2, L, D) / lnorm ;
237	if (err > 1e-6)
238	    error ('!') ;
239	end
240
241	% now test resymbol
242	LD4 = resymbol (LD3, S2) ;
243	[L,D] = ldlsplit (LD4) ;
244	err = ldl_normest (S2, L, D) / lnorm ;
245	if (err > 1e-6)
246	    error ('!') ;
247	end
248	fprintf ('after resymbol: %d\n', nnz (LD4)) ;
249
250	% compare resymbol with ldlchol
251	LD5 = ldlchol (S2) ;
252	if (nnz (LD5) ~= nnz (LD4))
253	    error ('!') ;
254	end
255
256        % revised June 30, 2020
257        if (nnz (GB_spones_mex (LD5) - GB_spones_mex (LD4)) ~= 0)
258            error ('!') ;
259        end
260
261%       if (nnz (spones (LD5) - spones (LD4)) ~= 0)
262%           % TODO: fails on ssget (878) because MATLAB changed spones(...).
263%           LD5 (262,246)
264%           LD4 (262,246)
265%           spones (LD5) - spones (LD4)
266%           error ('!') ;
267%       end
268
269	b = rand (n,2) ;
270	x = ldlsolve (LD4, b) ;
271	err1 = norm (S2*x-b,1) / norm (S,1) ;
272
273	fprintf ('CHOLMOD residual:  %6.1e\n', err1) ;
274
275	x = S2\b ;
276	err2 = norm (S2*x-b,1) / norm (S,1) ;
277	fprintf ('MATLAB  residual:  %6.1e\n', err2) ;
278
279	b = sprandn (n,3,0.4) ;
280	x = ldlsolve (LD4, b) ;
281	err1 = norm (S2*x-b,1) / norm (S,1) ;
282
283	fprintf ('CHOLMOD residual:  %6.1e (sparse b)\n', err1) ;
284
285	x = S2\b ;
286	err2 = norm (S2*x-b,1) / norm (S,1) ;
287	fprintf ('MATLAB  residual:  %6.1e (sparse b)\n', err2) ;
288
289        % ----------------------------------------------------------------------
290        % test the row delete
291        k = max (1, fix (n/2)) ;
292        tic
293        LD6 = ldlrowmod (LD,k) ;
294        t6 = toc ;
295	fprintf ('\nCHOLMOD time: ldlrowmod, delete    %10.4f  nnz(L) %d', ...
296	    t6, nnz (LD6)) ;
297
298        I = speye (n) ;
299        S2 = S ;
300        S2 (k,:) = I (k,:) ;
301        S2 (:,k) = I (:,k) ;
302
303	% check the factorization, LD6 has LDL' factorization of S2
304	[L,D] = ldlsplit (LD6) ;
305	err = ldl_normest (S2, L, D) / lnorm ;
306	if (err > 1e-6)
307	    error ('!') ;
308	end
309
310        % test the row add, by adding the row back in again
311        Ck = S (:,k) ;
312        S2 (:,k) = Ck ;
313        S2 (k,:) = Ck' ;
314
315        tic
316        LD7 = ldlrowmod (LD6,k,Ck) ;
317        t7 = toc ;
318	fprintf ('\nCHOLMOD time: ldlrowmod, add       %10.4f  nnz(L) %d', ...
319	    t7, nnz (LD7)) ;
320
321	% check the factorization, LD7 has LDL' factorization of S
322	[L,D] = ldlsplit (LD7) ;
323	err = ldl_normest (S, L, D) / lnorm ;
324	if (err > 1e-6)
325	    error ('!') ;
326	end
327
328        % ----------------------------------------------------------------------
329
330    % catch
331    %	fprintf ('failed\n') ;
332    % end
333
334    clear A S C L R LD LD2 LD3 D p q C1 C2 LD3 S2 LD4 b x LD5 I LDL6 LD7 Ck
335
336end
337
338fprintf ('test0 passed\n') ;
339