1function F = factorize (A,strategy,burble) 2%FACTORIZE an object-oriented method for solving linear systems 3% and least-squares problems, and for representing operations with the 4% inverse of a square matrix or the pseudo-inverse of a rectangular matrix. 5% 6% F = factorize(A) returns an object F that holds the factorization of A. 7% x=F\b then solves a linear system or a least-squares problem. S=inverse(F) 8% or S=inverse(A) returns a factorized representation of the inverse of A so 9% that inverse(A)*b is mathematically equivalent to pinv(A)*b, but the former 10% does not actually compute the inverse or pseudo-inverse of A. 11% 12% Example 13% 14% F = factorize(A) ; % LU, QR, or Cholesky factorization of A 15% x = F\b ; % solve A*x=b; same as x=A\b 16% S = inverse (F) ; % S represents the factorization of inv(A) 17% x = S*b ; % same as x = A\b. 18% E = A-B*inverse(D)*C % efficiently computes the Schur complement 19% E = A-B*inv(D)*C % bad method for computing the Schur complement 20% S = inverse(A) ; S(:,1) % compute just the first column of inv(A), 21% % without computing inv(A) 22% 23% F = factorize (A, strategy, burble) ; % optional 2nd and 3rd inputs 24% 25% A string can be specified as a second input parameter to select the strategy 26% used to factorize the matrix. The first two are meta-strategies: 27% 28% 'default' if rectangular 29% use QR for sparse A or A' (whichever is tall and thin); 30% use COD for full A 31% else 32% if symmetric 33% if positive real diagonal: try CHOL 34% else (or if CHOL fails): try LDL 35% end 36% if not yet factorized: try LU (fails if rank-deficient) 37% end 38% if all else fails, or if QR or LU report that the matrix 39% is singular (or nearly so): use COD 40% This strategy mimics backslash, except that backslash never 41% uses COD. Backslash also exploits other solvers, such as 42% specialized tridiagonal and banded solvers. 43% 44% 'symmetric' as 'default' above, but assumes A is symmetric without 45% checking, which is faster if you already know A is symmetric. 46% Uses tril(A) and assumes triu(A) is its transpose. Results 47% will be incorrect if A is not symmetric. If A is rectangular, 48% the 'default' strategy is used instead. 49% 50% 'unsymmetric' as 'default', but assumes A is unsymmetric. 51% 52% The next "strategies" just select a single method, listed in decreasing order 53% of generality and increasing order of speed and memory efficiency. All of 54% them except the SVD can exploit sparsity. 55% 56% 'svd' use SVD. Never fails ... unless it runs out of time or memory. 57% Coerces a sparse matrix A to full. 58% 59% 'cod' use COD. Almost as accurate as SVD, and much faster. Based 60% on dense or sparse QR with rank estimation. Handles 61% rank-deficient matrices, as long as it correctly estimates 62% the rank. If the rank is ill-defined, use the SVD instead. 63% Sparse COD requires the SPQR package to be installed 64% (see http://www.suitesparse.com). 65% 66% 'qr' use QR. Reports a warning if A is singular. 67% 68% 'lu' use LU. Fails if A is rectangular; warning if A singular. 69% 70% 'ldl' use LDL. Fails if A is rank-deficient or not symmetric, or if 71% A is sparse and complex. Uses tril(A) and assumes triu(A) 72% is the transpose of tril(A). 73% 74% 'chol' use CHOL. Fails if A is rank-deficient or not symmetric 75% positive definite. If A is sparse, it uses tril(A) and 76% assumes triu(A) is the transpose of tril(A). If A is dense, 77% triu(A) is used instead. 78% 79% A third input, burble, can be provided to tell this function to print what 80% methods it tries (if burble is nonzero). 81% 82% For a demo type "fdemo" in the Demo directory or see the Doc/ directory. 83% 84% See also inverse, slash, linsolve, spqr. 85 86% Copyright 2011-2012, Timothy A. Davis, http://www.suitesparse.com 87 88assert (ndims (A) == 2, 'Matrix must be 2D.') ; 89 90if (nargin < 2 || isempty (strategy)) 91 strategy = 'default' ; 92end 93if (nargin < 3) 94 burble = 0 ; 95end 96 97if (burble) 98 fprintf ('\nfactorize: strategy %s, A has size %d-by-%d, ', ... 99 strategy, size (A)) ; 100 if (issparse (A)) 101 fprintf ('sparse with %d nonzeros.\n', nnz (A)) ; 102 else 103 fprintf ('full.\n') ; 104 end 105end 106 107switch strategy 108 109 case 'default' 110 [F, me] = backslash_mimic (A, burble, 0) ; 111 112 case 'symmetric' 113 [F, me] = backslash_mimic (A, burble, 1) ; 114 115 case 'unsymmetric' 116 [F, me] = backslash_mimic (A, burble, 2) ; 117 118 case 'svd' 119 [F, me] = factorize_svd (A, burble) ; 120 121 case 'cod' 122 [F, me] = factorize_cod (A, burble) ; 123 124 case 'qr' 125 [F, me] = factorize_qr (A, burble, 0) ; 126 127 case 'lu' 128 % do not report a failure if the matrix is singular 129 [F, me] = factorize_lu (A, burble, 0) ; 130 131 case 'ldl' 132 [F, me] = factorize_ldl (A, burble) ; 133 134 case 'chol' 135 [F, me] = factorize_chol (A, burble) ; 136 137 otherwise 138 error ('FACTORIZE:invalidStrategy', 'unrecognized strategy.') ; 139 140end 141 142if (~isobject (F)) 143 throw (me) ; 144end 145 146%------------------------------------------------------------------------------- 147 148function [F, me] = backslash_mimic (A, burble, strategy) 149%BACKSLASH_MIMIC automatically select a method to factorize A. 150F = [ ] ; 151me = [ ] ; 152[m, n] = size (A) ; 153 154% If the following condition is true, then the QR, QRT, or LU factorizations 155% will report a failure if A is singular (or nearly so). This allows COD 156% or COD_SPARSE to then be used instead. COD_SPARSE relies on the SPQR 157% mexFunction in SuiteSparse, which might not be installed. In this case, 158% QR, QRT, and LU do not report failures for sparse matrices that are singular 159% (or nearly so), since there is no COD_SPARSE to fall back on. 160fail_if_singular = ~issparse (A) || (exist ('spqr') == 3) ; %#ok 161try_cod = true ; 162 163if (m ~= n) 164 if (issparse (A)) 165 % Use QR for the sparse rectangular case (ignore 'strategy' argument). 166 [F, me] = factorize_qr (A, burble, fail_if_singular) ; 167 else 168 % Use COD for the full rectangular case (ignore 'strategy' argument). 169 % If this fails, there's no reason to retry the COD below. If A has 170 % full rank, then COD is the same as QR with column pivoting (with the 171 % same cost in terms of run time and memory). Backslash in MATLAB uses 172 % QR with column pivoting alone, so this is just as fast as x=A\b in 173 % the full-rank case, but gives a more reliable result in the rank- 174 % deficient case. 175 try_cod = false ; 176 [F, me] = factorize_cod (A, burble) ; 177 end 178else 179 % square case: Cholesky, LDL, or LU factorization of A 180 switch strategy 181 case 0 182 is_symmetric = (nnz (A-A') == 0) ; 183 case 1 184 is_symmetric = true ; 185 case 2 186 is_symmetric = false ; 187 end 188 if (is_symmetric) 189 % A is symmetric (or assumed to be so) 190 d = diag (A) ; 191 if (all (d > 0) && nnz (imag (d)) == 0) 192 % try a Cholesky factorization 193 [F, me] = factorize_chol (A, burble) ; 194 end 195 if (~isobject (F) && (~issparse (A) || isreal (A))) 196 % try an LDL factorization. 197 % complex sparse LDL does not yet exist in MATLAB 198 [F, me] = factorize_ldl (A, burble) ; 199 end 200 end 201 if (~isobject (F)) 202 % use LU if Cholesky and/or LDL failed, or were skipped. 203 [F, me] = factorize_lu (A, burble, fail_if_singular) ; 204 end 205end 206if (~isobject (F) && try_cod) 207 % everything else failed, matrix is rank-deficient. Use COD 208 [F, me] = factorize_cod (A, burble) ; 209end 210 211 212%------------------------------------------------------------------------------- 213 214function [F, me] = factorize_qr (A, burble, fail_if_singular) 215% QR fails if the matrix is rank-deficient. 216F = [ ] ; 217me = [ ] ; 218try 219 [m, n] = size (A) ; 220 if (m >= n) 221 if (burble) 222 fprintf ('factorize: try QR of A ... ') ; 223 end 224 if (issparse (A)) 225 F = factorization_qr_sparse (A, fail_if_singular) ; 226 else 227 F = factorization_qr_dense (A, fail_if_singular) ; 228 end 229 else 230 if (burble) 231 fprintf ('factorize: try QR of A'' ... ') ; 232 end 233 if (issparse (A)) 234 F = factorization_qrt_sparse (A, fail_if_singular) ; 235 else 236 F = factorization_qrt_dense (A, fail_if_singular) ; 237 end 238 end 239 if (burble) 240 fprintf ('OK.\n') ; 241 end 242catch me 243 if (burble) 244 fprintf ('failed.\nfactorize: %s\n', me.message) ; 245 end 246end 247 248 249%------------------------------------------------------------------------------- 250 251function [F, me] = factorize_chol (A, burble) 252% LDL fails if the matrix is rectangular, rank-deficient, or not positive 253% definite. Only the lower triangular part of A is used. 254F = [ ] ; 255me = [ ] ; 256try 257 if (burble) 258 fprintf ('factorize: try CHOL ... ') ; 259 end 260 if (issparse (A)) 261 F = factorization_chol_sparse (A) ; 262 else 263 F = factorization_chol_dense (A) ; 264 end 265 if (burble) 266 fprintf ('OK.\n') ; 267 end 268catch me 269 if (burble) 270 fprintf ('failed.\nfactorize: %s\n', me.message) ; 271 end 272end 273 274 275%------------------------------------------------------------------------------- 276 277function [F, me] = factorize_ldl (A, burble) 278% LDL fails if the matrix is rectangular or rank-deficient. 279% As of MATLAB R2012a, ldl does not work for complex sparse matrices. 280% Only the lower triangular part of A is used. 281F = [ ] ; 282me = [ ] ; 283try 284 if (burble) 285 fprintf ('factorize: try LDL ... ') ; 286 end 287 if (issparse (A)) 288 F = factorization_ldl_sparse (A) ; 289 else 290 F = factorization_ldl_dense (A) ; 291 end 292 if (burble) 293 fprintf ('OK.\n') ; 294 end 295catch me 296 if (burble) 297 fprintf ('failed.\nfactorize: %s\n', me.message) ; 298 end 299end 300 301 302%------------------------------------------------------------------------------- 303 304function [F, me] = factorize_lu (A, burble, fail_if_singular) 305% LU fails if the matrix is rectangular or rank-deficient. 306F = [ ] ; 307me = [ ] ; 308try 309 if (burble) 310 fprintf ('factorize: try LU ... ') ; 311 end 312 if (issparse (A)) 313 F = factorization_lu_sparse (A, fail_if_singular) ; 314 else 315 F = factorization_lu_dense (A, fail_if_singular) ; 316 end 317 if (burble) 318 fprintf ('OK.\n') ; 319 end 320catch me 321 if (burble) 322 fprintf ('failed.\nfactorize: %s\n', me.message) ; 323 end 324end 325 326 327%------------------------------------------------------------------------------- 328 329function [F, me] = factorize_cod (A, burble) 330% COD only fails when it runs out of memory. 331F = [ ] ; 332me = [ ] ; 333try 334 if (burble) 335 fprintf ('factorize: try COD ... ') ; 336 end 337 if (issparse (A)) 338 F = factorization_cod_sparse (A) ; 339 else 340 F = factorization_cod_dense (A) ; 341 end 342 if (burble) 343 fprintf ('OK.\n') ; 344 end 345catch me 346 if (burble) 347 fprintf ('failed.\nfactorize: %s\n', me.message) ; 348 end 349end 350 351%------------------------------------------------------------------------------- 352 353function [F, me] = factorize_svd (A, burble) 354% SVD only fails when it runs out of memory. 355F = [ ] ; 356me = [ ] ; 357try 358 if (burble) 359 fprintf ('factorize: try SVD ... ') ; 360 end 361 F = factorization_svd (A) ; 362 if (burble) 363 fprintf ('OK.\n') ; 364 end 365catch me 366 if (burble) 367 fprintf ('failed.\nfactorize: %s\n', me.message) ; 368 end 369end 370