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