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