dbgheap.c бросает исключение нарушения прав доступа

Question

я разработал следующий фрагмент кода, который работает хорошо:

#include "header.h" 

int test_DGGEV_11_14a(){ 
    const int n=3; 
    double a[n][n]={1,1,1,2,3,4,3,5,2}; 
    double b[n][n]={-10,-3,12,14,14,12,16,16,18}; 
    //double a[n][n]={1,7,3,2,9,12,5,22,7}; 
    //double b[n][n]={1,7,3,2,9,12,5,22,7}; 

    /*const int n=2; 
    double a[n][n]={1e-16,0,0,1e-15}; 
    double b[n][n]={1e-16,0,0,1e-15};*/ 

    lapack_int info; 
    double alphar[n]={0.0}; 
    double alphai[n]={0.0}; 
    double beta[n]={0.0}; 
    double vl[n][n]={0.0}; 
    double vr[n][n]={0.0}; 

    info=LAPACKE_dggev(LAPACK_ROW_MAJOR,'V','V',n,*a,n,*b,n,alphar,alphai,beta,*vl,n,*vr,n); 

    std::cout<<"right eigen vector (what we want):\n"; 
    for(int i=int(0);i<n;i++){ 
     for(int j=int(0);j<n;j++){ 
      printf("%1f ",vr[i][j]); 
     } 
     printf("\n"); 
    } 
    std::cout<<"left eigen vector:\n"; 
    for(int i=int(0);i<n;i++){ 
     for(int j=int(0);j<n;j++){ 
      printf("%1f ",vl[i][j]); 
     } 
     printf("\n"); 
    } 
    std::cout<<"eigen values:\n"; 
    for(int i=int(0);i<n;i++){ 
     if(beta[i]>DBL_MIN || beta[i]<-DBL_MIN){ 
      printf("%1f ",alphar[i]/beta[i]); 
      printf("\n"); 
     }else{ 
      printf("%1f ","beta is zero"); 
      printf("\n"); 
     } 
    } 
    return info; 
} 

Я изменял выше правильный код, чтобы использовать LAPACKE DGGEV процедуру для больших матриц, модифицированный код как показано ниже:

#include "header.h" 

int test_DGGEV_11_17a(){ 
    const int r=342; 
    const int c=342; 
    double**a=NULL;//stiffness 
    a=new double *[r]; 
    for(int i=int(0);i<r;i++) 
     a[i]=new double[c]; 
    readFile("Input_Files/OUTPUT_sub_2_stiffness.txt",a,r,c); 
    writeFile("Output_Files/K.txt",a,r,c);//to check if readFile was OK 

    double**b=NULL;//mass 
    b=new double*[r]; 
    for(int i=int(0);i<r;i++) 
     b[i]=new double[c]; 
    readFile("Input_Files/OUTPUT_sub_2_mass.txt",b,r,c); 
    writeFile("Output_Files/M.txt",b,r,c);//to check if readFile was OK 

    const int n=r;//r=c=n 
    lapack_int info=110; 
    double alphar[n]={0.0}; 
    double alphai[n]={0.0}; 
    double beta[n]={0.0}; 
    //double vl[n][n]={0.0};//generates stack overflow 
    double**vl=NULL; 
    vl=new double*[r]; 
    for(int i=int(0);i<r;i++) 
     vl[i]=new double[c]; 
    for(int i=int(0);i<r;i++) 
     for(int j=int(0);j<c;j++) 
      vl[i][j]=0.0; 
    //double vr[n][n]={0.0};//generates stack overflow 
    double**vr=NULL; 
    vr=new double*[r]; 
    for(int i=int(0);i<r;i++) 
     vr[i]=new double[c]; 
    for(int i=int(0);i<r;i++) 
     for(int j=int(0);j<c;j++) 
      vr[i][j]=0.0; 

    info=LAPACKE_dggev(LAPACK_ROW_MAJOR,'V','V',n,*a,n,*b,n,alphar,alphai,beta,*vl,n,*vr,n); 

    return info; 
} 

в приведенном выше модифицированного кода (для больших матриц), я должен выделить память из кучи, потому что в противном случае, стек получит переполнение. Проблема заключается в том, что, когда я выделить память из кучи по new я получаю следующее исключение, которое связано с кучи и происходит внутри dbgheap.c (функции Debug CRT Heap):

Кто-нибудь знает, почему это исключение происходит? возможно, это связано с тем, что DLL LAPACKE используют другую кучу для распределений ... Я не знаю.

---

РЕДАКТИРОВАТЬ:

трассировки стека заключается в следующем:

---

EDIT:

Окончательно решил проблему, заменив все 2D-массивы 1D-массивами. Следующий код - это исправленный код, который работает без ошибок. Для получения подробной информации об этом решении см. Ответ «Илья Кобелевский».

int test_DGGEV_11_18a(){ 
    const int r=342; 
    const int c=342; 
    double*a=NULL;//stiffness 
    a=new double [r*c]; 
    for(int i=int(0);i<r*c;i++) 
      a[i]=0.0; 
    readFile_1Darray("Input_Files/OUTPUT_sub_2_stiffness.txt",a,r,c); 
    writeFile_1Darray("Output_Files/K.txt",a,r,c);//to check if readFile was OK 

    double*b=NULL;//mass 
    b=new double[r*c]; 
    for(int i=int(0);i<r*c;i++) 
      b[i]=0.0; 
    readFile_1Darray("Input_Files/OUTPUT_sub_2_mass.txt",b,r,c); 
    writeFile_1Darray("Output_Files/M.txt",b,r,c);//to check if readFile was OK 

    const int n=r;//r=c=n 
    lapack_int info=110; 

    //double alphar[n]={0.0}; 
    double*alphar=NULL; 
    alphar=new double[n]; 
    for(int i=int(0);i<n;i++) 
     alphar[i]=0.0; 
    //double alphai[n]={0.0}; 
    double*alphai=NULL; 
    alphai=new double[n]; 
    for(int i=int(0);i<n;i++) 
     alphai[i]=0.0; 
    //double beta[n]={0.0}; 
    double*beta=NULL; 
    beta=new double[n]; 
    for(int i=int(0);i++;) 
     beta[i]=0.0; 
    //double vl[n][n]={0.0};//generates stack overflow 
    double*vl=NULL; 
    vl=new double[r*c]; 
    for(int i=int(0);i<r*c;i++) 
      vl[i]=0.0; 
    //double vr[n][n]={0.0};//generates stack overflow 
    double*vr=NULL; 
    vr=new double[r*c]; 
    for(int i=int(0);i<r*c;i++) 
      vr[i]=0.0; 

    info=LAPACKE_dggev(LAPACK_ROW_MAJOR,'V','V',n,a,n,b,n,alphar,alphai,beta,vl,n,vr,n); 
    std::cout<<"info returned by LAPACKE_dggev:\t"<<info<<'\n'; 

    double*eigValueReal=NULL; 
    eigValueReal=new double[n]; 
    for(int i=int(0);i<n;i++) 
     eigValueReal[i]=0.0; 
    for(int i=int(0);i<n;i++) 
     eigValueReal[i]=alphar[i]/beta[i]; 

    write1Darray("Output_Files/eigValueReal_LAPACKE_DGGEV.txt",eigValueReal,n); 
    write1Darray("Output_Files/beta.txt",beta,n); 
    writeFile_1Darray("Output_Files/eigVectorRight_LAPACKE_DGGEV.txt",vr,r,c); 

    delete[] a; 
    delete[] b; 
    delete[] alphar; 
    delete[] alphai; 
    delete[] beta; 
    delete[] vl; 
    delete[] vr; 
    delete[] eigValueReal; 

    return info; 
} 

Answer 1

Согласно documentation LAPACKE_dggev ожидает двойной * в качестве входного сигнала, так что все матрицы должны быть сохранены в виде линейных массивов.

Вместо

double**a=NULL;//stiffness 
a=new double *[r]; 
for(int i=int(0);i<r;i++) 
    a[i]=new double[c]; 

Вы должны использовать

double*a=new double[c*r];//stiffness 

с аналогичными изменениями для всех otehr матриц. Доступ к элементам матрицы можно получить как [i * c + j] для [i, j]. Теперь могут возникнуть другие проблемы с вашим кодом, но это очевидно, что может привести к ошибке, которую вы видите.