Actual source code: aijmatlab.c


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

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

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

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

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

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

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

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

 56:   if (mat->rmap->n < 0 && mat->cmap->n < 0) {
 57:     /* matrix has not yet had its size set */
 58:     PetscCall(MatSetSizes(mat, n, m, PETSC_DETERMINE, PETSC_DETERMINE));
 59:     PetscCall(MatSetUp(mat));
 60:   } else {
 61:     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);
 62:     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);
 63:   }
 64:   if (nz != aij->nz) {
 65:     /* number of nonzeros in matrix has changed, so need new data structure */
 66:     PetscCall(MatSeqXAIJFreeAIJ(mat, &aij->a, &aij->j, &aij->i));
 67:     aij->nz = nz;
 68:     PetscCall(PetscMalloc3(aij->nz, &aij->a, aij->nz, &aij->j, mat->rmap->n + 1, &aij->i));

 70:     aij->singlemalloc = 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: 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: 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: 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: PetscErrorCode MatFactorGetSolverType_seqaij_matlab(Mat A, MatSolverType *type)
169: {
170:   PetscFunctionBegin;
171:   *type = MATSOLVERMATLAB;
172:   PetscFunctionReturn(PETSC_SUCCESS);
173: }

175: 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: PETSC_EXTERN 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_EXTERN PetscErrorCode MatSolverTypeRegister_Matlab(void)
210: {
211:   PetscFunctionBegin;
212:   PetscCall(MatSolverTypeRegister(MATSOLVERMATLAB, MATSEQAIJ, MAT_FACTOR_LU, MatGetFactor_seqaij_matlab));
213:   PetscFunctionReturn(PETSC_SUCCESS);
214: }

216: PetscErrorCode MatView_Info_Matlab(Mat A, PetscViewer viewer)
217: {
218:   PetscFunctionBegin;
219:   PetscCall(PetscViewerASCIIPrintf(viewer, "MATLAB run parameters:  -- not written yet!\n"));
220:   PetscFunctionReturn(PETSC_SUCCESS);
221: }

223: PetscErrorCode MatView_Matlab(Mat A, PetscViewer viewer)
224: {
225:   PetscBool iascii;

227:   PetscFunctionBegin;
228:   PetscCall(MatView_SeqAIJ(A, viewer));
229:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
230:   if (iascii) {
231:     PetscViewerFormat format;

233:     PetscCall(PetscViewerGetFormat(viewer, &format));
234:     if (format == PETSC_VIEWER_ASCII_FACTOR_INFO) PetscCall(MatView_Info_Matlab(A, viewer));
235:   }
236:   PetscFunctionReturn(PETSC_SUCCESS);
237: }

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

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

244:   Options Database Key:
245: . -pc_factor_mat_solver_type matlab - selects MATLAB to do the sparse factorization

247:   Level: beginner

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