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