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