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 PETSCMAT_DLLEXPORT 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->n,B->m,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->m+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 PETSCMAT_DLLEXPORT 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(aij->singlemalloc,&aij->a,&aij->j,&aij->i);

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

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

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

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

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

     79:   return(0);
     80: }

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

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

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

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

    117: PetscErrorCode MatDestroy_Matlab(Mat A)
    118: {

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

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

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

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

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

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

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

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

    189: /* ---------------------------------------------------------------------------------*/
    192: PetscErrorCode MatSolve_Matlab_QR(Mat A,Vec b,Vec x)
    193: {
    195:   const char     *_A,*_b,*_x;

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

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

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

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

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

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

    250:   return(0);
    251: }

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

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

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

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

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

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

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

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

    334: PetscErrorCode PETSCMAT_DLLEXPORT MatConvert_SeqAIJ_Matlab(Mat A,MatType type,MatReuse reuse,Mat *newmat)
    335: {
    336:   /* This routine is only called to convert to MATMATLAB */
    337:   /* from MATSEQAIJ, so we will ignore 'MatType type'. */
    339:   Mat            B=*newmat;
    340:   Mat_Matlab     *lu;
    341:   PetscTruth     qr;

    344:   if (reuse == MAT_INITIAL_MATRIX) {
    345:     MatDuplicate(A,MAT_COPY_VALUES,&B);
    346:   }

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

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

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

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

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

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

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

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

    403:   Level: beginner

    405: .seealso: PCLU
    406: M*/

    411: PetscErrorCode PETSCMAT_DLLEXPORT MatCreate_Matlab(Mat A)
    412: {

    416:   PetscObjectChangeTypeName((PetscObject)A,MATMATLAB);
    417:   MatSetType(A,MATSEQAIJ);
    418:   MatConvert_SeqAIJ_Matlab(A,MATMATLAB,MAT_REUSE_MATRIX,&A);
    419:   return(0);
    420: }