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