error "Error in `. / test": double free or corruption (out): 0x00000000011c6f30" occurred in free (Hk). I can"t see what went wrong.
Code:
void linalg::sp_Davidson(const char uplo, const SparseMatrix_COO& a, double* diagValsA, double* eigVecs, double* eigVals, int nroot, int p, double tol, int MaxCycle, int MaxSpace, std::ostream& os)
{
int n = a.nRows(); // the size of matrix a
// workspaces
double* Vk = (double*)std::malloc(sizeof(double) * (n * MaxSpace));
double* Wk = (double*)std::malloc(sizeof(double) * (n * MaxSpace));
double* Hk = (double*)std::malloc(sizeof(double) * (MaxSpace * MaxSpace));
double* Sk = (double*)std::malloc(sizeof(double) * (MaxSpace)); // eigenvals of Hk
double* Sold = (double*)std::malloc(sizeof(double) * nroot);
double* Q1 = (double*)std::malloc(sizeof(double) * n);
double* Q2 = (double*)std::malloc(sizeof(double) * n);
// generate initial guess Vk
std::memset(Vk, 0, (sizeof(double) * n * p));
for(int i = 0; i < p; PPi) Vk[i * n + i] = 1.0;
// iterative diagonalization
int ispace = p; // current subspace size
int k = 0; // cycle count
while(k < MaxCycle)
{
os << "--> k: " << k << std::endl;
// current space size
os << "ispace: " << ispace << std::endl;
// MGS(Vk)
MGS(Vk, n, ispace);
// compute W_k = A V_k, size: [n, ispace] = [n, n] * [n, ispace]
sp_coosymm(uplo, a, Vk, ispace, Wk);
// compute H_k = V_K^T A V_k = V_K^T W_k
cblas_dgemm(CblasColMajor, CblasTrans, CblasNoTrans, ispace, ispace, n, 1.0, Vk, n, Wk, n, 0.0, Hk, ispace);
// Sk, Hk = eig(Hk)
LAPACKE_dsyevd(LAPACK_COL_MAJOR, "V", "U", ispace, Hk, ispace, Sk);
// compute eigenvecs eigVecs = Vk * Hk[:p], size: [n, p] = [n, ispace] * [ispace, p], Ritz approximations
cblas_dgemm(CblasColMajor, CblasNoTrans, CblasNoTrans, n, p, ispace, 1.0, Vk, n, Hk, ispace, 0.0, eigVecs, n);
double norm = 9999.0;
for(int j = 0; j < p; PPj)
{
// Q1 = Wk*Hk[j] - S[j]*Vk*Hk[j]
// Q1 = Wk*Hk[j]
cblas_dgemv(CblasColMajor, CblasNoTrans, n, ispace, 1.0, Wk, n, (Hk + j * ispace), 1, 0.0, Q1, 1);
// Q2 = Vk*Hk[j]
cblas_dgemv(CblasColMajor, CblasNoTrans, n, ispace, 1.0, Vk, n, (Hk + j * ispace), 1, 0.0, Q2, 1);
// Q1 = - Sk[j]*Q2 + Q1
cblas_daxpy(n, -Sk[j], Q2, 1, Q1, 1);
// Q1 = Q1 / (Sk[j] - A[j][j])
cblas_dscal(n, (Sk[j] - diagValsA[j]), Q1, 1);
// V[:,ispace] = q1
cblas_dcopy(n, Q1, 1, (Vk + n * ispace), 1);
PPispace;
if(ispace >= MaxSpace) // restart
{
ispace = p;
break;
}
}
// compute norm for n root
// Sold = -Sk + Sold
cblas_daxpy(nroot, -1.0, Sk, 1, Sold, 1);
norm = cblas_dnrm2(nroot, Sold, 1);
// Sold = Sk
cblas_dcopy(nroot, Sk, 1, Sold, 1);
if(norm < tol)
{
break;
}
PPk;
}
if(k < MaxCycle)
{
os << "------> convergence reached! <------" << std::endl;
os << "--> nloop: " << k << std::endl;
}
else
{
os << "------> unconverged! <------" << std::endl;
}
cblas_dcopy(nroot, Sold, 1, eigVals, 1);
std::free(Q2);
std::free(Q1);
std::free(Sold);
std::free(Sk);
std::free(Hk);
std::free(Wk);
std::free(Vk);
}because the MaxSpace is relatively large, there will be no out-of-bounds problems in the process. So what"s wrong with the free of the Hk array?
Thank you for the answer.