# 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 *i*th 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) 1n←rows[A] 2fori← 1ton3doπ[i] ←i4fork← 1ton5dop← 0 6fori←kton7do if|a| >_{ik}p8thenp← |a| 9_{ik}k' ←i10ifp= 0 11then error"singular matrix" 12 exchangeπ[k] ↔π[k'] 13fori← 1ton14doexchangea↔_{ki}a15_{k'i}fori←k1ton16doa←_{ik}a/_{ik}a17_{kk}forj←k1ton18doa←_{ij}a-_{ij}a_{ik}a_{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

*k*th and

*k*'th rows of

*A*in lines 13-14, thereby making the pivot element

*a*. (The entire rows are swapped because in the derivation of the method above, not only is

_{kk}*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."

*(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*^{3}), which is the same as that of LU-DECOMPOSITION. Thus, pivoting costs us at most a constant factor in time.