I have developed the following piece of code which works good:
#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;
}
I modified the above correct code, to use LAPACKE DGGEV routine for large matrices, the modified code is shown below:
#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;
}
In the above modified code (for large matrices), I have to allocate memory from heap because otherwise, stack would get overflow. The problem is that when I allocate memory from heap by new
I get the following exception which is related to heap and occurs inside dbgheap.c
(Debug CRT Heap Functions):
Does anybody know why this exception happens? maybe it is related to the fact that LAPACKE DLLs are using a different heap for allocations...I don't know.
EDIT:
the stack trace is this:
EDIT:
Finally solved the problem by replacing all the 2D arrays with 1D arrays. The following code is the corrected code which works without any error. Please see answer of "Ilya Kobelevskiy" for the details of this solution.
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;
}
According to documentation LAPACKE_dggev expects double* as an input, so all matrices need to be stored as linear arrays.
Instead of
You should use
With similar changes for all otehr matrices. Matrix elements can be accessed as a[i*c+j] for a[i,j]. Now, there might be other problems with your code, but this is an obvious one that may lead to the error you are seeing.