1
function [result, t] = linfactor (arg1, arg2)
2
%LINFACTOR factorize a matrix, or use the factors to solve Ax=b.
3
% Uses LU or CHOL to factorize A, or uses a previously computed factorization to
4
% solve a linear system. This function automatically selects an LU or Cholesky
5
% factorization, depending on the matrix. A better method would be for you to
6
% select it yourself. Note that mldivide uses a faster method for detecting
7
% whether or not A is a candidate for sparse Cholesky factorization (see spsym
8
% in the CHOLMOD package, for example).
11
% F = linfactor (A) ; % factorizes A into the object F
12
% x = linfactor (F,b) ; % uses F to solve Ax=b
15
% A second output is the time taken by the method, ignoring the overhead of
16
% determining which method to use. This makes for a fairer comparison between
17
% methods, since normally the user will know if the matrix is supposed to be
18
% symmetric positive definite or not, and whether or not the matrix is sparse.
19
% Also, the overhead here is much higher than mldivide or spsym.
21
% This function has its limitations:
23
% (1) determining whether or not the matrix is symmetric via nnz(A-A') is slow.
24
% mldivide (and spsym in CHOLMOD) do it much faster.
26
% (2) MATLAB really needs a sparse linsolve. See cs_lsolve, cs_ltsolve, and
27
% cs_usolve in CSparse, for example.
29
% (3) this function really needs to be written as a mexFunction.
31
% (4) the full power of mldivide is not brought to bear. For example, UMFPACK
32
% is not very fast for sparse tridiagonal matrices. It's about a factor of
33
% four slower than a specialized tridiagonal solver as used in mldivide.
35
% (5) permuting a sparse vector or matrix is slower in MATLAB than it should be;
36
% a built-in linfactor would reduce this overhead.
38
% (6) mldivide when using UMFPACK uses relaxed partial pivoting and then
39
% iterative refinement. This leads to sparser LU factors, and typically
40
% accurate results. linfactor uses sparse LU without iterative refinement.
42
% The primary purpose of this function is to answer The Perennially Asked
43
% Question (or The PAQ for short (*)): "Why not use x=inv(A)*b to solve Ax=b?
44
% How do I use LU or CHOL to solve Ax=b?" The full answer is below. The short
45
% answer to The PAQ (*) is "PAQ=LU ... ;-) ... never EVER use inv(A) to solve
48
% The secondary purpose of this function is to provide a prototype for some of
49
% the functionality of a true MATLAB built-in linfactor function.
51
% Finally, the third purpose of this function is that you might find it actually
52
% useful for production use, since its syntax is simpler than factorizing the
53
% matrix yourself and then using the factors to solve the system.
55
% See also lu, chol, mldivide, linsolve, umfpack, cholmod.
57
% Oh, did I tell you never to use inv(A) to solve Ax=b?
59
% Requires MATLAB 7.3 (R2006b) or later.
61
% Copyright 2007, Timothy A. Davis, University of Florida
62
% VERSION 1.1.0, Nov 1, 2007
64
if (nargin < 1 | nargin > 2 | nargout > 2) %#ok
65
error ('Usage: F=linfactor(A) or x=linfactor(F,b)') ;
70
%---------------------------------------------------------------------------
72
%---------------------------------------------------------------------------
77
error ('linfactor: A must be square') ;
82
% try sparse Cholesky (CHOLMOD): L*L' = P*A*P'
83
if (nnz (A-A') == 0 & all (diag (A) > 0)) %#ok
86
[L, g, PT] = chol (A, 'lower') ;
90
result.LT = L' ; % takes more memory, but solve is faster
91
result.P = PT' ; % ditto. Need a sparse linsolve here...
93
result.kind = 'sparse Cholesky: L*L'' = P*A*P''' ;
98
% matrix is symmetric, but not positive definite
99
% (or we ran out of memory)
103
% try sparse LU (UMFPACK, with row scaling): L*U = P*(R\A)*Q
105
[L, U, P, Q, R] = lu (A) ;
112
result.kind = 'sparse LU: L*U = P*(R\A)*Q where R is diagonal' ;
117
% try dense Cholesky (LAPACK): L*L' = A
118
if (nnz (A-A') == 0 & all (diag (A) > 0)) %#ok
121
L = chol (A, 'lower') ;
124
result.kind = 'dense Cholesky: L*L'' = A' ;
128
% matrix is symmetric, but not positive definite
129
% (or we ran out of memory)
133
% try dense LU (LAPACK): L*U = A(p,:)
135
[L, U, p] = lu (A, 'vector') ;
140
result.kind = 'dense LU: L*U = A(p,:)' ;
147
%---------------------------------------------------------------------------
148
% x = linfactor (F,b)
149
%---------------------------------------------------------------------------
156
% sparse Cholesky: MATLAB could use a sparse linsolve here ...
158
result = F.PT * (F.LT \ (F.L \ (F.P * b))) ;
163
% sparse LU: MATLAB could use a sparse linsolve here too ...
165
result = F.Q * (F.U \ (F.L \ (F.P * (F.R \ b)))) ;
170
% dense Cholesky: result = F.L' \ (F.L \ b) ;
173
upper.TRANSA = true ;
175
result = linsolve (F.L, linsolve (F.L, b, lower), upper) ;
180
% dense LU: result = F.U \ (F.L \ b (F.p,:)) ;
184
result = linsolve (F.U, linsolve (F.L, b (F.p,:), lower), upper) ;
added
removed
LINFACTOR/linfactor.m
