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 { 13 Mat_SeqAIJ *aij = (Mat_SeqAIJ*)B->data; 14 mwIndex *ii,*jj; 15 mxArray *mat; 16 PetscInt i; 17 18 PetscFunctionBegin; 19 mat = mxCreateSparse(B->cmap->n,B->rmap->n,aij->nz,mxREAL); 20 if (PetscArraycpy(mxGetPr(mat),aij->a,aij->nz)) return NULL; 21 /* MATLAB stores by column, not row so we pass in the transpose of the matrix */ 22 jj = mxGetIr(mat); 23 for (i=0; i<aij->nz; i++) jj[i] = aij->j[i]; 24 ii = mxGetJc(mat); 25 for (i=0; i<B->rmap->n+1; i++) ii[i] = aij->i[i]; 26 PetscFunctionReturn(mat); 27 } 28 29 PETSC_EXTERN PetscErrorCode MatlabEnginePut_SeqAIJ(PetscObject obj,void *mengine) 30 { 31 mxArray *mat; 32 33 PetscFunctionBegin; 34 mat = MatSeqAIJToMatlab((Mat)obj);PetscCheck(mat,PETSC_COMM_SELF,PETSC_ERR_LIB,"Cannot create MATLAB matrix"); 35 PetscCall(PetscObjectName(obj)); 36 engPutVariable((Engine*)mengine,obj->name,mat); 37 PetscFunctionReturn(0); 38 } 39 40 PETSC_EXTERN PetscErrorCode MatSeqAIJFromMatlab(mxArray *mmat,Mat mat) 41 { 42 PetscInt nz,n,m,*i,*j,k; 43 mwIndex nnz,nn,nm,*ii,*jj; 44 Mat_SeqAIJ *aij = (Mat_SeqAIJ*)mat->data; 45 46 PetscFunctionBegin; 47 nn = mxGetN(mmat); /* rows of transpose of matrix */ 48 nm = mxGetM(mmat); 49 nnz = (mxGetJc(mmat))[nn]; 50 ii = mxGetJc(mmat); 51 jj = mxGetIr(mmat); 52 n = (PetscInt) nn; 53 m = (PetscInt) nm; 54 nz = (PetscInt) nnz; 55 56 if (mat->rmap->n < 0 && mat->cmap->n < 0) { 57 /* matrix has not yet had its size set */ 58 PetscCall(MatSetSizes(mat,n,m,PETSC_DETERMINE,PETSC_DETERMINE)); 59 PetscCall(MatSetUp(mat)); 60 } else { 61 PetscCheckFalse(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); 62 PetscCheckFalse(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); 63 } 64 if (nz != aij->nz) { 65 /* number of nonzeros in matrix has changed, so need new data structure */ 66 PetscCall(MatSeqXAIJFreeAIJ(mat,&aij->a,&aij->j,&aij->i)); 67 aij->nz = nz; 68 PetscCall(PetscMalloc3(aij->nz,&aij->a,aij->nz,&aij->j,mat->rmap->n+1,&aij->i)); 69 70 aij->singlemalloc = 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(0); 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(0); 97 } 98 99 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(0); 117 } 118 119 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 sprintf(name,"_%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 sprintf(name,"_%s",_A); 151 PetscCall(PetscObjectSetName((PetscObject)F,name)); 152 PetscCall(PetscFree(name)); 153 154 F->ops->solve = MatSolve_Matlab; 155 } 156 PetscFunctionReturn(0); 157 } 158 159 PetscErrorCode MatLUFactorSymbolic_Matlab(Mat F,Mat A,IS r,IS c,const MatFactorInfo *info) 160 { 161 PetscFunctionBegin; 162 PetscCheckFalse(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(0); 166 } 167 168 PetscErrorCode MatFactorGetSolverType_seqaij_matlab(Mat A,MatSolverType *type) 169 { 170 PetscFunctionBegin; 171 *type = MATSOLVERMATLAB; 172 PetscFunctionReturn(0); 173 } 174 175 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 PetscFunctionReturn(0); 183 } 184 185 PETSC_EXTERN PetscErrorCode MatGetFactor_seqaij_matlab(Mat A,MatFactorType ftype,Mat *F) 186 { 187 PetscFunctionBegin; 188 PetscCheckFalse(A->cmap->N != A->rmap->N,PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"matrix must be square"); 189 PetscCall(MatCreate(PetscObjectComm((PetscObject)A),F)); 190 PetscCall(MatSetSizes(*F,A->rmap->n,A->cmap->n,A->rmap->n,A->cmap->n)); 191 PetscCall(PetscStrallocpy("matlab",&((PetscObject)*F)->type_name)); 192 PetscCall(MatSetUp(*F)); 193 194 (*F)->ops->destroy = MatDestroy_matlab; 195 (*F)->ops->getinfo = MatGetInfo_External; 196 (*F)->trivialsymbolic = PETSC_TRUE; 197 (*F)->ops->lufactorsymbolic = MatLUFactorSymbolic_Matlab; 198 (*F)->ops->ilufactorsymbolic = MatLUFactorSymbolic_Matlab; 199 200 PetscCall(PetscObjectComposeFunction((PetscObject)(*F),"MatFactorGetSolverType_C",MatFactorGetSolverType_seqaij_matlab)); 201 202 (*F)->factortype = ftype; 203 PetscCall(PetscFree((*F)->solvertype)); 204 PetscCall(PetscStrallocpy(MATSOLVERMATLAB,&(*F)->solvertype)); 205 PetscFunctionReturn(0); 206 } 207 208 PETSC_EXTERN PetscErrorCode MatSolverTypeRegister_Matlab(void) 209 { 210 PetscFunctionBegin; 211 PetscCall(MatSolverTypeRegister(MATSOLVERMATLAB,MATSEQAIJ, MAT_FACTOR_LU,MatGetFactor_seqaij_matlab)); 212 PetscFunctionReturn(0); 213 } 214 215 /* --------------------------------------------------------------------------------*/ 216 217 PetscErrorCode MatView_Info_Matlab(Mat A,PetscViewer viewer) 218 { 219 PetscFunctionBegin; 220 PetscCall(PetscViewerASCIIPrintf(viewer,"MATLAB run parameters: -- not written yet!\n")); 221 PetscFunctionReturn(0); 222 } 223 224 PetscErrorCode MatView_Matlab(Mat A,PetscViewer viewer) 225 { 226 PetscBool iascii; 227 228 PetscFunctionBegin; 229 PetscCall(MatView_SeqAIJ(A,viewer)); 230 PetscCall(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii)); 231 if (iascii) { 232 PetscViewerFormat format; 233 234 PetscCall(PetscViewerGetFormat(viewer,&format)); 235 if (format == PETSC_VIEWER_ASCII_FACTOR_INFO) PetscCall(MatView_Info_Matlab(A,viewer)); 236 } 237 PetscFunctionReturn(0); 238 } 239 240 /*MC 241 MATSOLVERMATLAB - "matlab" - Providing direct solver LU for sequential aij matrix via the external package MATLAB. 242 243 Works with MATSEQAIJ matrices. 244 245 Options Database Keys: 246 . -pc_factor_mat_solver_type matlab - selects MATLAB to do the sparse factorization 247 248 Notes: 249 You must ./configure with the options --with-matlab --with-matlab-engine 250 251 Level: beginner 252 253 .seealso: PCLU 254 255 .seealso: PCFactorSetMatSolverType(), MatSolverType 256 M*/ 257