1 2 /* 3 Provides an interface for the MATLAB engine sparse solver 4 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 */ 10 11 PETSC_EXTERN mxArray *MatSeqAIJToMatlab(Mat B) { 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 27 PETSC_EXTERN PetscErrorCode MatlabEnginePut_SeqAIJ(PetscObject obj, void *mengine) { 28 mxArray *mat; 29 30 PetscFunctionBegin; 31 mat = MatSeqAIJToMatlab((Mat)obj); 32 PetscCheck(mat, PETSC_COMM_SELF, PETSC_ERR_LIB, "Cannot create MATLAB matrix"); 33 PetscCall(PetscObjectName(obj)); 34 engPutVariable((Engine *)mengine, obj->name, mat); 35 PetscFunctionReturn(0); 36 } 37 38 PETSC_EXTERN PetscErrorCode MatSeqAIJFromMatlab(mxArray *mmat, Mat mat) { 39 PetscInt nz, n, m, *i, *j, k; 40 mwIndex nnz, nn, nm, *ii, *jj; 41 Mat_SeqAIJ *aij = (Mat_SeqAIJ *)mat->data; 42 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; 52 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 PetscCall(MatSeqXAIJFreeAIJ(mat, &aij->a, &aij->j, &aij->i)); 64 aij->nz = nz; 65 PetscCall(PetscMalloc3(aij->nz, &aij->a, aij->nz, &aij->j, mat->rmap->n + 1, &aij->i)); 66 67 aij->singlemalloc = PETSC_TRUE; 68 } 69 70 PetscCall(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]; 76 77 for (k = 0; k < mat->rmap->n; k++) aij->ilen[k] = aij->imax[k] = aij->i[k + 1] - aij->i[k]; 78 79 mat->nonzerostate++; /* since the nonzero structure can change anytime force the Inode information to always be rebuilt */ 80 PetscCall(MatAssemblyBegin(mat, MAT_FINAL_ASSEMBLY)); 81 PetscCall(MatAssemblyEnd(mat, MAT_FINAL_ASSEMBLY)); 82 PetscFunctionReturn(0); 83 } 84 85 PETSC_EXTERN PetscErrorCode MatlabEngineGet_SeqAIJ(PetscObject obj, void *mengine) { 86 Mat mat = (Mat)obj; 87 mxArray *mmat; 88 89 PetscFunctionBegin; 90 mmat = engGetVariable((Engine *)mengine, obj->name); 91 PetscCall(MatSeqAIJFromMatlab(mmat, mat)); 92 PetscFunctionReturn(0); 93 } 94 95 PetscErrorCode MatSolve_Matlab(Mat A, Vec b, Vec x) { 96 const char *_A, *_b, *_x; 97 98 PetscFunctionBegin; 99 /* make sure objects have names; use default if not */ 100 PetscCall(PetscObjectName((PetscObject)b)); 101 PetscCall(PetscObjectName((PetscObject)x)); 102 103 PetscCall(PetscObjectGetName((PetscObject)A, &_A)); 104 PetscCall(PetscObjectGetName((PetscObject)b, &_b)); 105 PetscCall(PetscObjectGetName((PetscObject)x, &_x)); 106 PetscCall(PetscMatlabEnginePut(PETSC_MATLAB_ENGINE_(PetscObjectComm((PetscObject)A)), (PetscObject)b)); 107 PetscCall(PetscMatlabEngineEvaluate(PETSC_MATLAB_ENGINE_(PetscObjectComm((PetscObject)A)), "%s = u%s\\(l%s\\(p%s*%s));", _x, _A, _A, _A, _b)); 108 PetscCall(PetscMatlabEngineEvaluate(PETSC_MATLAB_ENGINE_(PetscObjectComm((PetscObject)A)), "%s = 0;", _b)); 109 /* PetscCall(PetscMatlabEnginePrintOutput(PETSC_MATLAB_ENGINE_(PetscObjectComm((PetscObject)A)),stdout)); */ 110 PetscCall(PetscMatlabEngineGet(PETSC_MATLAB_ENGINE_(PetscObjectComm((PetscObject)A)), (PetscObject)x)); 111 PetscFunctionReturn(0); 112 } 113 114 PetscErrorCode MatLUFactorNumeric_Matlab(Mat F, Mat A, const MatFactorInfo *info) { 115 size_t len; 116 char *_A, *name; 117 PetscReal dtcol = info->dtcol; 118 119 PetscFunctionBegin; 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; 125 126 PetscCall(PetscMatlabEnginePut(PETSC_MATLAB_ENGINE_(PetscObjectComm((PetscObject)A)), (PetscObject)A)); 127 _A = ((PetscObject)A)->name; 128 PetscCall(PetscMatlabEngineEvaluate(PETSC_MATLAB_ENGINE_(PetscObjectComm((PetscObject)A)), "info_%s = struct('droptol',%g,'thresh',%g);", _A, info->dt, dtcol)); 129 PetscCall(PetscMatlabEngineEvaluate(PETSC_MATLAB_ENGINE_(PetscObjectComm((PetscObject)A)), "[l_%s,u_%s,p_%s] = luinc(%s',info_%s);", _A, _A, _A, _A, _A)); 130 PetscCall(PetscMatlabEngineEvaluate(PETSC_MATLAB_ENGINE_(PetscObjectComm((PetscObject)A)), "%s = 0;", _A)); 131 132 PetscCall(PetscStrlen(_A, &len)); 133 PetscCall(PetscMalloc1(len + 2, &name)); 134 sprintf(name, "_%s", _A); 135 PetscCall(PetscObjectSetName((PetscObject)F, name)); 136 PetscCall(PetscFree(name)); 137 } else { 138 PetscCall(PetscMatlabEnginePut(PETSC_MATLAB_ENGINE_(PetscObjectComm((PetscObject)A)), (PetscObject)A)); 139 _A = ((PetscObject)A)->name; 140 PetscCall(PetscMatlabEngineEvaluate(PETSC_MATLAB_ENGINE_(PetscObjectComm((PetscObject)A)), "[l_%s,u_%s,p_%s] = lu(%s',%g);", _A, _A, _A, _A, dtcol)); 141 PetscCall(PetscMatlabEngineEvaluate(PETSC_MATLAB_ENGINE_(PetscObjectComm((PetscObject)A)), "%s = 0;", _A)); 142 PetscCall(PetscStrlen(_A, &len)); 143 PetscCall(PetscMalloc1(len + 2, &name)); 144 sprintf(name, "_%s", _A); 145 PetscCall(PetscObjectSetName((PetscObject)F, name)); 146 PetscCall(PetscFree(name)); 147 148 F->ops->solve = MatSolve_Matlab; 149 } 150 PetscFunctionReturn(0); 151 } 152 153 PetscErrorCode MatLUFactorSymbolic_Matlab(Mat F, Mat A, IS r, IS c, const MatFactorInfo *info) { 154 PetscFunctionBegin; 155 PetscCheck(A->cmap->N == A->rmap->N, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "matrix must be square"); 156 F->ops->lufactornumeric = MatLUFactorNumeric_Matlab; 157 F->assembled = PETSC_TRUE; 158 PetscFunctionReturn(0); 159 } 160 161 PetscErrorCode MatFactorGetSolverType_seqaij_matlab(Mat A, MatSolverType *type) { 162 PetscFunctionBegin; 163 *type = MATSOLVERMATLAB; 164 PetscFunctionReturn(0); 165 } 166 167 PetscErrorCode MatDestroy_matlab(Mat A) { 168 const char *_A; 169 170 PetscFunctionBegin; 171 PetscCall(PetscObjectGetName((PetscObject)A, &_A)); 172 PetscCall(PetscMatlabEngineEvaluate(PETSC_MATLAB_ENGINE_(PetscObjectComm((PetscObject)A)), "delete %s l_%s u_%s;", _A, _A, _A)); 173 PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatFactorGetSolverType_C", NULL)); 174 PetscFunctionReturn(0); 175 } 176 177 PETSC_EXTERN PetscErrorCode MatGetFactor_seqaij_matlab(Mat A, MatFactorType ftype, Mat *F) { 178 PetscFunctionBegin; 179 PetscCheck(A->cmap->N == A->rmap->N, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "matrix must be square"); 180 PetscCall(MatCreate(PetscObjectComm((PetscObject)A), F)); 181 PetscCall(MatSetSizes(*F, A->rmap->n, A->cmap->n, A->rmap->n, A->cmap->n)); 182 PetscCall(PetscStrallocpy("matlab", &((PetscObject)*F)->type_name)); 183 PetscCall(MatSetUp(*F)); 184 185 (*F)->ops->destroy = MatDestroy_matlab; 186 (*F)->ops->getinfo = MatGetInfo_External; 187 (*F)->trivialsymbolic = PETSC_TRUE; 188 (*F)->ops->lufactorsymbolic = MatLUFactorSymbolic_Matlab; 189 (*F)->ops->ilufactorsymbolic = MatLUFactorSymbolic_Matlab; 190 191 PetscCall(PetscObjectComposeFunction((PetscObject)(*F), "MatFactorGetSolverType_C", MatFactorGetSolverType_seqaij_matlab)); 192 193 (*F)->factortype = ftype; 194 PetscCall(PetscFree((*F)->solvertype)); 195 PetscCall(PetscStrallocpy(MATSOLVERMATLAB, &(*F)->solvertype)); 196 PetscFunctionReturn(0); 197 } 198 199 PETSC_EXTERN PetscErrorCode MatSolverTypeRegister_Matlab(void) { 200 PetscFunctionBegin; 201 PetscCall(MatSolverTypeRegister(MATSOLVERMATLAB, MATSEQAIJ, MAT_FACTOR_LU, MatGetFactor_seqaij_matlab)); 202 PetscFunctionReturn(0); 203 } 204 205 /* --------------------------------------------------------------------------------*/ 206 207 PetscErrorCode MatView_Info_Matlab(Mat A, PetscViewer viewer) { 208 PetscFunctionBegin; 209 PetscCall(PetscViewerASCIIPrintf(viewer, "MATLAB run parameters: -- not written yet!\n")); 210 PetscFunctionReturn(0); 211 } 212 213 PetscErrorCode MatView_Matlab(Mat A, PetscViewer viewer) { 214 PetscBool iascii; 215 216 PetscFunctionBegin; 217 PetscCall(MatView_SeqAIJ(A, viewer)); 218 PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii)); 219 if (iascii) { 220 PetscViewerFormat format; 221 222 PetscCall(PetscViewerGetFormat(viewer, &format)); 223 if (format == PETSC_VIEWER_ASCII_FACTOR_INFO) PetscCall(MatView_Info_Matlab(A, viewer)); 224 } 225 PetscFunctionReturn(0); 226 } 227 228 /*MC 229 MATSOLVERMATLAB - "matlab" - Providing direct solver LU for `MATSEQAIJ` matrix via the external package MATLAB. 230 231 Options Database Keys: 232 . -pc_factor_mat_solver_type matlab - selects MATLAB to do the sparse factorization 233 234 Note: 235 You must ./configure with the options --with-matlab --with-matlab-engine 236 237 Level: beginner 238 239 .seealso: `PCLU`, `PCFactorSetMatSolverType()`, `MatSolverType` 240 M*/ 241