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