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