croutp#

Purpose#

Computes the Crout decomposition of a square matrix with partial (row) pivoting.

Format#

lup = croutp(x)#

Parameters:: x (NxN matrix) – square nonsingular matrix.
Returns:: lup ((N+1)xN matrix) – containing the lower ( $L$ ) and upper ( $U$ ) matrices of the Crout decomposition of a permuted x. The $N + 1$ row of the matrix lup gives the row order of the lup matrix. The matrix must be reordered prior to extracting the $L$ and $U$ matrices. Use lowmat() and upmat1() to extract the $L$ and $U$ matrices from the reordered lup matrix.

Examples#

This example illustrates a procedure for extracting $L$ and $U$ of the permuted x matrix. It continues by sorting the result of $L U$ to compare with the original matrix x.

X = { 1 2 -1,
      2 3 -2,
      1 -2 1 };

// Perform crout decomposition
lup = croutp(x);

If we view lup, we will see:

0000       0.50000       0.28571
lup = 2.0000        1.5000       -1.0000
0000       -3.5000      -0.57142
0000        3.0000        1.0000

/*
** This bottom row is the permutation index vector
** Calculate how many rows in 'lup'
*/
r = rows(lup);

/*
** Extract the index row and transpose it into a column
** vector
*/
index = lup[r, .]';

Viewing index will reveal:

        2
index = 3
        1

// Rearrange the rows of 'lup' based upon the index vector
z = lup[index, .];

// obtain L and U of permuted matrix X
L = lowmat(z);
U = upmat1(z);

/*
** Horizontally concatenate the index vector and the product
** of L*U then pass that result into the 'sortc' function
** which will sort this result based upon the first column
** (which is the index vector)
*/
q = sortc(index~(L*U), 1);

/*
** Remove the index vector, which we added by way of
** horizontal concatenation in the statement just above
*/
x2 = q[., 2:cols(q)];

Now at the end of this example, x2 is equal to x.