How to solve the linear system Ax = b in mex file when A is a sparse matrix

1 次查看（过去 30 天）
Rylan
Rylan 2022-1-8
评论： Rylan 2022-1-14
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
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.
Rylan
Rylan 2022-1-14
Well, thanks for you warmly reply @James Tursa. Actually, I have sort this out by calling MKL. And I've already confirmed that it can give the right answer. However, I run into another question. The thread control in mex function does not take effect. Since I need to use OPENMP to accelarete the computation in the mex function, the MKL should run in serial. I use the following thread control sentence in my program:
mkl_set_num_threads(1);
However, the thread number for each part of the mkl domain, i.e. MKL_DOMAIN_PARDISO, MKL_DOMAIN_BLAS, remain unchanged.
Do you have any idea about this? thanks.

请先登录，再进行评论。

请先登录，再回答此问题。

回答（0 个）

类别

在 Help Center 和 File Exchange 中查找有关 Sparse Matrices 的更多信息

标签

产品


版本

R2020b

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by