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