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 ierr = PetscMatlabEnginePut(PETSC_MATLAB_ENGINE_(((PetscObject)A)->comm),(PetscObject)A);CHKERRQ(ierr); 108 _A = ((PetscObject)A)->name; 109 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); 110 ierr = PetscMatlabEngineEvaluate(PETSC_MATLAB_ENGINE_(((PetscObject)A)->comm),"%s = 0;",_A);CHKERRQ(ierr); 111 ierr = PetscStrlen(_A,&len);CHKERRQ(ierr); 112 ierr = PetscMalloc((len+2)*sizeof(char),&name);CHKERRQ(ierr); 113 sprintf(name,"_%s",_A); 114 ierr = PetscObjectSetName((PetscObject)F,name);CHKERRQ(ierr); 115 ierr = PetscFree(name);CHKERRQ(ierr); 116 F->ops->solve = MatSolve_Matlab; 117 PetscFunctionReturn(0); 118 } 119 120 #undef __FUNCT__ 121 #define __FUNCT__ "MatLUFactorSymbolic_Matlab" 122 PetscErrorCode MatLUFactorSymbolic_Matlab(Mat F,Mat A,IS r,IS c,const MatFactorInfo *info) 123 { 124 PetscFunctionBegin; 125 if (A->cmap->N != A->rmap->N) SETERRQ(PETSC_ERR_ARG_SIZ,"matrix must be square"); 126 F->ops->lufactornumeric = MatLUFactorNumeric_Matlab; 127 PetscFunctionReturn(0); 128 } 129 130 EXTERN_C_BEGIN 131 #undef __FUNCT__ 132 #define __FUNCT__ "MatFactorGetSolverPackage_seqaij_matlab" 133 PetscErrorCode MatFactorGetSolverPackage_seqaij_matlab(Mat A,const MatSolverPackage *type) 134 { 135 PetscFunctionBegin; 136 *type = MAT_SOLVER_MATLAB; 137 PetscFunctionReturn(0); 138 } 139 EXTERN_C_END 140 141 #undef __FUNCT__ 142 #define __FUNCT__ "MatGetFactor_seqaij_matlab" 143 PetscErrorCode MatGetFactor_seqaij_matlab(Mat A,MatFactorType ftype,Mat *F) 144 { 145 PetscErrorCode ierr; 146 147 PetscFunctionBegin; 148 if (A->cmap->N != A->rmap->N) SETERRQ(PETSC_ERR_ARG_SIZ,"matrix must be square"); 149 ierr = MatCreate(((PetscObject)A)->comm,F);CHKERRQ(ierr); 150 ierr = MatSetSizes(*F,A->rmap->n,A->cmap->n,A->rmap->n,A->cmap->n);CHKERRQ(ierr); 151 ierr = MatSetType(*F,((PetscObject)A)->type_name);CHKERRQ(ierr); 152 ierr = MatSeqAIJSetPreallocation(*F,0,PETSC_NULL);CHKERRQ(ierr); 153 (*F)->ops->lufactorsymbolic = MatLUFactorSymbolic_Matlab; 154 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatFactorGetSolverPackage_C","MatFactorGetSolverPackage_seqaij_matlab",MatFactorGetSolverPackage_seqaij_matlab);CHKERRQ(ierr); 155 156 (*F)->factor = MAT_FACTOR_LU; 157 PetscFunctionReturn(0); 158 } 159 160 161 /* --------------------------------------------------------------------------------*/ 162 #undef __FUNCT__ 163 #define __FUNCT__ "MatILUDTFactor_Matlab" 164 PetscErrorCode MatILUDTFactor_Matlab(Mat A,IS isrow,IS iscol,const MatFactorInfo *info) 165 { 166 PetscErrorCode ierr; 167 size_t len; 168 char *_A,*name; 169 PetscReal dt,dtcol; 170 Mat F; 171 172 PetscFunctionBegin; 173 if (info->dt == PETSC_DEFAULT) dt = .005; 174 if (info->dtcol == PETSC_DEFAULT) dtcol = .01; 175 ierr = MatGetFactor(A,MAT_SOLVER_MATLAB,MAT_FACTOR_ILU,&F);CHKERRQ(ierr); 176 F->ops->solve = MatSolve_Matlab; 177 F->factor = MAT_FACTOR_LU; 178 ierr = PetscMatlabEnginePut(PETSC_MATLAB_ENGINE_(((PetscObject)A)->comm),(PetscObject)A);CHKERRQ(ierr); 179 _A = ((PetscObject)A)->name; 180 ierr = PetscMatlabEngineEvaluate(PETSC_MATLAB_ENGINE_(((PetscObject)A)->comm),"info_%s = struct('droptol',%g,'thresh',%g);",_A,dt,dtcol);CHKERRQ(ierr); 181 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); 182 ierr = PetscMatlabEngineEvaluate(PETSC_MATLAB_ENGINE_(((PetscObject)A)->comm),"%s = 0;",_A);CHKERRQ(ierr); 183 184 ierr = PetscStrlen(_A,&len);CHKERRQ(ierr); 185 ierr = PetscMalloc((len+2)*sizeof(char),&name);CHKERRQ(ierr); 186 sprintf(name,"_%s",_A); 187 ierr = PetscObjectSetName((PetscObject)F,name);CHKERRQ(ierr); 188 ierr = PetscFree(name);CHKERRQ(ierr); 189 PetscFunctionReturn(0); 190 } 191 192 #undef __FUNCT__ 193 #define __FUNCT__ "MatFactorInfo_Matlab" 194 PetscErrorCode MatFactorInfo_Matlab(Mat A,PetscViewer viewer) 195 { 196 PetscErrorCode ierr; 197 198 PetscFunctionBegin; 199 ierr = PetscViewerASCIIPrintf(viewer,"Matlab run parameters: -- not written yet!\n");CHKERRQ(ierr); 200 PetscFunctionReturn(0); 201 } 202 203 #undef __FUNCT__ 204 #define __FUNCT__ "MatView_Matlab" 205 PetscErrorCode MatView_Matlab(Mat A,PetscViewer viewer) 206 { 207 PetscErrorCode ierr; 208 PetscTruth iascii; 209 PetscViewerFormat format; 210 211 PetscFunctionBegin; 212 ierr = MatView_SeqAIJ(A,viewer);CHKERRQ(ierr); 213 ierr = PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_ASCII,&iascii);CHKERRQ(ierr); 214 if (iascii) { 215 ierr = PetscViewerGetFormat(viewer,&format);CHKERRQ(ierr); 216 if (format == PETSC_VIEWER_ASCII_FACTOR_INFO) { 217 ierr = MatFactorInfo_Matlab(A,viewer); 218 } 219 } 220 PetscFunctionReturn(0); 221 } 222 223 224 /*MC 225 MATMATLAB - MATMATLAB = "matlab" - A matrix type providing direct solvers (LU and QR) and drop tolerance 226 based ILU factorization (ILUDT) for sequential matrices via the external package Matlab. 227 228 If Matlab is instaled (see the manual for 229 instructions on how to declare the existence of external packages), 230 a matrix type can be constructed which invokes Matlab solvers. 231 After calling MatCreate(...,A), simply call MatSetType(A,MATMATLAB). 232 This matrix type is only supported for double precision real. 233 234 This matrix inherits from MATSEQAIJ. As a result, MatSeqAIJSetPreallocation is 235 supported for this matrix type. One can also call MatConvert for an inplace conversion to or from 236 the MATSEQAIJ type without data copy. 237 238 Options Database Keys: 239 + -mat_type matlab - sets the matrix type to "matlab" during a call to MatSetFromOptions() 240 - -mat_matlab_qr - sets the direct solver to be QR instead of LU 241 242 Level: beginner 243 244 .seealso: PCLU 245 M*/ 246 247