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