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 (*MatView)(Mat,PetscViewer); 23 int (*MatAssemblyEnd)(Mat,MatAssemblyType); 24 int (*MatLUFactorSymbolic)(Mat,IS,IS,MatFactorInfo*,Mat*); 25 int (*MatDestroy)(Mat); 26 27 /* Flag to clean up UMFPACK objects during Destroy */ 28 PetscTruth CleanUpUMFPACK; 29 } Mat_SeqAIJ_UMFPACK; 30 31 EXTERN_C_BEGIN 32 #undef __FUNCT__ 33 #define __FUNCT__ "MatConvert_UMFPACK_SeqAIJ" 34 int MatConvert_UMFPACK_SeqAIJ(Mat A,MatType type,Mat *newmat) { 35 /* This routine is only called to convert an unfactored PETSc-UMFPACK matrix */ 36 /* to its base PETSc type, so we will ignore 'MatType type'. */ 37 int ierr; 38 Mat B=*newmat; 39 Mat_SeqAIJ_UMFPACK *lu=(Mat_SeqAIJ_UMFPACK*)A->spptr; 40 41 PetscFunctionBegin; 42 if (B != A) { 43 /* This routine was inherited from SeqAIJ. */ 44 ierr = MatDuplicate(A,MAT_COPY_VALUES,&B);CHKERRQ(ierr); 45 } else { 46 /* Reset the original function pointers */ 47 A->ops->view = lu->MatView; 48 A->ops->assemblyend = lu->MatAssemblyEnd; 49 A->ops->lufactorsymbolic = lu->MatLUFactorSymbolic; 50 A->ops->destroy = lu->MatDestroy; 51 52 ierr = PetscFree(lu);CHKERRQ(ierr); 53 ierr = PetscObjectChangeTypeName((PetscObject)B,MATSEQAIJ);CHKERRQ(ierr); 54 } 55 *newmat = B; 56 PetscFunctionReturn(0); 57 } 58 EXTERN_C_END 59 60 #undef __FUNCT__ 61 #define __FUNCT__ "MatDestroy_SeqAIJ_UMFPACK" 62 int MatDestroy_SeqAIJ_UMFPACK(Mat A) 63 { 64 int ierr; 65 Mat_SeqAIJ_UMFPACK *lu = (Mat_SeqAIJ_UMFPACK*)A->spptr; 66 67 PetscFunctionBegin; 68 if (lu->CleanUpUMFPACK) { 69 umfpack_di_free_symbolic(&lu->Symbolic) ; 70 umfpack_di_free_numeric(&lu->Numeric) ; 71 ierr = PetscFree(lu->Wi);CHKERRQ(ierr); 72 ierr = PetscFree(lu->W);CHKERRQ(ierr); 73 74 if (lu->PetscMatOdering) { 75 ierr = PetscFree(lu->perm_c);CHKERRQ(ierr); 76 } 77 } 78 ierr = MatConvert_UMFPACK_SeqAIJ(A,MATSEQAIJ,&A);CHKERRQ(ierr); 79 ierr = (*A->ops->destroy)(A);CHKERRQ(ierr); 80 PetscFunctionReturn(0); 81 } 82 83 #undef __FUNCT__ 84 #define __FUNCT__ "MatSolve_SeqAIJ_UMFPACK" 85 int MatSolve_SeqAIJ_UMFPACK(Mat A,Vec b,Vec x) 86 { 87 Mat_SeqAIJ_UMFPACK *lu = (Mat_SeqAIJ_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 PetscFunctionReturn(0); 107 } 108 109 #undef __FUNCT__ 110 #define __FUNCT__ "MatLUFactorNumeric_SeqAIJ_UMFPACK" 111 int MatLUFactorNumeric_SeqAIJ_UMFPACK(Mat A,Mat *F) 112 { 113 Mat_SeqAIJ_UMFPACK *lu = (Mat_SeqAIJ_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_SeqAIJ_UMFPACK" 141 int MatLUFactorSymbolic_SeqAIJ_UMFPACK(Mat A,IS r,IS c,MatFactorInfo *info,Mat *F) 142 { 143 Mat B; 144 Mat_SeqAIJ *mat=(Mat_SeqAIJ*)A->data; 145 Mat_SeqAIJ_UMFPACK *lu; 146 int ierr,m=A->m,n=A->n,*ai=mat->i,*aj=mat->j,status,*ca; 147 PetscScalar *av=mat->a; 148 149 PetscFunctionBegin; 150 /* Create the factorization matrix F */ 151 ierr = MatCreate(A->comm,PETSC_DECIDE,PETSC_DECIDE,m,n,&B);CHKERRQ(ierr); 152 ierr = MatSetType(B,MATUMFPACK);CHKERRQ(ierr); 153 ierr = MatSeqAIJSetPreallocation(B,0,PETSC_NULL);CHKERRQ(ierr); 154 155 B->ops->lufactornumeric = MatLUFactorNumeric_SeqAIJ_UMFPACK; 156 B->ops->solve = MatSolve_SeqAIJ_UMFPACK; 157 B->factor = FACTOR_LU; 158 B->assembled = PETSC_TRUE; /* required by -sles_view */ 159 160 lu = (Mat_SeqAIJ_UMFPACK*)(B->spptr); 161 162 /* initializations */ 163 /* ------------------------------------------------*/ 164 /* get the default control parameters */ 165 umfpack_di_defaults (lu->Control) ; 166 lu->perm_c = PETSC_NULL; /* use defaul UMFPACK col permutation */ 167 168 ierr = PetscOptionsBegin(A->comm,A->prefix,"UMFPACK Options","Mat");CHKERRQ(ierr); 169 /* Control parameters used by reporting routiones */ 170 ierr = PetscOptionsReal("-mat_umfpack_prl","Control[UMFPACK_PRL]","None",lu->Control[UMFPACK_PRL],&lu->Control[UMFPACK_PRL],PETSC_NULL);CHKERRQ(ierr); 171 172 /* Control parameters for symbolic factorization */ 173 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); 174 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); 175 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); 176 177 /* Control parameters used by numeric factorization */ 178 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); 179 #if !defined(PETSC_HAVE_UMFPACK_41_OR_NEWER) 180 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); 181 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); 182 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); 183 #endif 184 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); 185 186 /* Control parameters used by solve */ 187 ierr = PetscOptionsReal("-mat_umfpack_irstep","Control[UMFPACK_IRSTEP]","None",lu->Control[UMFPACK_IRSTEP],&lu->Control[UMFPACK_IRSTEP],PETSC_NULL);CHKERRQ(ierr); 188 189 /* use Petsc mat ordering (notice size is for the transpose) */ 190 ierr = PetscOptionsHasName(PETSC_NULL,"-pc_lu_mat_ordering_type",&lu->PetscMatOdering);CHKERRQ(ierr); 191 if (lu->PetscMatOdering) { 192 ierr = ISGetIndices(c,&ca);CHKERRQ(ierr); 193 ierr = PetscMalloc(A->m*sizeof(int),&lu->perm_c);CHKERRQ(ierr); 194 ierr = PetscMemcpy(lu->perm_c,ca,A->m*sizeof(int));CHKERRQ(ierr); 195 ierr = ISRestoreIndices(c,&ca);CHKERRQ(ierr); 196 } 197 PetscOptionsEnd(); 198 199 /* print the control parameters */ 200 if( lu->Control[UMFPACK_PRL] > 1 ) umfpack_di_report_control (lu->Control); 201 202 /* symbolic factorization of A' */ 203 /* ---------------------------------------------------------------------- */ 204 #if defined(PETSC_HAVE_UMFPACK_41_OR_NEWER) 205 status = umfpack_di_qsymbolic(n,m,ai,aj,PETSC_NULL,lu->perm_c,&lu->Symbolic,lu->Control,lu->Info) ; 206 #else 207 status = umfpack_di_qsymbolic(n,m,ai,aj,lu->perm_c,&lu->Symbolic,lu->Control,lu->Info) ; 208 #endif 209 if (status < 0){ 210 umfpack_di_report_info(lu->Control, lu->Info) ; 211 umfpack_di_report_status(lu->Control, status) ; 212 SETERRQ(1,"umfpack_di_symbolic failed"); 213 } 214 /* report sumbolic factorization of A' when Control[PRL] > 3 */ 215 (void) umfpack_di_report_symbolic(lu->Symbolic, lu->Control) ; 216 217 lu->flg = DIFFERENT_NONZERO_PATTERN; 218 lu->ai = ai; 219 lu->aj = aj; 220 lu->av = av; 221 lu->CleanUpUMFPACK = PETSC_TRUE; 222 *F = B; 223 PetscFunctionReturn(0); 224 } 225 226 #undef __FUNCT__ 227 #define __FUNCT__ "MatAssemblyEnd_SeqAIJ_UMFPACK" 228 int MatAssemblyEnd_SeqAIJ_UMFPACK(Mat A,MatAssemblyType mode) { 229 int ierr; 230 Mat_SeqAIJ_UMFPACK *lu=(Mat_SeqAIJ_UMFPACK*)(A->spptr); 231 232 PetscFunctionBegin; 233 ierr = (*lu->MatAssemblyEnd)(A,mode);CHKERRQ(ierr); 234 lu->MatLUFactorSymbolic = A->ops->lufactorsymbolic; 235 A->ops->lufactorsymbolic = MatLUFactorSymbolic_SeqAIJ_UMFPACK; 236 PetscFunctionReturn(0); 237 } 238 239 /* used by -sles_view */ 240 #undef __FUNCT__ 241 #define __FUNCT__ "MatSeqAIJFactorInfo_UMFPACK" 242 int MatSeqAIJFactorInfo_UMFPACK(Mat A,PetscViewer viewer) 243 { 244 Mat_SeqAIJ_UMFPACK *lu= (Mat_SeqAIJ_UMFPACK*)A->spptr; 245 int ierr; 246 PetscFunctionBegin; 247 /* check if matrix is UMFPACK type */ 248 if (A->ops->solve != MatSolve_SeqAIJ_UMFPACK) PetscFunctionReturn(0); 249 250 ierr = PetscViewerASCIIPrintf(viewer,"UMFPACK run parameters:\n");CHKERRQ(ierr); 251 /* Control parameters used by reporting routiones */ 252 ierr = PetscViewerASCIIPrintf(viewer," Control[UMFPACK_PRL]: %g\n",lu->Control[UMFPACK_PRL]);CHKERRQ(ierr); 253 254 /* Control parameters used by symbolic factorization */ 255 ierr = PetscViewerASCIIPrintf(viewer," Control[UMFPACK_DENSE_COL]: %g\n",lu->Control[UMFPACK_DENSE_COL]);CHKERRQ(ierr); 256 ierr = PetscViewerASCIIPrintf(viewer," Control[UMFPACK_DENSE_ROW]: %g\n",lu->Control[UMFPACK_DENSE_ROW]);CHKERRQ(ierr); 257 ierr = PetscViewerASCIIPrintf(viewer," Control[UMFPACK_BLOCK_SIZE]: %g\n",lu->Control[UMFPACK_BLOCK_SIZE]);CHKERRQ(ierr); 258 259 /* Control parameters used by numeric factorization */ 260 ierr = PetscViewerASCIIPrintf(viewer," Control[UMFPACK_PIVOT_TOLERANCE]: %g\n",lu->Control[UMFPACK_PIVOT_TOLERANCE]);CHKERRQ(ierr); 261 #if !defined(PETSC_HAVE_UMFPACK_41_OR_NEWER) 262 ierr = PetscViewerASCIIPrintf(viewer," Control[UMFPACK_RELAXED_AMALGAMATION]: %g\n",lu->Control[UMFPACK_RELAXED_AMALGAMATION]);CHKERRQ(ierr); 263 ierr = PetscViewerASCIIPrintf(viewer," Control[UMFPACK_RELAXED2_AMALGAMATION]: %g\n",lu->Control[UMFPACK_RELAXED2_AMALGAMATION]);CHKERRQ(ierr); 264 ierr = PetscViewerASCIIPrintf(viewer," Control[UMFPACK_RELAXED3_AMALGAMATION]: %g\n",lu->Control[UMFPACK_RELAXED3_AMALGAMATION]);CHKERRQ(ierr); 265 #endif 266 ierr = PetscViewerASCIIPrintf(viewer," Control[UMFPACK_ALLOC_INIT]: %g\n",lu->Control[UMFPACK_ALLOC_INIT]);CHKERRQ(ierr); 267 268 /* Control parameters used by solve */ 269 ierr = PetscViewerASCIIPrintf(viewer," Control[UMFPACK_IRSTEP]: %g\n",lu->Control[UMFPACK_IRSTEP]);CHKERRQ(ierr); 270 271 /* mat ordering */ 272 if(!lu->PetscMatOdering) ierr = PetscViewerASCIIPrintf(viewer," UMFPACK default matrix ordering is used (not the PETSc matrix ordering) \n");CHKERRQ(ierr); 273 274 PetscFunctionReturn(0); 275 } 276 277 #undef __FUNCT__ 278 #define __FUNCT__ "MatView_SeqAIJ_UMFPACK" 279 int MatView_SeqAIJ_UMFPACK(Mat A,PetscViewer viewer) 280 { 281 int ierr; 282 PetscTruth isascii; 283 PetscViewerFormat format; 284 Mat_SeqAIJ_UMFPACK *lu=(Mat_SeqAIJ_UMFPACK*)(A->spptr); 285 286 PetscFunctionBegin; 287 ierr = (*lu->MatView)(A,viewer);CHKERRQ(ierr); 288 289 ierr = PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_ASCII,&isascii);CHKERRQ(ierr); 290 if (isascii) { 291 ierr = PetscViewerGetFormat(viewer,&format);CHKERRQ(ierr); 292 if (format == PETSC_VIEWER_ASCII_FACTOR_INFO) { 293 ierr = MatSeqAIJFactorInfo_UMFPACK(A,viewer);CHKERRQ(ierr); 294 } 295 } 296 PetscFunctionReturn(0); 297 } 298 299 EXTERN_C_BEGIN 300 #undef __FUNCT__ 301 #define __FUNCT__ "MatConvert_SeqAIJ_UMFPACK" 302 int MatConvert_SeqAIJ_UMFPACK(Mat A,MatType type,Mat *newmat) { 303 /* This routine is only called to convert to MATUMFPACK */ 304 /* from MATSEQAIJ, so we will ignore 'MatType type'. */ 305 int ierr; 306 Mat B=*newmat; 307 Mat_SeqAIJ_UMFPACK *lu; 308 309 PetscFunctionBegin; 310 if (B != A) { 311 ierr = MatDuplicate(A,MAT_COPY_VALUES,&B);CHKERRQ(ierr); 312 } 313 314 ierr = PetscNew(Mat_SeqAIJ_UMFPACK,&lu);CHKERRQ(ierr); 315 lu->MatView = A->ops->view; 316 lu->MatAssemblyEnd = A->ops->assemblyend; 317 lu->MatLUFactorSymbolic = A->ops->lufactorsymbolic; 318 lu->MatDestroy = A->ops->destroy; 319 lu->CleanUpUMFPACK = PETSC_FALSE; 320 321 B->spptr = (void*)lu; 322 B->ops->view = MatView_SeqAIJ_UMFPACK; 323 B->ops->assemblyend = MatAssemblyEnd_SeqAIJ_UMFPACK; 324 B->ops->lufactorsymbolic = MatLUFactorSymbolic_SeqAIJ_UMFPACK; 325 B->ops->destroy = MatDestroy_SeqAIJ_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 /*MC 339 MATUMFPACK - a matrix type providing direct solvers (LU) for sequential matrices 340 via the external package UMFPACK. 341 342 If UMFPACK is installed (see the manual for 343 instructions on how to declare the existence of external packages), 344 a matrix type can be constructed which invokes UMFPACK solvers. 345 After calling MatCreate(...,A), simply call MatSetType(A,UMFPACK). 346 This matrix type is only supported for double precision real. 347 348 This matrix inherits from MATSEQAIJ. As a result, MatSeqAIJSetPreallocation is 349 supported for this matrix type. 350 351 Consult UMFPACK documentation for more information about the Control parameters 352 which correspond to the options database keys below. 353 354 Options Database Keys: 355 + -mat_type umfpack - sets the matrix type to umfpack during a call to MatSetFromOptions() 356 . -mat_umfpack_prl - UMFPACK print level: Control[UMFPACK_PRL] 357 . -mat_umfpack_dense_col <alpha_c> - UMFPACK dense column threshold: Control[UMFPACK_DENSE_COL] 358 . -mat_umfpack_block_size <bs> - UMFPACK block size for BLAS-Level 3 calls: Control[UMFPACK_BLOCK_SIZE] 359 . -mat_umfpack_pivot_tolerance <delta> - UMFPACK partial pivot tolerance: Control[UMFPACK_PIVOT_TOLERANCE] 360 . -mat_umfpack_alloc_init <delta> - UMFPACK factorized matrix allocation modifier: Control[UMFPACK_ALLOC_INIT] 361 - -mat_umfpack_irstep <maxit> - UMFPACK maximum number of iterative refinement steps: Control[UMFPACK_IRSTEP] 362 363 Level: beginner 364 365 .seealso: PCLU 366 M*/ 367 368 EXTERN_C_BEGIN 369 #undef __FUNCT__ 370 #define __FUNCT__ "MatCreate_SeqAIJ_UMFPACK" 371 int MatCreate_SeqAIJ_UMFPACK(Mat A) { 372 int ierr; 373 374 PetscFunctionBegin; 375 ierr = MatSetType(A,MATSEQAIJ);CHKERRQ(ierr); 376 ierr = MatConvert_SeqAIJ_UMFPACK(A,MATUMFPACK,&A);CHKERRQ(ierr); 377 PetscFunctionReturn(0); 378 } 379 EXTERN_C_END 380