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 EXTERN_C_BEGIN 339 #undef __FUNCT__ 340 #define __FUNCT__ "MatCreate_SeqAIJ_UMFPACK" 341 int MatCreate_SeqAIJ_UMFPACK(Mat A) { 342 int ierr; 343 344 PetscFunctionBegin; 345 ierr = MatSetType(A,MATSEQAIJ);CHKERRQ(ierr); 346 ierr = MatConvert_SeqAIJ_UMFPACK(A,MATUMFPACK,&A);CHKERRQ(ierr); 347 PetscFunctionReturn(0); 348 } 349 EXTERN_C_END 350