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