Actual source code: aijmatlab.c
petsc-3.7.7 2017-09-25
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*/