LAPACK zlapmt_ freezes code [closed]

Question

I'm using LAPACK's zggev_ routine to solve some generalized eigenvalue problem. While it produces the correct results, I want the eigenvalues and according eigenvectors sorted by absolute value. For this I implemented a naive sorting algorithm which provides me with the permutation vector needed by the LAPACK routine zlapmt_ which should sort the eigenvectors.

When running the code it freezes before even entering the main() function. If the zlapmt_ call is commented out everything works fine.

Anyone experienced this issue as well? Or do I do something wrong in the code? I attached a MWE below

Thank you for your help!

#include <stdio.h>
#include <complex.h>
#include <math.h>

extern double dznrm2_( int *n, double complex *x, int *incx );
extern void zdrscl_( int *n, double *alpha, double complex *x, int *incx );
extern void zlapmt_( int *fwd, int *m, int *n, double complex *X, int *ldx, int *k );
extern void zggev_( char *jobvl, char *jobvr, int *n, double complex *A, int *lda, double complex *B, int *ldb, 
                double complex *alpha, double complex *beta, double complex *vl, int *ldvl, 
                double complex *vr, int *ldvr, double complex *work, int *lwork, double *rwork, int *info );

int main(void) {
  int m = 3, lwork = 2*m*m, pi[m], info;
  double rwork[8*m];
  char jobvl = 'N', jobvr = 'V';
  double complex theta[m], alpha[m], beta[m], ev[m*m], work[lwork];
  double complex A[9] = { 12, -14, 3, 4, -4, 13, 2,  0, -4 }; // matrices stored columnwise
  double complex B[9] = { 4, -9, 21, -2, 1, -6, 15, 17, 3};

  zggev_( &jobvl, &jobvr, &m, A, &m, B, &m, alpha, beta, NULL, &m, ev, &m, work, &lwork, rwork, &info );
  if( info ) printf("LAPACK error! Errorcode: %d\n", info);

  // obtain eigenvalues
   for ( int i=0; i<3; i++ ) theta[i] = alpha[i]/beta[i];

  // scale eigenvectors
  info = 1;
  for ( int i=0; i<m; i++ ) {
    rwork[0] = dznrm2_( &m, ev+i*m, &info );
    zdrscl_( &m, rwork, ev+i*m, &info );
  }

  // sort eigenvalues ascending in absolute value (implicit reordering)
  for( int i=0; i<m; i++ )
    pi[i] = i;
  for( int i=m; i>1; i-- ) {
    for (int j=0; j<i-1; j++ ) {
      if ( fabs(ev[pi[j]]) > fabs(ev[pi[j+1]]) ) {
        info = pi[j];
        pi[j] = pi[j+1];
        pi[j+1] = info;
      }
    }
  }
  info = 1;
  zlapmt_( &info, &m, &m, ev, &m, pi );

  // output
  for ( int i=0; i<3; i++ ) {
    printf("theta[%d] = %lf\nev[%d] = ",i, creal(theta[i]), i);
    for ( int j=0; j<3; j++ ) 
      printf("%lf ", creal(ev[i*m+j]));
    printf("\n\n");
  }

  return 0;
}

Answer 1

A couple of notes:

• By introducing printf statements throughout the code, you would be able to see that even when your zlapmt_ call is present, the code certainly enters the main function (it would be obscure if it does not).
• Right now, your sorting algorithm is certainly non-optimal. I would strongly suggest you use quicksort or something else. Preferably, not a custom implementation.
• In your sorting algorithm, you are using fabs to get an absolute value of the argument. In your case, it is double complex. According to the specification, you should use cabs in this case. I even got a compile error with fabs.

Regarding permutations:

Right now, you are using indexing from $0$ to $M-1$ in your permutation array pi. According to lapmt_ documentation, you indexing starts from $1$. So, changing to:

for( int i=0; i<m; i++ )
    pi[i] = i+1;

Should solve the problem.

Also, I think there is a problem with your output code in the end. You never sort your theta's; however, you output the sorted array ev. Make sure, you are doing what you actually intend to be doing there.