Is dlaswp functioning properly?

[lapakost]

Hey,

I'm wondering how to interpret the results produced by Lapack's dlaswp routine dlaswp.f, which implements pivoting of array elements (or, more generally, rows of a matrix) using the supplied pivot array.

When I send the routine a simple pivot array of [2, 1], I'm expecting the routine to swap the elements of the two-element data array, which send along. But this is not what the routing is doing: it in fact leaves the elements in their original places, or, rather, swaps them twice. This behavior seems evident from the code of the routine too (dlaswp.f v 3.12.0, lines 171-183).

Below is my simple test case that demonstrates this behavior on I64 Linux.

Thank you in advance for your comments!

Output:

original x:	 0.100000	 0.200000
pivotes:	 2	 1
pivoted x:	 0.100000	 0.200000

Compile: gcc tst.c -llapack -o tst

#include <stdio.h>
#include <stdlib.h>

extern int
dlaswp_(const int* N, double* A, const int* LDA,
	const int* K1, const int* K2, const int* IPIV, const int* INCX );

int main() {

  double x[] = {.1, .2};
  const int piv[] = {2, 1};

  static const int N = sizeof(x) / sizeof(x[0]);

  if (sizeof(x) / sizeof(x[0]) != sizeof(piv) / sizeof(piv[0])) {
    printf("Error: lengths of x and piv arrays mismatch, exiting.\n");
    exit(1);
  }
  
  printf("original x:");
  for (int i = 0; i < N; ++ i)
    printf("\t %f", x[i]);
  printf("\n");

  printf("pivotes:");
  for (int i = 0; i < N; ++ i)
    printf("\t %i", piv[i]);
  printf("\n");

  static const int ONE = 1;

  dlaswp_(&ONE, x, &N, &ONE, &N, piv, &ONE);

  printf("pivoted x:");
  for (int i = 0; i < N; ++ i)
    printf("\t %f", x[i]);
  printf("\n");

  return 0;
}

[langou]

The representation of the permutation vector in LAPACK LASWP is not the one you expect. The vector p=[ 2, 1 ] means that you will do two elementary permutations. First swap (1,2) and then swap (2,1). (So yes, in this case, at the end, you end up as you started. So the output of your test example is correct.) If you want to swap row 1 and row 2, you can use the vector p=[2,2].

Maybe using an "N" larger than 2 would make it easier to understand. Try N=9 and p=[3,2,4,9,2,7,8,8,3] then this'll apply the nine elementary permutations (1,3), then (2,2), then (3,4), then (4,9), then (5,2), then (6,7), then (7,8), then (8,8), then (9,3). (In this order if INCX = 1.) So that an array [0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9] will be permuted to [0.3 0.5 0.1 0.9 0.2 0.7 0.8 0.6 0.4].

[lapakost]

Thank you a lot for a prompt and detailed response! It would be nice to add something like that to the doc for this routine.

I stumbled upon this behavior when I tried to use dlaswp.f for pivoting, following semi-definite Cholesky factorization by dpstrf.f. The latter routine does populate the pivot array with [2,1], when pivoting between the two only rows and columns of the input matrix needs to be made. Hence my question: what is the recommended way of pivoting the right-hand side of a system of equations and its solution, following a semi-definite Cholesky factorization of its matrix by dpstrf()?

Below is my test case that produces [2,1] pivoting in the trivial case of a diagonal matrix with [1, 4] on the diagonal.

Thanks again!

Output:

info: 0, rank: 2
pivotes:	 2	 1
a: 
	 2	 0
	 0	 1

Program:

#include <stdio.h>
#include <stdlib.h>

extern int
dpstrf_(const char* uplo,
	const int* n,
	double* a,
	const int* lda,
	int* piv,
	int* rank,
	const double* tol,
	double* work,
	int* info,
	const int len );

int main() {

  static const int N = 2;

  double a[] = {1, 0,
		0, 4 };


  if (sizeof(a) / sizeof(a[0]) != N * N) {
    printf("Error: dimensionality of a mismatches, exiting.\n");
    exit(1);
  }

  const char U = 'U';
  int piv[N];
  int rank;
  const double tol = -1;
  double work[2 * N];
  int info;

  dpstrf_(&U, &N, a, &N, piv, &rank, &tol, work, &info, 1);

  printf("info: %i, rank: %i\n", info, rank);

  printf("pivotes:");
  for (int i = 0; i < N; ++ i)
    printf("\t %i", piv[i]);
  printf("\n");

  printf("a: \n");
  for (int i = 0; i < N; ++i) {
    for (int j = 0; j < N; ++ j)
      printf("\t %g", a[i + N * j]);
    printf("\n");
  }

  return 0;
}