1 2 /* 3 Provides an interface for the MATLAB engine sparse solver 4 5 */ 6 #include <../src/mat/impls/aij/seq/aij.h> 7 8 #include <engine.h> /* MATLAB include file */ 9 #include <mex.h> /* MATLAB include file */ 10 11 EXTERN_C_BEGIN 12 #undef __FUNCT__ 13 #define __FUNCT__ "MatSeqAIJToMatlab" 14 mxArray *MatSeqAIJToMatlab(Mat B) 15 { 16 PetscErrorCode ierr; 17 Mat_SeqAIJ *aij = (Mat_SeqAIJ*)B->data; 18 mwIndex *ii,*jj; 19 mxArray *mat; 20 PetscInt i; 21 22 PetscFunctionBegin; 23 mat = mxCreateSparse(B->cmap->n,B->rmap->n,aij->nz,mxREAL); 24 ierr = PetscMemcpy(mxGetPr(mat),aij->a,aij->nz*sizeof(PetscScalar)); 25 /* MATLAB stores by column, not row so we pass in the transpose of the matrix */ 26 jj = mxGetIr(mat); 27 for (i=0; i<aij->nz; i++) jj[i] = aij->j[i]; 28 ii = mxGetJc(mat); 29 for (i=0; i<B->rmap->n+1; i++) ii[i] = aij->i[i]; 30 PetscFunctionReturn(mat); 31 } 32 EXTERN_C_END 33 34 35 EXTERN_C_BEGIN 36 #undef __FUNCT__ 37 #define __FUNCT__ "MatlabEnginePut_SeqAIJ" 38 PetscErrorCode MatlabEnginePut_SeqAIJ(PetscObject obj,void *mengine) 39 { 40 PetscErrorCode ierr; 41 mxArray *mat; 42 43 PetscFunctionBegin; 44 mat = MatSeqAIJToMatlab((Mat)obj);if (!mat) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"Cannot create MATLAB matrix"); 45 ierr = PetscObjectName(obj);CHKERRQ(ierr); 46 engPutVariable((Engine*)mengine,obj->name,mat); 47 PetscFunctionReturn(0); 48 } 49 EXTERN_C_END 50 51 EXTERN_C_BEGIN 52 #undef __FUNCT__ 53 #define __FUNCT__ "MatSeqAIJFromMatlab" 54 /*@C 55 MatSeqAIJFromMatlab - Given a MATLAB sparse matrix, fills a SeqAIJ matrix with its transpose. 56 57 Not Collective 58 59 Input Parameters: 60 + mmat - a MATLAB sparse matris 61 - mat - a already created MATSEQAIJ 62 63 Level: intermediate 64 65 @*/ 66 PetscErrorCode MatSeqAIJFromMatlab(mxArray *mmat,Mat mat) 67 { 68 PetscErrorCode ierr; 69 PetscInt nz,n,m,*i,*j,k; 70 mwIndex nnz,nn,nm,*ii,*jj; 71 Mat_SeqAIJ *aij = (Mat_SeqAIJ*)mat->data; 72 73 PetscFunctionBegin; 74 nn = mxGetN(mmat); /* rows of transpose of matrix */ 75 nm = mxGetM(mmat); 76 nnz = (mxGetJc(mmat))[nn]; 77 ii = mxGetJc(mmat); 78 jj = mxGetIr(mmat); 79 n = (PetscInt) nn; 80 m = (PetscInt) nm; 81 nz = (PetscInt) nnz; 82 83 if (mat->rmap->n < 0 && mat->cmap->n < 0) { 84 /* matrix has not yet had its size set */ 85 ierr = MatSetSizes(mat,n,m,PETSC_DETERMINE,PETSC_DETERMINE);CHKERRQ(ierr); 86 ierr = MatSetUp(mat);CHKERRQ(ierr); 87 } else { 88 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); 89 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); 90 } 91 if (nz != aij->nz) { 92 /* number of nonzeros in matrix has changed, so need new data structure */ 93 ierr = MatSeqXAIJFreeAIJ(mat,&aij->a,&aij->j,&aij->i);CHKERRQ(ierr); 94 aij->nz = nz; 95 ierr = PetscMalloc3(aij->nz,PetscScalar,&aij->a,aij->nz,PetscInt,&aij->j,mat->rmap->n+1,PetscInt,&aij->i);CHKERRQ(ierr); 96 97 aij->singlemalloc = PETSC_TRUE; 98 } 99 100 ierr = PetscMemcpy(aij->a,mxGetPr(mmat),aij->nz*sizeof(PetscScalar));CHKERRQ(ierr); 101 /* MATLAB stores by column, not row so we pass in the transpose of the matrix */ 102 i = aij->i; 103 for (k=0; k<n+1; k++) i[k] = (PetscInt) ii[k]; 104 j = aij->j; 105 for (k=0; k<nz; k++) j[k] = (PetscInt) jj[k]; 106 107 for (k=0; k<mat->rmap->n; k++) aij->ilen[k] = aij->imax[k] = aij->i[k+1] - aij->i[k]; 108 109 ierr = MatAssemblyBegin(mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 110 ierr = MatAssemblyEnd(mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 111 PetscFunctionReturn(0); 112 } 113 EXTERN_C_END 114 115 116 EXTERN_C_BEGIN 117 #undef __FUNCT__ 118 #define __FUNCT__ "MatlabEngineGet_SeqAIJ" 119 PetscErrorCode MatlabEngineGet_SeqAIJ(PetscObject obj,void *mengine) 120 { 121 PetscErrorCode ierr; 122 Mat mat = (Mat)obj; 123 mxArray *mmat; 124 125 PetscFunctionBegin; 126 mmat = engGetVariable((Engine*)mengine,obj->name); 127 ierr = MatSeqAIJFromMatlab(mmat,mat);CHKERRQ(ierr); 128 PetscFunctionReturn(0); 129 } 130 EXTERN_C_END 131 132 #undef __FUNCT__ 133 #define __FUNCT__ "MatSolve_Matlab" 134 PetscErrorCode MatSolve_Matlab(Mat A,Vec b,Vec x) 135 { 136 PetscErrorCode ierr; 137 const char *_A,*_b,*_x; 138 139 PetscFunctionBegin; 140 /* make sure objects have names; use default if not */ 141 ierr = PetscObjectName((PetscObject)b);CHKERRQ(ierr); 142 ierr = PetscObjectName((PetscObject)x);CHKERRQ(ierr); 143 144 ierr = PetscObjectGetName((PetscObject)A,&_A);CHKERRQ(ierr); 145 ierr = PetscObjectGetName((PetscObject)b,&_b);CHKERRQ(ierr); 146 ierr = PetscObjectGetName((PetscObject)x,&_x);CHKERRQ(ierr); 147 ierr = PetscMatlabEnginePut(PETSC_MATLAB_ENGINE_(((PetscObject)A)->comm),(PetscObject)b);CHKERRQ(ierr); 148 ierr = PetscMatlabEngineEvaluate(PETSC_MATLAB_ENGINE_(((PetscObject)A)->comm),"%s = u%s\\(l%s\\(p%s*%s));",_x,_A,_A,_A,_b);CHKERRQ(ierr); 149 ierr = PetscMatlabEngineEvaluate(PETSC_MATLAB_ENGINE_(((PetscObject)A)->comm),"%s = 0;",_b);CHKERRQ(ierr); 150 /* ierr = PetscMatlabEnginePrintOutput(PETSC_MATLAB_ENGINE_(((PetscObject)A)->comm),stdout);CHKERRQ(ierr); */ 151 ierr = PetscMatlabEngineGet(PETSC_MATLAB_ENGINE_(((PetscObject)A)->comm),(PetscObject)x);CHKERRQ(ierr); 152 PetscFunctionReturn(0); 153 } 154 155 #undef __FUNCT__ 156 #define __FUNCT__ "MatLUFactorNumeric_Matlab" 157 PetscErrorCode MatLUFactorNumeric_Matlab(Mat F,Mat A,const MatFactorInfo *info) 158 { 159 PetscErrorCode ierr; 160 size_t len; 161 char *_A,*name; 162 PetscReal dtcol = info->dtcol; 163 164 PetscFunctionBegin; 165 if (F->factortype == MAT_FACTOR_ILU || info->dt > 0) { 166 if (info->dtcol == PETSC_DEFAULT) dtcol = .01; 167 F->ops->solve = MatSolve_Matlab; 168 F->factortype = MAT_FACTOR_LU; 169 170 ierr = PetscMatlabEnginePut(PETSC_MATLAB_ENGINE_(((PetscObject)A)->comm),(PetscObject)A);CHKERRQ(ierr); 171 _A = ((PetscObject)A)->name; 172 ierr = PetscMatlabEngineEvaluate(PETSC_MATLAB_ENGINE_(((PetscObject)A)->comm),"info_%s = struct('droptol',%g,'thresh',%g);",_A,info->dt,dtcol);CHKERRQ(ierr); 173 ierr = PetscMatlabEngineEvaluate(PETSC_MATLAB_ENGINE_(((PetscObject)A)->comm),"[l_%s,u_%s,p_%s] = luinc(%s',info_%s);",_A,_A,_A,_A,_A);CHKERRQ(ierr); 174 ierr = PetscMatlabEngineEvaluate(PETSC_MATLAB_ENGINE_(((PetscObject)A)->comm),"%s = 0;",_A);CHKERRQ(ierr); 175 176 ierr = PetscStrlen(_A,&len);CHKERRQ(ierr); 177 ierr = PetscMalloc((len+2)*sizeof(char),&name);CHKERRQ(ierr); 178 sprintf(name,"_%s",_A); 179 ierr = PetscObjectSetName((PetscObject)F,name);CHKERRQ(ierr); 180 ierr = PetscFree(name);CHKERRQ(ierr); 181 } else { 182 ierr = PetscMatlabEnginePut(PETSC_MATLAB_ENGINE_(((PetscObject)A)->comm),(PetscObject)A);CHKERRQ(ierr); 183 _A = ((PetscObject)A)->name; 184 ierr = PetscMatlabEngineEvaluate(PETSC_MATLAB_ENGINE_(((PetscObject)A)->comm),"[l_%s,u_%s,p_%s] = lu(%s',%g);",_A,_A,_A,_A,dtcol);CHKERRQ(ierr); 185 ierr = PetscMatlabEngineEvaluate(PETSC_MATLAB_ENGINE_(((PetscObject)A)->comm),"%s = 0;",_A);CHKERRQ(ierr); 186 ierr = PetscStrlen(_A,&len);CHKERRQ(ierr); 187 ierr = PetscMalloc((len+2)*sizeof(char),&name);CHKERRQ(ierr); 188 sprintf(name,"_%s",_A); 189 ierr = PetscObjectSetName((PetscObject)F,name);CHKERRQ(ierr); 190 ierr = PetscFree(name);CHKERRQ(ierr); 191 192 F->ops->solve = MatSolve_Matlab; 193 } 194 PetscFunctionReturn(0); 195 } 196 197 #undef __FUNCT__ 198 #define __FUNCT__ "MatLUFactorSymbolic_Matlab" 199 PetscErrorCode MatLUFactorSymbolic_Matlab(Mat F,Mat A,IS r,IS c,const MatFactorInfo *info) 200 { 201 PetscFunctionBegin; 202 if (A->cmap->N != A->rmap->N) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"matrix must be square"); 203 F->ops->lufactornumeric = MatLUFactorNumeric_Matlab; 204 F->assembled = PETSC_TRUE; 205 PetscFunctionReturn(0); 206 } 207 208 EXTERN_C_BEGIN 209 #undef __FUNCT__ 210 #define __FUNCT__ "MatFactorGetSolverPackage_seqaij_matlab" 211 PetscErrorCode MatFactorGetSolverPackage_seqaij_matlab(Mat A,const MatSolverPackage *type) 212 { 213 PetscFunctionBegin; 214 *type = MATSOLVERMATLAB; 215 PetscFunctionReturn(0); 216 } 217 EXTERN_C_END 218 219 EXTERN_C_BEGIN 220 #undef __FUNCT__ 221 #define __FUNCT__ "MatGetFactor_seqaij_matlab" 222 PetscErrorCode MatGetFactor_seqaij_matlab(Mat A,MatFactorType ftype,Mat *F) 223 { 224 PetscErrorCode ierr; 225 226 PetscFunctionBegin; 227 if (A->cmap->N != A->rmap->N) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"matrix must be square"); 228 ierr = MatCreate(((PetscObject)A)->comm,F);CHKERRQ(ierr); 229 ierr = MatSetSizes(*F,A->rmap->n,A->cmap->n,A->rmap->n,A->cmap->n);CHKERRQ(ierr); 230 ierr = MatSetType(*F,((PetscObject)A)->type_name);CHKERRQ(ierr); 231 ierr = MatSeqAIJSetPreallocation(*F,0,PETSC_NULL);CHKERRQ(ierr); 232 (*F)->ops->lufactorsymbolic = MatLUFactorSymbolic_Matlab; 233 (*F)->ops->ilufactorsymbolic = MatLUFactorSymbolic_Matlab; 234 235 ierr = PetscObjectComposeFunctionDynamic((PetscObject)(*F),"MatFactorGetSolverPackage_C","MatFactorGetSolverPackage_seqaij_matlab",MatFactorGetSolverPackage_seqaij_matlab);CHKERRQ(ierr); 236 237 (*F)->factortype = ftype; 238 PetscFunctionReturn(0); 239 } 240 EXTERN_C_END 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_type 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