我有代码进行数值优化(大型程序)。在 gdb 中,它以
Program received signal EXC_BAD_ACCESS, Could not access memory.
Reason: 13 at address: 0x0000000000000000
[Switching to process 8673 thread 0x1203]
0x00007fff95d2a9fb in small_malloc_from_free_list ()
在回溯期间,我被告知以下行包含错误,但无法找出问题所在。
177 U = calloc(cols*cols,sizeof(double));
包含该行的代码是
int fcn_invXX(double *X, int rows, int cols, double *invXX)
{
//---------------------------------------------------------------
// Declarations
//---------------------------------------------------------------
// Counters and info
int i, j, info;
// Upper triangular placeholder, U, and X'X
double *U, *XX;
// Indicator for upper triangular
char s;
// Alpha and beta for matrix multiplication - see dgemm documentation
double alpha, beta;
//---------------------------------------------------------------
// Memory Allocation
//---------------------------------------------------------------
U = calloc(cols*cols,sizeof(double));
XX = calloc(cols*cols,sizeof(double));
//---------------------------------------------------------------
// Compute inverse of X'X
//---------------------------------------------------------------
alpha = 1.0;
beta = 0.0;
// Get X'X
cblas_dgemm(CblasColMajor, CblasTrans, CblasNoTrans, cols, cols, rows, alpha, X, rows, X, rows, beta, XX, cols);
// CblasColMajor means that data are ordered in columns first, so if X is 4x4, X[0] and X[1] are first column
// CblasTrans means that the matrix should be transposed before multiplication. Here first matrix is transposed
// CblasNoTrans means that matrix should not be transposed. Here second matrix is not transposed
// X is the first matrix, and here, also the second
// XX stores the product
for (i = 0; i < cols*cols; i++)
{
U[i] = XX[i];
}
// Upper Cholesky factorization
s = 'U'; // upper
dpotrf_(&s, &cols, U, &cols, &info);
// info = LAPACKE_dpotrf(LAPACK_COL_MAJOR, s, cols, U, cols);
for (i = 0; i < cols*cols; i++)
{
invXX[i] = U[i];
}
// Get inverse
dpotri_(&s, &cols, invXX, &cols, &info);
// info = LAPACKE_dpotri(LAPACK_COL_MAJOR, s, cols, invXX, cols);
for (i = 0; i < cols; i++)
{
for (j = 0; j < cols; j++)
{
if (i > j)
{
invXX[i + j * cols] = invXX[j + i * cols];
}
}
}
//---------------------------------------------------------------
// Free Memory
//---------------------------------------------------------------
free(U);
free(XX);
return info;
}
希望你能帮忙。