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