## On the meaning of the permutation pivot indices, ipiv

Open discussion regarding features, bugs, issues, vendors, etc.

### On the meaning of the permutation pivot indices, ipiv

I am a newbie with CLAPACK with no special education in Algebra and no knowledge of Fortran; with the purpose of interfacing the library with a high-level language, I am trying to understand the meaning of the ipiv vector I get, for example, out of dgetrf.

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)`
its meaning is:
• 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.
If I were to search directly for swappings involving the remaining rows (5 and 6), I would have to scan the ipiv elements in search of 5 and 6 and take their index:
• Row 5 must be swapped with row 1.
• Row 6 must be swapped with row 2.
Correct?

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)`
if 4 is not found, 4 must be poked at index 4.
• Scan ipiv in search of 5 and poke its position 1 at index 5:
Code: Select all
`(4 5 2 3 0 1)`
if 5 is not found, 5 must be poked at index 5.
• 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)`
obtaining:
Code: Select all
`    0 1 2 3 4 5   ------------0 | 0 0 0 0 1 01 | 0 0 0 0 0 12 | 0 0 1 0 0 03 | 0 0 0 1 0 04 | 1 0 0 0 0 05 | 0 1 0 0 0 0`
this is what GNU Octave gives me.
Yes?

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);intmain (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);}voidprint_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");    }}`
marcomaggi

Posts: 7
Joined: Mon Feb 08, 2010 6:14 am

### Re: On the meaning of the permutation pivot indices, ipiv

It looks like the above is wrong. Here is another take on it from the documentation of the package I am writing.

The pivot indices stored in the vector represent the sequence of row-permutations to be applied to a matrix.

Let's say we have factored a matrix A, with 6 rows and 4 columns, using the LU decomposition; we have A = PLU, where P is the permutation matrix. Let's say the permutations are represented by the following array of pivot indices:
Code: Select all
`(5 6 3 4)               ;; Fortran convention`
we have to understand that rows and columns are indexed starting from 1, following the Fortran convention: the pivot index 5 is associated to row 1, the pivot index 6 is associated to row 2, the pivot index 3 is associated to row 3, ...

The meaning of the array is that the following sequence of operations must be applied to the matrix:
• Swap row 1 with row 5.
• Swap row 2 with row 6.
• Swap row 3 with row 3, that is: leave it untouched.
• Swap row 4 with row 4, that is: leave it untouched.
Let's say instead that the matrix A is square with dimension 4 and we got the following array of pivot indices:
Code: Select all
`(3 2 3 4)               ;; Fortran convention`
notice that 3 appears twice; its meaning is that the following sequence of operations must be applied to the matrix:
• Swap row 1 with row 3.
• Swap row 2 with row 2, that is: leave it untouched.
• Swap row 3 with row 3, that is: leave it untouched; we have already swapped the third row with the first, so we need to do nothing now.
• Swap row 4 with row 4, that is: leave it untouched.
Now given the pivots for the 6-by-4 matrix:
Code: Select all
`(5 6 3 4)               ;; Fortran convention`
to obtain the permutation matrix P which left-multiplies LU using the zero-based C language indexing:
• Subtract 1 from all the indices in ipiv:
Code: Select all
`(4 5 2 3)               ;; C convention`
• 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)           ;; C convention`
• Initialise the added positions with the position index itself:
Code: Select all
`(4 5 2 3 4 5)           ;; C convention`
• Scan the whole ipiv using i as counter, if i = ipiv(i) enter the following subloop:
• Scan ipiv using the counter j = 0, ..., i-1.
• If i = ipiv(j) set ipiv(i) := j.
we understand that the array is mutated to:
Code: Select all
`(4 5 2 3 0 1)           ;; C convention`
• 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)   ;; C convention(3, 3)  (4, 0)  (5, 1)`

obtaining:
Code: Select all
`    0 1 2 3 4 5   ------------0 | 0 0 0 0 1 01 | 0 0 0 0 0 12 | 0 0 1 0 0 03 | 0 0 0 1 0 04 | 1 0 0 0 0 05 | 0 1 0 0 0 0`

Following the same procedure for the 4-by-4 matrix having pivots array:
Code: Select all
`(3 2 3 4)               ;; Fortran convention`

decrementing by 1:
Code: Select all
`(2 1 2 3)               ;; C convention`

mutating it:
Code: Select all
`(2 1 0 3)               ;; C convention`

we get the following permutation matrix:
Code: Select all
`    0 1 2 3   --------0 | 0 0 1 01 | 0 1 0 02 | 1 0 0 03 | 0 0 0 1`

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);intmain (void){  {    static const int ipiv = { 5, 6, 3, 4 };    print_perm(ipiv, 4, 6);  }  {    static const int ipiv = { 3, 2, 3, 4 };    print_perm(ipiv, 4, 4);  }  exit(EXIT_SUCCESS);}voidprint_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;      if (full_ipiv[i] == i) {        for (j=0; j<i; ++j)          if (i == full_ipiv[j])            {              full_ipiv[i] = j;              break;            }      }    }  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("\nFull 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");    }}`
marcomaggi

Posts: 7
Joined: Mon Feb 08, 2010 6:14 am

### Re: On the meaning of the permutation pivot indices, ipiv

You do not need two nested loops but only one loop.
Julien.
Julien Langou

Posts: 835
Joined: Thu Dec 09, 2004 12:32 pm
Location: Denver, CO, USA

### Re: On the meaning of the permutation pivot indices, ipiv

Pivot column matrices are most difficult procedure. But the formula for calculating is very useful. Follow the rows and colmn according to sequence make it very easy.
harop23

Posts: 1
Joined: Tue Jun 10, 2014 3:24 am

### Re: On the meaning of the permutation pivot indices, ipiv

In LAPACK there are two sorts of index arrays: IPIV (eg DGETRF) and JPVT (eg DGEQPF).
IPIV is a sequence of swaps and JPVT is a column permutation.

I find in coding for performance I want to go the opposite way, that is, given JPVT find the
appropriate sequence of swap (IPIV) that does the same job.
Why?
because the swaps can be done in-place
because the number of swaps can be very much less than N
I coded this operation up and use it (although the same operation probably exists squirrelled
away some where in a LAPACK auxiliary routine) and find significant savings on large matrices
by swapping rather than permuting.
lawrence mulholland

Posts: 25
Joined: Mon Jun 11, 2012 6:33 am
Location: NAG Ltd, Oxford, UK

### Re: On the meaning of the permutation pivot indices, ipiv

Hi all,
I m currently doing some C++ wrapping of Lapack and I want to be sure to understand well the meaning of ipiv and jpiv
lawrence mulholland wrote:In LAPACK there are two sorts of index arrays: IPIV (eg DGETRF) and JPVT (eg DGEQPF).
IPIV is a sequence of swaps and JPVT is a column permutation.

Can you confirm me, that:

IPIV is the list of "transpositions" (id the list of swap operations)
JPIV is a "compressed" representation of a "permutation" matrix (and as a consequence can not be applied "inplace")

Thank you! VPixor

Posts: 1
Joined: Thu Oct 08, 2015 2:25 am