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