Hi. We are moving from MKL to Accelerate in order to accommodate the transition to apple silicon. However, we are occasionally getting crashes inside the Accelerate framework when factoring sparse matrices with the Cholesky decomposition.
In the following link you can find an Xcode project with a minimal reproducible example: https://drive.google.com/file/d/1rHmJZbA5yc4-68Z1-vm7g3IvIPRgbN_c/view?usp=share_link
I also pasted the code below.
First, we load a sparse matrix from the hard drive. Then we factor it. When factoring the matrix with LDLT, everything works fine. When factoring it with Cholesky decomposition, the SparseFactor sometimes returns -3, sometimes it returns -1, and sometimes it crashes. The code below tries the factoring 100 times and always crashes. The matrix is positive definite, I checked it in another tool. It's smallest (algebraically) eigenvalue is at around 23000. Other matrices, such as the one in the 'good' file attached in the link above, doesn't cause crashes.
Can someone please shed light onto what's causing the crash? What exactly in this matrix causes the crash, and is this a known issue in the Accelerate framework?
#include <iostream>
#include <fstream>
#include <vector>
#include <libgen.h>
#include <Accelerate/Accelerate.h>
int main(int argc, const char * argv[]) {
std::string filename = __FILE__;
std::string folder = dirname(const_cast<char*>(filename.c_str()));;
std::vector<std::string> files = {
folder + "/m_good.bin",
folder + "/m_bad.bin"
};
int GOOD = 0;
int BAD = 1;
std::string datafilename = files[BAD];
bool cholesky = true;
std::ifstream file(datafilename, std::ios::binary);
if (!file.is_open())
{
std::cout << " can't find file: " << datafilename << "\n";
return 1;
}
SparseMatrix_Double _A;
// Read the data into a buffer
file.read((char*)&_A, sizeof(_A));
std::vector<double> _values;
auto size = _values.size();
file.read((char*)&size, sizeof(size));
_values.resize(size);
file.read((char*)&_values[0], size * sizeof(double));
std::vector<long> _column_starts;
file.read((char*)&size, sizeof(size));
_column_starts.resize(size);
file.read((char*)&_column_starts[0], size * sizeof(long));
std::vector<int> _row_indices;
file.read((char*)&size, sizeof(size));
_row_indices.resize(size);
file.read((char*)&_row_indices[0], size * sizeof(int));
file.close();
_A.data = const_cast<double*>(&_values[0]);
_A.structure.columnStarts = &_column_starts[0];
_A.structure.rowIndices = &_row_indices[0];
for (int i = 0; i < 100; i++)
{
std::cout << "\nRun: " << i << " ";
auto _LLT = SparseFactor(cholesky ? SparseFactorizationCholesky : SparseFactorizationLDLT, _A);
std::cout << _LLT.status;
SparseCleanup(_LLT);
if (_LLT.status != SparseStatusOK)
{
std::cout << " Failed!!";
}
std::cout << std::endl;
// Close the file
}
return 0;
}