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 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 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 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 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 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 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 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 168 static PetscErrorCode MatFactorGetSolverType_seqaij_matlab(Mat A, MatSolverType *type) 169 { 170 PetscFunctionBegin; 171 *type = MATSOLVERMATLAB; 172 PetscFunctionReturn(PETSC_SUCCESS); 173 } 174 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 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 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