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