Actual source code: aijmatlab.c
1: #define PETSCMAT_DLL
3: /*
4: Provides an interface for the Matlab engine sparse solver
6: */
7: #include src/mat/impls/aij/seq/aij.h
9: #include "engine.h" /* Matlab include file */
10: #include "mex.h" /* Matlab include file */
12: typedef struct {
13: PetscErrorCode (*MatDuplicate)(Mat,MatDuplicateOption,Mat*);
14: PetscErrorCode (*MatView)(Mat,PetscViewer);
15: PetscErrorCode (*MatLUFactorSymbolic)(Mat,IS,IS,MatFactorInfo*,Mat*);
16: PetscErrorCode (*MatILUDTFactor)(Mat,IS,IS,MatFactorInfo*,Mat*);
17: PetscErrorCode (*MatDestroy)(Mat);
18: } Mat_Matlab;
24: PetscErrorCode MatMatlabEnginePut_Matlab(PetscObject obj,void *mengine)
25: {
27: Mat B = (Mat)obj;
28: mxArray *mat;
29: Mat_SeqAIJ *aij = (Mat_SeqAIJ*)B->data;
32: mat = mxCreateSparse(B->cmap.n,B->rmap.n,aij->nz,mxREAL);
33: PetscMemcpy(mxGetPr(mat),aij->a,aij->nz*sizeof(PetscScalar));
34: /* Matlab stores by column, not row so we pass in the transpose of the matrix */
35: PetscMemcpy(mxGetIr(mat),aij->j,aij->nz*sizeof(int));
36: PetscMemcpy(mxGetJc(mat),aij->i,(B->rmap.n+1)*sizeof(int));
38: /* Matlab indices start at 0 for sparse (what a surprise) */
39:
40: PetscObjectName(obj);
41: engPutVariable((Engine *)mengine,obj->name,mat);
42: return(0);
43: }
49: PetscErrorCode MatMatlabEngineGet_Matlab(PetscObject obj,void *mengine)
50: {
52: int ii;
53: Mat mat = (Mat)obj;
54: Mat_SeqAIJ *aij = (Mat_SeqAIJ*)mat->data;
55: mxArray *mmat;
58: MatSeqXAIJFreeAIJ(mat,&aij->a,&aij->j,&aij->i);
60: mmat = engGetVariable((Engine *)mengine,obj->name);
62: aij->nz = (mxGetJc(mmat))[mat->rmap.n];
63: PetscMalloc3(aij->nz,PetscScalar,&aij->a,aij->nz,PetscInt,&aij->j,mat->rmap.n+1,PetscInt,&aij->i);
64: aij->singlemalloc = PETSC_TRUE;
66: PetscMemcpy(aij->a,mxGetPr(mmat),aij->nz*sizeof(PetscScalar));
67: /* Matlab stores by column, not row so we pass in the transpose of the matrix */
68: PetscMemcpy(aij->j,mxGetIr(mmat),aij->nz*sizeof(int));
69: PetscMemcpy(aij->i,mxGetJc(mmat),(mat->rmap.n+1)*sizeof(int));
71: for (ii=0; ii<mat->rmap.n; ii++) {
72: aij->ilen[ii] = aij->imax[ii] = aij->i[ii+1] - aij->i[ii];
73: }
75: MatAssemblyBegin(mat,MAT_FINAL_ASSEMBLY);
76: MatAssemblyEnd(mat,MAT_FINAL_ASSEMBLY);
78: return(0);
79: }
85: PetscErrorCode MatConvert_Matlab_SeqAIJ(Mat A,MatType type,MatReuse reuse,Mat *newmat)
86: {
88: Mat B=*newmat;
89: Mat_Matlab *lu=(Mat_Matlab*)A->spptr;
92: if (reuse == MAT_INITIAL_MATRIX) {
93: MatDuplicate(A,MAT_COPY_VALUES,&B);
94: }
95: B->ops->duplicate = lu->MatDuplicate;
96: B->ops->view = lu->MatView;
97: B->ops->lufactorsymbolic = lu->MatLUFactorSymbolic;
98: B->ops->iludtfactor = lu->MatILUDTFactor;
99: B->ops->destroy = lu->MatDestroy;
100:
101: PetscFree(lu);
103: PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqaij_matlab_C","",PETSC_NULL);
104: PetscObjectComposeFunction((PetscObject)B,"MatConvert_matlab_seqaij_C","",PETSC_NULL);
105: PetscObjectComposeFunction((PetscObject)B,"PetscMatlabEnginePut_C","",PETSC_NULL);
106: PetscObjectComposeFunction((PetscObject)B,"PetscMatlabEngineGet_C","",PETSC_NULL);
108: PetscObjectChangeTypeName((PetscObject)B,MATSEQAIJ);
109: *newmat = B;
110: return(0);
111: }
116: PetscErrorCode MatDestroy_Matlab(Mat A)
117: {
121: MatConvert_Matlab_SeqAIJ(A,MATSEQAIJ,MAT_REUSE_MATRIX,&A);
122: (*A->ops->destroy)(A);
123: return(0);
124: }
128: PetscErrorCode MatSolve_Matlab(Mat A,Vec b,Vec x)
129: {
131: const char *_A,*_b,*_x;
134: /* make sure objects have names; use default if not */
135: PetscObjectName((PetscObject)b);
136: PetscObjectName((PetscObject)x);
138: PetscObjectGetName((PetscObject)A,&_A);
139: PetscObjectGetName((PetscObject)b,&_b);
140: PetscObjectGetName((PetscObject)x,&_x);
141: PetscMatlabEnginePut(PETSC_MATLAB_ENGINE_(A->comm),(PetscObject)b);
142: PetscMatlabEngineEvaluate(PETSC_MATLAB_ENGINE_(A->comm),"%s = u%s\\(l%s\\(p%s*%s));",_x,_A,_A,_A,_b);
143: PetscMatlabEngineEvaluate(PETSC_MATLAB_ENGINE_(A->comm),"%s = 0;",_b);
144: /* PetscMatlabEnginePrintOutput(PETSC_MATLAB_ENGINE_(A->comm),stdout); */
145: PetscMatlabEngineGet(PETSC_MATLAB_ENGINE_(A->comm),(PetscObject)x);
146: return(0);
147: }
151: PetscErrorCode MatLUFactorNumeric_Matlab(Mat A,MatFactorInfo *info,Mat *F)
152: {
154: size_t len;
155: char *_A,*name;
158: PetscMatlabEnginePut(PETSC_MATLAB_ENGINE_(A->comm),(PetscObject)A);
159: _A = A->name;
160: PetscMatlabEngineEvaluate(PETSC_MATLAB_ENGINE_(A->comm),"[l_%s,u_%s,p_%s] = lu(%s',%g);",_A,_A,_A,_A,info->dtcol);
161: PetscMatlabEngineEvaluate(PETSC_MATLAB_ENGINE_(A->comm),"%s = 0;",_A);
162: PetscStrlen(_A,&len);
163: PetscMalloc((len+2)*sizeof(char),&name);
164: sprintf(name,"_%s",_A);
165: PetscObjectSetName((PetscObject)*F,name);
166: PetscFree(name);
167: return(0);
168: }
172: PetscErrorCode MatLUFactorSymbolic_Matlab(Mat A,IS r,IS c,MatFactorInfo *info,Mat *F)
173: {
177: if (A->cmap.N != A->rmap.N) SETERRQ(PETSC_ERR_ARG_SIZ,"matrix must be square");
178: MatCreate(A->comm,F);
179: MatSetSizes(*F,A->rmap.n,A->cmap.n,A->rmap.n,A->cmap.n);
180: MatSetType(*F,A->type_name);
181: MatSeqAIJSetPreallocation(*F,0,PETSC_NULL);
182: (*F)->ops->solve = MatSolve_Matlab;
183: (*F)->ops->lufactornumeric = MatLUFactorNumeric_Matlab;
184: (*F)->factor = FACTOR_LU;
185: return(0);
186: }
188: /* ---------------------------------------------------------------------------------*/
191: PetscErrorCode MatSolve_Matlab_QR(Mat A,Vec b,Vec x)
192: {
194: const char *_A,*_b,*_x;
197: /* make sure objects have names; use default if not */
198: PetscObjectName((PetscObject)b);
199: PetscObjectName((PetscObject)x);
201: PetscObjectGetName((PetscObject)A,&_A);
202: PetscObjectGetName((PetscObject)b,&_b);
203: PetscObjectGetName((PetscObject)x,&_x);
204: PetscMatlabEnginePut(PETSC_MATLAB_ENGINE_(A->comm),(PetscObject)b);
205: PetscMatlabEngineEvaluate(PETSC_MATLAB_ENGINE_(A->comm),"%s = r%s\\(r%s'\\(%s*%s));",_x,_A,_A,_A+1,_b);
206: PetscMatlabEngineEvaluate(PETSC_MATLAB_ENGINE_(A->comm),"%s = 0;",_b);
207: /* PetscMatlabEnginePrintOutput(PETSC_MATLAB_ENGINE_(A->comm),stdout); */
208: PetscMatlabEngineGet(PETSC_MATLAB_ENGINE_(A->comm),(PetscObject)x);
209: return(0);
210: }
214: PetscErrorCode MatLUFactorNumeric_Matlab_QR(Mat A,MatFactorInfo *info,Mat *F)
215: {
217: size_t len;
218: char *_A,*name;
221: PetscMatlabEnginePut(PETSC_MATLAB_ENGINE_(A->comm),(PetscObject)A);
222: _A = A->name;
223: PetscMatlabEngineEvaluate(PETSC_MATLAB_ENGINE_(A->comm),"r_%s = qr(%s');",_A,_A);
224: PetscStrlen(_A,&len);
225: PetscMalloc((len+2)*sizeof(char),&name);
226: sprintf(name,"_%s",_A);
227: PetscObjectSetName((PetscObject)*F,name);
228: PetscFree(name);
229: return(0);
230: }
234: PetscErrorCode MatLUFactorSymbolic_Matlab_QR(Mat A,IS r,IS c,MatFactorInfo *info,Mat *F)
235: {
239: if (A->cmap.N != A->rmap.N) SETERRQ(PETSC_ERR_ARG_SIZ,"matrix must be square");
240: MatCreate(A->comm,F);
241: MatSetSizes(*F,A->rmap.n,A->cmap.n,A->rmap.n,A->cmap.n);
242: MatSetType(*F,A->type_name);
243: MatSeqAIJSetPreallocation(*F,0,PETSC_NULL);
244: (*F)->ops->solve = MatSolve_Matlab_QR;
245: (*F)->ops->lufactornumeric = MatLUFactorNumeric_Matlab_QR;
246: (*F)->factor = FACTOR_LU;
247: (*F)->assembled = PETSC_TRUE; /* required by -ksp_view */
249: return(0);
250: }
252: /* --------------------------------------------------------------------------------*/
255: PetscErrorCode MatILUDTFactor_Matlab(Mat A,IS isrow,IS iscol,MatFactorInfo *info,Mat *F)
256: {
258: size_t len;
259: char *_A,*name;
262: if (info->dt == PETSC_DEFAULT) info->dt = .005;
263: if (info->dtcol == PETSC_DEFAULT) info->dtcol = .01;
264: if (A->cmap.N != A->rmap.N) SETERRQ(PETSC_ERR_ARG_SIZ,"matrix must be square");
265: MatCreate(A->comm,F);
266: MatSetSizes(*F,A->rmap.n,A->cmap.n,A->rmap.n,A->cmap.n);
267: MatSetType(*F,A->type_name);
268: MatSeqAIJSetPreallocation(*F,0,PETSC_NULL);
269: (*F)->ops->solve = MatSolve_Matlab;
270: (*F)->factor = FACTOR_LU;
271: PetscMatlabEnginePut(PETSC_MATLAB_ENGINE_(A->comm),(PetscObject)A);
272: _A = A->name;
273: PetscMatlabEngineEvaluate(PETSC_MATLAB_ENGINE_(A->comm),"info_%s = struct('droptol',%g,'thresh',%g);",_A,info->dt,info->dtcol);
274: PetscMatlabEngineEvaluate(PETSC_MATLAB_ENGINE_(A->comm),"[l_%s,u_%s,p_%s] = luinc(%s',info_%s);",_A,_A,_A,_A,_A);
275: PetscMatlabEngineEvaluate(PETSC_MATLAB_ENGINE_(A->comm),"%s = 0;",_A);
277: PetscStrlen(_A,&len);
278: PetscMalloc((len+2)*sizeof(char),&name);
279: sprintf(name,"_%s",_A);
280: PetscObjectSetName((PetscObject)*F,name);
281: PetscFree(name);
282: return(0);
283: }
287: PetscErrorCode MatFactorInfo_Matlab(Mat A,PetscViewer viewer)
288: {
290:
292: PetscViewerASCIIPrintf(viewer,"Matlab run parameters: -- not written yet!\n");
293: return(0);
294: }
298: PetscErrorCode MatView_Matlab(Mat A,PetscViewer viewer)
299: {
300: PetscErrorCode ierr;
301: PetscTruth iascii;
302: PetscViewerFormat format;
303: Mat_Matlab *lu=(Mat_Matlab*)(A->spptr);
306: (*lu->MatView)(A,viewer);
307: PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_ASCII,&iascii);
308: if (iascii) {
309: PetscViewerGetFormat(viewer,&format);
310: if (format == PETSC_VIEWER_ASCII_FACTOR_INFO) {
311: MatFactorInfo_Matlab(A,viewer);
312: }
313: }
314: return(0);
315: }
319: PetscErrorCode MatDuplicate_Matlab(Mat A, MatDuplicateOption op, Mat *M)
320: {
322: Mat_Matlab *lu=(Mat_Matlab*)A->spptr;
325: (*lu->MatDuplicate)(A,op,M);
326: PetscMemcpy((*M)->spptr,lu,sizeof(Mat_Matlab));
327: return(0);
328: }
333: PetscErrorCode MatConvert_SeqAIJ_Matlab(Mat A,MatType type,MatReuse reuse,Mat *newmat)
334: {
336: Mat B=*newmat;
337: Mat_Matlab *lu;
338: PetscTruth qr;
341: if (reuse == MAT_INITIAL_MATRIX) {
342: MatDuplicate(A,MAT_COPY_VALUES,&B);
343: }
345: PetscNew(Mat_Matlab,&lu);
346: lu->MatDuplicate = A->ops->duplicate;
347: lu->MatView = A->ops->view;
348: lu->MatLUFactorSymbolic = A->ops->lufactorsymbolic;
349: lu->MatILUDTFactor = A->ops->iludtfactor;
350: lu->MatDestroy = A->ops->destroy;
352: B->spptr = (void*)lu;
353: B->ops->duplicate = MatDuplicate_Matlab;
354: B->ops->view = MatView_Matlab;
355: B->ops->lufactorsymbolic = MatLUFactorSymbolic_Matlab;
356: B->ops->iludtfactor = MatILUDTFactor_Matlab;
357: B->ops->destroy = MatDestroy_Matlab;
359: PetscOptionsHasName(A->prefix,"-mat_matlab_qr",&qr);
360: if (qr) {
361: B->ops->lufactorsymbolic = MatLUFactorSymbolic_Matlab_QR;
362: PetscInfo(0,"Using Matlab QR with iterative refinement for LU factorization and solves\n");
363: } else {
364: PetscInfo(0,"Using Matlab for LU factorizations and solves.\n");
365: }
366: PetscInfo(0,"Using Matlab for ILUDT factorizations and solves.\n");
368: PetscObjectComposeFunctionDynamic((PetscObject)B,"MatConvert_seqaij_matlab_C",
369: "MatConvert_SeqAIJ_Matlab",MatConvert_SeqAIJ_Matlab);
370: PetscObjectComposeFunctionDynamic((PetscObject)B,"MatConvert_matlab_seqaij_C",
371: "MatConvert_Matlab_SeqAIJ",MatConvert_Matlab_SeqAIJ);
372: PetscObjectComposeFunctionDynamic((PetscObject)B,"PetscMatlabEnginePut_C",
373: "MatMatlabEnginePut_Matlab",MatMatlabEnginePut_Matlab);
374: PetscObjectComposeFunctionDynamic((PetscObject)B,"PetscMatlabEngineGet_C",
375: "MatMatlabEngineGet_Matlab",MatMatlabEngineGet_Matlab);
376: PetscObjectChangeTypeName((PetscObject)B,MATMATLAB);
377: *newmat = B;
378: return(0);
379: }
382: /*MC
383: MATMATLAB - MATMATLAB = "matlab" - A matrix type providing direct solvers (LU and QR) and drop tolerance
384: based ILU factorization (ILUDT) for sequential matrices via the external package Matlab.
386: If Matlab is instaled (see the manual for
387: instructions on how to declare the existence of external packages),
388: a matrix type can be constructed which invokes Matlab solvers.
389: After calling MatCreate(...,A), simply call MatSetType(A,MATMATLAB).
390: This matrix type is only supported for double precision real.
392: This matrix inherits from MATSEQAIJ. As a result, MatSeqAIJSetPreallocation is
393: supported for this matrix type. One can also call MatConvert for an inplace conversion to or from
394: the MATSEQAIJ type without data copy.
396: Options Database Keys:
397: + -mat_type matlab - sets the matrix type to "matlab" during a call to MatSetFromOptions()
398: - -mat_matlab_qr - sets the direct solver to be QR instead of LU
400: Level: beginner
402: .seealso: PCLU
403: M*/
408: PetscErrorCode MatCreate_Matlab(Mat A)
409: {
413: MatSetType(A,MATSEQAIJ);
414: MatConvert_SeqAIJ_Matlab(A,MATMATLAB,MAT_REUSE_MATRIX,&A);
415: return(0);
416: }