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