|
5164
|
1 function x = umfpack_btf (A, b, Control) |
|
|
2 % UMFPACK_BTF |
|
|
3 % |
|
|
4 % x = umfpack_btf (A, b, Control) |
|
|
5 % |
|
|
6 % solve Ax=b by first permuting the matrix A to block triangular form via dmperm |
|
|
7 % and then using UMFPACK to factorize each diagonal block. Adjacent 1-by-1 |
|
|
8 % blocks are merged into a single upper triangular block, and solved via |
|
|
9 % MATLAB's \ operator. The Control parameter is optional (Type umfpack_details |
|
|
10 % and umfpack_report for details on its use). A must be square. |
|
|
11 % |
|
|
12 % See also: umfpack, umfpack_factorize, umfpack_details, dmperm |
|
|
13 |
|
|
14 % UMFPACK Version 4.3 (Jan. 16, 2004), Copyright (c) 2004 by Timothy A. |
|
|
15 % Davis. All Rights Reserved. Type umfpack_details for License. |
|
|
16 |
|
|
17 if (nargin < 2) |
|
|
18 help umfpack_btf |
|
|
19 error ('Usage: x = umfpack_btf (A, b, Control)') ; |
|
|
20 end |
|
|
21 |
|
|
22 [m n] = size (A) ; |
|
|
23 if (m ~= n) |
|
|
24 help umfpack_btf |
|
|
25 error ('umfpack_btf: A must be square') ; |
|
|
26 end |
|
|
27 [m1 n1] = size (b) ; |
|
|
28 if (m1 ~= n) |
|
|
29 help umfpack_btf |
|
|
30 error ('umfpack_btf: b has the wrong dimensions') ; |
|
|
31 end |
|
|
32 |
|
|
33 if (nargin < 3) |
|
|
34 Control = umfpack ; |
|
|
35 end |
|
|
36 |
|
|
37 %------------------------------------------------------------------------------- |
|
|
38 % find the block triangular form |
|
|
39 %------------------------------------------------------------------------------- |
|
|
40 |
|
|
41 [p,q,r] = dmperm (A) ; |
|
|
42 nblocks = length (r) - 1 ; |
|
|
43 |
|
|
44 %------------------------------------------------------------------------------- |
|
|
45 % solve the system |
|
|
46 %------------------------------------------------------------------------------- |
|
|
47 |
|
|
48 if (nblocks == 1 | sprank (A) < n) |
|
|
49 |
|
|
50 %--------------------------------------------------------------------------- |
|
|
51 % matrix is irreducible or structurally singular |
|
|
52 %--------------------------------------------------------------------------- |
|
|
53 |
|
|
54 x = umfpack_solve (A, '\\', b, Control) ; |
|
|
55 |
|
|
56 else |
|
|
57 |
|
|
58 %--------------------------------------------------------------------------- |
|
|
59 % A (p,q) is in block triangular form |
|
|
60 %--------------------------------------------------------------------------- |
|
|
61 |
|
|
62 b = b (p,:) ; |
|
|
63 A = A (p,q) ; |
|
|
64 x = zeros (size (b)) ; |
|
|
65 |
|
|
66 %--------------------------------------------------------------------------- |
|
|
67 % merge adjacent singletons into a single upper triangular block |
|
|
68 %--------------------------------------------------------------------------- |
|
|
69 |
|
|
70 [r, nblocks, is_triangular] = merge_singletons (r) ; |
|
|
71 |
|
|
72 %--------------------------------------------------------------------------- |
|
|
73 % solve the system: x (q) = A\b |
|
|
74 %--------------------------------------------------------------------------- |
|
|
75 |
|
|
76 for k = nblocks:-1:1 |
|
|
77 |
|
|
78 % get the kth block |
|
|
79 k1 = r (k) ; |
|
|
80 k2 = r (k+1) - 1 ; |
|
|
81 |
|
|
82 % solve the system |
|
|
83 x (k1:k2,:) = solver (A (k1:k2, k1:k2), b (k1:k2,:), ... |
|
|
84 is_triangular (k), Control) ; |
|
|
85 |
|
|
86 % off-diagonal block back substitution |
|
|
87 b (1:k1-1,:) = b (1:k1-1,:) - A (1:k1-1, k1:k2) * x (k1:k2,:) ; |
|
|
88 |
|
|
89 end |
|
|
90 |
|
|
91 x (q,:) = x ; |
|
|
92 |
|
|
93 end |
|
|
94 |
|
|
95 %------------------------------------------------------------------------------- |
|
|
96 % merge_singletons |
|
|
97 %------------------------------------------------------------------------------- |
|
|
98 |
|
|
99 function [r, nblocks, is_triangular] = merge_singletons (r) |
|
|
100 % |
|
|
101 % Given r from [p,q,r] = dmperm (A), where A is square, return a modified r that |
|
|
102 % reflects the merger of adjacent singletons into a single upper triangular |
|
|
103 % block. is_triangular (k) is 1 if the kth block is upper triangular. nblocks |
|
|
104 % is the number of new blocks. |
|
|
105 |
|
|
106 nblocks = length (r) - 1 ; |
|
|
107 bsize = r (2:nblocks+1) - r (1:nblocks) ; |
|
|
108 t = [0 (bsize == 1)] ; |
|
|
109 z = (t (1:nblocks) == 0 & t (2:nblocks+1) == 1) | t (2:nblocks+1) == 0 ; |
|
|
110 y = [(find (z)) nblocks+1] ; |
|
|
111 r = r (y) ; |
|
|
112 nblocks = length (y) - 1 ; |
|
|
113 is_triangular = y (2:nblocks+1) - y (1:nblocks) > 1 ; |
|
|
114 |
|
|
115 %------------------------------------------------------------------------------- |
|
|
116 % solve Ax=b, but check for small and/or triangular systems |
|
|
117 %------------------------------------------------------------------------------- |
|
|
118 |
|
|
119 function x = solver (A, b, is_triangular, Control) |
|
|
120 if (is_triangular) |
|
|
121 % back substitution only |
|
|
122 x = A \ b ; |
|
|
123 elseif (size (A,1) < 4) |
|
|
124 % a very small matrix, solve it as a dense linear system |
|
|
125 x = full (A) \ b ; |
|
|
126 else |
|
|
127 % solve it as a sparse linear system |
|
|
128 x = umfpack_solve (A, '\\', b, Control) ; |
|
|
129 end |
