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