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