1function umfpack_test (nmat) 2%UMFPACK_TEST for testing umfpack (requires ssget) 3% 4% Example: 5% umfpack_test 6% umfpack_test (100) % runs the first 100 matrices 7% See also umfpack 8 9% Copyright 1995-2007 by Timothy A. Davis. 10 11index = ssget ; 12 13f = find (index.nrows == index.ncols) ; 14[ignore, i] = sort (index.nrows (f)) ; 15f = f (i) ; 16 17if (nargin < 1) 18 nmat = length (f) ; 19else 20 nmat = min (nmat, length (f)) ; 21end 22nmat = max (nmat, 1) ; 23f = f (1:nmat) ; 24 25Control = umfpack ; 26Control.prl = 0 ; 27 28clf 29 30h = waitbar (0, 'UMFPACK test') ; 31 32for k = 1:nmat 33 34 i = f (k) ; 35 waitbar (k/nmat, h, 'UMFPACK test') ; 36 37% try 38 39 fprintf ('\nmatrix: %s %s %d\n', ... 40 index.Group{i}, index.Name{i}, index.nrows(i)) ; 41 42 Prob = ssget (i) ; 43 A = Prob.A ; 44 n = size (A,1) ; 45 46 b = rand (1,n) ; 47 c = b' ; 48 49 %----------------------------------------------------------------------- 50 % symbolic factorization 51 %----------------------------------------------------------------------- 52 53 [P1, Q1, Fr, Ch, Info] = umfpack (A, 'symbolic') ; 54 subplot (2,2,1) 55 spy (A) 56 title ('A') 57 58 subplot (2,2,2) 59 treeplot (Fr (1:end-1,2)') ; 60 title ('supercolumn etree') 61 62 %----------------------------------------------------------------------- 63 % P(R\A)Q = LU 64 %----------------------------------------------------------------------- 65 66 [L,U,P,Q,R,Info] = umfpack (A, struct ('details',1)) ; 67 err = lu_normest (P*(R\A)*Q, L, U) ; 68 fprintf ('norm est PR\\AQ-LU: %g relative: %g\n', ... 69 err, err / norm (A,1)) ; 70 71 subplot (2,2,3) 72 spy (P*A*Q) 73 title ('PAQ') ; 74 75 cs = Info.number_of_column_singletons ; % (57) ; 76 rs = Info.number_of_row_singletons ; % (58) ; 77 78 subplot (2,2,4) 79 hold off 80 try 81 spy (L+U) 82 catch 83 end 84 hold on 85 if (cs > 0) 86 plot ([0 cs n n 0] + .5, [0 cs cs 0 0]+.5, 'c') ; 87 end 88 if (rs > 0) 89 plot ([0 rs rs 0 0] + cs +.5, [cs cs+rs n n cs]+.5, 'r') ; 90 end 91 92 title ('LU factors') 93 drawnow 94 95 %----------------------------------------------------------------------- 96 % PAQ = LU 97 %----------------------------------------------------------------------- 98 99 [L,U,P,Q] = umfpack (A) ; 100 err = lu_normest (P*A*Q, L, U) ; 101 fprintf ('norm est PAQ-LU: %g relative: %g\n', ... 102 err, err / norm (A,1)) ; 103 104 %----------------------------------------------------------------------- 105 % solve 106 %----------------------------------------------------------------------- 107 108 x1 = b/A ; 109 y1 = A\c ; 110 m1 = norm (b-x1*A) ; 111 m2 = norm (A*y1-c) ; 112 113 % factor the transpose 114 Control.irstep = 2 ; 115 [x, info] = umfpack (A', '\', c, Control) ; 116 lunz0 = info.nnz_in_L_plus_U ; 117 r = norm (A'*x-c) ; 118 119 fprintf (':: %8.2e matlab: %8.2e %8.2e\n', r, m1, m2) ; 120 121 % factor the original matrix and solve xA=b 122 for ir = 0:4 123 Control.irstep = ir ; 124 [x, info] = umfpack (b, '/', A, Control) ; 125 r = norm (b-x*A) ; 126 if (ir == 0) 127 lunz1 = info.nnz_in_L_plus_U ; 128 end 129 fprintf ('%d: %8.2e %d\n', ir, r, info.iterative_refinement_steps) ; 130 end 131 132 % factor the original matrix and solve Ax=b 133 for ir = 0:4 134 Control.irstep = ir ; 135 [x, info] = umfpack (A, '\', c, Control) ; 136 r = norm (A*x-c) ; 137 fprintf ('%d: %8.2e %d\n', ir, r, info.iterative_refinement_steps) ; 138 end 139 140 fprintf (... 141 'lunz trans %12d no trans: %12d trans/notrans: %10.4f\n', ... 142 lunz0, lunz1, lunz0 / lunz1) ; 143 144 %----------------------------------------------------------------------- 145 % get the determinant 146 %----------------------------------------------------------------------- 147 148 det1 = det (A) ; 149 det2 = umfpack (A, 'det') ; 150 [det3 dexp3] = umfpack (A, 'det') ; 151 err = abs (det1-det2) ; 152 err3 = abs (det1 - (det3 * 10^dexp3)) ; 153 denom = det1 ; 154 if (denom == 0) 155 denom = 1 ; 156 end 157 err = err / denom ; 158 err3 = err3 / denom ; 159 fprintf ('det: %20.12e + (%20.12e)i MATLAB\n', ... 160 real(det1), imag(det1)) ; 161 fprintf ('det: %20.12e + (%20.12e)i umfpack\n', ... 162 real(det2), imag(det2)) ; 163 fprintf ('det: (%20.12e + (%20.12e)i) * 10^(%g) umfpack\n', ... 164 real(det3), imag(det3), dexp3) ; 165 fprintf ('diff %g %g\n', err, err3) ; 166 167% catch 168% % out-of-memory is OK, other errors are not 169% disp (lasterr) ; 170% if (isempty (strfind (lasterr, 'Out of memory'))) 171% error (lasterr) ; %#ok 172% else 173% fprintf ('test terminated early, but otherwise OK\n') ; 174% end 175% end 176 177end 178 179close (h) ; % close the waitbar 180