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