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