Actual source code: aijmatlab.c
1: /*
2: Provides an interface for the MATLAB engine sparse solver
4: */
5: #include <../src/mat/impls/aij/seq/aij.h>
6: #include <petscmatlab.h>
7: #include <engine.h> /* MATLAB include file */
8: #include <mex.h> /* MATLAB include file */
10: static mxArray *MatSeqAIJToMatlab(Mat B)
11: {
12: Mat_SeqAIJ *aij = (Mat_SeqAIJ *)B->data;
13: mwIndex *ii, *jj;
14: mxArray *mat = mxCreateSparse(B->cmap->n, B->rmap->n, aij->nz, mxREAL);
16: if (PetscArraycpy(mxGetPr(mat), aij->a, aij->nz)) return NULL;
17: /* MATLAB stores by column, not row so we pass in the transpose of the matrix */
18: jj = mxGetIr(mat);
19: for (PetscInt i = 0; i < aij->nz; i++) jj[i] = aij->j[i];
20: ii = mxGetJc(mat);
21: for (PetscInt i = 0; i < B->rmap->n + 1; i++) ii[i] = aij->i[i];
22: return mat;
23: }
25: PETSC_EXTERN PetscErrorCode MatlabEnginePut_SeqAIJ(PetscObject obj, void *mengine)
26: {
27: mxArray *mat;
29: PetscFunctionBegin;
30: mat = MatSeqAIJToMatlab((Mat)obj);
31: PetscCheck(mat, PETSC_COMM_SELF, PETSC_ERR_LIB, "Cannot create MATLAB matrix");
32: PetscCall(PetscObjectName(obj));
33: engPutVariable((Engine *)mengine, obj->name, mat);
34: PetscFunctionReturn(PETSC_SUCCESS);
35: }
37: static PetscErrorCode MatSeqAIJFromMatlab(mxArray *mmat, Mat mat)
38: {
39: PetscInt nz, n, m, *i, *j, k;
40: mwIndex nnz, nn, nm, *ii, *jj;
41: Mat_SeqAIJ *aij = (Mat_SeqAIJ *)mat->data;
43: PetscFunctionBegin;
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: PetscCall(MatSetSizes(mat, n, m, PETSC_DETERMINE, PETSC_DETERMINE));
56: PetscCall(MatSetUp(mat));
57: } else {
58: 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);
59: 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);
60: }
61: if (nz != aij->nz) {
62: /* number of nonzeros in matrix has changed, so need new data structure */
63: aij->nz = nz;
64: PetscCall(PetscShmgetAllocateArray(aij->nz, sizeof(PetscScalar), (void **)&aij->a));
65: PetscCall(PetscShmgetAllocateArray(aij->nz, sizeof(PetscInt), (void **)&aij->j));
66: PetscCall(PetscShmgetAllocateArray(mat->rmap->n + 1, sizeof(PetscInt), (void **)&aij->i));
67: aij->free_a = PETSC_TRUE;
68: aij->free_ij = PETSC_TRUE;
69: }
71: PetscCall(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: PetscCall(MatAssemblyBegin(mat, MAT_FINAL_ASSEMBLY));
82: PetscCall(MatAssemblyEnd(mat, MAT_FINAL_ASSEMBLY));
83: PetscFunctionReturn(PETSC_SUCCESS);
84: }
86: PETSC_EXTERN PetscErrorCode MatlabEngineGet_SeqAIJ(PetscObject obj, void *mengine)
87: {
88: Mat mat = (Mat)obj;
89: mxArray *mmat;
91: PetscFunctionBegin;
92: mmat = engGetVariable((Engine *)mengine, obj->name);
93: PetscCall(MatSeqAIJFromMatlab(mmat, mat));
94: PetscFunctionReturn(PETSC_SUCCESS);
95: }
97: static PetscErrorCode MatSolve_Matlab(Mat A, Vec b, Vec x)
98: {
99: const char *_A, *_b, *_x;
101: PetscFunctionBegin;
102: /* make sure objects have names; use default if not */
103: PetscCall(PetscObjectName((PetscObject)b));
104: PetscCall(PetscObjectName((PetscObject)x));
106: PetscCall(PetscObjectGetName((PetscObject)A, &_A));
107: PetscCall(PetscObjectGetName((PetscObject)b, &_b));
108: PetscCall(PetscObjectGetName((PetscObject)x, &_x));
109: PetscCall(PetscMatlabEnginePut(PETSC_MATLAB_ENGINE_(PetscObjectComm((PetscObject)A)), (PetscObject)b));
110: PetscCall(PetscMatlabEngineEvaluate(PETSC_MATLAB_ENGINE_(PetscObjectComm((PetscObject)A)), "%s = u%s\\(l%s\\(p%s*%s));", _x, _A, _A, _A, _b));
111: PetscCall(PetscMatlabEngineEvaluate(PETSC_MATLAB_ENGINE_(PetscObjectComm((PetscObject)A)), "%s = 0;", _b));
112: /* PetscCall(PetscMatlabEnginePrintOutput(PETSC_MATLAB_ENGINE_(PetscObjectComm((PetscObject)A)),stdout)); */
113: PetscCall(PetscMatlabEngineGet(PETSC_MATLAB_ENGINE_(PetscObjectComm((PetscObject)A)), (PetscObject)x));
114: PetscFunctionReturn(PETSC_SUCCESS);
115: }
117: static PetscErrorCode MatLUFactorNumeric_Matlab(Mat F, Mat A, const MatFactorInfo *info)
118: {
119: size_t len;
120: char *_A, *name;
121: PetscReal dtcol = info->dtcol;
123: PetscFunctionBegin;
124: if (F->factortype == MAT_FACTOR_ILU || info->dt > 0) {
125: /* the ILU form is not currently registered */
126: if (info->dtcol == PETSC_DEFAULT) dtcol = .01;
127: F->ops->solve = MatSolve_Matlab;
128: F->factortype = MAT_FACTOR_LU;
130: PetscCall(PetscMatlabEnginePut(PETSC_MATLAB_ENGINE_(PetscObjectComm((PetscObject)A)), (PetscObject)A));
131: _A = ((PetscObject)A)->name;
132: PetscCall(PetscMatlabEngineEvaluate(PETSC_MATLAB_ENGINE_(PetscObjectComm((PetscObject)A)), "info_%s = struct('droptol',%g,'thresh',%g);", _A, info->dt, dtcol));
133: PetscCall(PetscMatlabEngineEvaluate(PETSC_MATLAB_ENGINE_(PetscObjectComm((PetscObject)A)), "[l_%s,u_%s,p_%s] = luinc(%s',info_%s);", _A, _A, _A, _A, _A));
134: PetscCall(PetscMatlabEngineEvaluate(PETSC_MATLAB_ENGINE_(PetscObjectComm((PetscObject)A)), "%s = 0;", _A));
136: PetscCall(PetscStrlen(_A, &len));
137: PetscCall(PetscMalloc1(len + 2, &name));
138: PetscCall(PetscSNPrintf(name, len + 2, "_%s", _A));
139: PetscCall(PetscObjectSetName((PetscObject)F, name));
140: PetscCall(PetscFree(name));
141: } else {
142: PetscCall(PetscMatlabEnginePut(PETSC_MATLAB_ENGINE_(PetscObjectComm((PetscObject)A)), (PetscObject)A));
143: _A = ((PetscObject)A)->name;
144: PetscCall(PetscMatlabEngineEvaluate(PETSC_MATLAB_ENGINE_(PetscObjectComm((PetscObject)A)), "[l_%s,u_%s,p_%s] = lu(%s',%g);", _A, _A, _A, _A, dtcol));
145: PetscCall(PetscMatlabEngineEvaluate(PETSC_MATLAB_ENGINE_(PetscObjectComm((PetscObject)A)), "%s = 0;", _A));
146: PetscCall(PetscStrlen(_A, &len));
147: PetscCall(PetscMalloc1(len + 2, &name));
148: PetscCall(PetscSNPrintf(name, len + 2, "_%s", _A));
149: PetscCall(PetscObjectSetName((PetscObject)F, name));
150: PetscCall(PetscFree(name));
152: F->ops->solve = MatSolve_Matlab;
153: }
154: PetscFunctionReturn(PETSC_SUCCESS);
155: }
157: static PetscErrorCode MatLUFactorSymbolic_Matlab(Mat F, Mat A, IS r, IS c, const MatFactorInfo *info)
158: {
159: PetscFunctionBegin;
160: PetscCheck(A->cmap->N == A->rmap->N, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "matrix must be square");
161: F->ops->lufactornumeric = MatLUFactorNumeric_Matlab;
162: F->assembled = PETSC_TRUE;
163: PetscFunctionReturn(PETSC_SUCCESS);
164: }
166: static PetscErrorCode MatFactorGetSolverType_seqaij_matlab(Mat A, MatSolverType *type)
167: {
168: PetscFunctionBegin;
169: *type = MATSOLVERMATLAB;
170: PetscFunctionReturn(PETSC_SUCCESS);
171: }
173: static PetscErrorCode MatDestroy_matlab(Mat A)
174: {
175: const char *_A;
177: PetscFunctionBegin;
178: PetscCall(PetscObjectGetName((PetscObject)A, &_A));
179: PetscCall(PetscMatlabEngineEvaluate(PETSC_MATLAB_ENGINE_(PetscObjectComm((PetscObject)A)), "delete %s l_%s u_%s;", _A, _A, _A));
180: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatFactorGetSolverType_C", NULL));
181: PetscFunctionReturn(PETSC_SUCCESS);
182: }
184: static PetscErrorCode MatGetFactor_seqaij_matlab(Mat A, MatFactorType ftype, Mat *F)
185: {
186: PetscFunctionBegin;
187: PetscCheck(A->cmap->N == A->rmap->N, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "matrix must be square");
188: PetscCall(MatCreate(PetscObjectComm((PetscObject)A), F));
189: PetscCall(MatSetSizes(*F, A->rmap->n, A->cmap->n, A->rmap->n, A->cmap->n));
190: PetscCall(PetscStrallocpy("matlab", &((PetscObject)*F)->type_name));
191: PetscCall(MatSetUp(*F));
193: (*F)->ops->destroy = MatDestroy_matlab;
194: (*F)->ops->getinfo = MatGetInfo_External;
195: (*F)->trivialsymbolic = PETSC_TRUE;
196: (*F)->ops->lufactorsymbolic = MatLUFactorSymbolic_Matlab;
197: (*F)->ops->ilufactorsymbolic = MatLUFactorSymbolic_Matlab;
199: PetscCall(PetscObjectComposeFunction((PetscObject)*F, "MatFactorGetSolverType_C", MatFactorGetSolverType_seqaij_matlab));
201: (*F)->factortype = ftype;
202: PetscCall(PetscFree((*F)->solvertype));
203: PetscCall(PetscStrallocpy(MATSOLVERMATLAB, &(*F)->solvertype));
204: PetscFunctionReturn(PETSC_SUCCESS);
205: }
207: PETSC_INTERN PetscErrorCode MatSolverTypeRegister_Matlab(void)
208: {
209: PetscFunctionBegin;
210: PetscCall(MatSolverTypeRegister(MATSOLVERMATLAB, MATSEQAIJ, MAT_FACTOR_LU, MatGetFactor_seqaij_matlab));
211: PetscFunctionReturn(PETSC_SUCCESS);
212: }
214: /*MC
215: MATSOLVERMATLAB - "matlab" - Providing direct solver LU for `MATSEQAIJ` matrix via the external package MATLAB.
217: Use `./configure` with the options `--with-matlab` to install PETSc with this capability
219: Options Database Key:
220: . -pc_factor_mat_solver_type matlab - selects MATLAB to do the sparse factorization
222: Level: beginner
224: .seealso: [](ch_matrices), `Mat`, `PCLU`, `PCFactorSetMatSolverType()`, `MatSolverType`
225: M*/