Actual source code: aijmatlab.c

  1: #define PETSCMAT_DLL

  3: /* 
  4:         Provides an interface for the Matlab engine sparse solver

  6: */
 7:  #include src/mat/impls/aij/seq/aij.h

  9: #include "engine.h"   /* Matlab include file */
 10: #include "mex.h"      /* Matlab include file */

 12: typedef struct {
 13:   PetscErrorCode (*MatDuplicate)(Mat,MatDuplicateOption,Mat*);
 14:   PetscErrorCode (*MatView)(Mat,PetscViewer);
 15:   PetscErrorCode (*MatLUFactorSymbolic)(Mat,IS,IS,MatFactorInfo*,Mat*);
 16:   PetscErrorCode (*MatILUDTFactor)(Mat,IS,IS,MatFactorInfo*,Mat*);
 17:   PetscErrorCode (*MatDestroy)(Mat);
 18: } Mat_Matlab;


 24: PetscErrorCode  MatMatlabEnginePut_Matlab(PetscObject obj,void *mengine)
 25: {
 27:   Mat            B = (Mat)obj;
 28:   mxArray        *mat;
 29:   Mat_SeqAIJ     *aij = (Mat_SeqAIJ*)B->data;

 32:   mat  = mxCreateSparse(B->cmap.n,B->rmap.n,aij->nz,mxREAL);
 33:   PetscMemcpy(mxGetPr(mat),aij->a,aij->nz*sizeof(PetscScalar));
 34:   /* Matlab stores by column, not row so we pass in the transpose of the matrix */
 35:   PetscMemcpy(mxGetIr(mat),aij->j,aij->nz*sizeof(int));
 36:   PetscMemcpy(mxGetJc(mat),aij->i,(B->rmap.n+1)*sizeof(int));

 38:   /* Matlab indices start at 0 for sparse (what a surprise) */
 39: 
 40:   PetscObjectName(obj);
 41:   engPutVariable((Engine *)mengine,obj->name,mat);
 42:   return(0);
 43: }

 49: PetscErrorCode  MatMatlabEngineGet_Matlab(PetscObject obj,void *mengine)
 50: {
 52:   int            ii;
 53:   Mat            mat = (Mat)obj;
 54:   Mat_SeqAIJ     *aij = (Mat_SeqAIJ*)mat->data;
 55:   mxArray        *mmat;

 58:   MatSeqXAIJFreeAIJ(mat,&aij->a,&aij->j,&aij->i);

 60:   mmat = engGetVariable((Engine *)mengine,obj->name);

 62:   aij->nz           = (mxGetJc(mmat))[mat->rmap.n];
 63:   PetscMalloc3(aij->nz,PetscScalar,&aij->a,aij->nz,PetscInt,&aij->j,mat->rmap.n+1,PetscInt,&aij->i);
 64:   aij->singlemalloc = PETSC_TRUE;

 66:   PetscMemcpy(aij->a,mxGetPr(mmat),aij->nz*sizeof(PetscScalar));
 67:   /* Matlab stores by column, not row so we pass in the transpose of the matrix */
 68:   PetscMemcpy(aij->j,mxGetIr(mmat),aij->nz*sizeof(int));
 69:   PetscMemcpy(aij->i,mxGetJc(mmat),(mat->rmap.n+1)*sizeof(int));

 71:   for (ii=0; ii<mat->rmap.n; ii++) {
 72:     aij->ilen[ii] = aij->imax[ii] = aij->i[ii+1] - aij->i[ii];
 73:   }

 75:   MatAssemblyBegin(mat,MAT_FINAL_ASSEMBLY);
 76:   MatAssemblyEnd(mat,MAT_FINAL_ASSEMBLY);

 78:   return(0);
 79: }

 85: PetscErrorCode  MatConvert_Matlab_SeqAIJ(Mat A,MatType type,MatReuse reuse,Mat *newmat)
 86: {
 88:   Mat            B=*newmat;
 89:   Mat_Matlab    *lu=(Mat_Matlab*)A->spptr;

 92:   if (reuse == MAT_INITIAL_MATRIX) {
 93:     MatDuplicate(A,MAT_COPY_VALUES,&B);
 94:   }
 95:   B->ops->duplicate        = lu->MatDuplicate;
 96:   B->ops->view             = lu->MatView;
 97:   B->ops->lufactorsymbolic = lu->MatLUFactorSymbolic;
 98:   B->ops->iludtfactor      = lu->MatILUDTFactor;
 99:   B->ops->destroy          = lu->MatDestroy;
100: 
101:   PetscFree(lu);

103:   PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqaij_matlab_C","",PETSC_NULL);
104:   PetscObjectComposeFunction((PetscObject)B,"MatConvert_matlab_seqaij_C","",PETSC_NULL);
105:   PetscObjectComposeFunction((PetscObject)B,"PetscMatlabEnginePut_C","",PETSC_NULL);
106:   PetscObjectComposeFunction((PetscObject)B,"PetscMatlabEngineGet_C","",PETSC_NULL);

108:   PetscObjectChangeTypeName((PetscObject)B,MATSEQAIJ);
109:   *newmat = B;
110:   return(0);
111: }

116: PetscErrorCode MatDestroy_Matlab(Mat A)
117: {

121:   MatConvert_Matlab_SeqAIJ(A,MATSEQAIJ,MAT_REUSE_MATRIX,&A);
122:   (*A->ops->destroy)(A);
123:   return(0);
124: }

128: PetscErrorCode MatSolve_Matlab(Mat A,Vec b,Vec x)
129: {
131:   const char     *_A,*_b,*_x;

134:   /* make sure objects have names; use default if not */
135:   PetscObjectName((PetscObject)b);
136:   PetscObjectName((PetscObject)x);

138:   PetscObjectGetName((PetscObject)A,&_A);
139:   PetscObjectGetName((PetscObject)b,&_b);
140:   PetscObjectGetName((PetscObject)x,&_x);
141:   PetscMatlabEnginePut(PETSC_MATLAB_ENGINE_(A->comm),(PetscObject)b);
142:   PetscMatlabEngineEvaluate(PETSC_MATLAB_ENGINE_(A->comm),"%s = u%s\\(l%s\\(p%s*%s));",_x,_A,_A,_A,_b);
143:   PetscMatlabEngineEvaluate(PETSC_MATLAB_ENGINE_(A->comm),"%s = 0;",_b);
144:   /* PetscMatlabEnginePrintOutput(PETSC_MATLAB_ENGINE_(A->comm),stdout);  */
145:   PetscMatlabEngineGet(PETSC_MATLAB_ENGINE_(A->comm),(PetscObject)x);
146:   return(0);
147: }

151: PetscErrorCode MatLUFactorNumeric_Matlab(Mat A,MatFactorInfo *info,Mat *F)
152: {
154:   size_t         len;
155:   char           *_A,*name;

158:   PetscMatlabEnginePut(PETSC_MATLAB_ENGINE_(A->comm),(PetscObject)A);
159:   _A   = A->name;
160:   PetscMatlabEngineEvaluate(PETSC_MATLAB_ENGINE_(A->comm),"[l_%s,u_%s,p_%s] = lu(%s',%g);",_A,_A,_A,_A,info->dtcol);
161:   PetscMatlabEngineEvaluate(PETSC_MATLAB_ENGINE_(A->comm),"%s = 0;",_A);
162:   PetscStrlen(_A,&len);
163:   PetscMalloc((len+2)*sizeof(char),&name);
164:   sprintf(name,"_%s",_A);
165:   PetscObjectSetName((PetscObject)*F,name);
166:   PetscFree(name);
167:   return(0);
168: }

172: PetscErrorCode MatLUFactorSymbolic_Matlab(Mat A,IS r,IS c,MatFactorInfo *info,Mat *F)
173: {

177:   if (A->cmap.N != A->rmap.N) SETERRQ(PETSC_ERR_ARG_SIZ,"matrix must be square");
178:   MatCreate(A->comm,F);
179:   MatSetSizes(*F,A->rmap.n,A->cmap.n,A->rmap.n,A->cmap.n);
180:   MatSetType(*F,A->type_name);
181:   MatSeqAIJSetPreallocation(*F,0,PETSC_NULL);
182:   (*F)->ops->solve           = MatSolve_Matlab;
183:   (*F)->ops->lufactornumeric = MatLUFactorNumeric_Matlab;
184:   (*F)->factor               = FACTOR_LU;
185:   return(0);
186: }

188: /* ---------------------------------------------------------------------------------*/
191: PetscErrorCode MatSolve_Matlab_QR(Mat A,Vec b,Vec x)
192: {
194:   const char     *_A,*_b,*_x;

197:   /* make sure objects have names; use default if not */
198:   PetscObjectName((PetscObject)b);
199:   PetscObjectName((PetscObject)x);

201:   PetscObjectGetName((PetscObject)A,&_A);
202:   PetscObjectGetName((PetscObject)b,&_b);
203:   PetscObjectGetName((PetscObject)x,&_x);
204:   PetscMatlabEnginePut(PETSC_MATLAB_ENGINE_(A->comm),(PetscObject)b);
205:   PetscMatlabEngineEvaluate(PETSC_MATLAB_ENGINE_(A->comm),"%s = r%s\\(r%s'\\(%s*%s));",_x,_A,_A,_A+1,_b);
206:   PetscMatlabEngineEvaluate(PETSC_MATLAB_ENGINE_(A->comm),"%s = 0;",_b);
207:   /* PetscMatlabEnginePrintOutput(PETSC_MATLAB_ENGINE_(A->comm),stdout);  */
208:   PetscMatlabEngineGet(PETSC_MATLAB_ENGINE_(A->comm),(PetscObject)x);
209:   return(0);
210: }

214: PetscErrorCode MatLUFactorNumeric_Matlab_QR(Mat A,MatFactorInfo *info,Mat *F)
215: {
217:   size_t         len;
218:   char           *_A,*name;

221:   PetscMatlabEnginePut(PETSC_MATLAB_ENGINE_(A->comm),(PetscObject)A);
222:   _A   = A->name;
223:   PetscMatlabEngineEvaluate(PETSC_MATLAB_ENGINE_(A->comm),"r_%s = qr(%s');",_A,_A);
224:   PetscStrlen(_A,&len);
225:   PetscMalloc((len+2)*sizeof(char),&name);
226:   sprintf(name,"_%s",_A);
227:   PetscObjectSetName((PetscObject)*F,name);
228:   PetscFree(name);
229:   return(0);
230: }

234: PetscErrorCode MatLUFactorSymbolic_Matlab_QR(Mat A,IS r,IS c,MatFactorInfo *info,Mat *F)
235: {

239:   if (A->cmap.N != A->rmap.N) SETERRQ(PETSC_ERR_ARG_SIZ,"matrix must be square");
240:   MatCreate(A->comm,F);
241:   MatSetSizes(*F,A->rmap.n,A->cmap.n,A->rmap.n,A->cmap.n);
242:   MatSetType(*F,A->type_name);
243:   MatSeqAIJSetPreallocation(*F,0,PETSC_NULL);
244:   (*F)->ops->solve           = MatSolve_Matlab_QR;
245:   (*F)->ops->lufactornumeric = MatLUFactorNumeric_Matlab_QR;
246:   (*F)->factor               = FACTOR_LU;
247:   (*F)->assembled            = PETSC_TRUE;  /* required by -ksp_view */

249:   return(0);
250: }

252: /* --------------------------------------------------------------------------------*/
255: PetscErrorCode MatILUDTFactor_Matlab(Mat A,IS isrow,IS iscol,MatFactorInfo *info,Mat *F)
256: {
258:   size_t         len;
259:   char           *_A,*name;

262:   if (info->dt == PETSC_DEFAULT)      info->dt      = .005;
263:   if (info->dtcol == PETSC_DEFAULT)   info->dtcol   = .01;
264:   if (A->cmap.N != A->rmap.N) SETERRQ(PETSC_ERR_ARG_SIZ,"matrix must be square");
265:   MatCreate(A->comm,F);
266:   MatSetSizes(*F,A->rmap.n,A->cmap.n,A->rmap.n,A->cmap.n);
267:   MatSetType(*F,A->type_name);
268:   MatSeqAIJSetPreallocation(*F,0,PETSC_NULL);
269:   (*F)->ops->solve           = MatSolve_Matlab;
270:   (*F)->factor               = FACTOR_LU;
271:   PetscMatlabEnginePut(PETSC_MATLAB_ENGINE_(A->comm),(PetscObject)A);
272:   _A   = A->name;
273:   PetscMatlabEngineEvaluate(PETSC_MATLAB_ENGINE_(A->comm),"info_%s = struct('droptol',%g,'thresh',%g);",_A,info->dt,info->dtcol);
274:   PetscMatlabEngineEvaluate(PETSC_MATLAB_ENGINE_(A->comm),"[l_%s,u_%s,p_%s] = luinc(%s',info_%s);",_A,_A,_A,_A,_A);
275:   PetscMatlabEngineEvaluate(PETSC_MATLAB_ENGINE_(A->comm),"%s = 0;",_A);

277:   PetscStrlen(_A,&len);
278:   PetscMalloc((len+2)*sizeof(char),&name);
279:   sprintf(name,"_%s",_A);
280:   PetscObjectSetName((PetscObject)*F,name);
281:   PetscFree(name);
282:   return(0);
283: }

287: PetscErrorCode MatFactorInfo_Matlab(Mat A,PetscViewer viewer)
288: {
290: 
292:   PetscViewerASCIIPrintf(viewer,"Matlab run parameters:  -- not written yet!\n");
293:   return(0);
294: }

298: PetscErrorCode MatView_Matlab(Mat A,PetscViewer viewer)
299: {
300:   PetscErrorCode    ierr;
301:   PetscTruth        iascii;
302:   PetscViewerFormat format;
303:   Mat_Matlab        *lu=(Mat_Matlab*)(A->spptr);

306:   (*lu->MatView)(A,viewer);
307:   PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_ASCII,&iascii);
308:   if (iascii) {
309:     PetscViewerGetFormat(viewer,&format);
310:     if (format == PETSC_VIEWER_ASCII_FACTOR_INFO) {
311:       MatFactorInfo_Matlab(A,viewer);
312:     }
313:   }
314:   return(0);
315: }

319: PetscErrorCode MatDuplicate_Matlab(Mat A, MatDuplicateOption op, Mat *M)
320: {
322:   Mat_Matlab     *lu=(Mat_Matlab*)A->spptr;

325:   (*lu->MatDuplicate)(A,op,M);
326:   PetscMemcpy((*M)->spptr,lu,sizeof(Mat_Matlab));
327:   return(0);
328: }

333: PetscErrorCode  MatConvert_SeqAIJ_Matlab(Mat A,MatType type,MatReuse reuse,Mat *newmat)
334: {
336:   Mat            B=*newmat;
337:   Mat_Matlab     *lu;
338:   PetscTruth     qr;

341:   if (reuse == MAT_INITIAL_MATRIX) {
342:     MatDuplicate(A,MAT_COPY_VALUES,&B);
343:   }

345:   PetscNew(Mat_Matlab,&lu);
346:   lu->MatDuplicate         = A->ops->duplicate;
347:   lu->MatView              = A->ops->view;
348:   lu->MatLUFactorSymbolic  = A->ops->lufactorsymbolic;
349:   lu->MatILUDTFactor       = A->ops->iludtfactor;
350:   lu->MatDestroy           = A->ops->destroy;

352:   B->spptr                 = (void*)lu;
353:   B->ops->duplicate        = MatDuplicate_Matlab;
354:   B->ops->view             = MatView_Matlab;
355:   B->ops->lufactorsymbolic = MatLUFactorSymbolic_Matlab;
356:   B->ops->iludtfactor      = MatILUDTFactor_Matlab;
357:   B->ops->destroy          = MatDestroy_Matlab;

359:   PetscOptionsHasName(A->prefix,"-mat_matlab_qr",&qr);
360:   if (qr) {
361:     B->ops->lufactorsymbolic = MatLUFactorSymbolic_Matlab_QR;
362:     PetscInfo(0,"Using Matlab QR with iterative refinement for LU factorization and solves\n");
363:   } else {
364:     PetscInfo(0,"Using Matlab for LU factorizations and solves.\n");
365:   }
366:   PetscInfo(0,"Using Matlab for ILUDT factorizations and solves.\n");

368:   PetscObjectComposeFunctionDynamic((PetscObject)B,"MatConvert_seqaij_matlab_C",
369:                                            "MatConvert_SeqAIJ_Matlab",MatConvert_SeqAIJ_Matlab);
370:   PetscObjectComposeFunctionDynamic((PetscObject)B,"MatConvert_matlab_seqaij_C",
371:                                            "MatConvert_Matlab_SeqAIJ",MatConvert_Matlab_SeqAIJ);
372:   PetscObjectComposeFunctionDynamic((PetscObject)B,"PetscMatlabEnginePut_C",
373:                                            "MatMatlabEnginePut_Matlab",MatMatlabEnginePut_Matlab);
374:   PetscObjectComposeFunctionDynamic((PetscObject)B,"PetscMatlabEngineGet_C",
375:                                            "MatMatlabEngineGet_Matlab",MatMatlabEngineGet_Matlab);
376:   PetscObjectChangeTypeName((PetscObject)B,MATMATLAB);
377:   *newmat = B;
378:   return(0);
379: }

382: /*MC
383:   MATMATLAB - MATMATLAB = "matlab" - A matrix type providing direct solvers (LU and QR) and drop tolerance
384:   based ILU factorization (ILUDT) for sequential matrices via the external package Matlab.

386:   If Matlab is instaled (see the manual for
387:   instructions on how to declare the existence of external packages),
388:   a matrix type can be constructed which invokes Matlab solvers.
389:   After calling MatCreate(...,A), simply call MatSetType(A,MATMATLAB).
390:   This matrix type is only supported for double precision real.

392:   This matrix inherits from MATSEQAIJ.  As a result, MatSeqAIJSetPreallocation is 
393:   supported for this matrix type.  One can also call MatConvert for an inplace conversion to or from 
394:   the MATSEQAIJ type without data copy.

396:   Options Database Keys:
397: + -mat_type matlab - sets the matrix type to "matlab" during a call to MatSetFromOptions()
398: - -mat_matlab_qr   - sets the direct solver to be QR instead of LU

400:   Level: beginner

402: .seealso: PCLU
403: M*/

408: PetscErrorCode  MatCreate_Matlab(Mat A)
409: {

413:   MatSetType(A,MATSEQAIJ);
414:   MatConvert_SeqAIJ_Matlab(A,MATMATLAB,MAT_REUSE_MATRIX,&A);
415:   return(0);
416: }