1function [L,U,p] = LUModSimple( L, U, p, u, v, tau )
2%
3% Copyright (c) 2009-2014, Jack Poulson
4% All rights reserved.
5%
6% This file is part of Elemental and is under the BSD 2-Clause License,
7% which can be found in the LICENSE file in the root directory, or at
8% http://opensource.org/licenses/BSD-2-Clause
9%
10% Begin with an LU factorization with partial pivoting,
11%     A = P^T L U,
12% and turn it into a partially-pivoted LU factorization of
13%     A + u v',
14% say
15%     (A + u v') = P^T L ( U + w v'),
16% w = inv(L) P u.
17%
18% Please see subsection 2.1 from
19%     Peter Stange, Andreas Griewank, and Matthias Bollhofer,
20%     "On the efficient update of rectangular LU factorizations subject to
21%      low rank modifications"
22% which discusses the technique of Schwetlick and Kielbasinski described in
23% "Numerische Lineare Algebra".
24%
25[m,n] = size(L);
26minDim = min(m,n);
27if minDim ~= m,
28  error('It is assumed that height(L) <= width(L)');
29end
30
31% w := inv(L) P u
32w=u(p);
33w=L\w;
34
35% Maintain an external vector for the temporary subdiagonal of U
36uSub=zeros(minDim-1,1);
37
38% Reduce w to a multiple of e0
39for i=minDim-1:-1:1,
40  % Decide if we should pivot the i'th and i+1'th rows of w
41  rightTerm = abs(L(i+1,i)*w(i)+w(i+1));
42  if abs(w(i)) < tau*rightTerm,
43    % P := P_i P
44    p([i,i+1])=p([i+1,i]);
45
46    % Simultaneously perform
47    %   U := P_i U and
48    %   L := P_i L P_i^T
49    U([i,i+1],:) = U([i+1,i],:);
50    L([i,i+1],:) = L([i+1,i],:);
51    L(:,[i,i+1]) = L(:,[i+1,i]);
52
53    % Update
54    %     L := L T_{i,L}^{-1},
55    %     U := T_{i,L} U,
56    %     w := T_{i,L} P_i w,
57    % where T_{i,L} is the Gauss transform which zeros (P_i w)_{i+1}.
58    %
59    % More succinctly,
60    %     gamma    := w(i) / w(i+1),
61    %     L(:,i)   += gamma L(:,i+1),
62    %     U(i+1,:) -= gamma U(i,:),
63    %     w(i)     := w(i+1),
64    %     w(i+1)   := 0.
65    gamma = w(i) / w(i+1);
66    L(:,i) = L(:,i) + gamma*L(:,i+1);
67    U(i+1,:) = U(i+1,:) - gamma*U(i,:);
68    w(i) = w(i+1);
69    w(i+1) = 0;
70
71    % Force L back to *unit* lower-triangular form via the transform
72    %     L := L T_{i,U}^{-1} D^{-1},
73    % where D is diagonal and responsible for forcing L(i,i) and
74    % L(i+1,i+1) back to 1. The effect on L is:
75    %     eta       := L(i,i+1)/L(i,i),
76    %     L(:,i+1)  -= eta L(:,i),
77    %     delta_i   := L(i,i),
78    %     delta_ip1 := L(i+1,i+1),
79    %     L(:,i)   /= delta_i,
80    %     L(:,i+1) /= delta_ip1,
81    % while the effect on U is
82    %     U(i,:)   += eta U(i+1,:)
83    %     U(i,:)   *= delta_i,
84    %     U(i+1,:) *= delta_{i+1},
85    % while the effect on w is
86    %     w(i) *= delta_i w(i).
87    eta = L(i,i+1)/L(i,i);
88    L(:,i+1) = L(:,i+1) - eta*L(:,i);
89    delta_i   = L(i,i);
90    delta_ip1 = L(i+1,i+1);
91
92    L(:,i  ) = L(:,i  ) / delta_i;
93    L(:,i+1) = L(:,i+1) / delta_ip1;
94
95    U(i,:) = U(i,:) + eta*U(i+1,:);
96    U(i,:) = delta_i*U(i,:);
97    U(i+1,:) = delta_ip1*U(i+1,:);
98    w(i) = delta_i*w(i);
99  else
100    % Update
101    %     L := L T_{i,L}^{-1},
102    %     U := T_{i,L} U,
103    %     w := T_{i,L} w,
104    % where T_{i,L} is the Gauss transform which zeros w_{i+1}.
105    %
106    % More succinctly,
107    %     gamma    := w(i+1) / w(i),
108    %     L(:,i)   += gamma L(:,i+1),
109    %     U(i+1,:) -= gamma U(i,:),
110    %     w(i+1)   := 0.
111    gamma = w(i+1) / w(i);
112    L(:,i) = L(:,i) + gamma*L(:,i+1);
113    U(i+1,:) = U(i+1,:) - gamma*U(i,:);
114    w(i+1) = 0;
115  end
116end
117
118% Add the modified w v' into U
119U(1,:) = U(1,:) + w(1)*v';
120
121% Transform U from upper-Hessenberg to upper-triangular form
122for i=1:minDim-1,
123  % Decide if we should pivot the i'th and i+1'th rows U
124  rightTerm = abs(L(i+1,i)*U(i,i)+U(i+1,i));
125  if abs(U(i,i)) < tau*rightTerm,
126    % P := P_i P
127    p([i,i+1])=p([i+1,i]);
128
129    % Simultaneously perform
130    %   U := P_i U and
131    %   L := P_i L P_i^T
132    U([i,i+1],:) = U([i+1,i],:);
133    L([i,i+1],:) = L([i+1,i],:);
134    L(:,[i,i+1]) = L(:,[i+1,i]);
135
136    % Update
137    %     L := L T_{i,L}^{-1},
138    %     U := T_{i,L} U,
139    % where T_{i,L} is the Gauss transform which zeros U(i+1,i).
140    %
141    % More succinctly,
142    %     gamma    := U(i+1,i) / U(i,i),
143    %     L(:,i)   += gamma L(:,i+1),
144    %     U(i+1,:) -= gamma U(i,:),
145    gamma = U(i+1,i) / U(i,i);
146    L(:,i) = L(:,i) + gamma*L(:,i+1);
147    U(i+1,:) = U(i+1,:) - gamma*U(i,:);
148
149    % Force L back to *unit* lower-triangular form via the transform
150    %     L := L T_{i,U}^{-1} D^{-1},
151    % where D is diagonal and responsible for forcing L(i,i) and
152    % L(i+1,i+1) back to 1. The effect on L is:
153    %     eta       := L(i,i+1)/L(i,i),
154    %     L(:,i+1)  -= eta L(:,i),
155    %     delta_i   := L(i,i),
156    %     delta_ip1 := L(i+1,i+1),
157    %     L(:,i)   /= delta_i,
158    %     L(:,i+1) /= delta_ip1,
159    % while the effect on U is
160    %     U(i,:)   += eta U(i+1,:)
161    %     U(i,:)   *= delta(i),
162    %     U(i+1,:) *= delta(i+1).
163    eta = L(i,i+1)/L(i,i);
164    L(:,i+1) = L(:,i+1) - eta*L(:,i);
165    delta_i = L(i,i);
166    delta_ip1 = L(i+1,i+1);
167    L(:,i) = L(:,i) / delta_i;
168    L(:,i+1) = L(:,i+1) / delta_ip1;
169    U(i,:) = delta_i*(U(i,:) + eta*U(i+1,:));
170    U(i+1,:) = delta_ip1*U(i+1,:);
171  else
172    % Update
173    %     L := L T_{i,L}^{-1},
174    %     U := T_{i,L} U,
175    % where T_{i,L} is the Gauss transform which zeros U(i+1,i).
176    %
177    % More succinctly,
178    %     gamma    := U(i+1,i) / U(i,i),
179    %     L(:,i)   += gamma L(:,i+1),
180    %     U(i+1,:) -= gamma U(i,:).
181    gamma = U(i+1,i) / U(i,i);
182    L(:,i) = L(:,i) + gamma*L(:,i+1);
183    U(i+1,:) = U(i+1,:) - gamma*U(i,:);
184  end
185end
186