计算科学 - LAPACK zlapmt_冻结代码 - 吾爱随笔录

我正在使用 LAPACK 的 zggev_ 例程来解决一些广义特征值问题。虽然它产生了正确的结果，但我想要特征值和根据绝对值排序的特征向量。为此，我实现了一个简单的排序算法，它为我提供了 LAPACK 例程 zlapmt_ 所需的置换向量，它应该对特征向量进行排序。

运行代码时，它甚至在进入 main() 函数之前就冻结了。如果 zlapmt_ 调用被注释掉，一切正常。

有人也遇到过这个问题吗？还是我在代码中做错了什么？我在下面附上了一个 MWE

谢谢您的帮助！

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