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