Actual source code: aijmatlab.c

petsc-3.3-p7 2013-05-11
  2: /* 
  3:         Provides an interface for the MATLAB engine sparse solver

  5: */
  6: #include <../src/mat/impls/aij/seq/aij.h>

  8: #include <engine.h>   /* MATLAB include file */
  9: #include <mex.h>      /* MATLAB include file */

 11: EXTERN_C_BEGIN
 14: mxArray *MatSeqAIJToMatlab(Mat B)
 15: {
 17:   Mat_SeqAIJ     *aij = (Mat_SeqAIJ*)B->data;
 18:   mwIndex        *ii,*jj;
 19:   mxArray        *mat;
 20:   PetscInt       i;

 23:   mat  = mxCreateSparse(B->cmap->n,B->rmap->n,aij->nz,mxREAL);
 24:   PetscMemcpy(mxGetPr(mat),aij->a,aij->nz*sizeof(PetscScalar));
 25:   /* MATLAB stores by column, not row so we pass in the transpose of the matrix */
 26:   jj = mxGetIr(mat);
 27:   for (i=0; i<aij->nz; i++) jj[i] = aij->j[i];
 28:   ii = mxGetJc(mat);
 29:   for (i=0; i<B->rmap->n+1; i++) ii[i] = aij->i[i];
 30: 
 31:   PetscFunctionReturn(mat);
 32: }
 33: EXTERN_C_END


 36: EXTERN_C_BEGIN
 39: PetscErrorCode  MatlabEnginePut_SeqAIJ(PetscObject obj,void *mengine)
 40: {
 42:   mxArray        *mat;

 45:   mat = MatSeqAIJToMatlab((Mat)obj);if (!mat) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"Cannot create MATLAB matrix");
 46:   PetscObjectName(obj);
 47:   engPutVariable((Engine *)mengine,obj->name,mat);
 48:   return(0);
 49: }
 50: EXTERN_C_END

 52: EXTERN_C_BEGIN
 55: /*@C
 56:     MatSeqAIJFromMatlab - Given a MATLAB sparse matrix, fills a SeqAIJ matrix with its transpose.

 58:    Not Collective

 60:    Input Parameters:
 61: +     mmat - a MATLAB sparse matris
 62: -     mat - a already created MATSEQAIJ

 64:   Level: intermediate

 66: @*/
 67: PetscErrorCode  MatSeqAIJFromMatlab(mxArray *mmat,Mat mat)
 68: {
 70:   PetscInt       nz,n,m,*i,*j,k;
 71:   mwIndex        nnz,nn,nm,*ii,*jj;
 72:   Mat_SeqAIJ     *aij = (Mat_SeqAIJ*)mat->data;

 75:   nn  = mxGetN(mmat);   /* rows of transpose of matrix */
 76:   nm  = mxGetM(mmat);
 77:   nnz = (mxGetJc(mmat))[nn];
 78:   ii  = mxGetJc(mmat);
 79:   jj  = mxGetIr(mmat);
 80:   n   = (PetscInt) nn;
 81:   m   = (PetscInt) nm;
 82:   nz  = (PetscInt) nnz;

 84:   if (mat->rmap->n < 0 && mat->cmap->n < 0) {
 85:     /* matrix has not yet had its size set */
 86:     MatSetSizes(mat,n,m,PETSC_DETERMINE,PETSC_DETERMINE);
 87:     MatSetUp(mat);
 88:   } else {
 89:     if (mat->rmap->n != n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_SUP,"Cannot change size of PETSc matrix %D to %D",mat->rmap->n,n);
 90:     if (mat->cmap->n != m) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_SUP,"Cannot change size of PETSc matrix %D to %D",mat->cmap->n,m);
 91:   }
 92:   if (nz != aij->nz) {
 93:     /* number of nonzeros in matrix has changed, so need new data structure */
 94:     MatSeqXAIJFreeAIJ(mat,&aij->a,&aij->j,&aij->i);
 95:     aij->nz           = nz;
 96:     PetscMalloc3(aij->nz,PetscScalar,&aij->a,aij->nz,PetscInt,&aij->j,mat->rmap->n+1,PetscInt,&aij->i);
 97:     aij->singlemalloc = PETSC_TRUE;
 98:   }

100:   PetscMemcpy(aij->a,mxGetPr(mmat),aij->nz*sizeof(PetscScalar));
101:   /* MATLAB stores by column, not row so we pass in the transpose of the matrix */
102:   i = aij->i;
103:   for (k=0; k<n+1; k++) {
104:     i[k] = (PetscInt) ii[k];
105:   }
106:   j = aij->j;
107:   for (k=0; k<nz; k++) {
108:     j[k] = (PetscInt) jj[k];
109:   }

111:   for (k=0; k<mat->rmap->n; k++) {
112:     aij->ilen[k] = aij->imax[k] = aij->i[k+1] - aij->i[k];
113:   }

115:   MatAssemblyBegin(mat,MAT_FINAL_ASSEMBLY);
116:   MatAssemblyEnd(mat,MAT_FINAL_ASSEMBLY);
117:   return(0);
118: }
119: EXTERN_C_END


122: EXTERN_C_BEGIN
125: PetscErrorCode  MatlabEngineGet_SeqAIJ(PetscObject obj,void *mengine)
126: {
128:   Mat            mat = (Mat)obj;
129:   mxArray        *mmat;

132:   mmat = engGetVariable((Engine *)mengine,obj->name);
133:   MatSeqAIJFromMatlab(mmat,mat);
134:   return(0);
135: }
136: EXTERN_C_END

140: PetscErrorCode MatSolve_Matlab(Mat A,Vec b,Vec x)
141: {
143:   const char     *_A,*_b,*_x;

146:   /* make sure objects have names; use default if not */
147:   PetscObjectName((PetscObject)b);
148:   PetscObjectName((PetscObject)x);

150:   PetscObjectGetName((PetscObject)A,&_A);
151:   PetscObjectGetName((PetscObject)b,&_b);
152:   PetscObjectGetName((PetscObject)x,&_x);
153:   PetscMatlabEnginePut(PETSC_MATLAB_ENGINE_(((PetscObject)A)->comm),(PetscObject)b);
154:   PetscMatlabEngineEvaluate(PETSC_MATLAB_ENGINE_(((PetscObject)A)->comm),"%s = u%s\\(l%s\\(p%s*%s));",_x,_A,_A,_A,_b);
155:   PetscMatlabEngineEvaluate(PETSC_MATLAB_ENGINE_(((PetscObject)A)->comm),"%s = 0;",_b);
156:   /* PetscMatlabEnginePrintOutput(PETSC_MATLAB_ENGINE_(((PetscObject)A)->comm),stdout);  */
157:   PetscMatlabEngineGet(PETSC_MATLAB_ENGINE_(((PetscObject)A)->comm),(PetscObject)x);
158:   return(0);
159: }

163: PetscErrorCode MatLUFactorNumeric_Matlab(Mat F,Mat A,const MatFactorInfo *info)
164: {
166:   size_t         len;
167:   char           *_A,*name;
168:   PetscReal      dtcol = info->dtcol;

171:   if (F->factortype == MAT_FACTOR_ILU || info->dt > 0) {
172:     if (info->dtcol == PETSC_DEFAULT)  dtcol = .01;
173:     F->ops->solve           = MatSolve_Matlab;
174:     F->factortype           = MAT_FACTOR_LU;
175:     PetscMatlabEnginePut(PETSC_MATLAB_ENGINE_(((PetscObject)A)->comm),(PetscObject)A);
176:     _A   = ((PetscObject)A)->name;
177:     PetscMatlabEngineEvaluate(PETSC_MATLAB_ENGINE_(((PetscObject)A)->comm),"info_%s = struct('droptol',%g,'thresh',%g);",_A,info->dt,dtcol);
178:     PetscMatlabEngineEvaluate(PETSC_MATLAB_ENGINE_(((PetscObject)A)->comm),"[l_%s,u_%s,p_%s] = luinc(%s',info_%s);",_A,_A,_A,_A,_A);
179:     PetscMatlabEngineEvaluate(PETSC_MATLAB_ENGINE_(((PetscObject)A)->comm),"%s = 0;",_A);

181:     PetscStrlen(_A,&len);
182:     PetscMalloc((len+2)*sizeof(char),&name);
183:     sprintf(name,"_%s",_A);
184:     PetscObjectSetName((PetscObject)F,name);
185:     PetscFree(name);
186:   } else {
187:     PetscMatlabEnginePut(PETSC_MATLAB_ENGINE_(((PetscObject)A)->comm),(PetscObject)A);
188:     _A   = ((PetscObject)A)->name;
189:     PetscMatlabEngineEvaluate(PETSC_MATLAB_ENGINE_(((PetscObject)A)->comm),"[l_%s,u_%s,p_%s] = lu(%s',%g);",_A,_A,_A,_A,dtcol);
190:     PetscMatlabEngineEvaluate(PETSC_MATLAB_ENGINE_(((PetscObject)A)->comm),"%s = 0;",_A);
191:     PetscStrlen(_A,&len);
192:     PetscMalloc((len+2)*sizeof(char),&name);
193:     sprintf(name,"_%s",_A);
194:     PetscObjectSetName((PetscObject)F,name);
195:     PetscFree(name);
196:     F->ops->solve              = MatSolve_Matlab;
197:   }
198:   return(0);
199: }

203: PetscErrorCode MatLUFactorSymbolic_Matlab(Mat F,Mat A,IS r,IS c,const MatFactorInfo *info)
204: {
206:   if (A->cmap->N != A->rmap->N) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"matrix must be square");
207:   F->ops->lufactornumeric    = MatLUFactorNumeric_Matlab;
208:   F->assembled = PETSC_TRUE;
209:   return(0);
210: }

212: EXTERN_C_BEGIN
215: PetscErrorCode MatFactorGetSolverPackage_seqaij_matlab(Mat A,const MatSolverPackage *type)
216: {
218:   *type = MATSOLVERMATLAB;
219:   return(0);
220: }
221: EXTERN_C_END

223: EXTERN_C_BEGIN
226: PetscErrorCode MatGetFactor_seqaij_matlab(Mat A,MatFactorType ftype,Mat *F)
227: {

231:   if (A->cmap->N != A->rmap->N) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"matrix must be square");
232:   MatCreate(((PetscObject)A)->comm,F);
233:   MatSetSizes(*F,A->rmap->n,A->cmap->n,A->rmap->n,A->cmap->n);
234:   MatSetType(*F,((PetscObject)A)->type_name);
235:   MatSeqAIJSetPreallocation(*F,0,PETSC_NULL);
236:   (*F)->ops->lufactorsymbolic  = MatLUFactorSymbolic_Matlab;
237:   (*F)->ops->ilufactorsymbolic = MatLUFactorSymbolic_Matlab;
238:   PetscObjectComposeFunctionDynamic((PetscObject)(*F),"MatFactorGetSolverPackage_C","MatFactorGetSolverPackage_seqaij_matlab",MatFactorGetSolverPackage_seqaij_matlab);

240:   (*F)->factortype             = ftype;
241:   return(0);
242: }
243: EXTERN_C_END

245: /* --------------------------------------------------------------------------------*/

249: PetscErrorCode MatFactorInfo_Matlab(Mat A,PetscViewer viewer)
250: {
252: 
254:   PetscViewerASCIIPrintf(viewer,"MATLAB run parameters:  -- not written yet!\n");
255:   return(0);
256: }

260: PetscErrorCode MatView_Matlab(Mat A,PetscViewer viewer)
261: {
262:   PetscErrorCode    ierr;
263:   PetscBool         iascii;
264:   PetscViewerFormat format;

267:   MatView_SeqAIJ(A,viewer);
268:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
269:   if (iascii) {
270:     PetscViewerGetFormat(viewer,&format);
271:     if (format == PETSC_VIEWER_ASCII_FACTOR_INFO) {
272:       MatFactorInfo_Matlab(A,viewer);
273:     }
274:   }
275:   return(0);
276: }


279: /*MC
280:   MATSOLVERMATLAB - "matlab" - Providing direct solvers (LU and QR) and drop tolerance
281:   based ILU factorization (ILUDT) for sequential matrices via the external package MATLAB.


284:   Works with MATSEQAIJ matrices.

286:   Options Database Keys:
287: . -pc_factor_mat_solver_type matlab - selects MATLAB to do the sparse factorization


290:   Level: beginner

292: .seealso: PCLU

294: .seealso: PCFactorSetMatSolverPackage(), MatSolverPackage
295: M*/