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