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 B; 150 Mat_SeqAIJ *mat=(Mat_SeqAIJ*)A->data; 151 Mat_SeqAIJ_UMFPACK *lu; 152 int ierr,m=A->m,n=A->n,*ai=mat->i,*aj=mat->j,status,*ca; 153 PetscScalar *av=mat->a; 154 155 PetscFunctionBegin; 156 /* Create the factorization matrix F */ 157 ierr = MatCreate(A->comm,PETSC_DECIDE,PETSC_DECIDE,m,n,&B);CHKERRQ(ierr); 158 ierr = MatSetType(B,MATUMFPACK);CHKERRQ(ierr); 159 ierr = MatSeqAIJSetPreallocation(B,0,PETSC_NULL);CHKERRQ(ierr); 160 161 B->ops->lufactornumeric = MatLUFactorNumeric_SeqAIJ_UMFPACK; 162 B->ops->solve = MatSolve_SeqAIJ_UMFPACK; 163 B->factor = FACTOR_LU; 164 B->assembled = PETSC_TRUE; /* required by -sles_view */ 165 166 lu = (Mat_SeqAIJ_UMFPACK*)(B->spptr); 167 168 /* initializations */ 169 /* ------------------------------------------------*/ 170 /* get the default control parameters */ 171 umfpack_di_defaults (lu->Control) ; 172 lu->perm_c = PETSC_NULL; /* use defaul UMFPACK col permutation */ 173 174 ierr = PetscOptionsBegin(A->comm,A->prefix,"UMFPACK Options","Mat");CHKERRQ(ierr); 175 /* Control parameters used by reporting routiones */ 176 ierr = PetscOptionsReal("-mat_umfpack_prl","Control[UMFPACK_PRL]","None",lu->Control[UMFPACK_PRL],&lu->Control[UMFPACK_PRL],PETSC_NULL);CHKERRQ(ierr); 177 178 /* Control parameters for symbolic factorization */ 179 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); 180 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); 181 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); 182 183 /* Control parameters used by numeric factorization */ 184 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); 185 #if !defined(PETSC_HAVE_UMFPACK_41_OR_NEWER) 186 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); 187 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); 188 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); 189 #endif 190 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); 191 192 /* Control parameters used by solve */ 193 ierr = PetscOptionsReal("-mat_umfpack_irstep","Control[UMFPACK_IRSTEP]","None",lu->Control[UMFPACK_IRSTEP],&lu->Control[UMFPACK_IRSTEP],PETSC_NULL);CHKERRQ(ierr); 194 195 /* use Petsc mat ordering (notice size is for the transpose) */ 196 ierr = PetscOptionsHasName(PETSC_NULL,"-pc_lu_mat_ordering_type",&lu->PetscMatOdering);CHKERRQ(ierr); 197 if (lu->PetscMatOdering) { 198 ierr = ISGetIndices(c,&ca);CHKERRQ(ierr); 199 ierr = PetscMalloc(A->m*sizeof(int),&lu->perm_c);CHKERRQ(ierr); 200 ierr = PetscMemcpy(lu->perm_c,ca,A->m*sizeof(int));CHKERRQ(ierr); 201 ierr = ISRestoreIndices(c,&ca);CHKERRQ(ierr); 202 } 203 PetscOptionsEnd(); 204 205 /* print the control parameters */ 206 if( lu->Control[UMFPACK_PRL] > 1 ) umfpack_di_report_control (lu->Control); 207 208 /* symbolic factorization of A' */ 209 /* ---------------------------------------------------------------------- */ 210 #if defined(PETSC_HAVE_UMFPACK_41_OR_NEWER) 211 status = umfpack_di_qsymbolic(n,m,ai,aj,PETSC_NULL,lu->perm_c,&lu->Symbolic,lu->Control,lu->Info) ; 212 #else 213 status = umfpack_di_qsymbolic(n,m,ai,aj,lu->perm_c,&lu->Symbolic,lu->Control,lu->Info) ; 214 #endif 215 if (status < 0){ 216 umfpack_di_report_info(lu->Control, lu->Info) ; 217 umfpack_di_report_status(lu->Control, status) ; 218 SETERRQ(1,"umfpack_di_symbolic failed"); 219 } 220 /* report sumbolic factorization of A' when Control[PRL] > 3 */ 221 (void) umfpack_di_report_symbolic(lu->Symbolic, lu->Control) ; 222 223 lu->flg = DIFFERENT_NONZERO_PATTERN; 224 lu->ai = ai; 225 lu->aj = aj; 226 lu->av = av; 227 lu->CleanUpUMFPACK = PETSC_TRUE; 228 *F = B; 229 PetscFunctionReturn(0); 230 } 231 232 #undef __FUNCT__ 233 #define __FUNCT__ "MatUseUMFPACK_SeqAIJ" 234 int MatUseUMFPACK_SeqAIJ(Mat A) 235 { 236 PetscFunctionBegin; 237 A->ops->lufactorsymbolic = MatLUFactorSymbolic_SeqAIJ_UMFPACK; 238 PetscFunctionReturn(0); 239 } 240 241 /* used by -sles_view */ 242 #undef __FUNCT__ 243 #define __FUNCT__ "MatSeqAIJFactorInfo_UMFPACK" 244 int MatSeqAIJFactorInfo_UMFPACK(Mat A,PetscViewer viewer) 245 { 246 Mat_SeqAIJ_UMFPACK *lu= (Mat_SeqAIJ_UMFPACK*)A->spptr; 247 int ierr; 248 PetscFunctionBegin; 249 /* check if matrix is UMFPACK type */ 250 if (A->ops->solve != MatSolve_SeqAIJ_UMFPACK) PetscFunctionReturn(0); 251 252 ierr = PetscViewerASCIIPrintf(viewer,"UMFPACK run parameters:\n");CHKERRQ(ierr); 253 /* Control parameters used by reporting routiones */ 254 ierr = PetscViewerASCIIPrintf(viewer," Control[UMFPACK_PRL]: %g\n",lu->Control[UMFPACK_PRL]);CHKERRQ(ierr); 255 256 /* Control parameters used by symbolic factorization */ 257 ierr = PetscViewerASCIIPrintf(viewer," Control[UMFPACK_DENSE_COL]: %g\n",lu->Control[UMFPACK_DENSE_COL]);CHKERRQ(ierr); 258 ierr = PetscViewerASCIIPrintf(viewer," Control[UMFPACK_DENSE_ROW]: %g\n",lu->Control[UMFPACK_DENSE_ROW]);CHKERRQ(ierr); 259 ierr = PetscViewerASCIIPrintf(viewer," Control[UMFPACK_BLOCK_SIZE]: %g\n",lu->Control[UMFPACK_BLOCK_SIZE]);CHKERRQ(ierr); 260 261 /* Control parameters used by numeric factorization */ 262 ierr = PetscViewerASCIIPrintf(viewer," Control[UMFPACK_PIVOT_TOLERANCE]: %g\n",lu->Control[UMFPACK_PIVOT_TOLERANCE]);CHKERRQ(ierr); 263 #if !defined(PETSC_HAVE_UMFPACK_41_OR_NEWER) 264 ierr = PetscViewerASCIIPrintf(viewer," Control[UMFPACK_RELAXED_AMALGAMATION]: %g\n",lu->Control[UMFPACK_RELAXED_AMALGAMATION]);CHKERRQ(ierr); 265 ierr = PetscViewerASCIIPrintf(viewer," Control[UMFPACK_RELAXED2_AMALGAMATION]: %g\n",lu->Control[UMFPACK_RELAXED2_AMALGAMATION]);CHKERRQ(ierr); 266 ierr = PetscViewerASCIIPrintf(viewer," Control[UMFPACK_RELAXED3_AMALGAMATION]: %g\n",lu->Control[UMFPACK_RELAXED3_AMALGAMATION]);CHKERRQ(ierr); 267 #endif 268 ierr = PetscViewerASCIIPrintf(viewer," Control[UMFPACK_ALLOC_INIT]: %g\n",lu->Control[UMFPACK_ALLOC_INIT]);CHKERRQ(ierr); 269 270 /* Control parameters used by solve */ 271 ierr = PetscViewerASCIIPrintf(viewer," Control[UMFPACK_IRSTEP]: %g\n",lu->Control[UMFPACK_IRSTEP]);CHKERRQ(ierr); 272 273 /* mat ordering */ 274 if(!lu->PetscMatOdering) ierr = PetscViewerASCIIPrintf(viewer," UMFPACK default matrix ordering is used (not the PETSc matrix ordering) \n");CHKERRQ(ierr); 275 276 PetscFunctionReturn(0); 277 } 278 279 EXTERN_C_BEGIN 280 #undef __FUNCT__ 281 #define __FUNCT__ "MatCreate_SeqAIJ_UMFPACK" 282 int MatCreate_SeqAIJ_UMFPACK(Mat A) { 283 int ierr; 284 Mat_SeqAIJ_UMFPACK *lu; 285 286 PetscFunctionBegin; 287 ierr = MatSetType(A,MATSEQAIJ);CHKERRQ(ierr); 288 ierr = MatUseUMFPACK_SeqAIJ(A);CHKERRQ(ierr); 289 290 ierr = PetscNew(Mat_SeqAIJ_UMFPACK,&lu);CHKERRQ(ierr); 291 lu->MatView = A->ops->view; 292 lu->MatAssemblyEnd = A->ops->assemblyend; 293 lu->MatDestroy = A->ops->destroy; 294 lu->CleanUpUMFPACK = PETSC_FALSE; 295 A->spptr = (void*)lu; 296 A->ops->view = MatView_SeqAIJ_UMFPACK; 297 A->ops->assemblyend = MatAssemblyEnd_SeqAIJ_UMFPACK; 298 A->ops->destroy = MatDestroy_SeqAIJ_UMFPACK; 299 PetscFunctionReturn(0); 300 } 301 EXTERN_C_END 302