當我直接提供矩陣時，爲什麼這個LAPACK程序能夠正確工作，但是當我從文件中讀取時，不能正確工作？

Question

下面是一個用於對角矩陣A的LAPACK代碼，我以數組a的形式提供。這只是一個官方例子的輕微修改，似乎產生了正確的結果。這是不切實際的，因爲我必須直接提供數組。

#include <stdlib.h> 
#include <stdio.h> 
#include <fstream> 
#include <vector> 


/* DSYEV prototype */ 
extern "C"{ 
void dsyev(char* jobz, char* uplo, int* n, double* a, int* lda, 
       double* w, double* work, int* lwork, int* info); 
} 
/* Auxiliary routines prototypes */ 
extern "C"{ 
void print_matrix(char* desc, int m, int n, double* a, int lda); 
} 
/* Parameters */ 
#define N 5 
#define LDA N 

/* Main program */ 
int main() { 
     /* Locals */ 
     int n = N, lda = LDA, info, lwork; 
     double wkopt; 
     double* work; 
     /* Local arrays */ 
     double w[N]; 
     double a[LDA*N] = { 
      1.96, 0.00, 0.00, 0.00, 0.00, 
      -6.49, 3.80, 0.00, 0.00, 0.00, 
      -0.47, -6.39, 4.17, 0.00, 0.00, 
      -7.20, 1.50, -1.51, 5.70, 0.00, 
      -0.65, -6.34, 2.67, 1.80, -7.10 
     }; 
     /* Executable statements */ 
     printf(" DSYEV Example Program Results\n"); 
     /* Query and allocate the optimal workspace */ 
     lwork = -1; 
     dsyev("Vectors", "Upper", &n, a, &lda, w, &wkopt, &lwork, &info); 
     lwork = (int)wkopt; 
     work = (double*)malloc(lwork*sizeof(double)); 
     /* Solve eigenproblem */ 
     dsyev("Vectors", "Upper", &n, a, &lda, w, work, &lwork, &info); 
     /* Check for convergence */ 
     if(info > 0) { 
       printf("The algorithm failed to compute eigenvalues.\n"); 
       exit(1); 
     } 
     /* Print eigenvalues */ 
     print_matrix("Eigenvalues", 1, n, w, 1); 
     /* Print eigenvectors */ 
     print_matrix("Eigenvectors (stored columnwise)", n, n, a, lda); 
     /* Free workspace */ 
     free((void*)work); 
     exit(0); 
} /* End of DSYEV Example */ 

/* Auxiliary routine: printing a matrix */ 
void print_matrix(char* desc, int m, int n, double* a, int lda) { 
     int i, j; 
     printf("\n %s\n", desc); 
     for(i = 0; i < m; i++) { 
       for(j = 0; j < n; j++) printf(" %6.2f", a[i+j*lda]); 
       printf("\n"); 
     } 
} 

我僅要修改上面的代碼，這樣我可以從文件中讀取該陣列，而不是直接提供它的。爲此，我編寫了從文件peano_covariance.data中讀取數組的函數read_covariance。後面的數據文件的內容是：

1.96 0.00 0.00 0.00 0.00 
-6.49 3.80 0.00 0.00 0.00 
-0.47 -6.39 4.17 0.00 0.00 
-7.20 1.50 -1.51 5.70 0.00 
-0.65 -6.34 2.67 1.80 -7.10 

以下是我的嘗試，它會產生非常不正確的特徵值和特徵向量。

#include <stdlib.h> 
#include <stdio.h> 
#include <fstream> 
#include <vector> 


int read_covariance (std::vector<double> data) 
    { 
    double tmp; 

    std::ifstream fin("peano_covariance.data"); 

    while(fin >> tmp) 
    { 
     data.push_back(tmp); 
    } 

    return 0; 
} 

/* DSYEV prototype */ 
extern "C"{ 
void dsyev(char* jobz, char* uplo, int* n, double* a, int* lda, 
       double* w, double* work, int* lwork, int* info); 
} 
/* Auxiliary routines prototypes */ 
extern "C"{ 
void print_matrix(char* desc, int m, int n, double* a, int lda); 
} 
/* Parameters */ 
#define N 5 
#define LDA N 

/* Main program */ 
int main() { 
     /* Locals */ 
     std::vector<double> data; 
     int n = N, lda = LDA, info, lwork; 
     double wkopt; 
     double* work; 
     /* Local arrays */ 
     double w[N]; 
     double a[LDA*N]; 
     read_covariance(data); 

     std::copy(data.begin(), data.end(), a); 
     /* Executable statements */ 
     printf(" DSYEV Example Program Results\n"); 
     /* Query and allocate the optimal workspace */ 
     lwork = -1; 
     dsyev("Vectors", "Upper", &n, a, &lda, w, &wkopt, &lwork, &info); 
     lwork = (int)wkopt; 
     work = (double*)malloc(lwork*sizeof(double)); 
     /* Solve eigenproblem */ 
     dsyev("Vectors", "Upper", &n, a, &lda, w, work, &lwork, &info); 
     /* Check for convergence */ 
     if(info > 0) { 
       printf("The algorithm failed to compute eigenvalues.\n"); 
       exit(1); 
     } 
     /* Print eigenvalues */ 
     print_matrix("Eigenvalues", 1, n, w, 1); 
     /* Print eigenvectors */ 
     print_matrix("Eigenvectors (stored columnwise)", n, n, a, lda); 
     /* Free workspace */ 
     free((void*)work); 
     exit(0); 
} /* End of DSYEV Example */ 

/* Auxiliary routine: printing a matrix */ 
void print_matrix(char* desc, int m, int n, double* a, int lda) { 
     int i, j; 
     printf("\n %s\n", desc); 
     for(i = 0; i < m; i++) { 
       for(j = 0; j < n; j++) printf(" %e", a[i+j*lda]); 
       printf("\n"); 
     } 
} 

Answer 1

更換

int read_covariance (std::vector<double> data) 

與

int read_covariance (std::vector<double> & data) 

您要發送的數組的副本，而不是對它的引用。這是臨時副本，正在填充值。這是bg2b在他的評論中提到的。

就個人而言，雖然，我寧願喜歡寫東西

int read_covariance (const std::string & fname) 
{ 
    std::ifstream in(fname.c_str()); 
    double val; 
    std::vector<double> cov; 
    while(in >> val) cov.push_back(val); 
    return cov; 
} 

更妙的是使用合適的多維數組庫，而不是笨重的一維向量。有很多這樣的庫可用，我不確定哪一個是最好的（在C++標準庫中缺少一個好的多維數組類是我常常使用fortran的主要原因之一），但是ndarray的外觀有趣的 - 它旨在模仿Python的優秀numpy陣列模塊的功能。