Actual source code: aijmatlab.c

petsc-3.7.3 2016-08-01
Report Typos and Errors
  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 */

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

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

 34: PETSC_EXTERN PetscErrorCode MatlabEnginePut_SeqAIJ(PetscObject obj,void *mengine)
 35: {
 37:   mxArray        *mat;

 40:   mat  = MatSeqAIJToMatlab((Mat)obj);if (!mat) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"Cannot create MATLAB matrix");
 41:   PetscObjectName(obj);
 42:   engPutVariable((Engine*)mengine,obj->name,mat);
 43:   return(0);
 44: }

 48: /*@C
 49:     MatSeqAIJFromMatlab - Given a MATLAB sparse matrix, fills a SeqAIJ matrix with its transpose.

 51:    Not Collective

 53:    Input Parameters:
 54: +     mmat - a MATLAB sparse matris
 55: -     mat - an already created MATSEQAIJ

 57:   Level: intermediate

 59: @*/
 60: PETSC_EXTERN PetscErrorCode MatSeqAIJFromMatlab(mxArray *mmat,Mat mat)
 61: {
 63:   PetscInt       nz,n,m,*i,*j,k;
 64:   mwIndex        nnz,nn,nm,*ii,*jj;
 65:   Mat_SeqAIJ     *aij = (Mat_SeqAIJ*)mat->data;

 68:   nn  = mxGetN(mmat);   /* rows of transpose of matrix */
 69:   nm  = mxGetM(mmat);
 70:   nnz = (mxGetJc(mmat))[nn];
 71:   ii  = mxGetJc(mmat);
 72:   jj  = mxGetIr(mmat);
 73:   n   = (PetscInt) nn;
 74:   m   = (PetscInt) nm;
 75:   nz  = (PetscInt) nnz;

 77:   if (mat->rmap->n < 0 && mat->cmap->n < 0) {
 78:     /* matrix has not yet had its size set */
 79:     MatSetSizes(mat,n,m,PETSC_DETERMINE,PETSC_DETERMINE);
 80:     MatSetUp(mat);
 81:   } else {
 82:     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);
 83:     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);
 84:   }
 85:   if (nz != aij->nz) {
 86:     /* number of nonzeros in matrix has changed, so need new data structure */
 87:     MatSeqXAIJFreeAIJ(mat,&aij->a,&aij->j,&aij->i);
 88:     aij->nz = nz;
 89:     PetscMalloc3(aij->nz,&aij->a,aij->nz,&aij->j,mat->rmap->n+1,&aij->i);

 91:     aij->singlemalloc = PETSC_TRUE;
 92:   }

 94:   PetscMemcpy(aij->a,mxGetPr(mmat),aij->nz*sizeof(PetscScalar));
 95:   /* MATLAB stores by column, not row so we pass in the transpose of the matrix */
 96:   i = aij->i;
 97:   for (k=0; k<n+1; k++) i[k] = (PetscInt) ii[k];
 98:   j = aij->j;
 99:   for (k=0; k<nz; k++) j[k] = (PetscInt) jj[k];

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

103:   MatAssemblyBegin(mat,MAT_FINAL_ASSEMBLY);
104:   MatAssemblyEnd(mat,MAT_FINAL_ASSEMBLY);
105:   return(0);
106: }

110: PETSC_EXTERN PetscErrorCode  MatlabEngineGet_SeqAIJ(PetscObject obj,void *mengine)
111: {
113:   Mat            mat = (Mat)obj;
114:   mxArray        *mmat;

117:   mmat = engGetVariable((Engine*)mengine,obj->name);
118:   MatSeqAIJFromMatlab(mmat,mat);
119:   return(0);
120: }

124: PetscErrorCode MatSolve_Matlab(Mat A,Vec b,Vec x)
125: {
127:   const char     *_A,*_b,*_x;

130:   /* make sure objects have names; use default if not */
131:   PetscObjectName((PetscObject)b);
132:   PetscObjectName((PetscObject)x);

134:   PetscObjectGetName((PetscObject)A,&_A);
135:   PetscObjectGetName((PetscObject)b,&_b);
136:   PetscObjectGetName((PetscObject)x,&_x);
137:   PetscMatlabEnginePut(PETSC_MATLAB_ENGINE_(PetscObjectComm((PetscObject)A)),(PetscObject)b);
138:   PetscMatlabEngineEvaluate(PETSC_MATLAB_ENGINE_(PetscObjectComm((PetscObject)A)),"%s = u%s\\(l%s\\(p%s*%s));",_x,_A,_A,_A,_b);
139:   PetscMatlabEngineEvaluate(PETSC_MATLAB_ENGINE_(PetscObjectComm((PetscObject)A)),"%s = 0;",_b);
140:   /* PetscMatlabEnginePrintOutput(PETSC_MATLAB_ENGINE_(PetscObjectComm((PetscObject)A)),stdout);  */
141:   PetscMatlabEngineGet(PETSC_MATLAB_ENGINE_(PetscObjectComm((PetscObject)A)),(PetscObject)x);
142:   return(0);
143: }

147: PetscErrorCode MatLUFactorNumeric_Matlab(Mat F,Mat A,const MatFactorInfo *info)
148: {
150:   size_t         len;
151:   char           *_A,*name;
152:   PetscReal      dtcol = info->dtcol;

155:   if (F->factortype == MAT_FACTOR_ILU || info->dt > 0) {
156:     if (info->dtcol == PETSC_DEFAULT) dtcol = .01;
157:     F->ops->solve = MatSolve_Matlab;
158:     F->factortype = MAT_FACTOR_LU;

160:     PetscMatlabEnginePut(PETSC_MATLAB_ENGINE_(PetscObjectComm((PetscObject)A)),(PetscObject)A);
161:     _A   = ((PetscObject)A)->name;
162:     PetscMatlabEngineEvaluate(PETSC_MATLAB_ENGINE_(PetscObjectComm((PetscObject)A)),"info_%s = struct('droptol',%g,'thresh',%g);",_A,info->dt,dtcol);
163:     PetscMatlabEngineEvaluate(PETSC_MATLAB_ENGINE_(PetscObjectComm((PetscObject)A)),"[l_%s,u_%s,p_%s] = luinc(%s',info_%s);",_A,_A,_A,_A,_A);
164:     PetscMatlabEngineEvaluate(PETSC_MATLAB_ENGINE_(PetscObjectComm((PetscObject)A)),"%s = 0;",_A);

166:     PetscStrlen(_A,&len);
167:     PetscMalloc1(len+2,&name);
168:     sprintf(name,"_%s",_A);
169:     PetscObjectSetName((PetscObject)F,name);
170:     PetscFree(name);
171:   } else {
172:     PetscMatlabEnginePut(PETSC_MATLAB_ENGINE_(PetscObjectComm((PetscObject)A)),(PetscObject)A);
173:     _A   = ((PetscObject)A)->name;
174:     PetscMatlabEngineEvaluate(PETSC_MATLAB_ENGINE_(PetscObjectComm((PetscObject)A)),"[l_%s,u_%s,p_%s] = lu(%s',%g);",_A,_A,_A,_A,dtcol);
175:     PetscMatlabEngineEvaluate(PETSC_MATLAB_ENGINE_(PetscObjectComm((PetscObject)A)),"%s = 0;",_A);
176:     PetscStrlen(_A,&len);
177:     PetscMalloc1(len+2,&name);
178:     sprintf(name,"_%s",_A);
179:     PetscObjectSetName((PetscObject)F,name);
180:     PetscFree(name);

182:     F->ops->solve = MatSolve_Matlab;
183:   }
184:   return(0);
185: }

189: PetscErrorCode MatLUFactorSymbolic_Matlab(Mat F,Mat A,IS r,IS c,const MatFactorInfo *info)
190: {
192:   if (A->cmap->N != A->rmap->N) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"matrix must be square");
193:   F->ops->lufactornumeric = MatLUFactorNumeric_Matlab;
194:   F->assembled            = PETSC_TRUE;
195:   return(0);
196: }

200: PetscErrorCode MatFactorGetSolverPackage_seqaij_matlab(Mat A,const MatSolverPackage *type)
201: {
203:   *type = MATSOLVERMATLAB;
204:   return(0);
205: }

209: PETSC_EXTERN PetscErrorCode MatGetFactor_seqaij_matlab(Mat A,MatFactorType ftype,Mat *F)
210: {

214:   if (A->cmap->N != A->rmap->N) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"matrix must be square");
215:   MatCreate(PetscObjectComm((PetscObject)A),F);
216:   MatSetSizes(*F,A->rmap->n,A->cmap->n,A->rmap->n,A->cmap->n);
217:   MatSetType(*F,((PetscObject)A)->type_name);
218:   MatSeqAIJSetPreallocation(*F,0,NULL);
219:   (*F)->ops->lufactorsymbolic  = MatLUFactorSymbolic_Matlab;
220:   (*F)->ops->ilufactorsymbolic = MatLUFactorSymbolic_Matlab;

222:   PetscObjectComposeFunction((PetscObject)(*F),"MatFactorGetSolverPackage_C",MatFactorGetSolverPackage_seqaij_matlab);

224:   (*F)->factortype = ftype;
225:   PetscFree((*F)->solvertype);
226:   PetscStrallocpy(MATSOLVERMATLAB,&(*F)->solvertype);
227:   return(0);
228: }


233: PETSC_EXTERN PetscErrorCode MatSolverPackageRegister_Matlab(void)
234: {

238:   MatSolverPackageRegister(MATSOLVERMATLAB,MATSEQAIJ,        MAT_FACTOR_LU,MatGetFactor_seqaij_matlab);
239:   return(0);
240: }

242: /* --------------------------------------------------------------------------------*/

246: PetscErrorCode MatFactorInfo_Matlab(Mat A,PetscViewer viewer)
247: {

251:   PetscViewerASCIIPrintf(viewer,"MATLAB run parameters:  -- not written yet!\n");
252:   return(0);
253: }

257: PetscErrorCode MatView_Matlab(Mat A,PetscViewer viewer)
258: {
259:   PetscErrorCode    ierr;
260:   PetscBool         iascii;
261:   PetscViewerFormat format;

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


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


281:   Works with MATSEQAIJ matrices.

283:   Options Database Keys:
284: . -pc_factor_mat_solver_package matlab - selects MATLAB to do the sparse factorization


287:   Level: beginner

289: .seealso: PCLU

291: .seealso: PCFactorSetMatSolverPackage(), MatSolverPackage
292: M*/