Actual source code: aijmatlab.c
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: {
13: Mat_SeqAIJ *aij = (Mat_SeqAIJ*)B->data;
14: mwIndex *ii,*jj;
15: mxArray *mat;
16: PetscInt i;
18: mat = mxCreateSparse(B->cmap->n,B->rmap->n,aij->nz,mxREAL);
19: if (PetscArraycpy(mxGetPr(mat),aij->a,aij->nz)) return NULL;
20: /* MATLAB stores by column, not row so we pass in the transpose of the matrix */
21: jj = mxGetIr(mat);
22: for (i=0; i<aij->nz; i++) jj[i] = aij->j[i];
23: ii = mxGetJc(mat);
24: for (i=0; i<B->rmap->n+1; i++) ii[i] = aij->i[i];
25: return mat;
26: }
28: PETSC_EXTERN PetscErrorCode MatlabEnginePut_SeqAIJ(PetscObject obj,void *mengine)
29: {
30: mxArray *mat;
33: PetscObjectName(obj);
34: engPutVariable((Engine*)mengine,obj->name,mat);
35: return 0;
36: }
38: PETSC_EXTERN PetscErrorCode MatSeqAIJFromMatlab(mxArray *mmat,Mat mat)
39: {
40: PetscInt nz,n,m,*i,*j,k;
41: mwIndex nnz,nn,nm,*ii,*jj;
42: Mat_SeqAIJ *aij = (Mat_SeqAIJ*)mat->data;
44: nn = mxGetN(mmat); /* rows of transpose of matrix */
45: nm = mxGetM(mmat);
46: nnz = (mxGetJc(mmat))[nn];
47: ii = mxGetJc(mmat);
48: jj = mxGetIr(mmat);
49: n = (PetscInt) nn;
50: m = (PetscInt) nm;
51: nz = (PetscInt) nnz;
53: if (mat->rmap->n < 0 && mat->cmap->n < 0) {
54: /* matrix has not yet had its size set */
55: MatSetSizes(mat,n,m,PETSC_DETERMINE,PETSC_DETERMINE);
56: MatSetUp(mat);
57: } else {
60: }
61: if (nz != aij->nz) {
62: /* number of nonzeros in matrix has changed, so need new data structure */
63: MatSeqXAIJFreeAIJ(mat,&aij->a,&aij->j,&aij->i);
64: aij->nz = nz;
65: PetscMalloc3(aij->nz,&aij->a,aij->nz,&aij->j,mat->rmap->n+1,&aij->i);
67: aij->singlemalloc = PETSC_TRUE;
68: }
70: PetscArraycpy(aij->a,mxGetPr(mmat),aij->nz);
71: /* MATLAB stores by column, not row so we pass in the transpose of the matrix */
72: i = aij->i;
73: for (k=0; k<n+1; k++) i[k] = (PetscInt) ii[k];
74: j = aij->j;
75: for (k=0; k<nz; k++) j[k] = (PetscInt) jj[k];
77: for (k=0; k<mat->rmap->n; k++) aij->ilen[k] = aij->imax[k] = aij->i[k+1] - aij->i[k];
79: mat->nonzerostate++; /* since the nonzero structure can change anytime force the Inode information to always be rebuilt */
80: MatAssemblyBegin(mat,MAT_FINAL_ASSEMBLY);
81: MatAssemblyEnd(mat,MAT_FINAL_ASSEMBLY);
82: return 0;
83: }
85: PETSC_EXTERN PetscErrorCode MatlabEngineGet_SeqAIJ(PetscObject obj,void *mengine)
86: {
87: Mat mat = (Mat)obj;
88: mxArray *mmat;
90: mmat = engGetVariable((Engine*)mengine,obj->name);
91: MatSeqAIJFromMatlab(mmat,mat);
92: return 0;
93: }
95: PetscErrorCode MatSolve_Matlab(Mat A,Vec b,Vec x)
96: {
97: const char *_A,*_b,*_x;
99: /* make sure objects have names; use default if not */
100: PetscObjectName((PetscObject)b);
101: PetscObjectName((PetscObject)x);
103: PetscObjectGetName((PetscObject)A,&_A);
104: PetscObjectGetName((PetscObject)b,&_b);
105: PetscObjectGetName((PetscObject)x,&_x);
106: PetscMatlabEnginePut(PETSC_MATLAB_ENGINE_(PetscObjectComm((PetscObject)A)),(PetscObject)b);
107: PetscMatlabEngineEvaluate(PETSC_MATLAB_ENGINE_(PetscObjectComm((PetscObject)A)),"%s = u%s\\(l%s\\(p%s*%s));",_x,_A,_A,_A,_b);
108: PetscMatlabEngineEvaluate(PETSC_MATLAB_ENGINE_(PetscObjectComm((PetscObject)A)),"%s = 0;",_b);
109: /* PetscMatlabEnginePrintOutput(PETSC_MATLAB_ENGINE_(PetscObjectComm((PetscObject)A)),stdout); */
110: PetscMatlabEngineGet(PETSC_MATLAB_ENGINE_(PetscObjectComm((PetscObject)A)),(PetscObject)x);
111: return 0;
112: }
114: PetscErrorCode MatLUFactorNumeric_Matlab(Mat F,Mat A,const MatFactorInfo *info)
115: {
116: size_t len;
117: char *_A,*name;
118: PetscReal dtcol = info->dtcol;
120: if (F->factortype == MAT_FACTOR_ILU || info->dt > 0) {
121: /* the ILU form is not currently registered */
122: if (info->dtcol == PETSC_DEFAULT) dtcol = .01;
123: F->ops->solve = MatSolve_Matlab;
124: F->factortype = MAT_FACTOR_LU;
126: PetscMatlabEnginePut(PETSC_MATLAB_ENGINE_(PetscObjectComm((PetscObject)A)),(PetscObject)A);
127: _A = ((PetscObject)A)->name;
128: PetscMatlabEngineEvaluate(PETSC_MATLAB_ENGINE_(PetscObjectComm((PetscObject)A)),"info_%s = struct('droptol',%g,'thresh',%g);",_A,info->dt,dtcol);
129: PetscMatlabEngineEvaluate(PETSC_MATLAB_ENGINE_(PetscObjectComm((PetscObject)A)),"[l_%s,u_%s,p_%s] = luinc(%s',info_%s);",_A,_A,_A,_A,_A);
130: PetscMatlabEngineEvaluate(PETSC_MATLAB_ENGINE_(PetscObjectComm((PetscObject)A)),"%s = 0;",_A);
132: PetscStrlen(_A,&len);
133: PetscMalloc1(len+2,&name);
134: sprintf(name,"_%s",_A);
135: PetscObjectSetName((PetscObject)F,name);
136: PetscFree(name);
137: } else {
138: PetscMatlabEnginePut(PETSC_MATLAB_ENGINE_(PetscObjectComm((PetscObject)A)),(PetscObject)A);
139: _A = ((PetscObject)A)->name;
140: PetscMatlabEngineEvaluate(PETSC_MATLAB_ENGINE_(PetscObjectComm((PetscObject)A)),"[l_%s,u_%s,p_%s] = lu(%s',%g);",_A,_A,_A,_A,dtcol);
141: PetscMatlabEngineEvaluate(PETSC_MATLAB_ENGINE_(PetscObjectComm((PetscObject)A)),"%s = 0;",_A);
142: PetscStrlen(_A,&len);
143: PetscMalloc1(len+2,&name);
144: sprintf(name,"_%s",_A);
145: PetscObjectSetName((PetscObject)F,name);
146: PetscFree(name);
148: F->ops->solve = MatSolve_Matlab;
149: }
150: return 0;
151: }
153: PetscErrorCode MatLUFactorSymbolic_Matlab(Mat F,Mat A,IS r,IS c,const MatFactorInfo *info)
154: {
156: F->ops->lufactornumeric = MatLUFactorNumeric_Matlab;
157: F->assembled = PETSC_TRUE;
158: return 0;
159: }
161: PetscErrorCode MatFactorGetSolverType_seqaij_matlab(Mat A,MatSolverType *type)
162: {
163: *type = MATSOLVERMATLAB;
164: return 0;
165: }
167: PetscErrorCode MatDestroy_matlab(Mat A)
168: {
169: const char *_A;
171: PetscObjectGetName((PetscObject)A,&_A);
172: PetscMatlabEngineEvaluate(PETSC_MATLAB_ENGINE_(PetscObjectComm((PetscObject)A)),"delete %s l_%s u_%s;",_A,_A,_A);
173: return 0;
174: }
176: PETSC_EXTERN PetscErrorCode MatGetFactor_seqaij_matlab(Mat A,MatFactorType ftype,Mat *F)
177: {
179: MatCreate(PetscObjectComm((PetscObject)A),F);
180: MatSetSizes(*F,A->rmap->n,A->cmap->n,A->rmap->n,A->cmap->n);
181: PetscStrallocpy("matlab",&((PetscObject)*F)->type_name);
182: MatSetUp(*F);
184: (*F)->ops->destroy = MatDestroy_matlab;
185: (*F)->ops->getinfo = MatGetInfo_External;
186: (*F)->trivialsymbolic = PETSC_TRUE;
187: (*F)->ops->lufactorsymbolic = MatLUFactorSymbolic_Matlab;
188: (*F)->ops->ilufactorsymbolic = MatLUFactorSymbolic_Matlab;
190: PetscObjectComposeFunction((PetscObject)(*F),"MatFactorGetSolverType_C",MatFactorGetSolverType_seqaij_matlab);
192: (*F)->factortype = ftype;
193: PetscFree((*F)->solvertype);
194: PetscStrallocpy(MATSOLVERMATLAB,&(*F)->solvertype);
195: return 0;
196: }
198: PETSC_EXTERN PetscErrorCode MatSolverTypeRegister_Matlab(void)
199: {
200: MatSolverTypeRegister(MATSOLVERMATLAB,MATSEQAIJ, MAT_FACTOR_LU,MatGetFactor_seqaij_matlab);
201: return 0;
202: }
204: /* --------------------------------------------------------------------------------*/
206: PetscErrorCode MatView_Info_Matlab(Mat A,PetscViewer viewer)
207: {
208: PetscViewerASCIIPrintf(viewer,"MATLAB run parameters: -- not written yet!\n");
209: return 0;
210: }
212: PetscErrorCode MatView_Matlab(Mat A,PetscViewer viewer)
213: {
214: PetscBool iascii;
216: MatView_SeqAIJ(A,viewer);
217: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
218: if (iascii) {
219: PetscViewerFormat format;
221: PetscViewerGetFormat(viewer,&format);
222: if (format == PETSC_VIEWER_ASCII_FACTOR_INFO) MatView_Info_Matlab(A,viewer);
223: }
224: return 0;
225: }
227: /*MC
228: MATSOLVERMATLAB - "matlab" - Providing direct solver LU for sequential aij matrix via the external package MATLAB.
230: Works with MATSEQAIJ matrices.
232: Options Database Keys:
233: . -pc_factor_mat_solver_type matlab - selects MATLAB to do the sparse factorization
235: Notes:
236: You must ./configure with the options --with-matlab --with-matlab-engine
238: Level: beginner
240: .seealso: PCLU
242: .seealso: PCFactorSetMatSolverType(), MatSolverType
243: M*/