On the meaning of the permutation pivot indices, ipiv

[marcomaggi]

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 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


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

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


[marcomaggi]

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 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


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 0
1 | 0 1 0 0
2 | 1 0 0 0
3 | 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);

int
main (void)
{
  {
    static const int ipiv[4] = { 5, 6, 3, 4 };
    print_perm(ipiv, 4, 6);
  }
  {
    static const int ipiv[4] = { 3, 2, 3, 4 };
    print_perm(ipiv, 4, 4);
  }
  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;
      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");
    }
}


[Julien Langou]

You do not need two nested loops but only one loop.
Please see: viewtopic.php?f=2&t=874
Julien.

[harop23]

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.

[lawrence mulholland]

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.

[VPixor]

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

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!