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:   mat = MatSeqAIJToMatlab((Mat)obj);
 34:   PetscObjectName(obj);
 35:   engPutVariable((Engine *)mengine, obj->name, mat);
 36:   return 0;
 37: }

 39: PETSC_EXTERN 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:   nn  = mxGetN(mmat); /* rows of transpose of matrix */
 46:   nm  = mxGetM(mmat);
 47:   nnz = (mxGetJc(mmat))[nn];
 48:   ii  = mxGetJc(mmat);
 49:   jj  = mxGetIr(mmat);
 50:   n   = (PetscInt)nn;
 51:   m   = (PetscInt)nm;
 52:   nz  = (PetscInt)nnz;

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

 68:     aij->singlemalloc = PETSC_TRUE;
 69:   }

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

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

 80:   mat->nonzerostate++; /* since the nonzero structure can change anytime force the Inode information to always be rebuilt */
 81:   MatAssemblyBegin(mat, MAT_FINAL_ASSEMBLY);
 82:   MatAssemblyEnd(mat, MAT_FINAL_ASSEMBLY);
 83:   return 0;
 84: }

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

 91:   mmat = engGetVariable((Engine *)mengine, obj->name);
 92:   MatSeqAIJFromMatlab(mmat, mat);
 93:   return 0;
 94: }

 96: PetscErrorCode MatSolve_Matlab(Mat A, Vec b, Vec x)
 97: {
 98:   const char *_A, *_b, *_x;

100:   /* make sure objects have names; use default if not */
101:   PetscObjectName((PetscObject)b);
102:   PetscObjectName((PetscObject)x);

104:   PetscObjectGetName((PetscObject)A, &_A);
105:   PetscObjectGetName((PetscObject)b, &_b);
106:   PetscObjectGetName((PetscObject)x, &_x);
107:   PetscMatlabEnginePut(PETSC_MATLAB_ENGINE_(PetscObjectComm((PetscObject)A)), (PetscObject)b);
108:   PetscMatlabEngineEvaluate(PETSC_MATLAB_ENGINE_(PetscObjectComm((PetscObject)A)), "%s = u%s\\(l%s\\(p%s*%s));", _x, _A, _A, _A, _b);
109:   PetscMatlabEngineEvaluate(PETSC_MATLAB_ENGINE_(PetscObjectComm((PetscObject)A)), "%s = 0;", _b);
110:   /* PetscMatlabEnginePrintOutput(PETSC_MATLAB_ENGINE_(PetscObjectComm((PetscObject)A)),stdout);  */
111:   PetscMatlabEngineGet(PETSC_MATLAB_ENGINE_(PetscObjectComm((PetscObject)A)), (PetscObject)x);
112:   return 0;
113: }

115: PetscErrorCode MatLUFactorNumeric_Matlab(Mat F, Mat A, const MatFactorInfo *info)
116: {
117:   size_t    len;
118:   char     *_A, *name;
119:   PetscReal dtcol = info->dtcol;

121:   if (F->factortype == MAT_FACTOR_ILU || info->dt > 0) {
122:     /* the ILU form is not currently registered */
123:     if (info->dtcol == PETSC_DEFAULT) dtcol = .01;
124:     F->ops->solve = MatSolve_Matlab;
125:     F->factortype = MAT_FACTOR_LU;

127:     PetscMatlabEnginePut(PETSC_MATLAB_ENGINE_(PetscObjectComm((PetscObject)A)), (PetscObject)A);
128:     _A = ((PetscObject)A)->name;
129:     PetscMatlabEngineEvaluate(PETSC_MATLAB_ENGINE_(PetscObjectComm((PetscObject)A)), "info_%s = struct('droptol',%g,'thresh',%g);", _A, info->dt, dtcol);
130:     PetscMatlabEngineEvaluate(PETSC_MATLAB_ENGINE_(PetscObjectComm((PetscObject)A)), "[l_%s,u_%s,p_%s] = luinc(%s',info_%s);", _A, _A, _A, _A, _A);
131:     PetscMatlabEngineEvaluate(PETSC_MATLAB_ENGINE_(PetscObjectComm((PetscObject)A)), "%s = 0;", _A);

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

149:     F->ops->solve = MatSolve_Matlab;
150:   }
151:   return 0;
152: }

154: PetscErrorCode MatLUFactorSymbolic_Matlab(Mat F, Mat A, IS r, IS c, const MatFactorInfo *info)
155: {
157:   F->ops->lufactornumeric = MatLUFactorNumeric_Matlab;
158:   F->assembled            = PETSC_TRUE;
159:   return 0;
160: }

162: PetscErrorCode MatFactorGetSolverType_seqaij_matlab(Mat A, MatSolverType *type)
163: {
164:   *type = MATSOLVERMATLAB;
165:   return 0;
166: }

168: PetscErrorCode MatDestroy_matlab(Mat A)
169: {
170:   const char *_A;

172:   PetscObjectGetName((PetscObject)A, &_A);
173:   PetscMatlabEngineEvaluate(PETSC_MATLAB_ENGINE_(PetscObjectComm((PetscObject)A)), "delete %s l_%s u_%s;", _A, _A, _A);
174:   PetscObjectComposeFunction((PetscObject)A, "MatFactorGetSolverType_C", NULL);
175:   return 0;
176: }

178: PETSC_EXTERN PetscErrorCode MatGetFactor_seqaij_matlab(Mat A, MatFactorType ftype, Mat *F)
179: {
181:   MatCreate(PetscObjectComm((PetscObject)A), F);
182:   MatSetSizes(*F, A->rmap->n, A->cmap->n, A->rmap->n, A->cmap->n);
183:   PetscStrallocpy("matlab", &((PetscObject)*F)->type_name);
184:   MatSetUp(*F);

186:   (*F)->ops->destroy           = MatDestroy_matlab;
187:   (*F)->ops->getinfo           = MatGetInfo_External;
188:   (*F)->trivialsymbolic        = PETSC_TRUE;
189:   (*F)->ops->lufactorsymbolic  = MatLUFactorSymbolic_Matlab;
190:   (*F)->ops->ilufactorsymbolic = MatLUFactorSymbolic_Matlab;

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

194:   (*F)->factortype = ftype;
195:   PetscFree((*F)->solvertype);
196:   PetscStrallocpy(MATSOLVERMATLAB, &(*F)->solvertype);
197:   return 0;
198: }

200: PETSC_EXTERN PetscErrorCode MatSolverTypeRegister_Matlab(void)
201: {
202:   MatSolverTypeRegister(MATSOLVERMATLAB, MATSEQAIJ, MAT_FACTOR_LU, MatGetFactor_seqaij_matlab);
203:   return 0;
204: }

206: /* --------------------------------------------------------------------------------*/

208: PetscErrorCode MatView_Info_Matlab(Mat A, PetscViewer viewer)
209: {
210:   PetscViewerASCIIPrintf(viewer, "MATLAB run parameters:  -- not written yet!\n");
211:   return 0;
212: }

214: PetscErrorCode MatView_Matlab(Mat A, PetscViewer viewer)
215: {
216:   PetscBool iascii;

218:   MatView_SeqAIJ(A, viewer);
219:   PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii);
220:   if (iascii) {
221:     PetscViewerFormat format;

223:     PetscViewerGetFormat(viewer, &format);
224:     if (format == PETSC_VIEWER_ASCII_FACTOR_INFO) MatView_Info_Matlab(A, viewer);
225:   }
226:   return 0;
227: }

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

232:   Options Database Keys:
233: . -pc_factor_mat_solver_type matlab - selects MATLAB to do the sparse factorization

235:   Note:
236:     You must ./configure with the options --with-matlab

238:   Level: beginner

240: .seealso: `PCLU`, `PCFactorSetMatSolverType()`, `MatSolverType`
241: M*/