При входе в параметр DGEEV номер 9 имел незаконное значение

Question

Я пытаюсь в первый раз использовать LAPACK из C, чтобы диагонализировать матрицу, и я застрял.

Я пытался изменить этот пример 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 имел недопустимое значение **

EDIT: Ошибка возникает в вызове от dgeev охватывающие линии 48 до (включая) 53.

EDIT: Обратите внимание, что аргументы отличаются от приведенных здесь спецификаций http://www.netlib.org/lapack/explore-html/d9/d28/dgeev_8f.html в том смысле, что они были переведены на указатели. Это необходимо при использовании этих Fortran подпрограмм в C, как описано здесь: 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 определение

недействительным dgeev_ (символ * JOBVL, символ * JOBVR, Int * N, двойной * A, INT * LDA, двойной * WR, двойной * WI, двойной * 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]; 

и позвоните по телефону dgeev использованием WR и WI:

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