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 (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 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); 185 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); 186 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); 187 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); 188 189 /* Control parameters used by solve */ 190 ierr = PetscOptionsReal("-mat_umfpack_irstep","Control[UMFPACK_IRSTEP]","None",lu->Control[UMFPACK_IRSTEP],&lu->Control[UMFPACK_IRSTEP],PETSC_NULL);CHKERRQ(ierr); 191 192 /* use Petsc mat ordering (notice size is for the transpose) */ 193 ierr = PetscOptionsHasName(PETSC_NULL,"-pc_lu_mat_ordering_type",&lu->PetscMatOdering);CHKERRQ(ierr); 194 if (lu->PetscMatOdering) { 195 ierr = ISGetIndices(c,&ca);CHKERRQ(ierr); 196 ierr = PetscMalloc(A->m*sizeof(int),&lu->perm_c);CHKERRQ(ierr); 197 ierr = PetscMemcpy(lu->perm_c,ca,A->m*sizeof(int));CHKERRQ(ierr); 198 ierr = ISRestoreIndices(c,&ca);CHKERRQ(ierr); 199 } 200 PetscOptionsEnd(); 201 202 /* print the control parameters */ 203 if( lu->Control[UMFPACK_PRL] > 1 ) umfpack_di_report_control (lu->Control); 204 205 /* symbolic factorization of A' */ 206 /* ---------------------------------------------------------------------- */ 207 status = umfpack_di_qsymbolic(n,m,ai,aj,lu->perm_c,&lu->Symbolic,lu->Control,lu->Info) ; 208 if (status < 0){ 209 umfpack_di_report_info(lu->Control, lu->Info) ; 210 umfpack_di_report_status(lu->Control, status) ; 211 SETERRQ(1,"umfpack_di_symbolic failed"); 212 } 213 /* report sumbolic factorization of A' when Control[PRL] > 3 */ 214 (void) umfpack_di_report_symbolic(lu->Symbolic, lu->Control) ; 215 216 lu->flg = DIFFERENT_NONZERO_PATTERN; 217 lu->ai = ai; 218 lu->aj = aj; 219 lu->av = av; 220 PetscFunctionReturn(0); 221 } 222 223 #undef __FUNCT__ 224 #define __FUNCT__ "MatUseUMFPACK_SeqAIJ" 225 int MatUseUMFPACK_SeqAIJ(Mat A) 226 { 227 PetscFunctionBegin; 228 A->ops->lufactorsymbolic = MatLUFactorSymbolic_SeqAIJ_UMFPACK; 229 PetscFunctionReturn(0); 230 } 231 232 /* used by -sles_view */ 233 #undef __FUNCT__ 234 #define __FUNCT__ "MatSeqAIJFactorInfo_UMFPACK" 235 int MatSeqAIJFactorInfo_UMFPACK(Mat A,PetscViewer viewer) 236 { 237 Mat_SeqAIJ_UMFPACK *lu= (Mat_SeqAIJ_UMFPACK*)A->spptr; 238 int ierr; 239 PetscFunctionBegin; 240 /* check if matrix is UMFPACK type */ 241 if (A->ops->solve != MatSolve_SeqAIJ_UMFPACK) PetscFunctionReturn(0); 242 243 ierr = PetscViewerASCIIPrintf(viewer,"UMFPACK run parameters:\n");CHKERRQ(ierr); 244 /* Control parameters used by reporting routiones */ 245 ierr = PetscViewerASCIIPrintf(viewer," Control[UMFPACK_PRL]: %g\n",lu->Control[UMFPACK_PRL]);CHKERRQ(ierr); 246 247 /* Control parameters used by symbolic factorization */ 248 ierr = PetscViewerASCIIPrintf(viewer," Control[UMFPACK_DENSE_COL]: %g\n",lu->Control[UMFPACK_DENSE_COL]);CHKERRQ(ierr); 249 ierr = PetscViewerASCIIPrintf(viewer," Control[UMFPACK_DENSE_ROW]: %g\n",lu->Control[UMFPACK_DENSE_ROW]);CHKERRQ(ierr); 250 ierr = PetscViewerASCIIPrintf(viewer," Control[UMFPACK_BLOCK_SIZE]: %g\n",lu->Control[UMFPACK_BLOCK_SIZE]);CHKERRQ(ierr); 251 252 /* Control parameters used by numeric factorization */ 253 ierr = PetscViewerASCIIPrintf(viewer," Control[UMFPACK_PIVOT_TOLERANCE]: %g\n",lu->Control[UMFPACK_PIVOT_TOLERANCE]);CHKERRQ(ierr); 254 ierr = PetscViewerASCIIPrintf(viewer," Control[UMFPACK_RELAXED_AMALGAMATION]: %g\n",lu->Control[UMFPACK_RELAXED_AMALGAMATION]);CHKERRQ(ierr); 255 ierr = PetscViewerASCIIPrintf(viewer," Control[UMFPACK_RELAXED2_AMALGAMATION]: %g\n",lu->Control[UMFPACK_RELAXED2_AMALGAMATION]);CHKERRQ(ierr); 256 ierr = PetscViewerASCIIPrintf(viewer," Control[UMFPACK_RELAXED3_AMALGAMATION]: %g\n",lu->Control[UMFPACK_RELAXED3_AMALGAMATION]);CHKERRQ(ierr); 257 ierr = PetscViewerASCIIPrintf(viewer," Control[UMFPACK_ALLOC_INIT]: %g\n",lu->Control[UMFPACK_ALLOC_INIT]);CHKERRQ(ierr); 258 259 /* Control parameters used by solve */ 260 ierr = PetscViewerASCIIPrintf(viewer," Control[UMFPACK_IRSTEP]: %g\n",lu->Control[UMFPACK_IRSTEP]);CHKERRQ(ierr); 261 262 /* mat ordering */ 263 if(!lu->PetscMatOdering) ierr = PetscViewerASCIIPrintf(viewer," UMFPACK default matrix ordering is used (not the PETSc matrix ordering) \n");CHKERRQ(ierr); 264 265 PetscFunctionReturn(0); 266 } 267 268 EXTERN_C_BEGIN 269 #undef __FUNCT__ 270 #define __FUNCT__ "MatCreate_SeqAIJ_UMFPACK" 271 int MatCreate_SeqAIJ_UMFPACK(Mat A) { 272 int ierr; 273 Mat_SeqAIJ_UMFPACK *lu; 274 275 PetscFunctionBegin; 276 ierr = MatSetType(A,MATSEQAIJ);CHKERRQ(ierr); 277 ierr = MatUseUMFPACK_SeqAIJ(A);CHKERRQ(ierr); 278 279 ierr = PetscNew(Mat_SeqAIJ_UMFPACK,&lu);CHKERRQ(ierr); 280 lu->MatView = A->ops->view; 281 lu->MatAssemblyEnd = A->ops->assemblyend; 282 lu->MatDestroy = A->ops->destroy; 283 lu->CleanUpUMFPACK = PETSC_FALSE; 284 A->spptr = (void*)lu; 285 A->ops->view = MatView_SeqAIJ_UMFPACK; 286 A->ops->assemblyend = MatAssemblyEnd_SeqAIJ_UMFPACK; 287 A->ops->destroy = MatDestroy_SeqAIJ_UMFPACK; 288 PetscFunctionReturn(0); 289 } 290 EXTERN_C_END 291