1 /* 2 Provides an interface for the Matlab engine sparse solver 3 4 */ 5 #include "src/mat/impls/aij/seq/aij.h" 6 7 #include "engine.h" /* Matlab include file */ 8 #include "mex.h" /* Matlab include file */ 9 10 typedef struct { 11 PetscErrorCode (*MatDuplicate)(Mat,MatDuplicateOption,Mat*); 12 PetscErrorCode (*MatView)(Mat,PetscViewer); 13 PetscErrorCode (*MatLUFactorSymbolic)(Mat,IS,IS,MatFactorInfo*,Mat*); 14 PetscErrorCode (*MatILUDTFactor)(Mat,MatFactorInfo*,IS,IS,Mat*); 15 PetscErrorCode (*MatDestroy)(Mat); 16 } Mat_Matlab; 17 18 19 EXTERN_C_BEGIN 20 #undef __FUNCT__ 21 #define __FUNCT__ "MatMatlabEnginePut_Matlab" 22 PetscErrorCode MatMatlabEnginePut_Matlab(PetscObject obj,void *mengine) 23 { 24 PetscErrorCode ierr; 25 Mat B = (Mat)obj; 26 mxArray *mat; 27 Mat_SeqAIJ *aij = (Mat_SeqAIJ*)B->data; 28 29 PetscFunctionBegin; 30 mat = mxCreateSparse(B->n,B->m,aij->nz,mxREAL); 31 ierr = PetscMemcpy(mxGetPr(mat),aij->a,aij->nz*sizeof(PetscScalar));CHKERRQ(ierr); 32 /* Matlab stores by column, not row so we pass in the transpose of the matrix */ 33 ierr = PetscMemcpy(mxGetIr(mat),aij->j,aij->nz*sizeof(int));CHKERRQ(ierr); 34 ierr = PetscMemcpy(mxGetJc(mat),aij->i,(B->m+1)*sizeof(int));CHKERRQ(ierr); 35 36 /* Matlab indices start at 0 for sparse (what a surprise) */ 37 38 ierr = PetscObjectName(obj);CHKERRQ(ierr); 39 engPutVariable((Engine *)mengine,obj->name,mat); 40 PetscFunctionReturn(0); 41 } 42 EXTERN_C_END 43 44 EXTERN_C_BEGIN 45 #undef __FUNCT__ 46 #define __FUNCT__ "MatMatlabEngineGet_Matlab" 47 PetscErrorCode MatMatlabEngineGet_Matlab(PetscObject obj,void *mengine) 48 { 49 PetscErrorCode ierr; 50 int ii; 51 Mat mat = (Mat)obj; 52 Mat_SeqAIJ *aij = (Mat_SeqAIJ*)mat->data; 53 mxArray *mmat; 54 55 PetscFunctionBegin; 56 ierr = PetscFree(aij->a);CHKERRQ(ierr); 57 58 mmat = engGetVariable((Engine *)mengine,obj->name); 59 60 aij->nz = (mxGetJc(mmat))[mat->m]; 61 ierr = PetscMalloc(((size_t) aij->nz)*(sizeof(int)+sizeof(PetscScalar))+(mat->m+1)*sizeof(int),&aij->a);CHKERRQ(ierr); 62 aij->j = (int*)(aij->a + aij->nz); 63 aij->i = aij->j + aij->nz; 64 aij->singlemalloc = PETSC_TRUE; 65 aij->freedata = PETSC_TRUE; 66 67 ierr = PetscMemcpy(aij->a,mxGetPr(mmat),aij->nz*sizeof(PetscScalar));CHKERRQ(ierr); 68 /* Matlab stores by column, not row so we pass in the transpose of the matrix */ 69 ierr = PetscMemcpy(aij->j,mxGetIr(mmat),aij->nz*sizeof(int));CHKERRQ(ierr); 70 ierr = PetscMemcpy(aij->i,mxGetJc(mmat),(mat->m+1)*sizeof(int));CHKERRQ(ierr); 71 72 for (ii=0; ii<mat->m; ii++) { 73 aij->ilen[ii] = aij->imax[ii] = aij->i[ii+1] - aij->i[ii]; 74 } 75 76 ierr = MatAssemblyBegin(mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 77 ierr = MatAssemblyEnd(mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 78 79 PetscFunctionReturn(0); 80 } 81 EXTERN_C_END 82 83 EXTERN_C_BEGIN 84 #undef __FUNCT__ 85 #define __FUNCT__ "MatConvert_Matlab_SeqAIJ" 86 PetscErrorCode MatConvert_Matlab_SeqAIJ(Mat A,const MatType type,Mat *newmat) { 87 PetscErrorCode ierr; 88 Mat B=*newmat; 89 Mat_Matlab *lu=(Mat_Matlab*)A->spptr; 90 91 PetscFunctionBegin; 92 if (B != A) { 93 ierr = MatDuplicate(A,MAT_COPY_VALUES,&B);CHKERRQ(ierr); 94 } 95 A->ops->duplicate = lu->MatDuplicate; 96 A->ops->view = lu->MatView; 97 A->ops->lufactorsymbolic = lu->MatLUFactorSymbolic; 98 A->ops->iludtfactor = lu->MatILUDTFactor; 99 A->ops->destroy = lu->MatDestroy; 100 101 ierr = PetscFree(lu);CHKERRQ(ierr); 102 ierr = PetscObjectChangeTypeName((PetscObject)B,MATSEQAIJ);CHKERRQ(ierr); 103 *newmat = B; 104 PetscFunctionReturn(0); 105 } 106 EXTERN_C_END 107 108 #undef __FUNCT__ 109 #define __FUNCT__ "MatDestroy_Matlab" 110 PetscErrorCode MatDestroy_Matlab(Mat A) 111 { 112 PetscErrorCode ierr; 113 114 PetscFunctionBegin; 115 ierr = MatConvert_Matlab_SeqAIJ(A,MATSEQAIJ,&A);CHKERRQ(ierr); 116 ierr = (*A->ops->destroy)(A);CHKERRQ(ierr); 117 PetscFunctionReturn(0); 118 } 119 120 #undef __FUNCT__ 121 #define __FUNCT__ "MatSolve_Matlab" 122 PetscErrorCode MatSolve_Matlab(Mat A,Vec b,Vec x) 123 { 124 PetscErrorCode ierr; 125 char *_A,*_b,*_x; 126 127 PetscFunctionBegin; 128 /* make sure objects have names; use default if not */ 129 ierr = PetscObjectName((PetscObject)b);CHKERRQ(ierr); 130 ierr = PetscObjectName((PetscObject)x);CHKERRQ(ierr); 131 132 ierr = PetscObjectGetName((PetscObject)A,&_A);CHKERRQ(ierr); 133 ierr = PetscObjectGetName((PetscObject)b,&_b);CHKERRQ(ierr); 134 ierr = PetscObjectGetName((PetscObject)x,&_x);CHKERRQ(ierr); 135 ierr = PetscMatlabEnginePut(PETSC_MATLAB_ENGINE_(A->comm),(PetscObject)b);CHKERRQ(ierr); 136 ierr = PetscMatlabEngineEvaluate(PETSC_MATLAB_ENGINE_(A->comm),"%s = u%s\\(l%s\\(p%s*%s));",_x,_A,_A,_A,_b);CHKERRQ(ierr); 137 ierr = PetscMatlabEngineEvaluate(PETSC_MATLAB_ENGINE_(A->comm),"%s = 0;",_b);CHKERRQ(ierr); 138 /* ierr = PetscMatlabEnginePrintOutput(PETSC_MATLAB_ENGINE_(A->comm),stdout);CHKERRQ(ierr); */ 139 ierr = PetscMatlabEngineGet(PETSC_MATLAB_ENGINE_(A->comm),(PetscObject)x);CHKERRQ(ierr); 140 PetscFunctionReturn(0); 141 } 142 143 #undef __FUNCT__ 144 #define __FUNCT__ "MatLUFactorNumeric_Matlab" 145 PetscErrorCode MatLUFactorNumeric_Matlab(Mat A,Mat *F) 146 { 147 Mat_SeqAIJ *f = (Mat_SeqAIJ*)(*F)->data; 148 PetscErrorCode ierr; 149 size_t len; 150 char *_A,*name; 151 152 PetscFunctionBegin; 153 ierr = PetscMatlabEnginePut(PETSC_MATLAB_ENGINE_(A->comm),(PetscObject)A);CHKERRQ(ierr); 154 _A = A->name; 155 ierr = PetscMatlabEngineEvaluate(PETSC_MATLAB_ENGINE_(A->comm),"[l_%s,u_%s,p_%s] = lu(%s',%g);",_A,_A,_A,_A,f->lu_dtcol);CHKERRQ(ierr); 156 ierr = PetscMatlabEngineEvaluate(PETSC_MATLAB_ENGINE_(A->comm),"%s = 0;",_A);CHKERRQ(ierr); 157 ierr = PetscStrlen(_A,&len);CHKERRQ(ierr); 158 ierr = PetscMalloc((len+2)*sizeof(char),&name);CHKERRQ(ierr); 159 sprintf(name,"_%s",_A); 160 ierr = PetscObjectSetName((PetscObject)*F,name);CHKERRQ(ierr); 161 ierr = PetscFree(name);CHKERRQ(ierr); 162 PetscFunctionReturn(0); 163 } 164 165 #undef __FUNCT__ 166 #define __FUNCT__ "MatLUFactorSymbolic_Matlab" 167 PetscErrorCode MatLUFactorSymbolic_Matlab(Mat A,IS r,IS c,MatFactorInfo *info,Mat *F) 168 { 169 PetscErrorCode ierr; 170 Mat_SeqAIJ *f; 171 172 PetscFunctionBegin; 173 if (A->N != A->M) SETERRQ(PETSC_ERR_ARG_SIZ,"matrix must be square"); 174 ierr = MatCreate(A->comm,A->m,A->n,A->m,A->n,F);CHKERRQ(ierr); 175 ierr = MatSetType(*F,A->type_name);CHKERRQ(ierr); 176 ierr = MatSeqAIJSetPreallocation(*F,0,PETSC_NULL);CHKERRQ(ierr); 177 (*F)->ops->solve = MatSolve_Matlab; 178 (*F)->ops->lufactornumeric = MatLUFactorNumeric_Matlab; 179 (*F)->factor = FACTOR_LU; 180 f = (Mat_SeqAIJ*)(*F)->data; 181 f->lu_dtcol = info->dtcol; 182 PetscFunctionReturn(0); 183 } 184 185 /* ---------------------------------------------------------------------------------*/ 186 #undef __FUNCT__ 187 #define __FUNCT__ "MatSolve_Matlab_QR" 188 PetscErrorCode MatSolve_Matlab_QR(Mat A,Vec b,Vec x) 189 { 190 PetscErrorCode ierr; 191 char *_A,*_b,*_x; 192 193 PetscFunctionBegin; 194 /* make sure objects have names; use default if not */ 195 ierr = PetscObjectName((PetscObject)b);CHKERRQ(ierr); 196 ierr = PetscObjectName((PetscObject)x);CHKERRQ(ierr); 197 198 ierr = PetscObjectGetName((PetscObject)A,&_A);CHKERRQ(ierr); 199 ierr = PetscObjectGetName((PetscObject)b,&_b);CHKERRQ(ierr); 200 ierr = PetscObjectGetName((PetscObject)x,&_x);CHKERRQ(ierr); 201 ierr = PetscMatlabEnginePut(PETSC_MATLAB_ENGINE_(A->comm),(PetscObject)b);CHKERRQ(ierr); 202 ierr = PetscMatlabEngineEvaluate(PETSC_MATLAB_ENGINE_(A->comm),"%s = r%s\\(r%s'\\(%s*%s));",_x,_A,_A,_A+1,_b);CHKERRQ(ierr); 203 ierr = PetscMatlabEngineEvaluate(PETSC_MATLAB_ENGINE_(A->comm),"%s = 0;",_b);CHKERRQ(ierr); 204 /* ierr = PetscMatlabEnginePrintOutput(PETSC_MATLAB_ENGINE_(A->comm),stdout);CHKERRQ(ierr); */ 205 ierr = PetscMatlabEngineGet(PETSC_MATLAB_ENGINE_(A->comm),(PetscObject)x);CHKERRQ(ierr); 206 PetscFunctionReturn(0); 207 } 208 209 #undef __FUNCT__ 210 #define __FUNCT__ "MatLUFactorNumeric_Matlab_QR" 211 PetscErrorCode MatLUFactorNumeric_Matlab_QR(Mat A,Mat *F) 212 { 213 PetscErrorCode ierr; 214 size_t len; 215 char *_A,*name; 216 217 PetscFunctionBegin; 218 ierr = PetscMatlabEnginePut(PETSC_MATLAB_ENGINE_(A->comm),(PetscObject)A);CHKERRQ(ierr); 219 _A = A->name; 220 ierr = PetscMatlabEngineEvaluate(PETSC_MATLAB_ENGINE_(A->comm),"r_%s = qr(%s');",_A,_A);CHKERRQ(ierr); 221 ierr = PetscStrlen(_A,&len);CHKERRQ(ierr); 222 ierr = PetscMalloc((len+2)*sizeof(char),&name);CHKERRQ(ierr); 223 sprintf(name,"_%s",_A); 224 ierr = PetscObjectSetName((PetscObject)*F,name);CHKERRQ(ierr); 225 ierr = PetscFree(name);CHKERRQ(ierr); 226 PetscFunctionReturn(0); 227 } 228 229 #undef __FUNCT__ 230 #define __FUNCT__ "MatLUFactorSymbolic_Matlab_QR" 231 PetscErrorCode MatLUFactorSymbolic_Matlab_QR(Mat A,IS r,IS c,MatFactorInfo *info,Mat *F) 232 { 233 PetscErrorCode ierr; 234 235 PetscFunctionBegin; 236 if (A->N != A->M) SETERRQ(PETSC_ERR_ARG_SIZ,"matrix must be square"); 237 ierr = MatCreate(A->comm,A->m,A->n,A->m,A->n,F);CHKERRQ(ierr); 238 ierr = MatSetType(*F,A->type_name);CHKERRQ(ierr); 239 ierr = MatSeqAIJSetPreallocation(*F,0,PETSC_NULL);CHKERRQ(ierr); 240 (*F)->ops->solve = MatSolve_Matlab_QR; 241 (*F)->ops->lufactornumeric = MatLUFactorNumeric_Matlab_QR; 242 (*F)->factor = FACTOR_LU; 243 (*F)->assembled = PETSC_TRUE; /* required by -ksp_view */ 244 245 PetscFunctionReturn(0); 246 } 247 248 /* --------------------------------------------------------------------------------*/ 249 #undef __FUNCT__ 250 #define __FUNCT__ "MatILUDTFactor_Matlab" 251 PetscErrorCode MatILUDTFactor_Matlab(Mat A,MatFactorInfo *info,IS isrow,IS iscol,Mat *F) 252 { 253 PetscErrorCode ierr; 254 size_t len; 255 char *_A,*name; 256 257 PetscFunctionBegin; 258 if (info->dt == PETSC_DEFAULT) info->dt = .005; 259 if (info->dtcol == PETSC_DEFAULT) info->dtcol = .01; 260 if (A->N != A->M) SETERRQ(PETSC_ERR_ARG_SIZ,"matrix must be square"); 261 ierr = MatCreate(A->comm,A->m,A->n,A->m,A->n,F);CHKERRQ(ierr); 262 ierr = MatSetType(*F,A->type_name);CHKERRQ(ierr); 263 ierr = MatSeqAIJSetPreallocation(*F,0,PETSC_NULL);CHKERRQ(ierr); 264 (*F)->ops->solve = MatSolve_Matlab; 265 (*F)->factor = FACTOR_LU; 266 ierr = PetscMatlabEnginePut(PETSC_MATLAB_ENGINE_(A->comm),(PetscObject)A);CHKERRQ(ierr); 267 _A = A->name; 268 ierr = PetscMatlabEngineEvaluate(PETSC_MATLAB_ENGINE_(A->comm),"info_%s = struct('droptol',%g,'thresh',%g);",_A,info->dt,info->dtcol);CHKERRQ(ierr); 269 ierr = PetscMatlabEngineEvaluate(PETSC_MATLAB_ENGINE_(A->comm),"[l_%s,u_%s,p_%s] = luinc(%s',info_%s);",_A,_A,_A,_A,_A);CHKERRQ(ierr); 270 ierr = PetscMatlabEngineEvaluate(PETSC_MATLAB_ENGINE_(A->comm),"%s = 0;",_A);CHKERRQ(ierr); 271 272 ierr = PetscStrlen(_A,&len);CHKERRQ(ierr); 273 ierr = PetscMalloc((len+2)*sizeof(char),&name);CHKERRQ(ierr); 274 sprintf(name,"_%s",_A); 275 ierr = PetscObjectSetName((PetscObject)*F,name);CHKERRQ(ierr); 276 ierr = PetscFree(name);CHKERRQ(ierr); 277 PetscFunctionReturn(0); 278 } 279 280 #undef __FUNCT__ 281 #define __FUNCT__ "MatFactorInfo_Matlab" 282 PetscErrorCode MatFactorInfo_Matlab(Mat A,PetscViewer viewer) 283 { 284 PetscErrorCode ierr; 285 286 PetscFunctionBegin; 287 ierr = PetscViewerASCIIPrintf(viewer,"Matlab run parameters: -- not written yet!\n");CHKERRQ(ierr); 288 PetscFunctionReturn(0); 289 } 290 291 #undef __FUNCT__ 292 #define __FUNCT__ "MatView_Matlab" 293 PetscErrorCode MatView_Matlab(Mat A,PetscViewer viewer) { 294 PetscErrorCode ierr; 295 PetscTruth iascii; 296 PetscViewerFormat format; 297 Mat_Matlab *lu=(Mat_Matlab*)(A->spptr); 298 299 PetscFunctionBegin; 300 ierr = (*lu->MatView)(A,viewer);CHKERRQ(ierr); 301 ierr = PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_ASCII,&iascii);CHKERRQ(ierr); 302 if (iascii) { 303 ierr = PetscViewerGetFormat(viewer,&format);CHKERRQ(ierr); 304 if (format == PETSC_VIEWER_ASCII_FACTOR_INFO) { 305 ierr = MatFactorInfo_Matlab(A,viewer); 306 } 307 } 308 PetscFunctionReturn(0); 309 } 310 311 #undef __FUNCT__ 312 #define __FUNCT__ "MatDuplicate_Matlab" 313 PetscErrorCode MatDuplicate_Matlab(Mat A, MatDuplicateOption op, Mat *M) { 314 PetscErrorCode ierr; 315 Mat_Matlab *lu=(Mat_Matlab*)A->spptr; 316 317 PetscFunctionBegin; 318 ierr = (*lu->MatDuplicate)(A,op,M);CHKERRQ(ierr); 319 ierr = PetscMemcpy((*M)->spptr,lu,sizeof(Mat_Matlab));CHKERRQ(ierr); 320 PetscFunctionReturn(0); 321 } 322 323 EXTERN_C_BEGIN 324 #undef __FUNCT__ 325 #define __FUNCT__ "MatConvert_SeqAIJ_Matlab" 326 PetscErrorCode MatConvert_SeqAIJ_Matlab(Mat A,const MatType type,Mat *newmat) { 327 /* This routine is only called to convert to MATMATLAB */ 328 /* from MATSEQAIJ, so we will ignore 'MatType type'. */ 329 PetscErrorCode ierr; 330 Mat B=*newmat; 331 Mat_Matlab *lu; 332 PetscTruth qr; 333 334 PetscFunctionBegin; 335 if (B != A) { 336 ierr = MatDuplicate(A,MAT_COPY_VALUES,&B);CHKERRQ(ierr); 337 } 338 339 ierr = PetscNew(Mat_Matlab,&lu);CHKERRQ(ierr); 340 lu->MatDuplicate = A->ops->duplicate; 341 lu->MatView = A->ops->view; 342 lu->MatLUFactorSymbolic = A->ops->lufactorsymbolic; 343 lu->MatILUDTFactor = A->ops->iludtfactor; 344 lu->MatDestroy = A->ops->destroy; 345 346 B->spptr = (void*)lu; 347 B->ops->duplicate = MatDuplicate_Matlab; 348 B->ops->view = MatView_Matlab; 349 B->ops->lufactorsymbolic = MatLUFactorSymbolic_Matlab; 350 B->ops->iludtfactor = MatILUDTFactor_Matlab; 351 B->ops->destroy = MatDestroy_Matlab; 352 353 ierr = PetscOptionsHasName(A->prefix,"-mat_matlab_qr",&qr);CHKERRQ(ierr); 354 if (qr) { 355 B->ops->lufactorsymbolic = MatLUFactorSymbolic_Matlab_QR; 356 PetscLogInfo(0,"Using Matlab QR with iterative refinement for LU factorization and solves"); 357 } else { 358 PetscLogInfo(0,"Using Matlab for LU factorizations and solves."); 359 } 360 PetscLogInfo(0,"Using Matlab for ILUDT factorizations and solves."); 361 362 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatConvert_seqaij_matlab_C", 363 "MatConvert_SeqAIJ_Matlab",MatConvert_SeqAIJ_Matlab);CHKERRQ(ierr); 364 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatConvert_matlab_seqaij_C", 365 "MatConvert_Matlab_SeqAIJ",MatConvert_Matlab_SeqAIJ);CHKERRQ(ierr); 366 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"PetscMatlabEnginePut_C", 367 "MatMatlabEnginePut_Matlab",MatMatlabEnginePut_Matlab);CHKERRQ(ierr); 368 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"PetscMatlabEngineGet_C", 369 "MatMatlabEngineGet_Matlab",MatMatlabEngineGet_Matlab);CHKERRQ(ierr); 370 ierr = PetscObjectChangeTypeName((PetscObject)B,MATMATLAB);CHKERRQ(ierr); 371 *newmat = B; 372 PetscFunctionReturn(0); 373 } 374 EXTERN_C_END 375 376 /*MC 377 MATMATLAB - MATMATLAB = "matlab" - A matrix type providing direct solvers (LU and QR) and drop tolerance 378 based ILU factorization (ILUDT) for sequential matrices via the external package Matlab. 379 380 If Matlab is instaled (see the manual for 381 instructions on how to declare the existence of external packages), 382 a matrix type can be constructed which invokes Matlab solvers. 383 After calling MatCreate(...,A), simply call MatSetType(A,MATMATLAB). 384 This matrix type is only supported for double precision real. 385 386 This matrix inherits from MATSEQAIJ. As a result, MatSeqAIJSetPreallocation is 387 supported for this matrix type. One can also call MatConvert for an inplace conversion to or from 388 the MATSEQAIJ type without data copy. 389 390 Options Database Keys: 391 + -mat_type matlab - sets the matrix type to "matlab" during a call to MatSetFromOptions() 392 - -mat_matlab_qr - sets the direct solver to be QR instead of LU 393 394 Level: beginner 395 396 .seealso: PCLU 397 M*/ 398 399 EXTERN_C_BEGIN 400 #undef __FUNCT__ 401 #define __FUNCT__ "MatCreate_Matlab" 402 PetscErrorCode MatCreate_Matlab(Mat A) 403 { 404 PetscErrorCode ierr; 405 406 PetscFunctionBegin; 407 ierr = PetscObjectChangeTypeName((PetscObject)A,MATMATLAB);CHKERRQ(ierr); 408 ierr = MatSetType(A,MATSEQAIJ);CHKERRQ(ierr); 409 ierr = MatConvert_SeqAIJ_Matlab(A,MATMATLAB,&A);CHKERRQ(ierr); 410 PetscFunctionReturn(0); 411 } 412 EXTERN_C_END 413