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