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 if (info->dtcol == PETSC_DEFAULT) dtcol = .01; 157 F->ops->solve = MatSolve_Matlab; 158 F->factortype = MAT_FACTOR_LU; 159 160 ierr = PetscMatlabEnginePut(PETSC_MATLAB_ENGINE_(PetscObjectComm((PetscObject)A)),(PetscObject)A);CHKERRQ(ierr); 161 _A = ((PetscObject)A)->name; 162 ierr = PetscMatlabEngineEvaluate(PETSC_MATLAB_ENGINE_(PetscObjectComm((PetscObject)A)),"info_%s = struct('droptol',%g,'thresh',%g);",_A,info->dt,dtcol);CHKERRQ(ierr); 163 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); 164 ierr = PetscMatlabEngineEvaluate(PETSC_MATLAB_ENGINE_(PetscObjectComm((PetscObject)A)),"%s = 0;",_A);CHKERRQ(ierr); 165 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 } else { 172 ierr = PetscMatlabEnginePut(PETSC_MATLAB_ENGINE_(PetscObjectComm((PetscObject)A)),(PetscObject)A);CHKERRQ(ierr); 173 _A = ((PetscObject)A)->name; 174 ierr = PetscMatlabEngineEvaluate(PETSC_MATLAB_ENGINE_(PetscObjectComm((PetscObject)A)),"[l_%s,u_%s,p_%s] = lu(%s',%g);",_A,_A,_A,_A,dtcol);CHKERRQ(ierr); 175 ierr = PetscMatlabEngineEvaluate(PETSC_MATLAB_ENGINE_(PetscObjectComm((PetscObject)A)),"%s = 0;",_A);CHKERRQ(ierr); 176 ierr = PetscStrlen(_A,&len);CHKERRQ(ierr); 177 ierr = PetscMalloc1(len+2,&name);CHKERRQ(ierr); 178 sprintf(name,"_%s",_A); 179 ierr = PetscObjectSetName((PetscObject)F,name);CHKERRQ(ierr); 180 ierr = PetscFree(name);CHKERRQ(ierr); 181 182 F->ops->solve = MatSolve_Matlab; 183 } 184 PetscFunctionReturn(0); 185 } 186 187 #undef __FUNCT__ 188 #define __FUNCT__ "MatLUFactorSymbolic_Matlab" 189 PetscErrorCode MatLUFactorSymbolic_Matlab(Mat F,Mat A,IS r,IS c,const MatFactorInfo *info) 190 { 191 PetscFunctionBegin; 192 if (A->cmap->N != A->rmap->N) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"matrix must be square"); 193 F->ops->lufactornumeric = MatLUFactorNumeric_Matlab; 194 F->assembled = PETSC_TRUE; 195 PetscFunctionReturn(0); 196 } 197 198 #undef __FUNCT__ 199 #define __FUNCT__ "MatFactorGetSolverPackage_seqaij_matlab" 200 PetscErrorCode MatFactorGetSolverPackage_seqaij_matlab(Mat A,const MatSolverPackage *type) 201 { 202 PetscFunctionBegin; 203 *type = MATSOLVERMATLAB; 204 PetscFunctionReturn(0); 205 } 206 207 #undef __FUNCT__ 208 #define __FUNCT__ "MatGetFactor_seqaij_matlab" 209 PETSC_EXTERN PetscErrorCode MatGetFactor_seqaij_matlab(Mat A,MatFactorType ftype,Mat *F) 210 { 211 PetscErrorCode ierr; 212 213 PetscFunctionBegin; 214 if (A->cmap->N != A->rmap->N) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"matrix must be square"); 215 ierr = MatCreate(PetscObjectComm((PetscObject)A),F);CHKERRQ(ierr); 216 ierr = MatSetSizes(*F,A->rmap->n,A->cmap->n,A->rmap->n,A->cmap->n);CHKERRQ(ierr); 217 ierr = MatSetType(*F,((PetscObject)A)->type_name);CHKERRQ(ierr); 218 ierr = MatSeqAIJSetPreallocation(*F,0,NULL);CHKERRQ(ierr); 219 (*F)->ops->lufactorsymbolic = MatLUFactorSymbolic_Matlab; 220 (*F)->ops->ilufactorsymbolic = MatLUFactorSymbolic_Matlab; 221 222 ierr = PetscObjectComposeFunction((PetscObject)(*F),"MatFactorGetSolverPackage_C",MatFactorGetSolverPackage_seqaij_matlab);CHKERRQ(ierr); 223 224 (*F)->factortype = ftype; 225 ierr = PetscFree((*F)->solvertype);CHKERRQ(ierr); 226 ierr = PetscStrallocpy(MATSOLVERMATLAB,&(*F)->solvertype);CHKERRQ(ierr); 227 PetscFunctionReturn(0); 228 } 229 230 231 #undef __FUNCT__ 232 #define __FUNCT__ "MatSolverPackageRegister_Matlab" 233 PETSC_EXTERN PetscErrorCode MatSolverPackageRegister_Matlab(void) 234 { 235 PetscErrorCode ierr; 236 237 PetscFunctionBegin; 238 ierr = MatSolverPackageRegister(MATSOLVERMATLAB,MATSEQAIJ, MAT_FACTOR_LU,MatGetFactor_seqaij_matlab);CHKERRQ(ierr); 239 PetscFunctionReturn(0); 240 } 241 242 /* --------------------------------------------------------------------------------*/ 243 244 #undef __FUNCT__ 245 #define __FUNCT__ "MatFactorInfo_Matlab" 246 PetscErrorCode MatFactorInfo_Matlab(Mat A,PetscViewer viewer) 247 { 248 PetscErrorCode ierr; 249 250 PetscFunctionBegin; 251 ierr = PetscViewerASCIIPrintf(viewer,"MATLAB run parameters: -- not written yet!\n");CHKERRQ(ierr); 252 PetscFunctionReturn(0); 253 } 254 255 #undef __FUNCT__ 256 #define __FUNCT__ "MatView_Matlab" 257 PetscErrorCode MatView_Matlab(Mat A,PetscViewer viewer) 258 { 259 PetscErrorCode ierr; 260 PetscBool iascii; 261 PetscViewerFormat format; 262 263 PetscFunctionBegin; 264 ierr = MatView_SeqAIJ(A,viewer);CHKERRQ(ierr); 265 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr); 266 if (iascii) { 267 ierr = PetscViewerGetFormat(viewer,&format);CHKERRQ(ierr); 268 if (format == PETSC_VIEWER_ASCII_FACTOR_INFO) { 269 ierr = MatFactorInfo_Matlab(A,viewer); 270 } 271 } 272 PetscFunctionReturn(0); 273 } 274 275 276 /*MC 277 MATSOLVERMATLAB - "matlab" - Providing direct solvers (LU and QR) and drop tolerance 278 based ILU factorization (ILUDT) for sequential matrices via the external package MATLAB. 279 280 281 Works with MATSEQAIJ matrices. 282 283 Options Database Keys: 284 . -pc_factor_mat_solver_package matlab - selects MATLAB to do the sparse factorization 285 286 287 Level: beginner 288 289 .seealso: PCLU 290 291 .seealso: PCFactorSetMatSolverPackage(), MatSolverPackage 292 M*/ 293 294