Actual source code: aijmatlab.c

petsc-3.13.6 2020-09-29
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 */

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

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

 30: PETSC_EXTERN PetscErrorCode MatlabEnginePut_SeqAIJ(PetscObject obj,void *mengine)
 31: {
 33:   mxArray        *mat;

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

 42: PETSC_EXTERN PetscErrorCode MatSeqAIJFromMatlab(mxArray *mmat,Mat mat)
 43: {
 45:   PetscInt       nz,n,m,*i,*j,k;
 46:   mwIndex        nnz,nn,nm,*ii,*jj;
 47:   Mat_SeqAIJ     *aij = (Mat_SeqAIJ*)mat->data;

 50:   nn  = mxGetN(mmat);   /* rows of transpose of matrix */
 51:   nm  = mxGetM(mmat);
 52:   nnz = (mxGetJc(mmat))[nn];
 53:   ii  = mxGetJc(mmat);
 54:   jj  = mxGetIr(mmat);
 55:   n   = (PetscInt) nn;
 56:   m   = (PetscInt) nm;
 57:   nz  = (PetscInt) nnz;

 59:   if (mat->rmap->n < 0 && mat->cmap->n < 0) {
 60:     /* matrix has not yet had its size set */
 61:     MatSetSizes(mat,n,m,PETSC_DETERMINE,PETSC_DETERMINE);
 62:     MatSetUp(mat);
 63:   } else {
 64:     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);
 65:     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);
 66:   }
 67:   if (nz != aij->nz) {
 68:     /* number of nonzeros in matrix has changed, so need new data structure */
 69:     MatSeqXAIJFreeAIJ(mat,&aij->a,&aij->j,&aij->i);
 70:     aij->nz = nz;
 71:     PetscMalloc3(aij->nz,&aij->a,aij->nz,&aij->j,mat->rmap->n+1,&aij->i);

 73:     aij->singlemalloc = PETSC_TRUE;
 74:   }

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

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

 85:   mat->nonzerostate++; /* since the nonzero structure can change anytime force the Inode information to always be rebuilt */
 86:   MatAssemblyBegin(mat,MAT_FINAL_ASSEMBLY);
 87:   MatAssemblyEnd(mat,MAT_FINAL_ASSEMBLY);
 88:   return(0);
 89: }

 91: PETSC_EXTERN PetscErrorCode  MatlabEngineGet_SeqAIJ(PetscObject obj,void *mengine)
 92: {
 94:   Mat            mat = (Mat)obj;
 95:   mxArray        *mmat;

 98:   mmat = engGetVariable((Engine*)mengine,obj->name);
 99:   MatSeqAIJFromMatlab(mmat,mat);
100:   return(0);
101: }

103: PetscErrorCode MatSolve_Matlab(Mat A,Vec b,Vec x)
104: {
106:   const char     *_A,*_b,*_x;

109:   /* make sure objects have names; use default if not */
110:   PetscObjectName((PetscObject)b);
111:   PetscObjectName((PetscObject)x);

113:   PetscObjectGetName((PetscObject)A,&_A);
114:   PetscObjectGetName((PetscObject)b,&_b);
115:   PetscObjectGetName((PetscObject)x,&_x);
116:   PetscMatlabEnginePut(PETSC_MATLAB_ENGINE_(PetscObjectComm((PetscObject)A)),(PetscObject)b);
117:   PetscMatlabEngineEvaluate(PETSC_MATLAB_ENGINE_(PetscObjectComm((PetscObject)A)),"%s = u%s\\(l%s\\(p%s*%s));",_x,_A,_A,_A,_b);
118:   PetscMatlabEngineEvaluate(PETSC_MATLAB_ENGINE_(PetscObjectComm((PetscObject)A)),"%s = 0;",_b);
119:   /* PetscMatlabEnginePrintOutput(PETSC_MATLAB_ENGINE_(PetscObjectComm((PetscObject)A)),stdout);  */
120:   PetscMatlabEngineGet(PETSC_MATLAB_ENGINE_(PetscObjectComm((PetscObject)A)),(PetscObject)x);
121:   return(0);
122: }

124: PetscErrorCode MatLUFactorNumeric_Matlab(Mat F,Mat A,const MatFactorInfo *info)
125: {
127:   size_t         len;
128:   char           *_A,*name;
129:   PetscReal      dtcol = info->dtcol;

132:   if (F->factortype == MAT_FACTOR_ILU || info->dt > 0) {
133:     /* the ILU form is not currently registered */
134:     if (info->dtcol == PETSC_DEFAULT) dtcol = .01;
135:     F->ops->solve = MatSolve_Matlab;
136:     F->factortype = MAT_FACTOR_LU;

138:     PetscMatlabEnginePut(PETSC_MATLAB_ENGINE_(PetscObjectComm((PetscObject)A)),(PetscObject)A);
139:     _A   = ((PetscObject)A)->name;
140:     PetscMatlabEngineEvaluate(PETSC_MATLAB_ENGINE_(PetscObjectComm((PetscObject)A)),"info_%s = struct('droptol',%g,'thresh',%g);",_A,info->dt,dtcol);
141:     PetscMatlabEngineEvaluate(PETSC_MATLAB_ENGINE_(PetscObjectComm((PetscObject)A)),"[l_%s,u_%s,p_%s] = luinc(%s',info_%s);",_A,_A,_A,_A,_A);
142:     PetscMatlabEngineEvaluate(PETSC_MATLAB_ENGINE_(PetscObjectComm((PetscObject)A)),"%s = 0;",_A);

144:     PetscStrlen(_A,&len);
145:     PetscMalloc1(len+2,&name);
146:     sprintf(name,"_%s",_A);
147:     PetscObjectSetName((PetscObject)F,name);
148:     PetscFree(name);
149:   } else {
150:     PetscMatlabEnginePut(PETSC_MATLAB_ENGINE_(PetscObjectComm((PetscObject)A)),(PetscObject)A);
151:     _A   = ((PetscObject)A)->name;
152:     PetscMatlabEngineEvaluate(PETSC_MATLAB_ENGINE_(PetscObjectComm((PetscObject)A)),"[l_%s,u_%s,p_%s] = lu(%s',%g);",_A,_A,_A,_A,dtcol);
153:     PetscMatlabEngineEvaluate(PETSC_MATLAB_ENGINE_(PetscObjectComm((PetscObject)A)),"%s = 0;",_A);
154:     PetscStrlen(_A,&len);
155:     PetscMalloc1(len+2,&name);
156:     sprintf(name,"_%s",_A);
157:     PetscObjectSetName((PetscObject)F,name);
158:     PetscFree(name);

160:     F->ops->solve = MatSolve_Matlab;
161:   }
162:   return(0);
163: }

165: PetscErrorCode MatLUFactorSymbolic_Matlab(Mat F,Mat A,IS r,IS c,const MatFactorInfo *info)
166: {
168:   if (A->cmap->N != A->rmap->N) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"matrix must be square");
169:   F->ops->lufactornumeric = MatLUFactorNumeric_Matlab;
170:   F->assembled            = PETSC_TRUE;
171:   return(0);
172: }

174: PetscErrorCode MatFactorGetSolverType_seqaij_matlab(Mat A,MatSolverType *type)
175: {
177:   *type = MATSOLVERMATLAB;
178:   return(0);
179: }

181: PetscErrorCode MatDestroy_matlab(Mat A)
182: {
184:   const char     *_A;

187:   PetscObjectGetName((PetscObject)A,&_A);
188:   PetscMatlabEngineEvaluate(PETSC_MATLAB_ENGINE_(PetscObjectComm((PetscObject)A)),"delete %s l_%s u_%s;",_A,_A,_A);
189:   return(0);
190: }

192: PETSC_EXTERN PetscErrorCode MatGetFactor_seqaij_matlab(Mat A,MatFactorType ftype,Mat *F)
193: {

197:   if (A->cmap->N != A->rmap->N) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"matrix must be square");
198:   MatCreate(PetscObjectComm((PetscObject)A),F);
199:   MatSetSizes(*F,A->rmap->n,A->cmap->n,A->rmap->n,A->cmap->n);
200:   PetscStrallocpy("matlab",&((PetscObject)*F)->type_name);
201:   MatSetUp(*F);

203:   (*F)->ops->destroy           = MatDestroy_matlab;
204:   (*F)->ops->getinfo           = MatGetInfo_External;  
205:   (*F)->ops->lufactorsymbolic  = MatLUFactorSymbolic_Matlab;
206:   (*F)->ops->ilufactorsymbolic = MatLUFactorSymbolic_Matlab;

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

210:   (*F)->factortype = ftype;
211:   PetscFree((*F)->solvertype);
212:   PetscStrallocpy(MATSOLVERMATLAB,&(*F)->solvertype);
213:   return(0);
214: }


217: PETSC_EXTERN PetscErrorCode MatSolverTypeRegister_Matlab(void)
218: {

222:   MatSolverTypeRegister(MATSOLVERMATLAB,MATSEQAIJ,        MAT_FACTOR_LU,MatGetFactor_seqaij_matlab);
223:   return(0);
224: }

226: /* --------------------------------------------------------------------------------*/

228: PetscErrorCode MatView_Info_Matlab(Mat A,PetscViewer viewer)
229: {

233:   PetscViewerASCIIPrintf(viewer,"MATLAB run parameters:  -- not written yet!\n");
234:   return(0);
235: }

237: PetscErrorCode MatView_Matlab(Mat A,PetscViewer viewer)
238: {
239:   PetscErrorCode    ierr;
240:   PetscBool         iascii;
241:   PetscViewerFormat format;

244:   MatView_SeqAIJ(A,viewer);
245:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
246:   if (iascii) {
247:     PetscViewerGetFormat(viewer,&format);
248:     if (format == PETSC_VIEWER_ASCII_FACTOR_INFO) {
249:       MatView_Info_Matlab(A,viewer);
250:     }
251:   }
252:   return(0);
253: }


256: /*MC
257:   MATSOLVERMATLAB - "matlab" - Providing direct solver LU for sequential aij matrix via the external package MATLAB.


260:   Works with MATSEQAIJ matrices.

262:   Options Database Keys:
263: . -pc_factor_mat_solver_type matlab - selects MATLAB to do the sparse factorization

265:   Notes:
266:     You must ./configure with the options --with-matlab --with-matlab-engine

268:   Level: beginner

270: .seealso: PCLU

272: .seealso: PCFactorSetMatSolverType(), MatSolverType
273: M*/