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__ "MatlabEngineGet_SeqAIJ" 42 PetscErrorCode PETSCMAT_DLLEXPORT MatlabEngineGet_SeqAIJ(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 PetscReal dtcol = info->dtcol; 106 107 PetscFunctionBegin; 108 if (F->factor == MAT_FACTOR_ILU || info->dt > 0) { 109 if (info->dtcol == PETSC_DEFAULT) dtcol = .01; 110 F->ops->solve = MatSolve_Matlab; 111 F->factor = MAT_FACTOR_LU; 112 ierr = PetscMatlabEnginePut(PETSC_MATLAB_ENGINE_(((PetscObject)A)->comm),(PetscObject)A);CHKERRQ(ierr); 113 _A = ((PetscObject)A)->name; 114 ierr = PetscMatlabEngineEvaluate(PETSC_MATLAB_ENGINE_(((PetscObject)A)->comm),"info_%s = struct('droptol',%g,'thresh',%g);",_A,info->dt,dtcol);CHKERRQ(ierr); 115 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); 116 ierr = PetscMatlabEngineEvaluate(PETSC_MATLAB_ENGINE_(((PetscObject)A)->comm),"%s = 0;",_A);CHKERRQ(ierr); 117 118 ierr = PetscStrlen(_A,&len);CHKERRQ(ierr); 119 ierr = PetscMalloc((len+2)*sizeof(char),&name);CHKERRQ(ierr); 120 sprintf(name,"_%s",_A); 121 ierr = PetscObjectSetName((PetscObject)F,name);CHKERRQ(ierr); 122 ierr = PetscFree(name);CHKERRQ(ierr); 123 } else { 124 ierr = PetscMatlabEnginePut(PETSC_MATLAB_ENGINE_(((PetscObject)A)->comm),(PetscObject)A);CHKERRQ(ierr); 125 _A = ((PetscObject)A)->name; 126 ierr = PetscMatlabEngineEvaluate(PETSC_MATLAB_ENGINE_(((PetscObject)A)->comm),"[l_%s,u_%s,p_%s] = lu(%s',%g);",_A,_A,_A,_A,dtcol);CHKERRQ(ierr); 127 ierr = PetscMatlabEngineEvaluate(PETSC_MATLAB_ENGINE_(((PetscObject)A)->comm),"%s = 0;",_A);CHKERRQ(ierr); 128 ierr = PetscStrlen(_A,&len);CHKERRQ(ierr); 129 ierr = PetscMalloc((len+2)*sizeof(char),&name);CHKERRQ(ierr); 130 sprintf(name,"_%s",_A); 131 ierr = PetscObjectSetName((PetscObject)F,name);CHKERRQ(ierr); 132 ierr = PetscFree(name);CHKERRQ(ierr); 133 F->ops->solve = MatSolve_Matlab; 134 } 135 PetscFunctionReturn(0); 136 } 137 138 #undef __FUNCT__ 139 #define __FUNCT__ "MatLUFactorSymbolic_Matlab" 140 PetscErrorCode MatLUFactorSymbolic_Matlab(Mat F,Mat A,IS r,IS c,const MatFactorInfo *info) 141 { 142 PetscFunctionBegin; 143 if (A->cmap->N != A->rmap->N) SETERRQ(PETSC_ERR_ARG_SIZ,"matrix must be square"); 144 F->ops->lufactornumeric = MatLUFactorNumeric_Matlab; 145 F->assembled = PETSC_TRUE; 146 PetscFunctionReturn(0); 147 } 148 149 EXTERN_C_BEGIN 150 #undef __FUNCT__ 151 #define __FUNCT__ "MatFactorGetSolverPackage_seqaij_matlab" 152 PetscErrorCode MatFactorGetSolverPackage_seqaij_matlab(Mat A,const MatSolverPackage *type) 153 { 154 PetscFunctionBegin; 155 *type = MAT_SOLVER_MATLAB; 156 PetscFunctionReturn(0); 157 } 158 EXTERN_C_END 159 160 #undef __FUNCT__ 161 #define __FUNCT__ "MatGetFactor_seqaij_matlab" 162 PetscErrorCode MatGetFactor_seqaij_matlab(Mat A,MatFactorType ftype,Mat *F) 163 { 164 PetscErrorCode ierr; 165 166 PetscFunctionBegin; 167 if (A->cmap->N != A->rmap->N) SETERRQ(PETSC_ERR_ARG_SIZ,"matrix must be square"); 168 ierr = MatCreate(((PetscObject)A)->comm,F);CHKERRQ(ierr); 169 ierr = MatSetSizes(*F,A->rmap->n,A->cmap->n,A->rmap->n,A->cmap->n);CHKERRQ(ierr); 170 ierr = MatSetType(*F,((PetscObject)A)->type_name);CHKERRQ(ierr); 171 ierr = MatSeqAIJSetPreallocation(*F,0,PETSC_NULL);CHKERRQ(ierr); 172 (*F)->ops->lufactorsymbolic = MatLUFactorSymbolic_Matlab; 173 (*F)->ops->ilufactorsymbolic = MatLUFactorSymbolic_Matlab; 174 ierr = PetscObjectComposeFunctionDynamic((PetscObject)(*F),"MatFactorGetSolverPackage_C","MatFactorGetSolverPackage_seqaij_matlab",MatFactorGetSolverPackage_seqaij_matlab);CHKERRQ(ierr); 175 176 (*F)->factor = ftype; 177 PetscFunctionReturn(0); 178 } 179 180 181 /* --------------------------------------------------------------------------------*/ 182 183 #undef __FUNCT__ 184 #define __FUNCT__ "MatFactorInfo_Matlab" 185 PetscErrorCode MatFactorInfo_Matlab(Mat A,PetscViewer viewer) 186 { 187 PetscErrorCode ierr; 188 189 PetscFunctionBegin; 190 ierr = PetscViewerASCIIPrintf(viewer,"Matlab run parameters: -- not written yet!\n");CHKERRQ(ierr); 191 PetscFunctionReturn(0); 192 } 193 194 #undef __FUNCT__ 195 #define __FUNCT__ "MatView_Matlab" 196 PetscErrorCode MatView_Matlab(Mat A,PetscViewer viewer) 197 { 198 PetscErrorCode ierr; 199 PetscTruth iascii; 200 PetscViewerFormat format; 201 202 PetscFunctionBegin; 203 ierr = MatView_SeqAIJ(A,viewer);CHKERRQ(ierr); 204 ierr = PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_ASCII,&iascii);CHKERRQ(ierr); 205 if (iascii) { 206 ierr = PetscViewerGetFormat(viewer,&format);CHKERRQ(ierr); 207 if (format == PETSC_VIEWER_ASCII_FACTOR_INFO) { 208 ierr = MatFactorInfo_Matlab(A,viewer); 209 } 210 } 211 PetscFunctionReturn(0); 212 } 213 214 215 /*MC 216 MAT_SOLVER_MATLAB - "matlab" - Providing direct solvers (LU and QR) and drop tolerance 217 based ILU factorization (ILUDT) for sequential matrices via the external package Matlab. 218 219 220 Works with MATSEQAIJ matrices. 221 222 Options Database Keys: 223 . -pc_factor_mat_solver_type matlab - selects Matlab to do the sparse factorization 224 225 226 Level: beginner 227 228 .seealso: PCLU 229 230 .seealso: PCFactorSetMatSolverPackage(), MatSolverPackage 231 M*/ 232 233