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

4 次查看（过去 30 天）

Rylan 2022-1-8

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1625010-how-to-solve-the-linear-system-ax-b-in-mex-file-when-a-is-a-sparse-matrix

评论： 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 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 2022-1-14

在 MATLAB Online 中打开

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.

请先登录，再进行评论。

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