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