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 aij->singlemalloc = PETSC_TRUE; 97 } 98 99 ierr = PetscMemcpy(aij->a,mxGetPr(mmat),aij->nz*sizeof(PetscScalar));CHKERRQ(ierr); 100 /* MATLAB stores by column, not row so we pass in the transpose of the matrix */ 101 i = aij->i; 102 for (k=0; k<n+1; k++) { 103 i[k] = (PetscInt) ii[k]; 104 } 105 j = aij->j; 106 for (k=0; k<nz; k++) { 107 j[k] = (PetscInt) jj[k]; 108 } 109 110 for (k=0; k<mat->rmap->n; k++) { 111 aij->ilen[k] = aij->imax[k] = aij->i[k+1] - aij->i[k]; 112 } 113 114 ierr = MatAssemblyBegin(mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 115 ierr = MatAssemblyEnd(mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 116 PetscFunctionReturn(0); 117 } 118 EXTERN_C_END 119 120 121 EXTERN_C_BEGIN 122 #undef __FUNCT__ 123 #define __FUNCT__ "MatlabEngineGet_SeqAIJ" 124 PetscErrorCode MatlabEngineGet_SeqAIJ(PetscObject obj,void *mengine) 125 { 126 PetscErrorCode ierr; 127 Mat mat = (Mat)obj; 128 mxArray *mmat; 129 130 PetscFunctionBegin; 131 mmat = engGetVariable((Engine *)mengine,obj->name); 132 ierr = MatSeqAIJFromMatlab(mmat,mat);CHKERRQ(ierr); 133 PetscFunctionReturn(0); 134 } 135 EXTERN_C_END 136 137 #undef __FUNCT__ 138 #define __FUNCT__ "MatSolve_Matlab" 139 PetscErrorCode MatSolve_Matlab(Mat A,Vec b,Vec x) 140 { 141 PetscErrorCode ierr; 142 const char *_A,*_b,*_x; 143 144 PetscFunctionBegin; 145 /* make sure objects have names; use default if not */ 146 ierr = PetscObjectName((PetscObject)b);CHKERRQ(ierr); 147 ierr = PetscObjectName((PetscObject)x);CHKERRQ(ierr); 148 149 ierr = PetscObjectGetName((PetscObject)A,&_A);CHKERRQ(ierr); 150 ierr = PetscObjectGetName((PetscObject)b,&_b);CHKERRQ(ierr); 151 ierr = PetscObjectGetName((PetscObject)x,&_x);CHKERRQ(ierr); 152 ierr = PetscMatlabEnginePut(PETSC_MATLAB_ENGINE_(((PetscObject)A)->comm),(PetscObject)b);CHKERRQ(ierr); 153 ierr = PetscMatlabEngineEvaluate(PETSC_MATLAB_ENGINE_(((PetscObject)A)->comm),"%s = u%s\\(l%s\\(p%s*%s));",_x,_A,_A,_A,_b);CHKERRQ(ierr); 154 ierr = PetscMatlabEngineEvaluate(PETSC_MATLAB_ENGINE_(((PetscObject)A)->comm),"%s = 0;",_b);CHKERRQ(ierr); 155 /* ierr = PetscMatlabEnginePrintOutput(PETSC_MATLAB_ENGINE_(((PetscObject)A)->comm),stdout);CHKERRQ(ierr); */ 156 ierr = PetscMatlabEngineGet(PETSC_MATLAB_ENGINE_(((PetscObject)A)->comm),(PetscObject)x);CHKERRQ(ierr); 157 PetscFunctionReturn(0); 158 } 159 160 #undef __FUNCT__ 161 #define __FUNCT__ "MatLUFactorNumeric_Matlab" 162 PetscErrorCode MatLUFactorNumeric_Matlab(Mat F,Mat A,const MatFactorInfo *info) 163 { 164 PetscErrorCode ierr; 165 size_t len; 166 char *_A,*name; 167 PetscReal dtcol = info->dtcol; 168 169 PetscFunctionBegin; 170 if (F->factortype == MAT_FACTOR_ILU || info->dt > 0) { 171 if (info->dtcol == PETSC_DEFAULT) dtcol = .01; 172 F->ops->solve = MatSolve_Matlab; 173 F->factortype = MAT_FACTOR_LU; 174 ierr = PetscMatlabEnginePut(PETSC_MATLAB_ENGINE_(((PetscObject)A)->comm),(PetscObject)A);CHKERRQ(ierr); 175 _A = ((PetscObject)A)->name; 176 ierr = PetscMatlabEngineEvaluate(PETSC_MATLAB_ENGINE_(((PetscObject)A)->comm),"info_%s = struct('droptol',%g,'thresh',%g);",_A,info->dt,dtcol);CHKERRQ(ierr); 177 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); 178 ierr = PetscMatlabEngineEvaluate(PETSC_MATLAB_ENGINE_(((PetscObject)A)->comm),"%s = 0;",_A);CHKERRQ(ierr); 179 180 ierr = PetscStrlen(_A,&len);CHKERRQ(ierr); 181 ierr = PetscMalloc((len+2)*sizeof(char),&name);CHKERRQ(ierr); 182 sprintf(name,"_%s",_A); 183 ierr = PetscObjectSetName((PetscObject)F,name);CHKERRQ(ierr); 184 ierr = PetscFree(name);CHKERRQ(ierr); 185 } else { 186 ierr = PetscMatlabEnginePut(PETSC_MATLAB_ENGINE_(((PetscObject)A)->comm),(PetscObject)A);CHKERRQ(ierr); 187 _A = ((PetscObject)A)->name; 188 ierr = PetscMatlabEngineEvaluate(PETSC_MATLAB_ENGINE_(((PetscObject)A)->comm),"[l_%s,u_%s,p_%s] = lu(%s',%g);",_A,_A,_A,_A,dtcol);CHKERRQ(ierr); 189 ierr = PetscMatlabEngineEvaluate(PETSC_MATLAB_ENGINE_(((PetscObject)A)->comm),"%s = 0;",_A);CHKERRQ(ierr); 190 ierr = PetscStrlen(_A,&len);CHKERRQ(ierr); 191 ierr = PetscMalloc((len+2)*sizeof(char),&name);CHKERRQ(ierr); 192 sprintf(name,"_%s",_A); 193 ierr = PetscObjectSetName((PetscObject)F,name);CHKERRQ(ierr); 194 ierr = PetscFree(name);CHKERRQ(ierr); 195 F->ops->solve = MatSolve_Matlab; 196 } 197 PetscFunctionReturn(0); 198 } 199 200 #undef __FUNCT__ 201 #define __FUNCT__ "MatLUFactorSymbolic_Matlab" 202 PetscErrorCode MatLUFactorSymbolic_Matlab(Mat F,Mat A,IS r,IS c,const MatFactorInfo *info) 203 { 204 PetscFunctionBegin; 205 if (A->cmap->N != A->rmap->N) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"matrix must be square"); 206 F->ops->lufactornumeric = MatLUFactorNumeric_Matlab; 207 F->assembled = PETSC_TRUE; 208 PetscFunctionReturn(0); 209 } 210 211 EXTERN_C_BEGIN 212 #undef __FUNCT__ 213 #define __FUNCT__ "MatFactorGetSolverPackage_seqaij_matlab" 214 PetscErrorCode MatFactorGetSolverPackage_seqaij_matlab(Mat A,const MatSolverPackage *type) 215 { 216 PetscFunctionBegin; 217 *type = MATSOLVERMATLAB; 218 PetscFunctionReturn(0); 219 } 220 EXTERN_C_END 221 222 EXTERN_C_BEGIN 223 #undef __FUNCT__ 224 #define __FUNCT__ "MatGetFactor_seqaij_matlab" 225 PetscErrorCode MatGetFactor_seqaij_matlab(Mat A,MatFactorType ftype,Mat *F) 226 { 227 PetscErrorCode ierr; 228 229 PetscFunctionBegin; 230 if (A->cmap->N != A->rmap->N) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"matrix must be square"); 231 ierr = MatCreate(((PetscObject)A)->comm,F);CHKERRQ(ierr); 232 ierr = MatSetSizes(*F,A->rmap->n,A->cmap->n,A->rmap->n,A->cmap->n);CHKERRQ(ierr); 233 ierr = MatSetType(*F,((PetscObject)A)->type_name);CHKERRQ(ierr); 234 ierr = MatSeqAIJSetPreallocation(*F,0,PETSC_NULL);CHKERRQ(ierr); 235 (*F)->ops->lufactorsymbolic = MatLUFactorSymbolic_Matlab; 236 (*F)->ops->ilufactorsymbolic = MatLUFactorSymbolic_Matlab; 237 ierr = PetscObjectComposeFunctionDynamic((PetscObject)(*F),"MatFactorGetSolverPackage_C","MatFactorGetSolverPackage_seqaij_matlab",MatFactorGetSolverPackage_seqaij_matlab);CHKERRQ(ierr); 238 239 (*F)->factortype = ftype; 240 PetscFunctionReturn(0); 241 } 242 EXTERN_C_END 243 244 /* --------------------------------------------------------------------------------*/ 245 246 #undef __FUNCT__ 247 #define __FUNCT__ "MatFactorInfo_Matlab" 248 PetscErrorCode MatFactorInfo_Matlab(Mat A,PetscViewer viewer) 249 { 250 PetscErrorCode ierr; 251 252 PetscFunctionBegin; 253 ierr = PetscViewerASCIIPrintf(viewer,"MATLAB run parameters: -- not written yet!\n");CHKERRQ(ierr); 254 PetscFunctionReturn(0); 255 } 256 257 #undef __FUNCT__ 258 #define __FUNCT__ "MatView_Matlab" 259 PetscErrorCode MatView_Matlab(Mat A,PetscViewer viewer) 260 { 261 PetscErrorCode ierr; 262 PetscBool iascii; 263 PetscViewerFormat format; 264 265 PetscFunctionBegin; 266 ierr = MatView_SeqAIJ(A,viewer);CHKERRQ(ierr); 267 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr); 268 if (iascii) { 269 ierr = PetscViewerGetFormat(viewer,&format);CHKERRQ(ierr); 270 if (format == PETSC_VIEWER_ASCII_FACTOR_INFO) { 271 ierr = MatFactorInfo_Matlab(A,viewer); 272 } 273 } 274 PetscFunctionReturn(0); 275 } 276 277 278 /*MC 279 MATSOLVERMATLAB - "matlab" - Providing direct solvers (LU and QR) and drop tolerance 280 based ILU factorization (ILUDT) for sequential matrices via the external package MATLAB. 281 282 283 Works with MATSEQAIJ matrices. 284 285 Options Database Keys: 286 . -pc_factor_mat_solver_type matlab - selects MATLAB to do the sparse factorization 287 288 289 Level: beginner 290 291 .seealso: PCLU 292 293 .seealso: PCFactorSetMatSolverPackage(), MatSolverPackage 294 M*/ 295 296