How to solve the linear system Ax = b in mex file when A is a sparse matrix
4 次查看(过去 30 天)
显示 更早的评论
I have multiple linear systems Ax = b to be solved in mex function. When A is of small or medium size, A is stored as a full matrix. For this case, this system is solved with the following function.
void EquationSolve(double *dest, ptrdiff_t m, ptrdiff_t n, double *CoeMatrix, double dt, double ImplicitParam, int Nvar, int K3d, int Np)
{
double *TempCoeMatrix = malloc(m*m*sizeof(double));
ptrdiff_t *iPivot = malloc(m*sizeof(ptrdiff_t));
ptrdiff_t info;
ptrdiff_t p = m;
for (int var = 0; var < Nvar; var++){
memcpy(TempCoeMatrix, CoeMatrix + var*m*m, m*m*sizeof(double));
for (int col = 0; col < m; col++){
for (int row = 0; row < m; row++){
TempCoeMatrix[col*m + row] = -1 * dt*ImplicitParam*TempCoeMatrix[col*m + row];
}
TempCoeMatrix[col*m + col] += 1;
}
% This is the core of steps used to solve the linear system.
dgesv(&m, &n, TempCoeMatrix, &m, iPivot, dest + var*K3d*Np, &p, &info);
}
free(iPivot);
free(TempCoeMatrix);
}
However, if A is of large size, and much elements contained in A are zeros. Storing A as a full matrix is not wise. I want to convert A into a sparse one( this is easy to implement, since I already know the structure of A). But how to solve this system in mex function for such case.
I know, in matlab, solving this system can be easily implemented with x = A\b, no matter A is full or sparse.
2 个评论
James Tursa
2022-1-14
You would need to use a sparse library. I thought Tim Davis had such a library on the FEX, but I can't find it there anymore.
回答(0 个)
另请参阅
类别
在 Help Center 和 File Exchange 中查找有关 Operating on Diagonal Matrices 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!