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 PETSC_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 - a 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,PetscScalar,&aij->a,aij->nz,PetscInt,&aij->j,mat->rmap->n+1,PetscInt,&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 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 = PetscMalloc((len+2)*sizeof(char),&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 = PetscMalloc((len+2)*sizeof(char),&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 230 #undef __FUNCT__ 231 #define __FUNCT__ "MatFactorInfo_Matlab" 232 PetscErrorCode MatFactorInfo_Matlab(Mat A,PetscViewer viewer) 233 { 234 PetscErrorCode ierr; 235 236 PetscFunctionBegin; 237 ierr = PetscViewerASCIIPrintf(viewer,"MATLAB run parameters: -- not written yet!\n");CHKERRQ(ierr); 238 PetscFunctionReturn(0); 239 } 240 241 #undef __FUNCT__ 242 #define __FUNCT__ "MatView_Matlab" 243 PetscErrorCode MatView_Matlab(Mat A,PetscViewer viewer) 244 { 245 PetscErrorCode ierr; 246 PetscBool iascii; 247 PetscViewerFormat format; 248 249 PetscFunctionBegin; 250 ierr = MatView_SeqAIJ(A,viewer);CHKERRQ(ierr); 251 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr); 252 if (iascii) { 253 ierr = PetscViewerGetFormat(viewer,&format);CHKERRQ(ierr); 254 if (format == PETSC_VIEWER_ASCII_FACTOR_INFO) { 255 ierr = MatFactorInfo_Matlab(A,viewer); 256 } 257 } 258 PetscFunctionReturn(0); 259 } 260 261 262 /*MC 263 MATSOLVERMATLAB - "matlab" - Providing direct solvers (LU and QR) and drop tolerance 264 based ILU factorization (ILUDT) for sequential matrices via the external package MATLAB. 265 266 267 Works with MATSEQAIJ matrices. 268 269 Options Database Keys: 270 . -pc_factor_mat_solver_type matlab - selects MATLAB to do the sparse factorization 271 272 273 Level: beginner 274 275 .seealso: PCLU 276 277 .seealso: PCFactorSetMatSolverPackage(), MatSolverPackage 278 M*/ 279 280