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");
}
}