1 #define PETSCMAT_DLL 2 3 /* 4 Provides an interface to the UMFPACKv4.3 sparse solver 5 */ 6 7 #include "src/mat/impls/aij/seq/aij.h" 8 9 EXTERN_C_BEGIN 10 #include "umfpack.h" 11 EXTERN_C_END 12 13 typedef struct { 14 void *Symbolic, *Numeric; 15 double Info[UMFPACK_INFO], Control[UMFPACK_CONTROL],*W; 16 int *Wi,*ai,*aj,*perm_c; 17 PetscScalar *av; 18 MatStructure flg; 19 PetscTruth PetscMatOdering; 20 21 /* A few function pointers for inheritance */ 22 PetscErrorCode (*MatDuplicate)(Mat,MatDuplicateOption,Mat*); 23 PetscErrorCode (*MatView)(Mat,PetscViewer); 24 PetscErrorCode (*MatAssemblyEnd)(Mat,MatAssemblyType); 25 PetscErrorCode (*MatLUFactorSymbolic)(Mat,IS,IS,MatFactorInfo*,Mat*); 26 PetscErrorCode (*MatDestroy)(Mat); 27 28 /* Flag to clean up UMFPACK objects during Destroy */ 29 PetscTruth CleanUpUMFPACK; 30 } Mat_UMFPACK; 31 32 EXTERN PetscErrorCode MatDuplicate_UMFPACK(Mat,MatDuplicateOption,Mat*); 33 34 EXTERN_C_BEGIN 35 #undef __FUNCT__ 36 #define __FUNCT__ "MatConvert_UMFPACK_SeqAIJ" 37 PetscErrorCode PETSCMAT_DLLEXPORT MatConvert_UMFPACK_SeqAIJ(Mat A,MatType type,MatReuse reuse,Mat *newmat) 38 { 39 PetscErrorCode ierr; 40 Mat B=*newmat; 41 Mat_UMFPACK *lu=(Mat_UMFPACK*)A->spptr; 42 43 PetscFunctionBegin; 44 if (reuse == MAT_INITIAL_MATRIX) { 45 ierr = MatDuplicate(A,MAT_COPY_VALUES,&B);CHKERRQ(ierr); 46 } 47 /* Reset the original function pointers */ 48 B->ops->duplicate = lu->MatDuplicate; 49 B->ops->view = lu->MatView; 50 B->ops->assemblyend = lu->MatAssemblyEnd; 51 B->ops->lufactorsymbolic = lu->MatLUFactorSymbolic; 52 B->ops->destroy = lu->MatDestroy; 53 54 ierr = PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqaij_umfpack_C","",PETSC_NULL);CHKERRQ(ierr); 55 ierr = PetscObjectComposeFunction((PetscObject)B,"MatConvert_umfpack_seqaij_C","",PETSC_NULL);CHKERRQ(ierr); 56 57 ierr = PetscObjectChangeTypeName((PetscObject)B,MATSEQAIJ);CHKERRQ(ierr); 58 *newmat = B; 59 PetscFunctionReturn(0); 60 } 61 EXTERN_C_END 62 63 #undef __FUNCT__ 64 #define __FUNCT__ "MatDestroy_UMFPACK" 65 PetscErrorCode MatDestroy_UMFPACK(Mat A) 66 { 67 PetscErrorCode ierr; 68 Mat_UMFPACK *lu=(Mat_UMFPACK*)A->spptr; 69 70 PetscFunctionBegin; 71 if (lu->CleanUpUMFPACK) { 72 umfpack_di_free_symbolic(&lu->Symbolic); 73 umfpack_di_free_numeric(&lu->Numeric); 74 ierr = PetscFree(lu->Wi);CHKERRQ(ierr); 75 ierr = PetscFree(lu->W);CHKERRQ(ierr); 76 if (lu->PetscMatOdering) { 77 ierr = PetscFree(lu->perm_c);CHKERRQ(ierr); 78 } 79 } 80 ierr = MatConvert_UMFPACK_SeqAIJ(A,MATSEQAIJ,MAT_REUSE_MATRIX,&A);CHKERRQ(ierr); 81 ierr = (*A->ops->destroy)(A);CHKERRQ(ierr); 82 PetscFunctionReturn(0); 83 } 84 85 #undef __FUNCT__ 86 #define __FUNCT__ "MatSolve_UMFPACK" 87 PetscErrorCode MatSolve_UMFPACK(Mat A,Vec b,Vec x) 88 { 89 Mat_UMFPACK *lu = (Mat_UMFPACK*)A->spptr; 90 PetscScalar *av=lu->av,*ba,*xa; 91 PetscErrorCode ierr; 92 int *ai=lu->ai,*aj=lu->aj,status; 93 94 PetscFunctionBegin; 95 /* solve Ax = b by umfpack_di_wsolve */ 96 /* ----------------------------------*/ 97 ierr = VecGetArray(b,&ba); 98 ierr = VecGetArray(x,&xa); 99 100 status = umfpack_di_wsolve(UMFPACK_At,ai,aj,av,xa,ba,lu->Numeric,lu->Control,lu->Info,lu->Wi,lu->W); 101 umfpack_di_report_info(lu->Control, lu->Info); 102 if (status < 0){ 103 umfpack_di_report_status(lu->Control, status); 104 SETERRQ(PETSC_ERR_LIB,"umfpack_di_wsolve failed"); 105 } 106 107 ierr = VecRestoreArray(b,&ba); 108 ierr = VecRestoreArray(x,&xa); 109 110 PetscFunctionReturn(0); 111 } 112 113 #undef __FUNCT__ 114 #define __FUNCT__ "MatLUFactorNumeric_UMFPACK" 115 PetscErrorCode MatLUFactorNumeric_UMFPACK(Mat A,MatFactorInfo *info,Mat *F) 116 { 117 Mat_UMFPACK *lu=(Mat_UMFPACK*)(*F)->spptr; 118 PetscErrorCode ierr; 119 int *ai=lu->ai,*aj=lu->aj,m=A->rmap.n,status; 120 PetscScalar *av=lu->av; 121 122 PetscFunctionBegin; 123 /* numeric factorization of A' */ 124 /* ----------------------------*/ 125 if (lu->flg == SAME_NONZERO_PATTERN && lu->Numeric){ 126 umfpack_di_free_numeric(&lu->Numeric); 127 } 128 status = umfpack_di_numeric (ai,aj,av,lu->Symbolic,&lu->Numeric,lu->Control,lu->Info); 129 if (status < 0) SETERRQ(PETSC_ERR_LIB,"umfpack_di_numeric failed"); 130 /* report numeric factorization of A' when Control[PRL] > 3 */ 131 (void) umfpack_di_report_numeric (lu->Numeric, lu->Control); 132 133 if (lu->flg == DIFFERENT_NONZERO_PATTERN){ /* first numeric factorization */ 134 /* allocate working space to be used by Solve */ 135 ierr = PetscMalloc(m * sizeof(int), &lu->Wi);CHKERRQ(ierr); 136 ierr = PetscMalloc(5*m * sizeof(double), &lu->W);CHKERRQ(ierr); 137 } 138 lu->flg = SAME_NONZERO_PATTERN; 139 lu->CleanUpUMFPACK = PETSC_TRUE; 140 PetscFunctionReturn(0); 141 } 142 143 /* 144 Note the r permutation is ignored 145 */ 146 #undef __FUNCT__ 147 #define __FUNCT__ "MatLUFactorSymbolic_UMFPACK" 148 PetscErrorCode MatLUFactorSymbolic_UMFPACK(Mat A,IS r,IS c,MatFactorInfo *info,Mat *F) 149 { 150 Mat B; 151 Mat_SeqAIJ *mat=(Mat_SeqAIJ*)A->data; 152 Mat_UMFPACK *lu; 153 PetscErrorCode ierr; 154 int m=A->rmap.n,n=A->cmap.n,*ai=mat->i,*aj=mat->j,status,*ra,idx; 155 PetscScalar *av=mat->a; 156 const char *strategy[]={"AUTO","UNSYMMETRIC","SYMMETRIC","2BY2"}, 157 *scale[]={"NONE","SUM","MAX"}; 158 PetscTruth flg; 159 160 PetscFunctionBegin; 161 /* Create the factorization matrix F */ 162 ierr = MatCreate(A->comm,&B);CHKERRQ(ierr); 163 ierr = MatSetSizes(B,PETSC_DECIDE,PETSC_DECIDE,m,n);CHKERRQ(ierr); 164 ierr = MatSetType(B,A->type_name);CHKERRQ(ierr); 165 ierr = MatSeqAIJSetPreallocation(B,0,PETSC_NULL);CHKERRQ(ierr); 166 167 B->ops->lufactornumeric = MatLUFactorNumeric_UMFPACK; 168 B->ops->solve = MatSolve_UMFPACK; 169 B->factor = FACTOR_LU; 170 B->assembled = PETSC_TRUE; /* required by -ksp_view */ 171 172 lu = (Mat_UMFPACK*)(B->spptr); 173 174 /* initializations */ 175 /* ------------------------------------------------*/ 176 /* get the default control parameters */ 177 umfpack_di_defaults (lu->Control); 178 lu->perm_c = PETSC_NULL; /* use defaul UMFPACK col permutation */ 179 lu->Control[UMFPACK_IRSTEP] = 0; /* max num of iterative refinement steps to attempt */ 180 181 ierr = PetscOptionsBegin(A->comm,A->prefix,"UMFPACK Options","Mat");CHKERRQ(ierr); 182 /* Control parameters used by reporting routiones */ 183 ierr = PetscOptionsReal("-mat_umfpack_prl","Control[UMFPACK_PRL]","None",lu->Control[UMFPACK_PRL],&lu->Control[UMFPACK_PRL],PETSC_NULL);CHKERRQ(ierr); 184 185 /* Control parameters for symbolic factorization */ 186 ierr = PetscOptionsEList("-mat_umfpack_strategy","ordering and pivoting strategy","None",strategy,4,strategy[0],&idx,&flg);CHKERRQ(ierr); 187 if (flg) { 188 switch (idx){ 189 case 0: lu->Control[UMFPACK_STRATEGY] = UMFPACK_STRATEGY_AUTO; break; 190 case 1: lu->Control[UMFPACK_STRATEGY] = UMFPACK_STRATEGY_UNSYMMETRIC; break; 191 case 2: lu->Control[UMFPACK_STRATEGY] = UMFPACK_STRATEGY_SYMMETRIC; break; 192 case 3: lu->Control[UMFPACK_STRATEGY] = UMFPACK_STRATEGY_2BY2; break; 193 } 194 } 195 ierr = PetscOptionsReal("-mat_umfpack_dense_col","Control[UMFPACK_DENSE_COL]","None",lu->Control[UMFPACK_DENSE_COL],&lu->Control[UMFPACK_DENSE_COL],PETSC_NULL);CHKERRQ(ierr); 196 ierr = PetscOptionsReal("-mat_umfpack_dense_row","Control[UMFPACK_DENSE_ROW]","None",lu->Control[UMFPACK_DENSE_ROW],&lu->Control[UMFPACK_DENSE_ROW],PETSC_NULL);CHKERRQ(ierr); 197 ierr = PetscOptionsReal("-mat_umfpack_amd_dense","Control[UMFPACK_AMD_DENSE]","None",lu->Control[UMFPACK_AMD_DENSE],&lu->Control[UMFPACK_AMD_DENSE],PETSC_NULL);CHKERRQ(ierr); 198 ierr = PetscOptionsReal("-mat_umfpack_block_size","Control[UMFPACK_BLOCK_SIZE]","None",lu->Control[UMFPACK_BLOCK_SIZE],&lu->Control[UMFPACK_BLOCK_SIZE],PETSC_NULL);CHKERRQ(ierr); 199 ierr = PetscOptionsReal("-mat_umfpack_2by2_tolerance","Control[UMFPACK_2BY2_TOLERANCE]","None",lu->Control[UMFPACK_2BY2_TOLERANCE],&lu->Control[UMFPACK_2BY2_TOLERANCE],PETSC_NULL);CHKERRQ(ierr); 200 ierr = PetscOptionsReal("-mat_umfpack_fixq","Control[UMFPACK_FIXQ]","None",lu->Control[UMFPACK_FIXQ],&lu->Control[UMFPACK_FIXQ],PETSC_NULL);CHKERRQ(ierr); 201 ierr = PetscOptionsReal("-mat_umfpack_aggressive","Control[UMFPACK_AGGRESSIVE]","None",lu->Control[UMFPACK_AGGRESSIVE],&lu->Control[UMFPACK_AGGRESSIVE],PETSC_NULL);CHKERRQ(ierr); 202 203 /* Control parameters used by numeric factorization */ 204 ierr = PetscOptionsReal("-mat_umfpack_pivot_tolerance","Control[UMFPACK_PIVOT_TOLERANCE]","None",lu->Control[UMFPACK_PIVOT_TOLERANCE],&lu->Control[UMFPACK_PIVOT_TOLERANCE],PETSC_NULL);CHKERRQ(ierr); 205 ierr = PetscOptionsReal("-mat_umfpack_sym_pivot_tolerance","Control[UMFPACK_SYM_PIVOT_TOLERANCE]","None",lu->Control[UMFPACK_SYM_PIVOT_TOLERANCE],&lu->Control[UMFPACK_SYM_PIVOT_TOLERANCE],PETSC_NULL);CHKERRQ(ierr); 206 ierr = PetscOptionsEList("-mat_umfpack_scale","Control[UMFPACK_SCALE]","None",scale,3,scale[0],&idx,&flg);CHKERRQ(ierr); 207 if (flg) { 208 switch (idx){ 209 case 0: lu->Control[UMFPACK_SCALE] = UMFPACK_SCALE_NONE; break; 210 case 1: lu->Control[UMFPACK_SCALE] = UMFPACK_SCALE_SUM; break; 211 case 2: lu->Control[UMFPACK_SCALE] = UMFPACK_SCALE_MAX; break; 212 } 213 } 214 ierr = PetscOptionsReal("-mat_umfpack_alloc_init","Control[UMFPACK_ALLOC_INIT]","None",lu->Control[UMFPACK_ALLOC_INIT],&lu->Control[UMFPACK_ALLOC_INIT],PETSC_NULL);CHKERRQ(ierr); 215 ierr = PetscOptionsReal("-mat_umfpack_front_alloc_init","Control[UMFPACK_FRONT_ALLOC_INIT]","None",lu->Control[UMFPACK_FRONT_ALLOC_INIT],&lu->Control[UMFPACK_ALLOC_INIT],PETSC_NULL);CHKERRQ(ierr); 216 ierr = PetscOptionsReal("-mat_umfpack_droptol","Control[UMFPACK_DROPTOL]","None",lu->Control[UMFPACK_DROPTOL],&lu->Control[UMFPACK_DROPTOL],PETSC_NULL);CHKERRQ(ierr); 217 218 /* Control parameters used by solve */ 219 ierr = PetscOptionsReal("-mat_umfpack_irstep","Control[UMFPACK_IRSTEP]","None",lu->Control[UMFPACK_IRSTEP],&lu->Control[UMFPACK_IRSTEP],PETSC_NULL);CHKERRQ(ierr); 220 221 /* use Petsc mat ordering (note: size is for the transpose, and PETSc r = Umfpack perm_c) */ 222 ierr = PetscOptionsHasName(PETSC_NULL,"-pc_factor_mat_ordering_type",&lu->PetscMatOdering);CHKERRQ(ierr); 223 if (lu->PetscMatOdering) { 224 ierr = ISGetIndices(r,&ra);CHKERRQ(ierr); 225 ierr = PetscMalloc(A->rmap.n*sizeof(int),&lu->perm_c);CHKERRQ(ierr); 226 ierr = PetscMemcpy(lu->perm_c,ra,A->rmap.n*sizeof(int));CHKERRQ(ierr); 227 ierr = ISRestoreIndices(r,&ra);CHKERRQ(ierr); 228 } 229 PetscOptionsEnd(); 230 231 /* print the control parameters */ 232 if( lu->Control[UMFPACK_PRL] > 1 ) umfpack_di_report_control (lu->Control); 233 234 /* symbolic factorization of A' */ 235 /* ---------------------------------------------------------------------- */ 236 if (lu->PetscMatOdering) { /* use Petsc row ordering */ 237 status = umfpack_di_qsymbolic(n,m,ai,aj,av,lu->perm_c,&lu->Symbolic,lu->Control,lu->Info); 238 } else { /* use Umfpack col ordering */ 239 status = umfpack_di_symbolic(n,m,ai,aj,av,&lu->Symbolic,lu->Control,lu->Info); 240 } 241 if (status < 0){ 242 umfpack_di_report_info(lu->Control, lu->Info); 243 umfpack_di_report_status(lu->Control, status); 244 SETERRQ(PETSC_ERR_LIB,"umfpack_di_symbolic failed"); 245 } 246 /* report sumbolic factorization of A' when Control[PRL] > 3 */ 247 (void) umfpack_di_report_symbolic(lu->Symbolic, lu->Control); 248 249 lu->flg = DIFFERENT_NONZERO_PATTERN; 250 lu->ai = ai; 251 lu->aj = aj; 252 lu->av = av; 253 lu->CleanUpUMFPACK = PETSC_TRUE; 254 *F = B; 255 PetscFunctionReturn(0); 256 } 257 258 #undef __FUNCT__ 259 #define __FUNCT__ "MatAssemblyEnd_UMFPACK" 260 PetscErrorCode MatAssemblyEnd_UMFPACK(Mat A,MatAssemblyType mode) 261 { 262 PetscErrorCode ierr; 263 Mat_UMFPACK *lu=(Mat_UMFPACK*)(A->spptr); 264 265 PetscFunctionBegin; 266 ierr = (*lu->MatAssemblyEnd)(A,mode);CHKERRQ(ierr); 267 lu->MatLUFactorSymbolic = A->ops->lufactorsymbolic; 268 A->ops->lufactorsymbolic = MatLUFactorSymbolic_UMFPACK; 269 PetscFunctionReturn(0); 270 } 271 272 /* used by -ksp_view */ 273 #undef __FUNCT__ 274 #define __FUNCT__ "MatFactorInfo_UMFPACK" 275 PetscErrorCode MatFactorInfo_UMFPACK(Mat A,PetscViewer viewer) 276 { 277 Mat_UMFPACK *lu= (Mat_UMFPACK*)A->spptr; 278 PetscErrorCode ierr; 279 280 PetscFunctionBegin; 281 /* check if matrix is UMFPACK type */ 282 if (A->ops->solve != MatSolve_UMFPACK) PetscFunctionReturn(0); 283 284 ierr = PetscViewerASCIIPrintf(viewer,"UMFPACK run parameters:\n");CHKERRQ(ierr); 285 /* Control parameters used by reporting routiones */ 286 ierr = PetscViewerASCIIPrintf(viewer," Control[UMFPACK_PRL]: %g\n",lu->Control[UMFPACK_PRL]);CHKERRQ(ierr); 287 288 /* Control parameters used by symbolic factorization */ 289 ierr = PetscViewerASCIIPrintf(viewer," Control[UMFPACK_STRATEGY]: %g\n",lu->Control[UMFPACK_STRATEGY]);CHKERRQ(ierr); 290 ierr = PetscViewerASCIIPrintf(viewer," Control[UMFPACK_DENSE_COL]: %g\n",lu->Control[UMFPACK_DENSE_COL]);CHKERRQ(ierr); 291 ierr = PetscViewerASCIIPrintf(viewer," Control[UMFPACK_DENSE_ROW]: %g\n",lu->Control[UMFPACK_DENSE_ROW]);CHKERRQ(ierr); 292 ierr = PetscViewerASCIIPrintf(viewer," Control[UMFPACK_AMD_DENSE]: %g\n",lu->Control[UMFPACK_AMD_DENSE]);CHKERRQ(ierr); 293 ierr = PetscViewerASCIIPrintf(viewer," Control[UMFPACK_BLOCK_SIZE]: %g\n",lu->Control[UMFPACK_BLOCK_SIZE]);CHKERRQ(ierr); 294 ierr = PetscViewerASCIIPrintf(viewer," Control[UMFPACK_2BY2_TOLERANCE]: %g\n",lu->Control[UMFPACK_2BY2_TOLERANCE]);CHKERRQ(ierr); 295 ierr = PetscViewerASCIIPrintf(viewer," Control[UMFPACK_FIXQ]: %g\n",lu->Control[UMFPACK_FIXQ]);CHKERRQ(ierr); 296 ierr = PetscViewerASCIIPrintf(viewer," Control[UMFPACK_AGGRESSIVE]: %g\n",lu->Control[UMFPACK_AGGRESSIVE]);CHKERRQ(ierr); 297 298 /* Control parameters used by numeric factorization */ 299 ierr = PetscViewerASCIIPrintf(viewer," Control[UMFPACK_PIVOT_TOLERANCE]: %g\n",lu->Control[UMFPACK_PIVOT_TOLERANCE]);CHKERRQ(ierr); 300 ierr = PetscViewerASCIIPrintf(viewer," Control[UMFPACK_SYM_PIVOT_TOLERANCE]: %g\n",lu->Control[UMFPACK_SYM_PIVOT_TOLERANCE]);CHKERRQ(ierr); 301 ierr = PetscViewerASCIIPrintf(viewer," Control[UMFPACK_SCALE]: %g\n",lu->Control[UMFPACK_SCALE]);CHKERRQ(ierr); 302 ierr = PetscViewerASCIIPrintf(viewer," Control[UMFPACK_ALLOC_INIT]: %g\n",lu->Control[UMFPACK_ALLOC_INIT]);CHKERRQ(ierr); 303 ierr = PetscViewerASCIIPrintf(viewer," Control[UMFPACK_DROPTOL]: %g\n",lu->Control[UMFPACK_DROPTOL]);CHKERRQ(ierr); 304 305 /* Control parameters used by solve */ 306 ierr = PetscViewerASCIIPrintf(viewer," Control[UMFPACK_IRSTEP]: %g\n",lu->Control[UMFPACK_IRSTEP]);CHKERRQ(ierr); 307 308 /* mat ordering */ 309 if(!lu->PetscMatOdering) ierr = PetscViewerASCIIPrintf(viewer," UMFPACK default matrix ordering is used (not the PETSc matrix ordering) \n");CHKERRQ(ierr); 310 311 PetscFunctionReturn(0); 312 } 313 314 #undef __FUNCT__ 315 #define __FUNCT__ "MatView_UMFPACK" 316 PetscErrorCode MatView_UMFPACK(Mat A,PetscViewer viewer) 317 { 318 PetscErrorCode ierr; 319 PetscTruth iascii; 320 PetscViewerFormat format; 321 Mat_UMFPACK *lu=(Mat_UMFPACK*)(A->spptr); 322 323 PetscFunctionBegin; 324 ierr = (*lu->MatView)(A,viewer);CHKERRQ(ierr); 325 326 ierr = PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_ASCII,&iascii);CHKERRQ(ierr); 327 if (iascii) { 328 ierr = PetscViewerGetFormat(viewer,&format);CHKERRQ(ierr); 329 if (format == PETSC_VIEWER_ASCII_INFO) { 330 ierr = MatFactorInfo_UMFPACK(A,viewer);CHKERRQ(ierr); 331 } 332 } 333 PetscFunctionReturn(0); 334 } 335 336 EXTERN_C_BEGIN 337 #undef __FUNCT__ 338 #define __FUNCT__ "MatConvert_SeqAIJ_UMFPACK" 339 PetscErrorCode PETSCMAT_DLLEXPORT MatConvert_SeqAIJ_UMFPACK(Mat A,MatType type,MatReuse reuse,Mat *newmat) 340 { 341 PetscErrorCode ierr; 342 Mat B=*newmat; 343 Mat_UMFPACK *lu; 344 345 PetscFunctionBegin; 346 if (reuse == MAT_INITIAL_MATRIX) { 347 ierr = MatDuplicate(A,MAT_COPY_VALUES,&B);CHKERRQ(ierr); 348 } 349 350 ierr = PetscNew(Mat_UMFPACK,&lu);CHKERRQ(ierr); 351 lu->MatDuplicate = A->ops->duplicate; 352 lu->MatView = A->ops->view; 353 lu->MatAssemblyEnd = A->ops->assemblyend; 354 lu->MatLUFactorSymbolic = A->ops->lufactorsymbolic; 355 lu->MatDestroy = A->ops->destroy; 356 lu->CleanUpUMFPACK = PETSC_FALSE; 357 358 B->spptr = (void*)lu; 359 B->ops->duplicate = MatDuplicate_UMFPACK; 360 B->ops->view = MatView_UMFPACK; 361 B->ops->assemblyend = MatAssemblyEnd_UMFPACK; 362 B->ops->lufactorsymbolic = MatLUFactorSymbolic_UMFPACK; 363 B->ops->destroy = MatDestroy_UMFPACK; 364 365 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatConvert_seqaij_umfpack_C", 366 "MatConvert_SeqAIJ_UMFPACK",MatConvert_SeqAIJ_UMFPACK);CHKERRQ(ierr); 367 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatConvert_umfpack_seqaij_C", 368 "MatConvert_UMFPACK_SeqAIJ",MatConvert_UMFPACK_SeqAIJ);CHKERRQ(ierr); 369 ierr = PetscInfo(0,"Using UMFPACK for SeqAIJ LU factorization and solves.\n");CHKERRQ(ierr); 370 ierr = PetscObjectChangeTypeName((PetscObject)B,MATUMFPACK);CHKERRQ(ierr); 371 *newmat = B; 372 PetscFunctionReturn(0); 373 } 374 EXTERN_C_END 375 376 #undef __FUNCT__ 377 #define __FUNCT__ "MatDuplicate_UMFPACK" 378 PetscErrorCode MatDuplicate_UMFPACK(Mat A, MatDuplicateOption op, Mat *M) 379 { 380 PetscErrorCode ierr; 381 Mat_UMFPACK *lu=(Mat_UMFPACK*)A->spptr; 382 383 PetscFunctionBegin; 384 ierr = (*lu->MatDuplicate)(A,op,M);CHKERRQ(ierr); 385 ierr = PetscMemcpy((*M)->spptr,lu,sizeof(Mat_UMFPACK));CHKERRQ(ierr); 386 PetscFunctionReturn(0); 387 } 388 389 /*MC 390 MATUMFPACK - MATUMFPACK = "umfpack" - A matrix type providing direct solvers (LU) for sequential matrices 391 via the external package UMFPACK. 392 393 If UMFPACK is installed (see the manual for 394 instructions on how to declare the existence of external packages), 395 a matrix type can be constructed which invokes UMFPACK solvers. 396 After calling MatCreate(...,A), simply call MatSetType(A,UMFPACK). 397 This matrix type is only supported for double precision real. 398 399 This matrix inherits from MATSEQAIJ. As a result, MatSeqAIJSetPreallocation is 400 supported for this matrix type. One can also call MatConvert for an inplace conversion to or from 401 the MATSEQAIJ type without data copy. 402 403 Consult UMFPACK documentation for more information about the Control parameters 404 which correspond to the options database keys below. 405 406 Options Database Keys: 407 + -mat_type umfpack - sets the matrix type to "umfpack" during a call to MatSetFromOptions() 408 . -mat_umfpack_prl - UMFPACK print level: Control[UMFPACK_PRL] 409 . -mat_umfpack_dense_col <alpha_c> - UMFPACK dense column threshold: Control[UMFPACK_DENSE_COL] 410 . -mat_umfpack_block_size <bs> - UMFPACK block size for BLAS-Level 3 calls: Control[UMFPACK_BLOCK_SIZE] 411 . -mat_umfpack_pivot_tolerance <delta> - UMFPACK partial pivot tolerance: Control[UMFPACK_PIVOT_TOLERANCE] 412 . -mat_umfpack_alloc_init <delta> - UMFPACK factorized matrix allocation modifier: Control[UMFPACK_ALLOC_INIT] 413 - -mat_umfpack_irstep <maxit> - UMFPACK maximum number of iterative refinement steps: Control[UMFPACK_IRSTEP] 414 415 Level: beginner 416 417 .seealso: PCLU 418 M*/ 419 420 EXTERN_C_BEGIN 421 #undef __FUNCT__ 422 #define __FUNCT__ "MatCreate_UMFPACK" 423 PetscErrorCode PETSCMAT_DLLEXPORT MatCreate_UMFPACK(Mat A) 424 { 425 PetscErrorCode ierr; 426 427 PetscFunctionBegin; 428 ierr = MatSetType(A,MATSEQAIJ);CHKERRQ(ierr); 429 ierr = MatConvert_SeqAIJ_UMFPACK(A,MATUMFPACK,MAT_REUSE_MATRIX,&A);CHKERRQ(ierr); 430 PetscFunctionReturn(0); 431 } 432 EXTERN_C_END 433