1 /*
2 Provides an interface for the MATLAB engine sparse solver
3
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 */
9
MatSeqAIJToMatlab(Mat B)10 static mxArray *MatSeqAIJToMatlab(Mat B)
11 {
12 Mat_SeqAIJ *aij = (Mat_SeqAIJ *)B->data;
13 mwIndex *ii, *jj;
14 mxArray *mat;
15 PetscInt i;
16
17 mat = mxCreateSparse(B->cmap->n, B->rmap->n, aij->nz, mxREAL);
18 if (PetscArraycpy(mxGetPr(mat), aij->a, aij->nz)) return NULL;
19 /* MATLAB stores by column, not row so we pass in the transpose of the matrix */
20 jj = mxGetIr(mat);
21 for (i = 0; i < aij->nz; i++) jj[i] = aij->j[i];
22 ii = mxGetJc(mat);
23 for (i = 0; i < B->rmap->n + 1; i++) ii[i] = aij->i[i];
24 return mat;
25 }
26
MatlabEnginePut_SeqAIJ(PetscObject obj,void * mengine)27 PETSC_EXTERN PetscErrorCode MatlabEnginePut_SeqAIJ(PetscObject obj, void *mengine)
28 {
29 mxArray *mat;
30
31 PetscFunctionBegin;
32 mat = MatSeqAIJToMatlab((Mat)obj);
33 PetscCheck(mat, PETSC_COMM_SELF, PETSC_ERR_LIB, "Cannot create MATLAB matrix");
34 PetscCall(PetscObjectName(obj));
35 engPutVariable((Engine *)mengine, obj->name, mat);
36 PetscFunctionReturn(PETSC_SUCCESS);
37 }
38
MatSeqAIJFromMatlab(mxArray * mmat,Mat mat)39 static 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;
44
45 PetscFunctionBegin;
46 nn = mxGetN(mmat); /* rows of transpose of matrix */
47 nm = mxGetM(mmat);
48 nnz = (mxGetJc(mmat))[nn];
49 ii = mxGetJc(mmat);
50 jj = mxGetIr(mmat);
51 n = (PetscInt)nn;
52 m = (PetscInt)nm;
53 nz = (PetscInt)nnz;
54
55 if (mat->rmap->n < 0 && mat->cmap->n < 0) {
56 /* matrix has not yet had its size set */
57 PetscCall(MatSetSizes(mat, n, m, PETSC_DETERMINE, PETSC_DETERMINE));
58 PetscCall(MatSetUp(mat));
59 } else {
60 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);
61 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);
62 }
63 if (nz != aij->nz) {
64 /* number of nonzeros in matrix has changed, so need new data structure */
65 aij->nz = nz;
66 PetscCall(PetscShmgetAllocateArray(aij->nz, sizeof(PetscScalar), (void **)&aij->a));
67 PetscCall(PetscShmgetAllocateArray(aij->nz, sizeof(PetscInt), (void **)&aij->j));
68 PetscCall(PetscShmgetAllocateArray(mat->rmap->n + 1, sizeof(PetscInt), (void **)&aij->i));
69 aij->free_a = PETSC_TRUE;
70 aij->free_ij = PETSC_TRUE;
71 }
72
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];
79
80 for (k = 0; k < mat->rmap->n; k++) aij->ilen[k] = aij->imax[k] = aij->i[k + 1] - aij->i[k];
81
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 }
87
MatlabEngineGet_SeqAIJ(PetscObject obj,void * mengine)88 PETSC_EXTERN PetscErrorCode MatlabEngineGet_SeqAIJ(PetscObject obj, void *mengine)
89 {
90 Mat mat = (Mat)obj;
91 mxArray *mmat;
92
93 PetscFunctionBegin;
94 mmat = engGetVariable((Engine *)mengine, obj->name);
95 PetscCall(MatSeqAIJFromMatlab(mmat, mat));
96 PetscFunctionReturn(PETSC_SUCCESS);
97 }
98
MatSolve_Matlab(Mat A,Vec b,Vec x)99 static PetscErrorCode MatSolve_Matlab(Mat A, Vec b, Vec x)
100 {
101 const char *_A, *_b, *_x;
102
103 PetscFunctionBegin;
104 /* make sure objects have names; use default if not */
105 PetscCall(PetscObjectName((PetscObject)b));
106 PetscCall(PetscObjectName((PetscObject)x));
107
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 }
118
MatLUFactorNumeric_Matlab(Mat F,Mat A,const MatFactorInfo * info)119 static PetscErrorCode MatLUFactorNumeric_Matlab(Mat F, Mat A, const MatFactorInfo *info)
120 {
121 size_t len;
122 char *_A, *name;
123 PetscReal dtcol = info->dtcol;
124
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;
131
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));
137
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));
153
154 F->ops->solve = MatSolve_Matlab;
155 }
156 PetscFunctionReturn(PETSC_SUCCESS);
157 }
158
MatLUFactorSymbolic_Matlab(Mat F,Mat A,IS r,IS c,const MatFactorInfo * info)159 static 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 }
167
MatFactorGetSolverType_seqaij_matlab(Mat A,MatSolverType * type)168 static PetscErrorCode MatFactorGetSolverType_seqaij_matlab(Mat A, MatSolverType *type)
169 {
170 PetscFunctionBegin;
171 *type = MATSOLVERMATLAB;
172 PetscFunctionReturn(PETSC_SUCCESS);
173 }
174
MatDestroy_matlab(Mat A)175 static PetscErrorCode MatDestroy_matlab(Mat A)
176 {
177 const char *_A;
178
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 }
185
MatGetFactor_seqaij_matlab(Mat A,MatFactorType ftype,Mat * F)186 static 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));
194
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;
200
201 PetscCall(PetscObjectComposeFunction((PetscObject)*F, "MatFactorGetSolverType_C", MatFactorGetSolverType_seqaij_matlab));
202
203 (*F)->factortype = ftype;
204 PetscCall(PetscFree((*F)->solvertype));
205 PetscCall(PetscStrallocpy(MATSOLVERMATLAB, &(*F)->solvertype));
206 PetscFunctionReturn(PETSC_SUCCESS);
207 }
208
MatSolverTypeRegister_Matlab(void)209 PETSC_INTERN PetscErrorCode MatSolverTypeRegister_Matlab(void)
210 {
211 PetscFunctionBegin;
212 PetscCall(MatSolverTypeRegister(MATSOLVERMATLAB, MATSEQAIJ, MAT_FACTOR_LU, MatGetFactor_seqaij_matlab));
213 PetscFunctionReturn(PETSC_SUCCESS);
214 }
215
216 /*MC
217 MATSOLVERMATLAB - "matlab" - Providing direct solver LU for `MATSEQAIJ` matrix via the external package MATLAB.
218
219 Use `./configure` with the options `--with-matlab` to install PETSc with this capability
220
221 Options Database Key:
222 . -pc_factor_mat_solver_type matlab - selects MATLAB to do the sparse factorization
223
224 Level: beginner
225
226 .seealso: [](ch_matrices), `Mat`, `PCLU`, `PCFactorSetMatSolverType()`, `MatSolverType`
227 M*/
228