1function err = test_factorize (A, strategy) 2%TEST_FACTORIZE test the accuracy of the factorization object 3% 4% Example 5% test_factorize (A) ; % where A is square or rectangular, sparse or dense 6% test_factorize (A, strategy) ; % forces a particular strategy; 7% % works only if the matrix is compatible. 8% 9% See also test_all, factorize, inverse, mldivide 10 11% Copyright 2011-2012, Timothy A. Davis, http://www.suitesparse.com 12 13reset_rand ; 14if (nargin < 1) 15 A = rand (100) ; 16end 17 18if (nargin < 2) 19 strategy = '' ; 20end 21 22% do not check the sparsity of the result when using the SVD 23spcheck = ~(strcmp (strategy, 'svd')) ; 24 25err = 0 ; 26if (strcmp (strategy, 'ldl') && issparse (A) && ~isreal(A)) 27 % do not test ldl on sparse complex matrices 28 return ; 29end 30 31[m, n] = size (A) ; 32if (min (m,n) > 0) 33 anorm = norm (A,1) ; 34else 35 anorm = 1 ; 36end 37 38is_symmetric = ((m == n) && (nnz (A-A') == 0)) ; 39 40for nrhs = 1:3 41 42 for bsparse = 0:1 43 44 b = rand (m,nrhs) ; 45 if (bsparse) 46 b = sparse (b) ; 47 end 48 49 %----------------------------------------------------------------------- 50 % test backslash and related methods 51 %----------------------------------------------------------------------- 52 53 for a = [1 pi (pi+1i)] 54 55 % method 0: 56 x = (a*A)\b ; 57 err = check_resid (err, anorm, a*A, x, b, spcheck) ; 58 59 % method 1: 60 S = inverse (A)/a ; 61 x = S*b ; 62 err = check_resid (err, anorm, a*A, x, b, spcheck) ; 63 64 % method 3: 65 F = testfac (A, strategy) ; 66 S = inverse (F)/a ; 67 x = S*b ; 68 err = check_resid (err, anorm, a*A, x, b, spcheck) ; 69 70 % method 4: 71 F = a*testfac (A, strategy) ; 72 x = F\b ; 73 err = check_resid (err, anorm, a*A, x, b, spcheck) ; 74 75 % method 5: 76 S = inverse (F) ; 77 x = S*b ; 78 err = check_resid (err, anorm, a*A, x, b, spcheck) ; 79 80 % method 6: 81 if (m == n) 82 [L, U, p] = lu (A, 'vector') ; 83 x = a * (U \ (L \ (b (p,:)))) ; 84 err = check_resid (err, anorm, A/a, x, b, spcheck) ; 85 86 % method 7: 87 if (is_symmetric) 88 F = a*factorize (A, 'symmetric') ; 89 x = F\b ; 90 err = check_resid (err, anorm, a*A, x, b, spcheck) ; 91 else 92 F = a*factorize (A, 'unsymmetric') ; 93 x = F\b ; 94 err = check_resid (err, anorm, a*A, x, b, spcheck) ; 95 end 96 end 97 98 end 99 100 clear S F 101 102 %------------------------------------------------------------------ 103 % test transposed backslash and related methods 104 %------------------------------------------------------------------ 105 106 b = rand (n,nrhs) ; 107 if (bsparse) 108 b = sparse (b) ; 109 end 110 111 for a = [1 pi (pi+1i)] 112 113 % method 0: 114 x = (a*A)'\b ; 115 err = check_resid (err, anorm, (a*A)', x, b, spcheck) ; 116 117 % method 1: 118 S = inverse (A) / a ; 119 x = S'*b ; 120 err = check_resid (err, anorm, (a*A)', x, b, spcheck) ; 121 122 % method 2: 123 F = a*testfac (A, strategy) ; 124 x = F'\b ; 125 err = check_resid (err, anorm, (a*A)', x, b, spcheck) ; 126 127 % method 3: 128 S = inverse (F') ; 129 x = S*b ; 130 err = check_resid (err, anorm, (a*A)', x, b, spcheck) ; 131 end 132 133 clear S F 134 135 %------------------------------------------------------------------ 136 % test mtimes, times, plus, minus, rdivide, and ldivide 137 %------------------------------------------------------------------ 138 139 for a = [1 pi pi+2i] 140 141 F = a*testfac (A, strategy) ; 142 S = inverse (F) ; 143 B = a*A ; % F is the factorization of a*A 144 145 % test mtimes 146 d = rand (n,1) ; 147 x = F*d ; 148 y = B*d ; 149 z = S\d ; 150 err = check_error (err, norm (mtx(x)-mtx(y),1)) ; 151 err = check_error (err, norm (mtx(x)-mtx(z),1)) ; 152 153 % test mtimes transpose 154 d = rand (m,1) ; 155 x = F'*d ; 156 y = B'*d ; 157 z = S'\d ; 158 err = check_error (err, norm (mtx(x)-mtx(y),1)) ; 159 err = check_error (err, norm (mtx(x)-mtx(z),1)) ; 160 161 % test for scalars 162 for s = [1 42 3-2i] 163 E = s\B - double (s\F) ; err = check_error (err, norm (E,1)) ; 164 E = B/s - double (F/s) ; err = check_error (err, norm (E,1)) ; 165% E = B.*s - F.*s ; err = check_error (err, norm (E,1)) ; 166% E = s.*B - s.*F ; err = check_error (err, norm (E,1)) ; 167% E = s.\B - s.\F ; err = check_error (err, norm (E,1)) ; 168% E = B./s - F./s ; err = check_error (err, norm (E,1)) ; 169 E = B*s - double (F*s) ; err = check_error (err, norm (E,1)) ; 170 E = s*B - double (s*F) ; err = check_error (err, norm (E,1)) ; 171 end 172 173 end 174 175 clear S F C 176 177 %------------------------------------------------------------------ 178 % test slash and related methods 179 %------------------------------------------------------------------ 180 181 b = rand (nrhs,n) ; 182 if (bsparse) 183 b = sparse (b) ; 184 end 185 186 % method 0: 187 x = b/A ; 188 err = check_resid (err, anorm, A', x', b', spcheck) ; 189 190 % method 1: 191 S = inverse (A) ; 192 x = b*S ; 193 err = check_resid (err, anorm, A', x', b', spcheck) ; 194 195 % method 4: 196 F = testfac (A, strategy) ; 197 x = b/F ; 198 err = check_resid (err, anorm, A', x', b', spcheck) ; 199 200 % method 5: 201 S = inverse (F) ; 202 x = b*S ; 203 err = check_resid (err, anorm, A', x', b', spcheck) ; 204 205 % method 6: 206 if (m == n) 207 [L, U, p] = lu (A, 'vector') ; 208 x = (b / U) / L ; x (:,p) = x ; 209 err = check_resid (err, anorm, A', x', b', spcheck) ; 210 end 211 212 %------------------------------------------------------------------ 213 % test transpose slash and related methods 214 %------------------------------------------------------------------ 215 216 b = rand (nrhs,m) ; 217 if (bsparse) 218 b = sparse (b) ; 219 end 220 221 % method 0: 222 x = b/A' ; 223 err = check_resid (err, anorm, A, x', b', spcheck) ; 224 225 % method 1: 226 S = inverse (A)' ; 227 x = b*S ; 228 err = check_resid (err, anorm, A, x', b', spcheck) ; 229 230 % method 4: 231 F = testfac (A, strategy)' ; 232 x = b/F ; 233 err = check_resid (err, anorm, A, x', b', spcheck) ; 234 235 % method 5: 236 S = inverse (F) ; 237 x = b*S ; 238 err = check_resid (err, anorm, A, x', b', spcheck) ; 239 240 %------------------------------------------------------------------ 241 % test double 242 %------------------------------------------------------------------ 243 244 Y = double (inverse (A)) ; 245 if (m == n) 246 Z = inv (A) ; 247 else 248 Z = pinv (full (A)) ; 249 end 250 e = norm (Y-Z,1) ; 251 if (n > 0) 252 e = e / norm (Z,1) ; 253 end 254 err = check_error (e, err) ; 255 256 %------------------------------------------------------------------ 257 % test subsref 258 %------------------------------------------------------------------ 259 260 F = testfac (A, strategy) ; 261 Y = inverse (A) ; 262 if (numel (A) > 1) 263 if (F (end) ~= A (end)) 264 error ('factorization subsref error') ; 265 end 266 if (F.A (end) ~= A (end)) 267 error ('factorization subsref error') ; 268 end 269 end 270 if (n > 0) 271 if (F (1,1) ~= A (1,1)) 272 error ('factorization subsref error') ; 273 end 274 if (F.A (1,1) ~= A (1,1)) 275 error ('factorization subsref error') ; 276 end 277 e = abs (Y (1,1) - Z (1,1)) ; 278 err = check_error (e,err) ; 279 if (m > 1 && n > 1) 280 e = norm (Y (1:2,1:2) - Z (1:2,1:2), 1) ; 281 err = check_error (e,err) ; 282 end 283 end 284 if (m > 3 && n > 1) 285 if (any (F (2:end, 1:2) - A (2:end, 1:2))) 286 error ('factorization subsref error') ; 287 end 288 if (any (F (2:4, :) - A (2:4, :))) 289 error ('factorization subsref error') ; 290 end 291 if (any (F (:, 1:2) - A (:, 1:2))) 292 error ('factorization subsref error') ; 293 end 294 end 295 296 %------------------------------------------------------------------ 297 % test transposed subsref 298 %------------------------------------------------------------------ 299 300 FT = F' ; 301 YT = Y' ; 302 AT = A' ; 303 ZT = Z' ; 304 if (numel (AT) > 1) 305 if (FT (end) ~= AT (end)) 306 error ('factorization subsref error') ; 307 end 308 if (F.A (end) ~= A (end)) 309 error ('factorization subsref error') ; 310 end 311 end 312 313 if (n > 0) 314 if (FT (1,1) ~= AT (1,1)) 315 error ('factorization subsref error') ; 316 end 317 if (F.A (1,1) ~= A (1,1)) 318 error ('factorization subsref error') ; 319 end 320 e = abs (YT (1,1) - ZT (1,1)) ; 321 err = check_error (e,err) ; 322 if (m > 1 && n > 1) 323 e = norm (YT (1:2,1:2) - ZT (1:2,1:2), 1) ; 324 err = check_error (e,err) ; 325 end 326 end 327 328 if (m > 1 && n > 3) 329 if (any (FT (2:end, 1:2) - AT (2:end, 1:2))) 330 error ('factorization subsref error') ; 331 end 332 if (any (FT (2:4, :) - AT (2:4, :))) 333 error ('factorization subsref error') ; 334 end 335 if (any (FT (:, 1:2) - AT (:, 1:2))) 336 error ('factorization subsref error') ; 337 end 338 end 339 340 %------------------------------------------------------------------ 341 % test update/downdate 342 %------------------------------------------------------------------ 343 344 if (isa (F, 'factorization_chol_dense')) 345 w = rand (n,1) ; 346 b = rand (n,1) ; 347 % update 348 G = cholupdate (F,w) ; 349 x = G\b ; 350 err = check_resid (err, anorm, A+w*w', x, b, spcheck) ; 351 % downdate 352 G = cholupdate (G,w,'-') ; 353 x = G\b ; 354 err = check_resid (err, anorm, A, x, b, spcheck) ; 355 clear G 356 end 357 358 %------------------------------------------------------------------ 359 % test size 360 %------------------------------------------------------------------ 361 362 [m1, n1] = size (F) ; 363 [m, n] = size (A) ; 364 if (m1 ~= m || n1 ~= n) 365 error ('size error') ; 366 end 367 [m1, n1] = size (Y) ; 368 if (m1 ~= n || n1 ~= m) 369 error ('pinv size error') ; 370 end 371 if (size (Y,1) ~= n || size (Y,2) ~= m) 372 error ('pinv size error') ; 373 end 374 if (size (F,1) ~= m || size (F,2) ~= n) 375 error ('size error') ; 376 end 377 378 %------------------------------------------------------------------ 379 % test mtimes 380 %------------------------------------------------------------------ 381 382 clear S F Y 383 384 for a = [1 pi pi+2i] 385 386 F = a * testfac (A, strategy) ; 387 S = inverse (F) ; 388 d = rand (1,m) ; 389 x = d*F ; 390 y = d*(A*a) ; 391 z = d/S ; 392 393 err = check_error (err, norm (mtx(x)-mtx(y),1)) ; 394 err = check_error (err, norm (mtx(x)-mtx(z),1)) ; 395 396 d = rand (1,n) ; 397 x = d*F' ; 398 y = d*(A*a)' ; 399 z = d/S' ; 400 401 err = check_error (err, norm (mtx(x)-mtx(y),1)) ; 402 err = check_error (err, norm (mtx(x)-mtx(z),1)) ; 403 404 end 405 406 %------------------------------------------------------------------ 407 % test inverse 408 %------------------------------------------------------------------ 409 410 Y = double (inverse (inverse (A))) ; 411 e = norm (A-Y,1) ; 412 if (e > 0) 413 error ('inverse error') ; 414 end 415 416 %------------------------------------------------------------------ 417 % test mldivide and mrdivide with a matrix b, transpose, etc 418 %------------------------------------------------------------------ 419 420 if (max (m,n) < 100) 421 F = testfac (A, strategy) ; 422 B = rand (m,n) ; 423 err = check_error (err, norm (B\A - mtx(B\F), 1) / anorm) ; 424 err = check_error (err, norm (A/B - mtx(F/B), 1) / anorm) ; 425 err = check_error (err, norm (B*pinv(full(A))-mtx(B/F), 1) / anorm); 426 end 427 end 428end 429 430fprintf ('.') ; 431 432%-------------------------------------------------------------------------- 433 434function err = check_resid (err, anorm, A, x, b, spcheck) 435[m, n] = size (A) ; 436 437x = mtx (x) ; 438 439if (m <= n) 440 e = norm (A*x-b,1) / (anorm + norm (x,1)) ; 441else 442 e = norm (A'*(A*x)-A'*b,1) / (anorm + norm (x,1)) ; 443end 444 445if (min (m,n) > 1 && spcheck) 446 if (issparse (A) && issparse (b)) 447 if (~issparse (x)) 448 error ('x must be sparse') ; 449 end 450 else 451 if (issparse (x)) 452 error ('x must be full') ; 453 end 454 end 455end 456 457err = check_error (e, err) ; 458 459%-------------------------------------------------------------------------- 460 461function x = mtx (x) 462% make sure that x is a matrix. It might be a factorization. 463if (isobject (x)) 464 x = double (x) ; 465end 466 467%-------------------------------------------------------------------------- 468 469function err = check_error (err1, err2) 470err = max (err1, err2) ; 471if (err > 1e-8) 472 fprintf ('error: %8.3e\n', full (err)) ; 473 error ('error too high') ; 474end 475err = full (err) ; 476 477%-------------------------------------------------------------------------- 478 479function F = testfac (A, strategy) 480if (isempty (strategy)) 481 F = factorize (A) ; 482else 483 F = factorize (A, strategy) ; 484end 485