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 PetscErrorCode ierr; 14 Mat_SeqAIJ *aij = (Mat_SeqAIJ*)B->data; 15 mwIndex *ii,*jj; 16 mxArray *mat; 17 PetscInt i; 18 19 PetscFunctionBegin; 20 mat = mxCreateSparse(B->cmap->n,B->rmap->n,aij->nz,mxREAL); 21 ierr = PetscArraycpy(mxGetPr(mat),aij->a,aij->nz);if (ierr) return NULL; 22 /* MATLAB stores by column, not row so we pass in the transpose of the matrix */ 23 jj = mxGetIr(mat); 24 for (i=0; i<aij->nz; i++) jj[i] = aij->j[i]; 25 ii = mxGetJc(mat); 26 for (i=0; i<B->rmap->n+1; i++) ii[i] = aij->i[i]; 27 PetscFunctionReturn(mat); 28 } 29 30 PETSC_EXTERN PetscErrorCode MatlabEnginePut_SeqAIJ(PetscObject obj,void *mengine) 31 { 32 mxArray *mat; 33 34 PetscFunctionBegin; 35 mat = MatSeqAIJToMatlab((Mat)obj);PetscCheckFalse(!mat,PETSC_COMM_SELF,PETSC_ERR_LIB,"Cannot create MATLAB matrix"); 36 CHKERRQ(PetscObjectName(obj)); 37 engPutVariable((Engine*)mengine,obj->name,mat); 38 PetscFunctionReturn(0); 39 } 40 41 PETSC_EXTERN PetscErrorCode MatSeqAIJFromMatlab(mxArray *mmat,Mat mat) 42 { 43 PetscInt nz,n,m,*i,*j,k; 44 mwIndex nnz,nn,nm,*ii,*jj; 45 Mat_SeqAIJ *aij = (Mat_SeqAIJ*)mat->data; 46 47 PetscFunctionBegin; 48 nn = mxGetN(mmat); /* rows of transpose of matrix */ 49 nm = mxGetM(mmat); 50 nnz = (mxGetJc(mmat))[nn]; 51 ii = mxGetJc(mmat); 52 jj = mxGetIr(mmat); 53 n = (PetscInt) nn; 54 m = (PetscInt) nm; 55 nz = (PetscInt) nnz; 56 57 if (mat->rmap->n < 0 && mat->cmap->n < 0) { 58 /* matrix has not yet had its size set */ 59 CHKERRQ(MatSetSizes(mat,n,m,PETSC_DETERMINE,PETSC_DETERMINE)); 60 CHKERRQ(MatSetUp(mat)); 61 } else { 62 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); 63 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); 64 } 65 if (nz != aij->nz) { 66 /* number of nonzeros in matrix has changed, so need new data structure */ 67 CHKERRQ(MatSeqXAIJFreeAIJ(mat,&aij->a,&aij->j,&aij->i)); 68 aij->nz = nz; 69 CHKERRQ(PetscMalloc3(aij->nz,&aij->a,aij->nz,&aij->j,mat->rmap->n+1,&aij->i)); 70 71 aij->singlemalloc = PETSC_TRUE; 72 } 73 74 CHKERRQ(PetscArraycpy(aij->a,mxGetPr(mmat),aij->nz)); 75 /* MATLAB stores by column, not row so we pass in the transpose of the matrix */ 76 i = aij->i; 77 for (k=0; k<n+1; k++) i[k] = (PetscInt) ii[k]; 78 j = aij->j; 79 for (k=0; k<nz; k++) j[k] = (PetscInt) jj[k]; 80 81 for (k=0; k<mat->rmap->n; k++) aij->ilen[k] = aij->imax[k] = aij->i[k+1] - aij->i[k]; 82 83 mat->nonzerostate++; /* since the nonzero structure can change anytime force the Inode information to always be rebuilt */ 84 CHKERRQ(MatAssemblyBegin(mat,MAT_FINAL_ASSEMBLY)); 85 CHKERRQ(MatAssemblyEnd(mat,MAT_FINAL_ASSEMBLY)); 86 PetscFunctionReturn(0); 87 } 88 89 PETSC_EXTERN PetscErrorCode MatlabEngineGet_SeqAIJ(PetscObject obj,void *mengine) 90 { 91 Mat mat = (Mat)obj; 92 mxArray *mmat; 93 94 PetscFunctionBegin; 95 mmat = engGetVariable((Engine*)mengine,obj->name); 96 CHKERRQ(MatSeqAIJFromMatlab(mmat,mat)); 97 PetscFunctionReturn(0); 98 } 99 100 PetscErrorCode MatSolve_Matlab(Mat A,Vec b,Vec x) 101 { 102 const char *_A,*_b,*_x; 103 104 PetscFunctionBegin; 105 /* make sure objects have names; use default if not */ 106 CHKERRQ(PetscObjectName((PetscObject)b)); 107 CHKERRQ(PetscObjectName((PetscObject)x)); 108 109 CHKERRQ(PetscObjectGetName((PetscObject)A,&_A)); 110 CHKERRQ(PetscObjectGetName((PetscObject)b,&_b)); 111 CHKERRQ(PetscObjectGetName((PetscObject)x,&_x)); 112 CHKERRQ(PetscMatlabEnginePut(PETSC_MATLAB_ENGINE_(PetscObjectComm((PetscObject)A)),(PetscObject)b)); 113 CHKERRQ(PetscMatlabEngineEvaluate(PETSC_MATLAB_ENGINE_(PetscObjectComm((PetscObject)A)),"%s = u%s\\(l%s\\(p%s*%s));",_x,_A,_A,_A,_b)); 114 CHKERRQ(PetscMatlabEngineEvaluate(PETSC_MATLAB_ENGINE_(PetscObjectComm((PetscObject)A)),"%s = 0;",_b)); 115 /* CHKERRQ(PetscMatlabEnginePrintOutput(PETSC_MATLAB_ENGINE_(PetscObjectComm((PetscObject)A)),stdout)); */ 116 CHKERRQ(PetscMatlabEngineGet(PETSC_MATLAB_ENGINE_(PetscObjectComm((PetscObject)A)),(PetscObject)x)); 117 PetscFunctionReturn(0); 118 } 119 120 PetscErrorCode MatLUFactorNumeric_Matlab(Mat F,Mat A,const MatFactorInfo *info) 121 { 122 size_t len; 123 char *_A,*name; 124 PetscReal dtcol = info->dtcol; 125 126 PetscFunctionBegin; 127 if (F->factortype == MAT_FACTOR_ILU || info->dt > 0) { 128 /* the ILU form is not currently registered */ 129 if (info->dtcol == PETSC_DEFAULT) dtcol = .01; 130 F->ops->solve = MatSolve_Matlab; 131 F->factortype = MAT_FACTOR_LU; 132 133 CHKERRQ(PetscMatlabEnginePut(PETSC_MATLAB_ENGINE_(PetscObjectComm((PetscObject)A)),(PetscObject)A)); 134 _A = ((PetscObject)A)->name; 135 CHKERRQ(PetscMatlabEngineEvaluate(PETSC_MATLAB_ENGINE_(PetscObjectComm((PetscObject)A)),"info_%s = struct('droptol',%g,'thresh',%g);",_A,info->dt,dtcol)); 136 CHKERRQ(PetscMatlabEngineEvaluate(PETSC_MATLAB_ENGINE_(PetscObjectComm((PetscObject)A)),"[l_%s,u_%s,p_%s] = luinc(%s',info_%s);",_A,_A,_A,_A,_A)); 137 CHKERRQ(PetscMatlabEngineEvaluate(PETSC_MATLAB_ENGINE_(PetscObjectComm((PetscObject)A)),"%s = 0;",_A)); 138 139 CHKERRQ(PetscStrlen(_A,&len)); 140 CHKERRQ(PetscMalloc1(len+2,&name)); 141 sprintf(name,"_%s",_A); 142 CHKERRQ(PetscObjectSetName((PetscObject)F,name)); 143 CHKERRQ(PetscFree(name)); 144 } else { 145 CHKERRQ(PetscMatlabEnginePut(PETSC_MATLAB_ENGINE_(PetscObjectComm((PetscObject)A)),(PetscObject)A)); 146 _A = ((PetscObject)A)->name; 147 CHKERRQ(PetscMatlabEngineEvaluate(PETSC_MATLAB_ENGINE_(PetscObjectComm((PetscObject)A)),"[l_%s,u_%s,p_%s] = lu(%s',%g);",_A,_A,_A,_A,dtcol)); 148 CHKERRQ(PetscMatlabEngineEvaluate(PETSC_MATLAB_ENGINE_(PetscObjectComm((PetscObject)A)),"%s = 0;",_A)); 149 CHKERRQ(PetscStrlen(_A,&len)); 150 CHKERRQ(PetscMalloc1(len+2,&name)); 151 sprintf(name,"_%s",_A); 152 CHKERRQ(PetscObjectSetName((PetscObject)F,name)); 153 CHKERRQ(PetscFree(name)); 154 155 F->ops->solve = MatSolve_Matlab; 156 } 157 PetscFunctionReturn(0); 158 } 159 160 PetscErrorCode MatLUFactorSymbolic_Matlab(Mat F,Mat A,IS r,IS c,const MatFactorInfo *info) 161 { 162 PetscFunctionBegin; 163 PetscCheckFalse(A->cmap->N != A->rmap->N,PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"matrix must be square"); 164 F->ops->lufactornumeric = MatLUFactorNumeric_Matlab; 165 F->assembled = PETSC_TRUE; 166 PetscFunctionReturn(0); 167 } 168 169 PetscErrorCode MatFactorGetSolverType_seqaij_matlab(Mat A,MatSolverType *type) 170 { 171 PetscFunctionBegin; 172 *type = MATSOLVERMATLAB; 173 PetscFunctionReturn(0); 174 } 175 176 PetscErrorCode MatDestroy_matlab(Mat A) 177 { 178 const char *_A; 179 180 PetscFunctionBegin; 181 CHKERRQ(PetscObjectGetName((PetscObject)A,&_A)); 182 CHKERRQ(PetscMatlabEngineEvaluate(PETSC_MATLAB_ENGINE_(PetscObjectComm((PetscObject)A)),"delete %s l_%s u_%s;",_A,_A,_A)); 183 PetscFunctionReturn(0); 184 } 185 186 PETSC_EXTERN PetscErrorCode MatGetFactor_seqaij_matlab(Mat A,MatFactorType ftype,Mat *F) 187 { 188 PetscFunctionBegin; 189 PetscCheckFalse(A->cmap->N != A->rmap->N,PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"matrix must be square"); 190 CHKERRQ(MatCreate(PetscObjectComm((PetscObject)A),F)); 191 CHKERRQ(MatSetSizes(*F,A->rmap->n,A->cmap->n,A->rmap->n,A->cmap->n)); 192 CHKERRQ(PetscStrallocpy("matlab",&((PetscObject)*F)->type_name)); 193 CHKERRQ(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 CHKERRQ(PetscObjectComposeFunction((PetscObject)(*F),"MatFactorGetSolverType_C",MatFactorGetSolverType_seqaij_matlab)); 202 203 (*F)->factortype = ftype; 204 CHKERRQ(PetscFree((*F)->solvertype)); 205 CHKERRQ(PetscStrallocpy(MATSOLVERMATLAB,&(*F)->solvertype)); 206 PetscFunctionReturn(0); 207 } 208 209 PETSC_EXTERN PetscErrorCode MatSolverTypeRegister_Matlab(void) 210 { 211 PetscFunctionBegin; 212 CHKERRQ(MatSolverTypeRegister(MATSOLVERMATLAB,MATSEQAIJ, MAT_FACTOR_LU,MatGetFactor_seqaij_matlab)); 213 PetscFunctionReturn(0); 214 } 215 216 /* --------------------------------------------------------------------------------*/ 217 218 PetscErrorCode MatView_Info_Matlab(Mat A,PetscViewer viewer) 219 { 220 PetscFunctionBegin; 221 CHKERRQ(PetscViewerASCIIPrintf(viewer,"MATLAB run parameters: -- not written yet!\n")); 222 PetscFunctionReturn(0); 223 } 224 225 PetscErrorCode MatView_Matlab(Mat A,PetscViewer viewer) 226 { 227 PetscBool iascii; 228 229 PetscFunctionBegin; 230 CHKERRQ(MatView_SeqAIJ(A,viewer)); 231 CHKERRQ(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii)); 232 if (iascii) { 233 PetscViewerFormat format; 234 235 CHKERRQ(PetscViewerGetFormat(viewer,&format)); 236 if (format == PETSC_VIEWER_ASCII_FACTOR_INFO) CHKERRQ(MatView_Info_Matlab(A,viewer)); 237 } 238 PetscFunctionReturn(0); 239 } 240 241 /*MC 242 MATSOLVERMATLAB - "matlab" - Providing direct solver LU for sequential aij matrix via the external package MATLAB. 243 244 Works with MATSEQAIJ matrices. 245 246 Options Database Keys: 247 . -pc_factor_mat_solver_type matlab - selects MATLAB to do the sparse factorization 248 249 Notes: 250 You must ./configure with the options --with-matlab --with-matlab-engine 251 252 Level: beginner 253 254 .seealso: PCLU 255 256 .seealso: PCFactorSetMatSolverType(), MatSolverType 257 M*/ 258