Solving Systems Of Linear Equations

which is a permutation matrix, since it is the product of two permutation matrices. We now have

yielding the LUP decomposition. Because L′ is unit lower-triangular, so is L, and because U′ is upper-triangular, so is U.

Notice that in this derivation, unlike the one for LU decomposition, both the column vector v/a_k1 and the Schur complement A′ - vw^T/a_k1 must be multiplied by the permutation matrix P′.

Like LU-DECOMPOSITION, our pseudocode for LUP decomposition replaces the recursion with an iteration loop. As an improvement over a direct implementation of the recursion, we dynamically maintain the permutation matrix P as an array π, where π[i] = j means that the ith row of P contains a 1 in column j. We also implement the code to compute L and U "in place" in the matrix A. Thus, when the procedure terminates,

	LUP-DECOMPOSITION(A)
 1 n ← rows[A]
 2 for i ← 1 to n
 3      do π[i] ← i
 4 for k ← 1 to n
 5      do p ← 0
 6         for i ← k to n
 7             do if |a_ik| > p
 8                   then p ← |a_ik|
 9                        k' ← i
10         if p = 0
11            then error "singular matrix"
12         exchange π[k] ↔ π[k']
13         for i ← 1 to n
14              do exchange a_ki ↔ a_k'i
15         for i ← k   1 to n
16              do a_ik ← a_ik/a_kk
17                 for j ← k   1 to n
18                     do a_ij ← a_ij - a_ika_kj

Figure 28.2 illustrates how LUP-DECOMPOSITION factors a matrix. The array π is initialized by lines 2-3 to represent the identity permutation. The outer for loop beginning in line 4 implements the recursion. Each time through the outer loop, lines 5-9 determine the element a_k'k with largest absolute value of those in the current first column (column k) of the (n - k 1) × (n - k 1) matrix whose LU decomposition must be found. If all elements in the current first column are zero, lines 10-11 report that the matrix is singular. To pivot, we exchange π[k'] with π[k] in line 12 and exchange the kth and k'th rows of A in lines 13-14, thereby making the pivot element a_kk. (The entire rows are swapped because in the derivation of the method above, not only is A′ - vw^T/a_k1 multiplied by P′, but so is v/a_k1.) Finally, the Schur complement is computed by lines 15-18 in much the same way as it is computed by lines 4-9 of LU-DECOMPOSITION, except that here the operation is written to work "in place."

Figure 28.2: The operation of LUP-DECOMPOSITION. (a) The input matrix A with the identity permutation of the rows on the left. The first step of the algorithm determines that the element 5 in the black circle in the third row is the pivot for the first column. (b) Rows 1 and 3 are swapped and the permutation is updated. The shaded column and row represent v and w^T. (c) The vector v is replaced by v/5, and the lower right of the matrix is updated with the Schur complement. Lines divide the matrix into three regions: elements of U (above), elements of L (left), and elements of the Schur complement (lower right). (d)-(f) The second step. (g)-(i) The third step. No further changes occur on the fourth, and final, step. (j) The LUP decomposition PA = LU.

Because of its triply nested loop structure, LUP-DECOMPOSITION has a running time of Θ(n³), which is the same as that of LU-DECOMPOSITION. Thus, pivoting costs us at most a constant factor in time.