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;
32: mat = MatSeqAIJToMatlab((Mat)obj);
34: PetscObjectName(obj);
35: engPutVariable((Engine *)mengine, obj->name, mat);
36: return 0;
37: }
39: PETSC_EXTERN PetscErrorCode MatSeqAIJFromMatlab(mxArray *mmat, Mat mat)
40: {
41: PetscInt nz, n, m, *i, *j, k;
42: mwIndex nnz, nn, nm, *ii, *jj;
43: Mat_SeqAIJ *aij = (Mat_SeqAIJ *)mat->data;
45: nn = mxGetN(mmat); /* rows of transpose of matrix */
46: nm = mxGetM(mmat);
47: nnz = (mxGetJc(mmat))[nn];
48: ii = mxGetJc(mmat);
49: jj = mxGetIr(mmat);
50: n = (PetscInt)nn;
51: m = (PetscInt)nm;
52: nz = (PetscInt)nnz;
54: if (mat->rmap->n < 0 && mat->cmap->n < 0) {
55: /* matrix has not yet had its size set */
56: MatSetSizes(mat, n, m, PETSC_DETERMINE, PETSC_DETERMINE);
57: MatSetUp(mat);
58: } else {
61: }
62: if (nz != aij->nz) {
63: /* number of nonzeros in matrix has changed, so need new data structure */
64: MatSeqXAIJFreeAIJ(mat, &aij->a, &aij->j, &aij->i);
65: aij->nz = nz;
66: PetscMalloc3(aij->nz, &aij->a, aij->nz, &aij->j, mat->rmap->n + 1, &aij->i);
68: aij->singlemalloc = PETSC_TRUE;
69: }
71: PetscArraycpy(aij->a, mxGetPr(mmat), aij->nz);
72: /* MATLAB stores by column, not row so we pass in the transpose of the matrix */
73: i = aij->i;
74: for (k = 0; k < n + 1; k++) i[k] = (PetscInt)ii[k];
75: j = aij->j;
76: for (k = 0; k < nz; k++) j[k] = (PetscInt)jj[k];
78: for (k = 0; k < mat->rmap->n; k++) aij->ilen[k] = aij->imax[k] = aij->i[k + 1] - aij->i[k];
80: mat->nonzerostate++; /* since the nonzero structure can change anytime force the Inode information to always be rebuilt */
81: MatAssemblyBegin(mat, MAT_FINAL_ASSEMBLY);
82: MatAssemblyEnd(mat, MAT_FINAL_ASSEMBLY);
83: return 0;
84: }
86: PETSC_EXTERN PetscErrorCode MatlabEngineGet_SeqAIJ(PetscObject obj, void *mengine)
87: {
88: Mat mat = (Mat)obj;
89: mxArray *mmat;
91: mmat = engGetVariable((Engine *)mengine, obj->name);
92: MatSeqAIJFromMatlab(mmat, mat);
93: return 0;
94: }
96: PetscErrorCode MatSolve_Matlab(Mat A, Vec b, Vec x)
97: {
98: const char *_A, *_b, *_x;
100: /* make sure objects have names; use default if not */
101: PetscObjectName((PetscObject)b);
102: PetscObjectName((PetscObject)x);
104: PetscObjectGetName((PetscObject)A, &_A);
105: PetscObjectGetName((PetscObject)b, &_b);
106: PetscObjectGetName((PetscObject)x, &_x);
107: PetscMatlabEnginePut(PETSC_MATLAB_ENGINE_(PetscObjectComm((PetscObject)A)), (PetscObject)b);
108: PetscMatlabEngineEvaluate(PETSC_MATLAB_ENGINE_(PetscObjectComm((PetscObject)A)), "%s = u%s\\(l%s\\(p%s*%s));", _x, _A, _A, _A, _b);
109: PetscMatlabEngineEvaluate(PETSC_MATLAB_ENGINE_(PetscObjectComm((PetscObject)A)), "%s = 0;", _b);
110: /* PetscMatlabEnginePrintOutput(PETSC_MATLAB_ENGINE_(PetscObjectComm((PetscObject)A)),stdout); */
111: PetscMatlabEngineGet(PETSC_MATLAB_ENGINE_(PetscObjectComm((PetscObject)A)), (PetscObject)x);
112: return 0;
113: }
115: PetscErrorCode MatLUFactorNumeric_Matlab(Mat F, Mat A, const MatFactorInfo *info)
116: {
117: size_t len;
118: char *_A, *name;
119: PetscReal dtcol = info->dtcol;
121: if (F->factortype == MAT_FACTOR_ILU || info->dt > 0) {
122: /* the ILU form is not currently registered */
123: if (info->dtcol == PETSC_DEFAULT) dtcol = .01;
124: F->ops->solve = MatSolve_Matlab;
125: F->factortype = MAT_FACTOR_LU;
127: PetscMatlabEnginePut(PETSC_MATLAB_ENGINE_(PetscObjectComm((PetscObject)A)), (PetscObject)A);
128: _A = ((PetscObject)A)->name;
129: PetscMatlabEngineEvaluate(PETSC_MATLAB_ENGINE_(PetscObjectComm((PetscObject)A)), "info_%s = struct('droptol',%g,'thresh',%g);", _A, info->dt, dtcol);
130: PetscMatlabEngineEvaluate(PETSC_MATLAB_ENGINE_(PetscObjectComm((PetscObject)A)), "[l_%s,u_%s,p_%s] = luinc(%s',info_%s);", _A, _A, _A, _A, _A);
131: PetscMatlabEngineEvaluate(PETSC_MATLAB_ENGINE_(PetscObjectComm((PetscObject)A)), "%s = 0;", _A);
133: PetscStrlen(_A, &len);
134: PetscMalloc1(len + 2, &name);
135: sprintf(name, "_%s", _A);
136: PetscObjectSetName((PetscObject)F, name);
137: PetscFree(name);
138: } else {
139: PetscMatlabEnginePut(PETSC_MATLAB_ENGINE_(PetscObjectComm((PetscObject)A)), (PetscObject)A);
140: _A = ((PetscObject)A)->name;
141: PetscMatlabEngineEvaluate(PETSC_MATLAB_ENGINE_(PetscObjectComm((PetscObject)A)), "[l_%s,u_%s,p_%s] = lu(%s',%g);", _A, _A, _A, _A, dtcol);
142: PetscMatlabEngineEvaluate(PETSC_MATLAB_ENGINE_(PetscObjectComm((PetscObject)A)), "%s = 0;", _A);
143: PetscStrlen(_A, &len);
144: PetscMalloc1(len + 2, &name);
145: sprintf(name, "_%s", _A);
146: PetscObjectSetName((PetscObject)F, name);
147: PetscFree(name);
149: F->ops->solve = MatSolve_Matlab;
150: }
151: return 0;
152: }
154: PetscErrorCode MatLUFactorSymbolic_Matlab(Mat F, Mat A, IS r, IS c, const MatFactorInfo *info)
155: {
157: F->ops->lufactornumeric = MatLUFactorNumeric_Matlab;
158: F->assembled = PETSC_TRUE;
159: return 0;
160: }
162: PetscErrorCode MatFactorGetSolverType_seqaij_matlab(Mat A, MatSolverType *type)
163: {
164: *type = MATSOLVERMATLAB;
165: return 0;
166: }
168: PetscErrorCode MatDestroy_matlab(Mat A)
169: {
170: const char *_A;
172: PetscObjectGetName((PetscObject)A, &_A);
173: PetscMatlabEngineEvaluate(PETSC_MATLAB_ENGINE_(PetscObjectComm((PetscObject)A)), "delete %s l_%s u_%s;", _A, _A, _A);
174: PetscObjectComposeFunction((PetscObject)A, "MatFactorGetSolverType_C", NULL);
175: return 0;
176: }
178: PETSC_EXTERN PetscErrorCode MatGetFactor_seqaij_matlab(Mat A, MatFactorType ftype, Mat *F)
179: {
181: MatCreate(PetscObjectComm((PetscObject)A), F);
182: MatSetSizes(*F, A->rmap->n, A->cmap->n, A->rmap->n, A->cmap->n);
183: PetscStrallocpy("matlab", &((PetscObject)*F)->type_name);
184: MatSetUp(*F);
186: (*F)->ops->destroy = MatDestroy_matlab;
187: (*F)->ops->getinfo = MatGetInfo_External;
188: (*F)->trivialsymbolic = PETSC_TRUE;
189: (*F)->ops->lufactorsymbolic = MatLUFactorSymbolic_Matlab;
190: (*F)->ops->ilufactorsymbolic = MatLUFactorSymbolic_Matlab;
192: PetscObjectComposeFunction((PetscObject)(*F), "MatFactorGetSolverType_C", MatFactorGetSolverType_seqaij_matlab);
194: (*F)->factortype = ftype;
195: PetscFree((*F)->solvertype);
196: PetscStrallocpy(MATSOLVERMATLAB, &(*F)->solvertype);
197: return 0;
198: }
200: PETSC_EXTERN PetscErrorCode MatSolverTypeRegister_Matlab(void)
201: {
202: MatSolverTypeRegister(MATSOLVERMATLAB, MATSEQAIJ, MAT_FACTOR_LU, MatGetFactor_seqaij_matlab);
203: return 0;
204: }
206: /* --------------------------------------------------------------------------------*/
208: PetscErrorCode MatView_Info_Matlab(Mat A, PetscViewer viewer)
209: {
210: PetscViewerASCIIPrintf(viewer, "MATLAB run parameters: -- not written yet!\n");
211: return 0;
212: }
214: PetscErrorCode MatView_Matlab(Mat A, PetscViewer viewer)
215: {
216: PetscBool iascii;
218: MatView_SeqAIJ(A, viewer);
219: PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii);
220: if (iascii) {
221: PetscViewerFormat format;
223: PetscViewerGetFormat(viewer, &format);
224: if (format == PETSC_VIEWER_ASCII_FACTOR_INFO) MatView_Info_Matlab(A, viewer);
225: }
226: return 0;
227: }
229: /*MC
230: MATSOLVERMATLAB - "matlab" - Providing direct solver LU for `MATSEQAIJ` matrix via the external package MATLAB.
232: Options Database Keys:
233: . -pc_factor_mat_solver_type matlab - selects MATLAB to do the sparse factorization
235: Note:
236: You must ./configure with the options --with-matlab
238: Level: beginner
240: .seealso: `PCLU`, `PCFactorSetMatSolverType()`, `MatSolverType`
241: M*/