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 (*MatDestroy)(Mat); 25 26 /* Flag to clean up UMFPACK objects during Destroy */ 27 PetscTruth CleanUpUMFPACK; 28 } Mat_SeqAIJ_UMFPACK; 29 30 #undef __FUNCT__ 31 #define __FUNCT__ "MatDestroy_SeqAIJ_UMFPACK" 32 int MatDestroy_SeqAIJ_UMFPACK(Mat A) 33 { 34 Mat_SeqAIJ_UMFPACK *lu = (Mat_SeqAIJ_UMFPACK*)A->spptr; 35 int ierr,(*destroy)(Mat); 36 37 PetscFunctionBegin; 38 if (lu->CleanUpUMFPACK) { 39 umfpack_di_free_symbolic(&lu->Symbolic) ; 40 umfpack_di_free_numeric(&lu->Numeric) ; 41 ierr = PetscFree(lu->Wi);CHKERRQ(ierr); 42 ierr = PetscFree(lu->W);CHKERRQ(ierr); 43 44 if (lu->PetscMatOdering) { 45 ierr = PetscFree(lu->perm_c);CHKERRQ(ierr); 46 } 47 } 48 49 destroy = lu->MatDestroy; 50 ierr = PetscFree(lu);CHKERRQ(ierr); 51 ierr = (*destroy)(A);CHKERRQ(ierr); 52 PetscFunctionReturn(0); 53 } 54 55 #undef __FUNCT__ 56 #define __FUNCT__ "MatView_SeqAIJ_UMFPACK" 57 int MatView_SeqAIJ_UMFPACK(Mat A,PetscViewer viewer) 58 { 59 int ierr; 60 PetscTruth isascii; 61 PetscViewerFormat format; 62 Mat_SeqAIJ_UMFPACK *lu=(Mat_SeqAIJ_UMFPACK*)(A->spptr); 63 64 PetscFunctionBegin; 65 ierr = (*lu->MatView)(A,viewer);CHKERRQ(ierr); 66 67 ierr = PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_ASCII,&isascii);CHKERRQ(ierr); 68 if (isascii) { 69 ierr = PetscViewerGetFormat(viewer,&format);CHKERRQ(ierr); 70 if (format == PETSC_VIEWER_ASCII_FACTOR_INFO) { 71 ierr = MatSeqAIJFactorInfo_UMFPACK(A,viewer);CHKERRQ(ierr); 72 } 73 } 74 PetscFunctionReturn(0); 75 } 76 77 #undef __FUNCT__ 78 #define __FUNCT__ "MatAssemblyEnd_SeqAIJ_UMFPACK" 79 int MatAssemblyEnd_SeqAIJ_UMFPACK(Mat A,MatAssemblyType mode) { 80 int ierr; 81 Mat_SeqAIJ_UMFPACK *lu=(Mat_SeqAIJ_UMFPACK*)(A->spptr); 82 83 PetscFunctionBegin; 84 ierr = (*lu->MatAssemblyEnd)(A,mode);CHKERRQ(ierr); 85 ierr = MatUseUMFPACK_SeqAIJ(A);CHKERRQ(ierr); 86 PetscFunctionReturn(0); 87 } 88 89 #undef __FUNCT__ 90 #define __FUNCT__ "MatSolve_SeqAIJ_UMFPACK" 91 int MatSolve_SeqAIJ_UMFPACK(Mat A,Vec b,Vec x) 92 { 93 Mat_SeqAIJ_UMFPACK *lu = (Mat_SeqAIJ_UMFPACK*)A->spptr; 94 PetscScalar *av=lu->av,*ba,*xa; 95 int ierr,*ai=lu->ai,*aj=lu->aj,status; 96 97 PetscFunctionBegin; 98 /* solve Ax = b by umfpack_di_wsolve */ 99 /* ----------------------------------*/ 100 ierr = VecGetArray(b,&ba); 101 ierr = VecGetArray(x,&xa); 102 103 status = umfpack_di_wsolve(UMFPACK_At,ai,aj,av,xa,ba,lu->Numeric,lu->Control,lu->Info,lu->Wi,lu->W); 104 umfpack_di_report_info(lu->Control, lu->Info); 105 if (status < 0){ 106 umfpack_di_report_status(lu->Control, status) ; 107 SETERRQ(1,"umfpack_di_wsolve failed") ; 108 } 109 110 ierr = VecRestoreArray(b,&ba); 111 ierr = VecRestoreArray(x,&xa); 112 PetscFunctionReturn(0); 113 } 114 115 #undef __FUNCT__ 116 #define __FUNCT__ "MatLUFactorNumeric_SeqAIJ_UMFPACK" 117 int MatLUFactorNumeric_SeqAIJ_UMFPACK(Mat A,Mat *F) 118 { 119 Mat_SeqAIJ_UMFPACK *lu = (Mat_SeqAIJ_UMFPACK*)(*F)->spptr; 120 int *ai=lu->ai,*aj=lu->aj,m=A->m,status,ierr; 121 PetscScalar *av=lu->av; 122 123 PetscFunctionBegin; 124 /* numeric factorization of A' */ 125 /* ----------------------------*/ 126 status = umfpack_di_numeric (ai,aj,av,lu->Symbolic,&lu->Numeric,lu->Control,lu->Info) ; 127 if (status < 0) SETERRQ(1,"umfpack_di_numeric failed"); 128 /* report numeric factorization of A' when Control[PRL] > 3 */ 129 (void) umfpack_di_report_numeric (lu->Numeric, lu->Control) ; 130 131 if (lu->flg == DIFFERENT_NONZERO_PATTERN){ /* first numeric factorization */ 132 /* allocate working space to be used by Solve */ 133 ierr = PetscMalloc(m * sizeof(int), &lu->Wi);CHKERRQ(ierr); 134 ierr = PetscMalloc(5*m * sizeof(double), &lu->W);CHKERRQ(ierr); 135 136 lu->flg = SAME_NONZERO_PATTERN; 137 } 138 139 PetscFunctionReturn(0); 140 } 141 142 /* 143 Note the r permutation is ignored 144 */ 145 #undef __FUNCT__ 146 #define __FUNCT__ "MatLUFactorSymbolic_SeqAIJ_UMFPACK" 147 int MatLUFactorSymbolic_SeqAIJ_UMFPACK(Mat A,IS r,IS c,MatFactorInfo *info,Mat *F) 148 { 149 Mat_SeqAIJ *mat=(Mat_SeqAIJ*)A->data; 150 Mat_SeqAIJ_UMFPACK *lu; 151 int ierr,m=A->m,n=A->n,*ai=mat->i,*aj=mat->j,status,*ca; 152 PetscScalar *av=mat->a; 153 154 PetscFunctionBegin; 155 /* Create the factorization matrix F */ 156 ierr = MatCreateSeqAIJ(A->comm,m,n,PETSC_NULL,PETSC_NULL,F);CHKERRQ(ierr); 157 158 (*F)->ops->lufactornumeric = MatLUFactorNumeric_SeqAIJ_UMFPACK; 159 (*F)->ops->solve = MatSolve_SeqAIJ_UMFPACK; 160 (*F)->ops->destroy = MatDestroy_SeqAIJ_UMFPACK; 161 (*F)->factor = FACTOR_LU; 162 (*F)->assembled = PETSC_TRUE; /* required by -sles_view */ 163 164 ierr = PetscNew(Mat_SeqAIJ_UMFPACK,&lu);CHKERRQ(ierr); 165 (*F)->spptr = (void*)lu; 166 167 /* initializations */ 168 /* ------------------------------------------------*/ 169 /* get the default control parameters */ 170 umfpack_di_defaults (lu->Control) ; 171 lu->perm_c = PETSC_NULL; /* use defaul UMFPACK col permutation */ 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 = PetscOptionsReal("-mat_umfpack_dense_col","Control[UMFPACK_DENSE_COL]","None",lu->Control[UMFPACK_DENSE_COL],&lu->Control[UMFPACK_DENSE_COL],PETSC_NULL);CHKERRQ(ierr); 179 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); 180 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); 181 182 /* Control parameters used by numeric factorization */ 183 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); 184 #if !defined(PETSC_HAVE_UMFPACK_42_OR_NEWER) 185 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); 186 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); 187 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); 188 #endif 189 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); 190 191 /* Control parameters used by solve */ 192 ierr = PetscOptionsReal("-mat_umfpack_irstep","Control[UMFPACK_IRSTEP]","None",lu->Control[UMFPACK_IRSTEP],&lu->Control[UMFPACK_IRSTEP],PETSC_NULL);CHKERRQ(ierr); 193 194 /* use Petsc mat ordering (notice size is for the transpose) */ 195 ierr = PetscOptionsHasName(PETSC_NULL,"-pc_lu_mat_ordering_type",&lu->PetscMatOdering);CHKERRQ(ierr); 196 if (lu->PetscMatOdering) { 197 ierr = ISGetIndices(c,&ca);CHKERRQ(ierr); 198 ierr = PetscMalloc(A->m*sizeof(int),&lu->perm_c);CHKERRQ(ierr); 199 ierr = PetscMemcpy(lu->perm_c,ca,A->m*sizeof(int));CHKERRQ(ierr); 200 ierr = ISRestoreIndices(c,&ca);CHKERRQ(ierr); 201 } 202 PetscOptionsEnd(); 203 204 /* print the control parameters */ 205 if( lu->Control[UMFPACK_PRL] > 1 ) umfpack_di_report_control (lu->Control); 206 207 /* symbolic factorization of A' */ 208 /* ---------------------------------------------------------------------- */ 209 status = umfpack_di_qsymbolic(n,m,ai,aj,lu->perm_c,&lu->Symbolic,lu->Control,lu->Info) ; 210 if (status < 0){ 211 umfpack_di_report_info(lu->Control, lu->Info) ; 212 umfpack_di_report_status(lu->Control, status) ; 213 SETERRQ(1,"umfpack_di_symbolic failed"); 214 } 215 /* report sumbolic factorization of A' when Control[PRL] > 3 */ 216 (void) umfpack_di_report_symbolic(lu->Symbolic, lu->Control) ; 217 218 lu->flg = DIFFERENT_NONZERO_PATTERN; 219 lu->ai = ai; 220 lu->aj = aj; 221 lu->av = av; 222 PetscFunctionReturn(0); 223 } 224 225 #undef __FUNCT__ 226 #define __FUNCT__ "MatUseUMFPACK_SeqAIJ" 227 int MatUseUMFPACK_SeqAIJ(Mat A) 228 { 229 PetscFunctionBegin; 230 A->ops->lufactorsymbolic = MatLUFactorSymbolic_SeqAIJ_UMFPACK; 231 PetscFunctionReturn(0); 232 } 233 234 /* used by -sles_view */ 235 #undef __FUNCT__ 236 #define __FUNCT__ "MatSeqAIJFactorInfo_UMFPACK" 237 int MatSeqAIJFactorInfo_UMFPACK(Mat A,PetscViewer viewer) 238 { 239 Mat_SeqAIJ_UMFPACK *lu= (Mat_SeqAIJ_UMFPACK*)A->spptr; 240 int ierr; 241 PetscFunctionBegin; 242 /* check if matrix is UMFPACK type */ 243 if (A->ops->solve != MatSolve_SeqAIJ_UMFPACK) PetscFunctionReturn(0); 244 245 ierr = PetscViewerASCIIPrintf(viewer,"UMFPACK run parameters:\n");CHKERRQ(ierr); 246 /* Control parameters used by reporting routiones */ 247 ierr = PetscViewerASCIIPrintf(viewer," Control[UMFPACK_PRL]: %g\n",lu->Control[UMFPACK_PRL]);CHKERRQ(ierr); 248 249 /* Control parameters used by symbolic factorization */ 250 ierr = PetscViewerASCIIPrintf(viewer," Control[UMFPACK_DENSE_COL]: %g\n",lu->Control[UMFPACK_DENSE_COL]);CHKERRQ(ierr); 251 ierr = PetscViewerASCIIPrintf(viewer," Control[UMFPACK_DENSE_ROW]: %g\n",lu->Control[UMFPACK_DENSE_ROW]);CHKERRQ(ierr); 252 ierr = PetscViewerASCIIPrintf(viewer," Control[UMFPACK_BLOCK_SIZE]: %g\n",lu->Control[UMFPACK_BLOCK_SIZE]);CHKERRQ(ierr); 253 254 /* Control parameters used by numeric factorization */ 255 ierr = PetscViewerASCIIPrintf(viewer," Control[UMFPACK_PIVOT_TOLERANCE]: %g\n",lu->Control[UMFPACK_PIVOT_TOLERANCE]);CHKERRQ(ierr); 256 #if !defined(PETSC_HAVE_UMFPACK_42_OR_NEWER) 257 ierr = PetscViewerASCIIPrintf(viewer," Control[UMFPACK_RELAXED_AMALGAMATION]: %g\n",lu->Control[UMFPACK_RELAXED_AMALGAMATION]);CHKERRQ(ierr); 258 ierr = PetscViewerASCIIPrintf(viewer," Control[UMFPACK_RELAXED2_AMALGAMATION]: %g\n",lu->Control[UMFPACK_RELAXED2_AMALGAMATION]);CHKERRQ(ierr); 259 ierr = PetscViewerASCIIPrintf(viewer," Control[UMFPACK_RELAXED3_AMALGAMATION]: %g\n",lu->Control[UMFPACK_RELAXED3_AMALGAMATION]);CHKERRQ(ierr); 260 #endif 261 ierr = PetscViewerASCIIPrintf(viewer," Control[UMFPACK_ALLOC_INIT]: %g\n",lu->Control[UMFPACK_ALLOC_INIT]);CHKERRQ(ierr); 262 263 /* Control parameters used by solve */ 264 ierr = PetscViewerASCIIPrintf(viewer," Control[UMFPACK_IRSTEP]: %g\n",lu->Control[UMFPACK_IRSTEP]);CHKERRQ(ierr); 265 266 /* mat ordering */ 267 if(!lu->PetscMatOdering) ierr = PetscViewerASCIIPrintf(viewer," UMFPACK default matrix ordering is used (not the PETSc matrix ordering) \n");CHKERRQ(ierr); 268 269 PetscFunctionReturn(0); 270 } 271 272 EXTERN_C_BEGIN 273 #undef __FUNCT__ 274 #define __FUNCT__ "MatCreate_SeqAIJ_UMFPACK" 275 int MatCreate_SeqAIJ_UMFPACK(Mat A) { 276 int ierr; 277 Mat_SeqAIJ_UMFPACK *lu; 278 279 PetscFunctionBegin; 280 ierr = MatSetType(A,MATSEQAIJ);CHKERRQ(ierr); 281 ierr = MatUseUMFPACK_SeqAIJ(A);CHKERRQ(ierr); 282 283 ierr = PetscNew(Mat_SeqAIJ_UMFPACK,&lu);CHKERRQ(ierr); 284 lu->MatView = A->ops->view; 285 lu->MatAssemblyEnd = A->ops->assemblyend; 286 lu->MatDestroy = A->ops->destroy; 287 lu->CleanUpUMFPACK = PETSC_FALSE; 288 A->spptr = (void*)lu; 289 A->ops->view = MatView_SeqAIJ_UMFPACK; 290 A->ops->assemblyend = MatAssemblyEnd_SeqAIJ_UMFPACK; 291 A->ops->destroy = MatDestroy_SeqAIJ_UMFPACK; 292 PetscFunctionReturn(0); 293 } 294 EXTERN_C_END 295