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