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