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