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