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