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