At an intermediate stage of LU factorization, the 2D array A that originally contained the matrix A has a mix of L, U, and partially processed parts of A. Because of that, mixing the terminology of array versus matrix is almost unavoidable. At the start, array A contains the matrix A. At the end, it is a data structure representing the matrices L and U. All LU factorization algorithms for dense matrices work in-place, so in general assume that A refers to the array.

Another common nomenclature: at an intermediate stage of LU factorization, the parts of the L and U factors that have already been computed are called the reduced part. The rest of the array, usually its trailing principle submatrix, is called the unreduced part.

When using any of the LU factorization methods (rank-1 update or matrix-vector product) of LU factorization, dividing by A(k,k) in the scaling operation will cause problems if the absolute value |A(k,k)| is small. To avoid this, algorithms search for a large value in the unreduced part of the array A, then swap rows and sometimes columns to move it to the (k,k) diagonal position.

Full pivoting means bringing the largest element in absolute value into that position by swapping rows and columns in the matrix to put the largest element in the (k,k) diagonal position. Partial pivoting means searching for the largest only among the entries in the k-th column, and so only require row swaps. Full pivoting is rarely used because it requires O(n³) searching and data movement operations¹. Also, LU factorization is "usually stable" when partial pivoting is used. Some classes of problems (e.g., two point boundary value problems in ODEs) can make it fail, but most often partial pivoting suffices. In general when the term "pivoting" is used, it refers to partial pivoting, not full pivoting. After pivoting is done to move a larger (or largest) value into the (k,k) diagonal position, it is called the "pivot" or "pivot entry" or "pivot element". In any case it is the diagonal element A(k,k), chosen in an attempt to avoid dividing by a too-small value when scaling the subdiagonal part of column k.

for k = 1:n
  p = index of max element in |A(k:n,k)|
  piv(k) = p
  swap rows k and p of A: A(k,:) ↔ A(p,:)
  A(k+1:n,k) = A(k+1:n,k)/A(k,k)
  A(k+1:n,k+1:n) = A(k+1:n,k+1:n) − A(k+1:n,k) * A(k,k+1:n)
end for 

Obviously, the three lines for implementing pivoting in the algorithms above are simple, and simple to implement. Here are some tips:

• If k = p, don't call a function to swap a row with itself (well, duh!). Since swapping rows constitutes data movement without any computation occurring, it has a load/store ratio of ∞. This is the worst possible ratio of memory references to flops, so should be avoided when possible.

• The goal of pivoting is not to get the largest possible element into the pivot position, but to avoid dividing by too-small of an entry. So the choice of the pivot element can be relaxed. Markovitz pivoting with parameter s searches for the subdiagonal element with largest absolute value, but only swaps the rows if |A(k,p)| > s*|A(k,k)|. The value s = 10 is a common choice. In other words, if the current diagonal element A(k,k) is within a factor of 10 of the largest possible pivot, use it as the pivot.

• Always choosing the largest subdiagonal element as the pivot can be seen as setting the "Markovitz parameter" s = 1. Never pivoting can be seen as setting s = ∞. This is a case where actually using one of the IEEE special floating point values can make the coding simpler.

• Choosing the Markovitz parameter s too large can lead to numerical instability, which is what pivoting is supposed to avoid. Part of the issue is just how much time is spent in swapping rows for a give value of s. At this point, you have the skills needed to determine that ... and it would be foolish not to insert the necessary timing calls in advance.

• Pivoting does mean that LU factorization is not completely in-place, using only the original array A for everything. However, it introduces only an integer vector of length n and so is minor compared to copying or replicating the original data for A which requires n² words of storage.

• The BLAS level-1 function dswap() will swap two vectors in-place, and was in large part defined to be used for LU factorization pivoting. This also explains why dswap has the signature dswap(n, x, incx, y, incy): the entries of the vectors x and y are rows of the 2D array A. When A is in column-major order, the strides incx and incy are greater than 1. What should those strides be?

• The Matlab idiom for swapping rows k and p in an array A is A([k p], :) = A([p k], :) which is more straightforward than the usual loop with a temporary approach required by C/C++ or Java. Just beware when using Fortran: it requires commas separate the two row indices: A([k, p], :) = A([p, k], :) . Using brackets [k, p] was introduced in Fortran 2008; before that the notation had to use (\k, p \) which is still valid but just plain ugly.