在进入 DGEEV 参数号 9 时具有非法值

Question

我第一次尝试使用 C 中的 LAPACK 对矩阵进行对角化，但我被卡住了。

我一直在尝试将这个示例 http://rcabreral.blogspot.co.uk/2010/05/eigenvalues-clapack.html 从 zgeev 修改为 dgeev。我看过DGEEV 输入参数，http://www.netlib.org/lapack/explore-html/d9/d28/dgeev_8f.html 但似乎我不太了解。

因此，下面的代码产生：

**** 在进入 DGEEV 参数编号 9 时具有非法值**

编辑：错误发生在dgeev 的调用中，跨越第 48 行到（包括）第 53 行。

编辑：请注意，参数与此处的规范不同 http://www.netlib.org/lapack/explore-html/d9/d28/dgeev_8f.html 因为它们已被翻译成指针。在 C 中使用这些 Fortran 例程时，这是必要的，如下所述： http://www.physics.orst.edu/~rubin/nacphy/lapack/cprogp.html

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

void dgeTranspose( double *Transposed, double *M ,int n) {    
  int i,j;
  for(i=0;i<n;i++)
    for(j=0;j<n;j++)
      Transposed[i+n*j] = M[i*n+j];
}
//.........................................................................
//  MatrixComplexEigensystem: computes the eigenvectors and eigenValues of input matrix A
//  The eigenvectors are stored in columns
//.........................................................................
void MatrixComplexEigensystem( double *eigenvectorsVR, double *eigenvaluesW, double *A, int N){
  int i;
  double *AT = (double *) malloc( N*N*sizeof(double ) );

  dgeTranspose( AT, A , N);

  char JOBVL ='N';   // Compute Right eigenvectors
  char JOBVR ='V';   // Do not compute Left eigenvectors
  double VL[1];

  int LDVL = 1;
  int LDVR = N;
  int LWORK = 4*N; 

  double *WORK =  (double *)malloc( LWORK*sizeof(double));
  double *RWORK = (double *)malloc( 2*N*sizeof(double));
  int INFO;

  dgeev_( &JOBVL, &JOBVR, &N, AT ,  &N , eigenvaluesW ,   
       VL, &LDVL,
       eigenvectorsVR, &LDVR, 
       WORK, 
       &LWORK, RWORK, &INFO );

  dgeTranspose( AT, eigenvectorsVR , N);

  for(i=0;i<N*N;i++) eigenvectorsVR[i]=AT[i];

  free(WORK);
  free(RWORK);
  free(AT);
}

int main(){
  int i,j;
  const int N = 3;
  double A[] = { 1.+I , 2. ,  3 , 4. , 5.+I , 6. , 7., 8., 9. + I};
  double eigenVectors[N*N];
  double eigenValues[N];

  MatrixComplexEigensystem( eigenVectors, eigenValues, A, N);

  printf("\nEigenvectors\n");
  for(i=0;i<N;i++){
    for(j=0;j<N;j++) printf("%e", eigenVectors[i*N + j]);
    printf("\n");
  }
  printf("\nEigenvalues \n");
  for(i=0;i<N;i++) printf("%e",  eigenValues[i] );
  printf("\n------------------------------------------------------------\n"); 

  return 0;
}

Answer 1

您不能直接从zgeev 移植到dgeev。 zgeev 得到一个复数矩阵并计算复数特征值。而dgeev 得到一个实矩阵并计算复杂的特征值。为了保持一致，LAPACK 使用WR 和WI，分别用于每个特征值的实部和虚部。

所以注意dgeev的定义是

void dgeev_(char* JOBVL, char* JOBVR, int* N, double* A, int* LDA, double* WR, double* WI, double* VL, int* LDVL, double* VR, int* LDVR, double* WORK, int* LWORK, int* INFO);

我对您的示例的建议是删除：

#include <complex.h>

从双精度矩阵中删除I：

double A[] = { 1. , 2. ,  3 , 4. , 5. , 6. , 7., 8., 9.};

然后将特征值向量的大小加倍：

double eigenValues[2*N];

并使用WR 和WI 调用dgeev：

double *eigenvaluesWR = eigenvaluesW;
double *eigenvaluesWI = eigenvaluesW+N;
dgeev_(&JOBVL, &JOBVR, &N, AT, &N, 
    eigenvaluesWR, eigenvaluesWI, 
    VL, &LDVL, 
    eigenvectorsVR, &LDVR,
    WORK, &LWORK, &INFO);