|
5164
|
1 % UMFPACK_TEST: test UMFPACK solve: b/A, A\b with iterative refinement |
|
|
2 % Requires the UFsparse package for downloading matrices from the UF |
|
|
3 % sparse matrix library. |
|
|
4 % |
|
|
5 % UMFPACK Version 4.4, Copyright (c) 2005 by Timothy A. Davis. |
|
|
6 % All Rights Reserved. Type umfpack_details for License. |
|
|
7 |
|
|
8 index = UFget ; |
|
|
9 |
|
|
10 f = find (index.nrows == index.ncols) ; |
|
|
11 [ignore, i] = sort (index.nrows (f)) ; |
|
|
12 f = f (i) ; |
|
|
13 |
|
|
14 Control = umfpack ; |
|
|
15 Control (1) = 0 ; |
|
|
16 |
|
|
17 warning ('off', 'all') ; |
|
|
18 figure (1) |
|
|
19 clf |
|
|
20 |
|
|
21 |
|
|
22 for i = f |
|
|
23 |
|
|
24 fprintf ('\nmatrix: %s %s %d\n', index.Group{i}, index.Name{i}, index.nrows(i)) ; |
|
|
25 |
|
|
26 Prob = UFget (i) ; |
|
|
27 A = Prob.A ; |
|
|
28 n = size (A,1) ; |
|
|
29 |
|
|
30 b = rand (1,n) ; |
|
|
31 c = b' ; |
|
|
32 |
|
|
33 try |
|
|
34 |
|
|
35 %----------------------------------------------------------------------- |
|
|
36 % symbolic factorization |
|
|
37 %----------------------------------------------------------------------- |
|
|
38 |
|
|
39 [P1, Q1, Fr, Ch, Info] = umfpack (A, 'symbolic') ; |
|
|
40 subplot (2,2,1) |
|
|
41 spy (A) |
|
|
42 title ('A') |
|
|
43 |
|
|
44 subplot (2,2,2) |
|
|
45 treeplot (Fr (1:end-1,2)') ; |
|
|
46 title ('supercolumn etree') |
|
|
47 |
|
|
48 %----------------------------------------------------------------------- |
|
|
49 % P(R\A)Q = LU |
|
|
50 %----------------------------------------------------------------------- |
|
|
51 |
|
|
52 [L,U,P,Q,R,Info] = umfpack (A) ; |
|
|
53 err = lu_normest (P*(R\A)*Q, L, U) ; |
|
|
54 fprintf ('norm est PR\\AQ-LU: %g relative: %g\n', ... |
|
|
55 err, err / norm (A,1)) ; |
|
|
56 |
|
|
57 subplot (2,2,3) |
|
|
58 spy (P*A*Q) |
|
|
59 title ('PAQ') ; |
|
|
60 |
|
|
61 cs = Info (57) ; |
|
|
62 rs = Info (58) ; |
|
|
63 |
|
|
64 subplot (2,2,4) |
|
|
65 hold off |
|
|
66 spy (L|U) |
|
|
67 hold on |
|
|
68 if (cs > 0) |
|
|
69 plot ([0 cs n n 0] + .5, [0 cs cs 0 0]+.5, 'c') ; |
|
|
70 end |
|
|
71 if (rs > 0) |
|
|
72 plot ([0 rs rs 0 0] + cs +.5, [cs cs+rs n n cs]+.5, 'r') ; |
|
|
73 end |
|
|
74 |
|
|
75 title ('LU factors') |
|
|
76 drawnow |
|
|
77 |
|
|
78 %----------------------------------------------------------------------- |
|
|
79 % PAQ = LU |
|
|
80 %----------------------------------------------------------------------- |
|
|
81 |
|
|
82 [L,U,P,Q] = umfpack (A) ; |
|
|
83 err = lu_normest (P*A*Q, L, U) ; |
|
|
84 fprintf ('norm est PAQ-LU: %g relative: %g\n', ... |
|
|
85 err, err / norm (A,1)) ; |
|
|
86 |
|
|
87 %----------------------------------------------------------------------- |
|
|
88 % solve |
|
|
89 %----------------------------------------------------------------------- |
|
|
90 |
|
|
91 x1 = b/A ; |
|
|
92 y1 = A\c ; |
|
|
93 m1 = norm (b-x1*A) ; |
|
|
94 m2 = norm (A*y1-c) ; |
|
|
95 |
|
|
96 % factor the transpose |
|
|
97 Control (8) = 2 ; |
|
|
98 [x, info] = umfpack (A', '\', c, Control) ; |
|
|
99 lunz0 = info (44) + info (45) - info (67) ; |
|
|
100 r = norm (A'*x-c) ; |
|
|
101 |
|
|
102 fprintf (':: %8.2e matlab: %8.2e %8.2e\n', r, m1, m2) ; |
|
|
103 |
|
|
104 % factor the original matrix and solve xA=b |
|
|
105 for ir = 0:4 |
|
|
106 Control (8) = ir ; |
|
|
107 [x, info] = umfpack (b, '/', A, Control) ; |
|
|
108 r = norm (b-x*A) ; |
|
|
109 if (ir == 0) |
|
|
110 lunz1 = info (44) + info (45) - info (67) ; |
|
|
111 end |
|
|
112 fprintf ('%d: %8.2e %d %d\n', ir, r, info (81), info (82)) ; |
|
|
113 end |
|
|
114 |
|
|
115 % factor the original matrix and solve Ax=b |
|
|
116 for ir = 0:4 |
|
|
117 Control (8) = ir ; |
|
|
118 [x, info] = umfpack (A, '\', c, Control) ; |
|
|
119 r = norm (A*x-c) ; |
|
|
120 fprintf ('%d: %8.2e %d %d\n', ir, r, info (81), info (82)) ; |
|
|
121 end |
|
|
122 |
|
|
123 fprintf ('lunz trans %12d no trans: %12d trans/notrans: %10.4f\n', ... |
|
|
124 lunz0, lunz1, lunz0 / lunz1) ; |
|
|
125 |
|
|
126 %----------------------------------------------------------------------- |
|
|
127 % get the determinant |
|
|
128 %----------------------------------------------------------------------- |
|
|
129 |
|
|
130 det1 = det (A) ; |
|
|
131 det2 = umfpack (A, 'det') ; |
|
|
132 [det3 dexp3] = umfpack (A, 'det') ; |
|
|
133 err = abs (det1-det2) ; |
|
|
134 err3 = abs (det1 - (det3 * 10^dexp3)) ; |
|
|
135 denom = det1 ; |
|
|
136 if (denom == 0) |
|
|
137 denom = 1 ; |
|
|
138 end |
|
|
139 err = err / denom ; |
|
|
140 err3 = err3 / denom ; |
|
|
141 fprintf ('det: %24.16e + (%24.16e)i MATLAB\n', real(det1), imag(det1)) ; |
|
|
142 fprintf ('det: %24.16e + (%24.16e)i umfpack\n',real(det2), imag(det2)) ; |
|
|
143 fprintf ('det: (%24.16e + (%24.16e)i) * 10^(%g) umfpack\n', real(det3), imag(det3), dexp3) ; |
|
|
144 fprintf ('diff %g %g\n', err, err3) ; |
|
|
145 |
|
|
146 catch |
|
|
147 fprintf ('failed\n') ; |
|
|
148 end |
|
|
149 |
|
|
150 % pause |
|
|
151 |
|
|
152 end |
