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*/