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

Postby marcomaggi » Wed Feb 10, 2010 4:19 am

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
 
Posts: 6
Joined: Mon Feb 08, 2010 6:14 am

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

Postby marcomaggi » Wed Feb 10, 2010 7:24 am

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");
    }
}
marcomaggi
 
Posts: 6
Joined: Mon Feb 08, 2010 6:14 am

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

Postby Julien Langou » Wed Feb 10, 2010 9:23 pm

You do not need two nested loops but only one loop.
Please see: viewtopic.php?f=2&t=874
Julien.
Julien Langou
 
Posts: 747
Joined: Thu Dec 09, 2004 12:32 pm
Location: Denver, CO, USA

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

Postby harop23 » Wed Jun 11, 2014 1:36 am

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

Postby lawrence mulholland » Wed Jul 02, 2014 7:55 am

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: 18
Joined: Mon Jun 11, 2012 6:33 am
Location: NAG Ltd, Oxford, UK


Return to User Discussion

Who is online

Users browsing this forum: Yahoo [Bot] and 2 guests