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
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
|
function [L,U,prow,pcol] = superlu(A,psparse)
% SUPERLU : Supernodal LU factorization
%
% Executive summary:
%
% [L,U,p] = superlu(A) is like [L,U,P] = lu(A), but faster.
% [L,U,prow,pcol] = superlu(A) preorders the columns of A by min degree,
% yielding A(prow,pcol) = L*U.
%
% Details and options:
%
% With one input and two or three outputs, SUPERLU has the same effect as LU,
% except that the pivoting permutation is returned as a vector, not a matrix:
%
% [L,U,p] = superlu(A) returns unit lower triangular L, upper triangular U,
% and permutation vector p with A(p,:) = L*U.
% [L,U] = superlu(A) returns permuted triangular L and upper triangular U
% with A = L*U.
%
% With a second input, the columns of A are permuted before factoring:
%
% [L,U,prow] = superlu(A,psparse) returns triangular L and U and permutation
% prow with A(prow,psparse) = L*U.
% [L,U] = superlu(A,psparse) returns permuted triangular L and triangular U
% with A(:,psparse) = L*U.
% Here psparse will normally be a user-supplied permutation matrix or vector
% to be applied to the columns of A for sparsity. COLAMD is one way to get
% such a permutation; see below to make SUPERLU compute it automatically.
% (If psparse is a permutation matrix, the matrix factored is A*psparse'.)
%
% With a fourth output, a column permutation is computed and applied:
%
% [L,U,prow,pcol] = superlu(A,psparse) returns triangular L and U and
% permutations prow and pcol with A(prow,pcol) = L*U.
% Here psparse is a user-supplied column permutation for sparsity,
% and the matrix factored is A(:,psparse) (or A*psparse' if the
% input is a permutation matrix). Output pcol is a permutation
% that first performs psparse, then postorders the etree of the
% column intersection graph of A. The postorder does not affect
% sparsity, but makes supernodes in L consecutive.
% [L,U,prow,pcol] = superlu(A,0) is the same as ... = superlu(A,I); it does
% not permute for sparsity but it does postorder the etree.
% [L,U,prow,pcol] = superlu(A) is the same as ... = superlu(A,colamd(A));
% it uses column minimum degree to permute columns for sparsity,
% then postorders the etree and factors.
%
% This m-file calls the mex-file MEXSUPERLU to do the work.
%
% See also COLAMD, LUSOLVE.
%
% John Gilbert, 6 April 1995.
% Copyright (c) 1995 by Xerox Corporation. All rights reserved.
% HELP COPYRIGHT for complete copyright and licensing notice.
% Note on permutation matrices and vectors:
%
% Names beginning with p are permutation vectors;
% names beginning with P are the corresponding permutation matrices.
%
% Thus A(pfoo,pbar) is the same as Pfoo*A*Pbar'.
%
% We don't actually form any permutation matrices except Psparse,
% but matrix notation is easier to untangle in the comments.
[m,n] = size(A);
if m ~= n
error('matrix must be square.');
end;
if n == 0
L = []; U = []; prow = []; pcol = [];
return;
end;
% If necessary, compute the column sparsity permutation.
if nargin < 2
if nargout >= 4
psparse = colamd(A);
else
psparse = 1:n;
end;
end;
if max(size(psparse)) <= 1
psparse = 1:n;
end;
% Compute the permutation-matrix version of psparse,
% which fits the internal data structures of mexsuperlu.
if min(size(psparse)) == 1
Psparse = sparse(1:n,psparse,1,n,n);
else
Psparse = psparse;
psparse = Psparse*[1:n]';
end;
% Make sure the matrices are sparse.
if ~issparse(A)
A = sparse(A);
end;
if ~issparse(Psparse)
Psparse = sparse(Psparse);
end;
% The output permutations from the mex-file are dense permutation vectors.
[L,U,prowInv,pcolInv] = mexsuperlu(A,Psparse);
prow = zeros(1,n);
prow(prowInv) = 1:n;
pcol = zeros(1,n);
pcol(pcolInv) = 1:n;
% -- Added 12/14/2011 --
% The row indices of L & U matrices are not sorted from SuperLU, but
% Matlab requires the matrices to be sorted.
% Now, we do double-transpose to achieve sorting. This acts like bucket sort,
% should be faster than find / sparse combination, which probably uses
% quicksort.
% [i,j,s] = find(L); L = sparse(i,j,s,m,n);
L1 = L'; L = L1';
U1 = U'; U = U1';
% ----
% We now have
%
% Prow*A*Psparse'*Post' = L*U (1)
% Pcol' = Psparse'*Post'
%
% (though we actually have the vectors, not the matrices, and
% we haven't computed Post explicitly from Pcol and Psparse).
% Now we figure out what the user really wanted returned,
% and rewrite (1) accordingly.
if nargout == 4
% Return both row and column permutations.
% This is what we've already got.
% (1) becomes Prow*A*Pcol' = L*U.
elseif nargout == 3
% Return row permutation only. Fold the postorder perm
% but not the sparsity perm into it, and into L and U.
% This preserves triangularity of L and U.
% (1) becomes (Post'*Prow) * A * Psparse' = (Post'*L*Post) * (Post'*U*Post).
postInv = pcolInv(psparse);
prow = prow(postInv);
L = L(postInv,postInv);
U = U(postInv,postInv);
else % nargout <= 2
% Return no permutations. Fold the postorder perm but
% not the sparsity perm into L and U. Fold the pivoting
% perm into L, which destroys its triangularity.
% (1) becomes A * Psparse' = (Prow'*L*Post) * (Post'*U*Post).
postInv = pcolInv(psparse);
L = L(prowInv,postInv);
U = U(postInv,postInv);
end;
|
