Actual source code: aijmatlab.c

      1: /*
      2:         Provides an interface for the MATLAB engine sparse solver

      4: */
      5: #include <../src/mat/impls/aij/seq/aij.h>
      6: #include <petscmatlab.h>
      7: #include <engine.h> /* MATLAB include file */
      8: #include <mex.h>    /* MATLAB include file */

     10: static mxArray *MatSeqAIJToMatlab(Mat B)
     11: {
     12:   Mat_SeqAIJ *aij = (Mat_SeqAIJ *)B->data;
     13:   mwIndex    *ii, *jj;
     14:   mxArray    *mat;
     15:   PetscInt    i;

     17:   mat = mxCreateSparse(B->cmap->n, B->rmap->n, aij->nz, mxREAL);
     18:   if (PetscArraycpy(mxGetPr(mat), aij->a, aij->nz)) return NULL;
     19:   /* MATLAB stores by column, not row so we pass in the transpose of the matrix */
     20:   jj = mxGetIr(mat);
     21:   for (i = 0; i < aij->nz; i++) jj[i] = aij->j[i];
     22:   ii = mxGetJc(mat);
     23:   for (i = 0; i < B->rmap->n + 1; i++) ii[i] = aij->i[i];
     24:   return mat;
     25: }

     27: PETSC_EXTERN PetscErrorCode MatlabEnginePut_SeqAIJ(PetscObject obj, void *mengine)
     28: {
     29:   mxArray *mat;

     31:   PetscFunctionBegin;
     32:   mat = MatSeqAIJToMatlab((Mat)obj);
     33:   PetscCheck(mat, PETSC_COMM_SELF, PETSC_ERR_LIB, "Cannot create MATLAB matrix");
     34:   PetscCall(PetscObjectName(obj));
     35:   engPutVariable((Engine *)mengine, obj->name, mat);
     36:   PetscFunctionReturn(PETSC_SUCCESS);
     37: }

     39: static PetscErrorCode MatSeqAIJFromMatlab(mxArray *mmat, Mat mat)
     40: {
     41:   PetscInt    nz, n, m, *i, *j, k;
     42:   mwIndex     nnz, nn, nm, *ii, *jj;
     43:   Mat_SeqAIJ *aij = (Mat_SeqAIJ *)mat->data;

     45:   PetscFunctionBegin;
     46:   nn  = mxGetN(mmat); /* rows of transpose of matrix */
     47:   nm  = mxGetM(mmat);
     48:   nnz = (mxGetJc(mmat))[nn];
     49:   ii  = mxGetJc(mmat);
     50:   jj  = mxGetIr(mmat);
     51:   n   = (PetscInt)nn;
     52:   m   = (PetscInt)nm;
     53:   nz  = (PetscInt)nnz;

     55:   if (mat->rmap->n < 0 && mat->cmap->n < 0) {
     56:     /* matrix has not yet had its size set */
     57:     PetscCall(MatSetSizes(mat, n, m, PETSC_DETERMINE, PETSC_DETERMINE));
     58:     PetscCall(MatSetUp(mat));
     59:   } else {
     60:     PetscCheck(mat->rmap->n == n, PETSC_COMM_SELF, PETSC_ERR_SUP, "Cannot change size of PETSc matrix %" PetscInt_FMT " to %" PetscInt_FMT, mat->rmap->n, n);
     61:     PetscCheck(mat->cmap->n == m, PETSC_COMM_SELF, PETSC_ERR_SUP, "Cannot change size of PETSc matrix %" PetscInt_FMT " to %" PetscInt_FMT, mat->cmap->n, m);
     62:   }
     63:   if (nz != aij->nz) {
     64:     /* number of nonzeros in matrix has changed, so need new data structure */
     65:     aij->nz = nz;
     66:     PetscCall(PetscShmgetAllocateArray(aij->nz, sizeof(PetscScalar), (void **)&aij->a));
     67:     PetscCall(PetscShmgetAllocateArray(aij->nz, sizeof(PetscInt), (void **)&aij->j));
     68:     PetscCall(PetscShmgetAllocateArray(mat->rmap->n + 1, sizeof(PetscInt), (void **)&aij->i));
     69:     aij->free_a  = PETSC_TRUE;
     70:     aij->free_ij = PETSC_TRUE;
     71:   }

     73:   PetscCall(PetscArraycpy(aij->a, mxGetPr(mmat), aij->nz));
     74:   /* MATLAB stores by column, not row so we pass in the transpose of the matrix */
     75:   i = aij->i;
     76:   for (k = 0; k < n + 1; k++) i[k] = (PetscInt)ii[k];
     77:   j = aij->j;
     78:   for (k = 0; k < nz; k++) j[k] = (PetscInt)jj[k];

     80:   for (k = 0; k < mat->rmap->n; k++) aij->ilen[k] = aij->imax[k] = aij->i[k + 1] - aij->i[k];

     82:   mat->nonzerostate++; /* since the nonzero structure can change anytime force the Inode information to always be rebuilt */
     83:   PetscCall(MatAssemblyBegin(mat, MAT_FINAL_ASSEMBLY));
     84:   PetscCall(MatAssemblyEnd(mat, MAT_FINAL_ASSEMBLY));
     85:   PetscFunctionReturn(PETSC_SUCCESS);
     86: }

     88: PETSC_EXTERN PetscErrorCode MatlabEngineGet_SeqAIJ(PetscObject obj, void *mengine)
     89: {
     90:   Mat      mat = (Mat)obj;
     91:   mxArray *mmat;

     93:   PetscFunctionBegin;
     94:   mmat = engGetVariable((Engine *)mengine, obj->name);
     95:   PetscCall(MatSeqAIJFromMatlab(mmat, mat));
     96:   PetscFunctionReturn(PETSC_SUCCESS);
     97: }

     99: static PetscErrorCode MatSolve_Matlab(Mat A, Vec b, Vec x)
    100: {
    101:   const char *_A, *_b, *_x;

    103:   PetscFunctionBegin;
    104:   /* make sure objects have names; use default if not */
    105:   PetscCall(PetscObjectName((PetscObject)b));
    106:   PetscCall(PetscObjectName((PetscObject)x));

    108:   PetscCall(PetscObjectGetName((PetscObject)A, &_A));
    109:   PetscCall(PetscObjectGetName((PetscObject)b, &_b));
    110:   PetscCall(PetscObjectGetName((PetscObject)x, &_x));
    111:   PetscCall(PetscMatlabEnginePut(PETSC_MATLAB_ENGINE_(PetscObjectComm((PetscObject)A)), (PetscObject)b));
    112:   PetscCall(PetscMatlabEngineEvaluate(PETSC_MATLAB_ENGINE_(PetscObjectComm((PetscObject)A)), "%s = u%s\\(l%s\\(p%s*%s));", _x, _A, _A, _A, _b));
    113:   PetscCall(PetscMatlabEngineEvaluate(PETSC_MATLAB_ENGINE_(PetscObjectComm((PetscObject)A)), "%s = 0;", _b));
    114:   /* PetscCall(PetscMatlabEnginePrintOutput(PETSC_MATLAB_ENGINE_(PetscObjectComm((PetscObject)A)),stdout));  */
    115:   PetscCall(PetscMatlabEngineGet(PETSC_MATLAB_ENGINE_(PetscObjectComm((PetscObject)A)), (PetscObject)x));
    116:   PetscFunctionReturn(PETSC_SUCCESS);
    117: }

    119: static PetscErrorCode MatLUFactorNumeric_Matlab(Mat F, Mat A, const MatFactorInfo *info)
    120: {
    121:   size_t    len;
    122:   char     *_A, *name;
    123:   PetscReal dtcol = info->dtcol;

    125:   PetscFunctionBegin;
    126:   if (F->factortype == MAT_FACTOR_ILU || info->dt > 0) {
    127:     /* the ILU form is not currently registered */
    128:     if (info->dtcol == PETSC_DEFAULT) dtcol = .01;
    129:     F->ops->solve = MatSolve_Matlab;
    130:     F->factortype = MAT_FACTOR_LU;

    132:     PetscCall(PetscMatlabEnginePut(PETSC_MATLAB_ENGINE_(PetscObjectComm((PetscObject)A)), (PetscObject)A));
    133:     _A = ((PetscObject)A)->name;
    134:     PetscCall(PetscMatlabEngineEvaluate(PETSC_MATLAB_ENGINE_(PetscObjectComm((PetscObject)A)), "info_%s = struct('droptol',%g,'thresh',%g);", _A, info->dt, dtcol));
    135:     PetscCall(PetscMatlabEngineEvaluate(PETSC_MATLAB_ENGINE_(PetscObjectComm((PetscObject)A)), "[l_%s,u_%s,p_%s] = luinc(%s',info_%s);", _A, _A, _A, _A, _A));
    136:     PetscCall(PetscMatlabEngineEvaluate(PETSC_MATLAB_ENGINE_(PetscObjectComm((PetscObject)A)), "%s = 0;", _A));

    138:     PetscCall(PetscStrlen(_A, &len));
    139:     PetscCall(PetscMalloc1(len + 2, &name));
    140:     PetscCall(PetscSNPrintf(name, len + 2, "_%s", _A));
    141:     PetscCall(PetscObjectSetName((PetscObject)F, name));
    142:     PetscCall(PetscFree(name));
    143:   } else {
    144:     PetscCall(PetscMatlabEnginePut(PETSC_MATLAB_ENGINE_(PetscObjectComm((PetscObject)A)), (PetscObject)A));
    145:     _A = ((PetscObject)A)->name;
    146:     PetscCall(PetscMatlabEngineEvaluate(PETSC_MATLAB_ENGINE_(PetscObjectComm((PetscObject)A)), "[l_%s,u_%s,p_%s] = lu(%s',%g);", _A, _A, _A, _A, dtcol));
    147:     PetscCall(PetscMatlabEngineEvaluate(PETSC_MATLAB_ENGINE_(PetscObjectComm((PetscObject)A)), "%s = 0;", _A));
    148:     PetscCall(PetscStrlen(_A, &len));
    149:     PetscCall(PetscMalloc1(len + 2, &name));
    150:     PetscCall(PetscSNPrintf(name, len + 2, "_%s", _A));
    151:     PetscCall(PetscObjectSetName((PetscObject)F, name));
    152:     PetscCall(PetscFree(name));

    154:     F->ops->solve = MatSolve_Matlab;
    155:   }
    156:   PetscFunctionReturn(PETSC_SUCCESS);
    157: }

    159: static PetscErrorCode MatLUFactorSymbolic_Matlab(Mat F, Mat A, IS r, IS c, const MatFactorInfo *info)
    160: {
    161:   PetscFunctionBegin;
    162:   PetscCheck(A->cmap->N == A->rmap->N, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "matrix must be square");
    163:   F->ops->lufactornumeric = MatLUFactorNumeric_Matlab;
    164:   F->assembled            = PETSC_TRUE;
    165:   PetscFunctionReturn(PETSC_SUCCESS);
    166: }

    168: static PetscErrorCode MatFactorGetSolverType_seqaij_matlab(Mat A, MatSolverType *type)
    169: {
    170:   PetscFunctionBegin;
    171:   *type = MATSOLVERMATLAB;
    172:   PetscFunctionReturn(PETSC_SUCCESS);
    173: }

    175: static PetscErrorCode MatDestroy_matlab(Mat A)
    176: {
    177:   const char *_A;

    179:   PetscFunctionBegin;
    180:   PetscCall(PetscObjectGetName((PetscObject)A, &_A));
    181:   PetscCall(PetscMatlabEngineEvaluate(PETSC_MATLAB_ENGINE_(PetscObjectComm((PetscObject)A)), "delete %s l_%s u_%s;", _A, _A, _A));
    182:   PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatFactorGetSolverType_C", NULL));
    183:   PetscFunctionReturn(PETSC_SUCCESS);
    184: }

    186: static PetscErrorCode MatGetFactor_seqaij_matlab(Mat A, MatFactorType ftype, Mat *F)
    187: {
    188:   PetscFunctionBegin;
    189:   PetscCheck(A->cmap->N == A->rmap->N, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "matrix must be square");
    190:   PetscCall(MatCreate(PetscObjectComm((PetscObject)A), F));
    191:   PetscCall(MatSetSizes(*F, A->rmap->n, A->cmap->n, A->rmap->n, A->cmap->n));
    192:   PetscCall(PetscStrallocpy("matlab", &((PetscObject)*F)->type_name));
    193:   PetscCall(MatSetUp(*F));

    195:   (*F)->ops->destroy           = MatDestroy_matlab;
    196:   (*F)->ops->getinfo           = MatGetInfo_External;
    197:   (*F)->trivialsymbolic        = PETSC_TRUE;
    198:   (*F)->ops->lufactorsymbolic  = MatLUFactorSymbolic_Matlab;
    199:   (*F)->ops->ilufactorsymbolic = MatLUFactorSymbolic_Matlab;

    201:   PetscCall(PetscObjectComposeFunction((PetscObject)*F, "MatFactorGetSolverType_C", MatFactorGetSolverType_seqaij_matlab));

    203:   (*F)->factortype = ftype;
    204:   PetscCall(PetscFree((*F)->solvertype));
    205:   PetscCall(PetscStrallocpy(MATSOLVERMATLAB, &(*F)->solvertype));
    206:   PetscFunctionReturn(PETSC_SUCCESS);
    207: }

    209: PETSC_INTERN PetscErrorCode MatSolverTypeRegister_Matlab(void)
    210: {
    211:   PetscFunctionBegin;
    212:   PetscCall(MatSolverTypeRegister(MATSOLVERMATLAB, MATSEQAIJ, MAT_FACTOR_LU, MatGetFactor_seqaij_matlab));
    213:   PetscFunctionReturn(PETSC_SUCCESS);
    214: }

    216: /*MC
    217:   MATSOLVERMATLAB - "matlab" - Providing direct solver LU for `MATSEQAIJ` matrix via the external package MATLAB.

    219:   Use `./configure` with the options `--with-matlab` to install PETSc with this capability

    221:   Options Database Key:
    222: . -pc_factor_mat_solver_type matlab - selects MATLAB to do the sparse factorization

    224:   Level: beginner

    226: .seealso: [](ch_matrices), `Mat`, `PCLU`, `PCFactorSetMatSolverType()`, `MatSolverType`
    227: M*/