It is my understanding that rows and columns are indexed starting from 1, so when I get, after factoring a 6-rows by 4-columns matrix, the ipiv:

- Code: Select all
`(5 6 3 4)`

- Row 1 must be swapped with row 5.
- Row 2 must be swapped with row 6.
- Row 3 must be swapped with row 3, that is it has to be left untouched.
- Row 4 must be swapped with row 4, that is it has to be left untouched.

- Row 5 must be swapped with row 1.
- Row 6 must be swapped with row 2.

Now to obtain the permutation matrix P which left multiplies LU using the zero-based C language indexing:

- Subtract 1 from all the indices:
- Code: Select all
`(4 5 2 3)`

- Reallocate ipiv so that it has a number of elements equal to the number of rows in the original A matrix:
- Code: Select all
`(4 5 2 3 x x)`

- Scan ipiv in search of 4 and poke its position 0 at index 4:
- Code: Select all
`(4 5 2 3 0 x)`

- Scan ipiv in search of 5 and poke its position 1 at index 5:
- Code: Select all
`(4 5 2 3 0 1)`

- Build a zero-filled matrix 6-by-6 (the dimension must be equal to the number of rows of the original matrix A), and poke 1 to the positions:
- Code: Select all
`(0, 4) (1, 5) (2, 2) (3, 3) (4, 0) (5, 1)`

- Code: Select all
`0 1 2 3 4 5`

------------

0 | 0 0 0 0 1 0

1 | 0 0 0 0 0 1

2 | 0 0 1 0 0 0

3 | 0 0 0 1 0 0

4 | 1 0 0 0 0 0

5 | 0 1 0 0 0 0

The following is a C program printing the permutation matrix from a given ipiv:

- Code: Select all
`#include <stdio.h>`

#include <stdlib.h>

#include <string.h>

static void print_perm (const int * ipiv, size_t number_of_indices, int number_of_rows);

int

main (void)

{

#define NUMBER_OF_INDICES 4

#define NUMBER_OF_ROWS 6

static const int ipiv[NUMBER_OF_INDICES] = { 5, 6, 3, 4 };

print_perm(ipiv, NUMBER_OF_INDICES, NUMBER_OF_ROWS);

exit(EXIT_SUCCESS);

}

void

print_perm (const int * ipiv, size_t number_of_indices, int number_of_rows)

{

int full_ipiv[number_of_rows];

int P[number_of_rows][number_of_rows];

size_t i, j;

for (i=0; i<number_of_indices; ++i)

full_ipiv[i] = ipiv[i] - 1;

for (; i<number_of_rows; ++i)

{

full_ipiv[i] = i;

for (j=0; j<number_of_indices; ++j)

if (i == full_ipiv[j])

{

full_ipiv[i] = j;

break;

}

}

memset(&P, '\0', number_of_rows * number_of_rows * sizeof(int));

for (i=0; i<number_of_rows; ++i)

P[i][full_ipiv[i]] = 1;

printf("Full ipiv: ");

for (i=0; i<number_of_rows; ++i)

printf("%d ", full_ipiv[i]);

printf("\n\n");

printf("Permutation matrix:\n");

for (i=0; i<number_of_rows; ++i)

{

for (j=0; j<number_of_rows; ++j)

printf("%d ", P[i][j]);

printf("\n");

}

}