c - 在进入 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但似乎我还不够了解。

因此，下面的代码产生:

**** On entry to DGEEV parameter number 9 had an illegal value**

编辑:错误发生在 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>

从 double 矩阵中删除I:

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

然后将特征值 vector 的大小加倍:

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