Actual source code: aijmatlab.c
petsc-3.9.4 2018-09-11
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: PetscMemcpy(mxGetPr(mat),aij->a,aij->nz*sizeof(PetscScalar));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: /*@C
43: MatSeqAIJFromMatlab - Given a MATLAB sparse matrix, fills a SeqAIJ matrix with its transpose.
45: Not Collective
47: Input Parameters:
48: + mmat - a MATLAB sparse matris
49: - mat - an already created MATSEQAIJ
51: Level: intermediate
53: @*/
54: PETSC_EXTERN PetscErrorCode MatSeqAIJFromMatlab(mxArray *mmat,Mat mat)
55: {
57: PetscInt nz,n,m,*i,*j,k;
58: mwIndex nnz,nn,nm,*ii,*jj;
59: Mat_SeqAIJ *aij = (Mat_SeqAIJ*)mat->data;
62: nn = mxGetN(mmat); /* rows of transpose of matrix */
63: nm = mxGetM(mmat);
64: nnz = (mxGetJc(mmat))[nn];
65: ii = mxGetJc(mmat);
66: jj = mxGetIr(mmat);
67: n = (PetscInt) nn;
68: m = (PetscInt) nm;
69: nz = (PetscInt) nnz;
71: if (mat->rmap->n < 0 && mat->cmap->n < 0) {
72: /* matrix has not yet had its size set */
73: MatSetSizes(mat,n,m,PETSC_DETERMINE,PETSC_DETERMINE);
74: MatSetUp(mat);
75: } else {
76: 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);
77: 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);
78: }
79: if (nz != aij->nz) {
80: /* number of nonzeros in matrix has changed, so need new data structure */
81: MatSeqXAIJFreeAIJ(mat,&aij->a,&aij->j,&aij->i);
82: aij->nz = nz;
83: PetscMalloc3(aij->nz,&aij->a,aij->nz,&aij->j,mat->rmap->n+1,&aij->i);
85: aij->singlemalloc = PETSC_TRUE;
86: }
88: PetscMemcpy(aij->a,mxGetPr(mmat),aij->nz*sizeof(PetscScalar));
89: /* MATLAB stores by column, not row so we pass in the transpose of the matrix */
90: i = aij->i;
91: for (k=0; k<n+1; k++) i[k] = (PetscInt) ii[k];
92: j = aij->j;
93: for (k=0; k<nz; k++) j[k] = (PetscInt) jj[k];
95: for (k=0; k<mat->rmap->n; k++) aij->ilen[k] = aij->imax[k] = aij->i[k+1] - aij->i[k];
97: MatAssemblyBegin(mat,MAT_FINAL_ASSEMBLY);
98: MatAssemblyEnd(mat,MAT_FINAL_ASSEMBLY);
99: return(0);
100: }
102: PETSC_EXTERN PetscErrorCode MatlabEngineGet_SeqAIJ(PetscObject obj,void *mengine)
103: {
105: Mat mat = (Mat)obj;
106: mxArray *mmat;
109: mmat = engGetVariable((Engine*)mengine,obj->name);
110: MatSeqAIJFromMatlab(mmat,mat);
111: return(0);
112: }
114: PetscErrorCode MatSolve_Matlab(Mat A,Vec b,Vec x)
115: {
117: const char *_A,*_b,*_x;
120: /* make sure objects have names; use default if not */
121: PetscObjectName((PetscObject)b);
122: PetscObjectName((PetscObject)x);
124: PetscObjectGetName((PetscObject)A,&_A);
125: PetscObjectGetName((PetscObject)b,&_b);
126: PetscObjectGetName((PetscObject)x,&_x);
127: PetscMatlabEnginePut(PETSC_MATLAB_ENGINE_(PetscObjectComm((PetscObject)A)),(PetscObject)b);
128: PetscMatlabEngineEvaluate(PETSC_MATLAB_ENGINE_(PetscObjectComm((PetscObject)A)),"%s = u%s\\(l%s\\(p%s*%s));",_x,_A,_A,_A,_b);
129: PetscMatlabEngineEvaluate(PETSC_MATLAB_ENGINE_(PetscObjectComm((PetscObject)A)),"%s = 0;",_b);
130: /* PetscMatlabEnginePrintOutput(PETSC_MATLAB_ENGINE_(PetscObjectComm((PetscObject)A)),stdout); */
131: PetscMatlabEngineGet(PETSC_MATLAB_ENGINE_(PetscObjectComm((PetscObject)A)),(PetscObject)x);
132: return(0);
133: }
135: PetscErrorCode MatLUFactorNumeric_Matlab(Mat F,Mat A,const MatFactorInfo *info)
136: {
138: size_t len;
139: char *_A,*name;
140: PetscReal dtcol = info->dtcol;
143: if (F->factortype == MAT_FACTOR_ILU || info->dt > 0) {
144: /* the ILU form is not currently registered */
145: if (info->dtcol == PETSC_DEFAULT) dtcol = .01;
146: F->ops->solve = MatSolve_Matlab;
147: F->factortype = MAT_FACTOR_LU;
149: PetscMatlabEnginePut(PETSC_MATLAB_ENGINE_(PetscObjectComm((PetscObject)A)),(PetscObject)A);
150: _A = ((PetscObject)A)->name;
151: PetscMatlabEngineEvaluate(PETSC_MATLAB_ENGINE_(PetscObjectComm((PetscObject)A)),"info_%s = struct('droptol',%g,'thresh',%g);",_A,info->dt,dtcol);
152: PetscMatlabEngineEvaluate(PETSC_MATLAB_ENGINE_(PetscObjectComm((PetscObject)A)),"[l_%s,u_%s,p_%s] = luinc(%s',info_%s);",_A,_A,_A,_A,_A);
153: PetscMatlabEngineEvaluate(PETSC_MATLAB_ENGINE_(PetscObjectComm((PetscObject)A)),"%s = 0;",_A);
155: PetscStrlen(_A,&len);
156: PetscMalloc1(len+2,&name);
157: sprintf(name,"_%s",_A);
158: PetscObjectSetName((PetscObject)F,name);
159: PetscFree(name);
160: } else {
161: PetscMatlabEnginePut(PETSC_MATLAB_ENGINE_(PetscObjectComm((PetscObject)A)),(PetscObject)A);
162: _A = ((PetscObject)A)->name;
163: PetscMatlabEngineEvaluate(PETSC_MATLAB_ENGINE_(PetscObjectComm((PetscObject)A)),"[l_%s,u_%s,p_%s] = lu(%s',%g);",_A,_A,_A,_A,dtcol);
164: PetscMatlabEngineEvaluate(PETSC_MATLAB_ENGINE_(PetscObjectComm((PetscObject)A)),"%s = 0;",_A);
165: PetscStrlen(_A,&len);
166: PetscMalloc1(len+2,&name);
167: sprintf(name,"_%s",_A);
168: PetscObjectSetName((PetscObject)F,name);
169: PetscFree(name);
171: F->ops->solve = MatSolve_Matlab;
172: }
173: return(0);
174: }
176: PetscErrorCode MatLUFactorSymbolic_Matlab(Mat F,Mat A,IS r,IS c,const MatFactorInfo *info)
177: {
179: if (A->cmap->N != A->rmap->N) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"matrix must be square");
180: F->ops->lufactornumeric = MatLUFactorNumeric_Matlab;
181: F->assembled = PETSC_TRUE;
182: return(0);
183: }
185: PetscErrorCode MatFactorGetSolverType_seqaij_matlab(Mat A,MatSolverType *type)
186: {
188: *type = MATSOLVERMATLAB;
189: return(0);
190: }
192: PetscErrorCode MatDestroy_matlab(Mat A)
193: {
195: const char *_A;
198: PetscObjectGetName((PetscObject)A,&_A);
199: PetscMatlabEngineEvaluate(PETSC_MATLAB_ENGINE_(PetscObjectComm((PetscObject)A)),"delete %s l_%s u_%s;",_A,_A,_A);
200: return(0);
201: }
203: PETSC_EXTERN PetscErrorCode MatGetFactor_seqaij_matlab(Mat A,MatFactorType ftype,Mat *F)
204: {
208: if (A->cmap->N != A->rmap->N) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"matrix must be square");
209: MatCreate(PetscObjectComm((PetscObject)A),F);
210: MatSetSizes(*F,A->rmap->n,A->cmap->n,A->rmap->n,A->cmap->n);
211: PetscStrallocpy("matlab",&((PetscObject)*F)->type_name);
212: MatSetUp(*F);
214: (*F)->ops->destroy = MatDestroy_matlab;
215: (*F)->ops->getinfo = MatGetInfo_External;
216: (*F)->ops->lufactorsymbolic = MatLUFactorSymbolic_Matlab;
217: (*F)->ops->ilufactorsymbolic = MatLUFactorSymbolic_Matlab;
219: PetscObjectComposeFunction((PetscObject)(*F),"MatFactorGetSolverType_C",MatFactorGetSolverType_seqaij_matlab);
221: (*F)->factortype = ftype;
222: PetscFree((*F)->solvertype);
223: PetscStrallocpy(MATSOLVERMATLAB,&(*F)->solvertype);
224: return(0);
225: }
228: PETSC_EXTERN PetscErrorCode MatSolverTypeRegister_Matlab(void)
229: {
233: MatSolverTypeRegister(MATSOLVERMATLAB,MATSEQAIJ, MAT_FACTOR_LU,MatGetFactor_seqaij_matlab);
234: return(0);
235: }
237: /* --------------------------------------------------------------------------------*/
239: PetscErrorCode MatView_Info_Matlab(Mat A,PetscViewer viewer)
240: {
244: PetscViewerASCIIPrintf(viewer,"MATLAB run parameters: -- not written yet!\n");
245: return(0);
246: }
248: PetscErrorCode MatView_Matlab(Mat A,PetscViewer viewer)
249: {
250: PetscErrorCode ierr;
251: PetscBool iascii;
252: PetscViewerFormat format;
255: MatView_SeqAIJ(A,viewer);
256: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
257: if (iascii) {
258: PetscViewerGetFormat(viewer,&format);
259: if (format == PETSC_VIEWER_ASCII_FACTOR_INFO) {
260: MatView_Info_Matlab(A,viewer);
261: }
262: }
263: return(0);
264: }
267: /*MC
268: MATSOLVERMATLAB - "matlab" - Providing direct solver LU for sequential aij matrix via the external package MATLAB.
271: Works with MATSEQAIJ matrices.
273: Options Database Keys:
274: . -pc_factor_mat_solver_type matlab - selects MATLAB to do the sparse factorization
276: Notes: You must ./configure with the options --with-matlab --with-matlab-engine
278: Level: beginner
280: .seealso: PCLU
282: .seealso: PCFactorSetMatSolverType(), MatSolverType
283: M*/