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