1function test7 2%TEST7 test cs_lu 3% 4% Example: 5% test7 6% See also: testall 7 8% Copyright 2006-2012, Timothy A. Davis, http://www.suitesparse.com 9 10index = ssget ; 11[ignore f] = sort (max (index.nrows, index.ncols)) ; 12f = f (1:100) ; 13 14clf 15 16maxerr1 = 0 ; 17maxerr2 = 0 ; 18 19for i = f 20 Prob = ssget (i) ; 21 disp (Prob) ; 22 A = Prob.A ; 23 24 [m n] = size (A) ; 25 if (m ~= n) 26 continue 27 end 28 29 for cmplex = 0:1 30 31 if (cmplex) 32 A = A + norm(A,1) * sprand (A) / 3 ; 33 end 34 35 [L,U,P] = lu (A) ; 36 37 udiag = full (diag (U)) ; 38 umin = min (abs (udiag)) ; 39 fprintf ('umin %g\n', umin) ; 40 41 if (umin > 1e-14) 42 43 [L2,U2,p] = cs_lu (A) ; 44 45 subplot (3,4,1) ; spy (A) ; 46 subplot (3,4,2) ; spy (A(p,:)) ; 47 subplot (3,4,3) ; spy (L2) ; 48 subplot (3,4,4) ; spy (U2) ; 49 50 err1 = norm (L*U-P*A,1) / norm (A,1) ; 51 err2 = norm (L2*U2-A(p,:),1) / norm (A,1) ; 52 fprintf ('err %g %g\n', err1, err2) ; 53 54 if (err1 > 1e-10 | err2 > 1e-10) %#ok 55 error ('!') ; 56 end 57 maxerr1 = max (maxerr1, err1) ; 58 maxerr2 = max (maxerr2, err2) ; 59 60 end 61 62 q = colamd (A) ; 63 64 [L,U,P] = lu (A (:,q)) ; 65 66 udiag = full (diag (U)) ; 67 umin = min (abs (udiag)) ; 68 fprintf ('umin %g with q\n', umin) ; 69 70 if (umin > 1e-14) 71 72 [L2,U2,p,q2] = cs_lu (A) ; 73 74 subplot (3,4,5) ; spy (A) ; 75 subplot (3,4,6) ; spy (A(p,q2)) ; 76 subplot (3,4,7) ; spy (L2) ; 77 subplot (3,4,8) ; spy (U2) ; 78 79 err1 = norm (L*U-P*A(:,q),1) / norm (A,1) ; 80 err2 = norm (L2*U2-A(p,q2),1) / norm (A,1) ; 81 fprintf ('err %g %g\n', err1, err2) ; 82 83 if (err1 > 1e-10 | err2 > 1e-10) %#ok 84 error ('!') ; 85 end 86 maxerr1 = max (maxerr1, err1) ; 87 maxerr2 = max (maxerr2, err2) ; 88 end 89 90 91 try 92 q = amd (A) ; 93 catch 94 q = symamd (A) ; 95 end 96 97 tol = 0.01 ; 98 99 [L,U,P] = lu (A (q,q), tol) ; 100 101 udiag = full (diag (U)) ; 102 umin = min (abs (udiag)) ; 103 fprintf ('umin %g with amd q\n', umin) ; 104 105 if (umin > 1e-14) 106 107 [L2,U2,p,q2] = cs_lu (A,tol) ; 108 109 subplot (3,4,9) ; spy (A) ; 110 subplot (3,4,10) ; spy (A(p,q2)) ; 111 subplot (3,4,11) ; spy (L2) ; 112 subplot (3,4,12) ; spy (U2) ; 113 114 err1 = norm (L*U-P*A(q,q),1) / norm (A,1) ; 115 err2 = norm (L2*U2-A(p,q2),1) / norm (A,1) ; 116 lbig = full (max (max (abs (L2)))) ; 117 fprintf ('err %g %g lbig %g\n', err1, err2, lbig) ; 118 if (lbig > 1/tol) 119 error ('L!') ; 120 end 121 122 if (err1 > 1e-10 | err2 > 1e-10) %#ok 123 error ('!') ; 124 end 125 maxerr1 = max (maxerr1, err1) ; 126 maxerr2 = max (maxerr2, err2) ; 127 end 128 129 drawnow 130 % pause 131 132 end 133end 134 135fprintf ('maxerr %g %g\n', maxerr1, maxerr2) ; 136