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