Strange bug in LUP Decomposition

I have written LU decomposition with partial pivot, but for some matrices, the entry L[[1,1]] isn't equal to 1 like it should be, throwing off the entire decomposition. Any ideas what is causing this?

(* LU DECOMP w/ pivot *)

Clear[U, L, A, j, k, a, b]
lupdecomp[A_?MatrixQ] :=
Module[{P, max, index, x, i, n, y, dim, d, temp, j, k, a, U, L},
dim = Dimensions[A];
n = Length[A];
U = A;
L = IdentityMatrix[n];
P = IdentityMatrix[n];

For[k = 1, k < n, k++,
max = Max[U[[k ;;, k]]];
index = Flatten[Position[U[[k ;;, k]], max]][[1]];
U[[{k, index}]] = U[[{index, k}]];
L[[{k, index}, 1 ;; k - 1]] = L[[{index, k}, 1 ;; k - 1]];
P[[{k, index}]] = P[[{index, k}]];
For[j = k + 1, j <=  n, j++,
L[[j, k]] = U[[j, k]]/U[[k, k]];
For[a = k, a <= n, a++,
U[[j, a]] = U[[j, a]] - L[[j, k]]*U[[k, a]];
];

];
];
Return[{P, L, U}];
];


For example,

b= {{4.40669, 3.43566, 9.61987}, {3.13933, 7.78681, 8.60175}, {9.29858, 4.73565, 3.41126}}

Yields

{P,L,U}={{{0, 1, 0}, {0, 0, 1}, {1, 0, 0}}, {{0.337614, 0, 0}, {1, 1, 0}, {0.47391, 0.251578, 1}}, {{4.44089*^(-16), 6.18799, 7.45006}, {9.29858, 4.73565, 3.41126}, {0., 0., 7.14504}}}

Which has an entry L[[1,1]] not equal to 1. However, for other matrices, such as

b= {{7.01497, 9.75886, 4.12399}, {3.91071, 2.29947, 6.53455}, {8.28621, 7.69486, 1.44275}}

{P,L,U}={{{0, 0, 1}, {0, 1, 0}, {1, 0, 0}}, {{1, 0, 0}, {0.471954, 1, 0}, {0.846584, -2.43555, 1}}, {{8.28621, 7.69486, 1.44275}, {0., -1.33215, 5.85364}, {-8.88178*^(-16), 0., 17.1594}}}

Which is a proper decomposition, as LU=Pb

Any ideas what I"m doing wrong?

The problem is in these two lines:

max = Max[U[[k ;;, k]]];
index = Flatten[Position[U[[k ;;, k]], max]][[1]];


You are looking for the pivoting row, so you have to add k-1 to index. In fact, you can have this easier by loading Needs["LinearAlgebraBLAS"] and by using

index = IAMAX[U[[k ;;, k]]] + (k - 1);


Here is also a somewhat more optimized version of the algorithm. It (i) gets rid of two For loops in favor of fast matrix operations and (ii) ensures that packed arrays are used if A is a matrix of machine precision reals.

Needs["LinearAlgebraBLAS"]

lupdecomp[A_?MatrixQ] := Module[{n, P, U, L, index, k},
n = Length[A];
U = DeveloperToPackedArray[A];
L = IdentityMatrix[n, WorkingPrecision -> Precision[A]];
P = IdentityMatrix[n, WorkingPrecision -> Precision[A]];
Do[
index = IAMAX[U[[k ;;, k]]] + (k - 1);
U[[{k, index}]] = U[[{index, k}]];
L[[{k, index}, 1 ;; k - 1]] = L[[{index, k}, 1 ;; k - 1]];
P[[{k, index}]] = P[[{index, k}]];
L[[k + 1 ;;, k]] = U[[k + 1 ;;, k]]/U[[k, k]];
U[[k + 1 ;;, k ;;]] = U[[k + 1 ;;, k ;;]] - KroneckerProduct[L[[k + 1 ;;, k]], U[[k, k ;;]]];
, {k, 1, n}];
{P, L, U}
];


Here is a test:

n = 15;
A = RandomReal[{-1, 1}, {n, n}];
{P, L, U} = lupdecomp[A];
Max[Abs[P\[Transpose].L.U - A]]


4.44089*10^-16

• This is because 1) we start the indexing with 1 and 2) we are finding the position of the max element in terms of numbers after k, so to get the index from the start of the matrix, we have to add (k-1). If the indexing started at 0, we would add k, correct? Commented Apr 4, 2019 at 20:16
• Right. If you could start indexing with 0 which you cannot in Mathematica. Commented Apr 4, 2019 at 20:23
• I was thinking more along the lines of "If i had to write this in c++, it'd be +k and not +(k-1)". If i remember correctly C++ indexing begins with 0`. I can't believe I didn't spot that, thanks as always @Henrik! Commented Apr 4, 2019 at 20:25
• You're welcome. Commented Apr 4, 2019 at 20:26