1 #define PETSCKSP_DLL 2 3 #include "../src/ksp/pc/impls/factor/factor.h" /*I "petscpc.h" I*/ 4 5 /* ------------------------------------------------------------------------------------------*/ 6 7 EXTERN_C_BEGIN 8 #undef __FUNCT__ 9 #define __FUNCT__ "PCFactorSetZeroPivot_Factor" 10 PetscErrorCode PETSCKSP_DLLEXPORT PCFactorSetZeroPivot_Factor(PC pc,PetscReal z) 11 { 12 PC_Factor *ilu = (PC_Factor*)pc->data; 13 14 PetscFunctionBegin; 15 ilu->info.zeropivot = z; 16 PetscFunctionReturn(0); 17 } 18 EXTERN_C_END 19 20 EXTERN_C_BEGIN 21 #undef __FUNCT__ 22 #define __FUNCT__ "PCFactorSetShiftType_Factor" 23 PetscErrorCode PETSCKSP_DLLEXPORT PCFactorSetShiftType_Factor(PC pc,MatFactorShiftType shifttype) 24 { 25 PC_Factor *dir = (PC_Factor*)pc->data; 26 27 PetscFunctionBegin; 28 if (shifttype == (MatFactorShiftType)PETSC_DECIDE){ 29 dir->info.shifttype = (PetscReal) MAT_SHIFT_NONE; 30 } else { 31 dir->info.shifttype = (PetscReal) shifttype; 32 if (shifttype == MAT_SHIFT_NONZERO && dir->info.shiftamount == 0.0){ 33 dir->info.shiftamount = 1.e-12; /* set default amount if user has not called PCFactorSetShiftAmount() yet */ 34 } 35 } 36 PetscFunctionReturn(0); 37 } 38 39 #undef __FUNCT__ 40 #define __FUNCT__ "PCFactorSetShiftAmount_Factor" 41 PetscErrorCode PETSCKSP_DLLEXPORT PCFactorSetShiftAmount_Factor(PC pc,PetscReal shiftamount) 42 { 43 PC_Factor *dir = (PC_Factor*)pc->data; 44 45 PetscFunctionBegin; 46 if (shiftamount == (PetscReal) PETSC_DECIDE){ 47 dir->info.shiftamount = 1.e-12; 48 } else { 49 dir->info.shiftamount = shiftamount; 50 } 51 PetscFunctionReturn(0); 52 } 53 EXTERN_C_END 54 55 EXTERN_C_BEGIN 56 #undef __FUNCT__ 57 #define __FUNCT__ "PCFactorSetDropTolerance_Factor" 58 PetscErrorCode PETSCKSP_DLLEXPORT PCFactorSetDropTolerance_Factor(PC pc,PetscReal dt,PetscReal dtcol,PetscInt dtcount) 59 { 60 PC_Factor *ilu = (PC_Factor*)pc->data; 61 62 PetscFunctionBegin; 63 if (pc->setupcalled && (!ilu->info.usedt || ((PC_Factor*)ilu)->info.dt != dt || ((PC_Factor*)ilu)->info.dtcol != dtcol || ((PC_Factor*)ilu)->info.dtcount != dtcount)) { 64 SETERRQ(((PetscObject)pc)->comm,PETSC_ERR_ARG_WRONGSTATE,"Cannot change tolerance after use"); 65 } 66 ilu->info.usedt = PETSC_TRUE; 67 ilu->info.dt = dt; 68 ilu->info.dtcol = dtcol; 69 ilu->info.dtcount = dtcount; 70 ilu->info.fill = PETSC_DEFAULT; 71 PetscFunctionReturn(0); 72 } 73 EXTERN_C_END 74 75 EXTERN_C_BEGIN 76 #undef __FUNCT__ 77 #define __FUNCT__ "PCFactorSetFill_Factor" 78 PetscErrorCode PETSCKSP_DLLEXPORT PCFactorSetFill_Factor(PC pc,PetscReal fill) 79 { 80 PC_Factor *dir = (PC_Factor*)pc->data; 81 82 PetscFunctionBegin; 83 dir->info.fill = fill; 84 PetscFunctionReturn(0); 85 } 86 EXTERN_C_END 87 88 EXTERN_C_BEGIN 89 #undef __FUNCT__ 90 #define __FUNCT__ "PCFactorSetMatOrderingType_Factor" 91 PetscErrorCode PETSCKSP_DLLEXPORT PCFactorSetMatOrderingType_Factor(PC pc,const MatOrderingType ordering) 92 { 93 PC_Factor *dir = (PC_Factor*)pc->data; 94 PetscErrorCode ierr; 95 PetscBool flg; 96 97 PetscFunctionBegin; 98 if (!pc->setupcalled) { 99 ierr = PetscFree(dir->ordering);CHKERRQ(ierr); 100 ierr = PetscStrallocpy(ordering,&dir->ordering);CHKERRQ(ierr); 101 } else { 102 ierr = PetscStrcmp(dir->ordering,ordering,&flg);CHKERRQ(ierr); 103 if (!flg) SETERRQ(((PetscObject)pc)->comm,PETSC_ERR_ARG_WRONGSTATE,"Cannot change ordering after use"); 104 } 105 PetscFunctionReturn(0); 106 } 107 EXTERN_C_END 108 109 EXTERN_C_BEGIN 110 #undef __FUNCT__ 111 #define __FUNCT__ "PCFactorSetLevels_Factor" 112 PetscErrorCode PETSCKSP_DLLEXPORT PCFactorSetLevels_Factor(PC pc,PetscInt levels) 113 { 114 PC_Factor *ilu = (PC_Factor*)pc->data; 115 116 PetscFunctionBegin; 117 if (!pc->setupcalled) { 118 ilu->info.levels = levels; 119 } else if (ilu->info.usedt || ilu->info.levels != levels) SETERRQ(((PetscObject)pc)->comm,PETSC_ERR_ARG_WRONGSTATE,"Cannot change levels after use"); 120 PetscFunctionReturn(0); 121 } 122 EXTERN_C_END 123 124 EXTERN_C_BEGIN 125 #undef __FUNCT__ 126 #define __FUNCT__ "PCFactorSetAllowDiagonalFill_Factor" 127 PetscErrorCode PETSCKSP_DLLEXPORT PCFactorSetAllowDiagonalFill_Factor(PC pc) 128 { 129 PC_Factor *dir = (PC_Factor*)pc->data; 130 131 PetscFunctionBegin; 132 dir->info.diagonal_fill = 1; 133 PetscFunctionReturn(0); 134 } 135 EXTERN_C_END 136 137 138 /* ------------------------------------------------------------------------------------------*/ 139 140 EXTERN_C_BEGIN 141 #undef __FUNCT__ 142 #define __FUNCT__ "PCFactorSetPivotInBlocks_Factor" 143 PetscErrorCode PETSCKSP_DLLEXPORT PCFactorSetPivotInBlocks_Factor(PC pc,PetscBool pivot) 144 { 145 PC_Factor *dir = (PC_Factor*)pc->data; 146 147 PetscFunctionBegin; 148 dir->info.pivotinblocks = pivot ? 1.0 : 0.0; 149 PetscFunctionReturn(0); 150 } 151 EXTERN_C_END 152 153 #undef __FUNCT__ 154 #define __FUNCT__ "PCFactorGetMatrix_Factor" 155 PetscErrorCode PETSCKSP_DLLEXPORT PCFactorGetMatrix_Factor(PC pc,Mat *mat) 156 { 157 PC_Factor *ilu = (PC_Factor*)pc->data; 158 159 PetscFunctionBegin; 160 if (!ilu->fact) SETERRQ(((PetscObject)pc)->comm,PETSC_ERR_ORDER,"Matrix not yet factored; call after KSPSetUp() or PCSetUp()"); 161 *mat = ilu->fact; 162 PetscFunctionReturn(0); 163 } 164 165 EXTERN_C_BEGIN 166 #undef __FUNCT__ 167 #define __FUNCT__ "PCFactorSetMatSolverPackage_Factor" 168 PetscErrorCode PETSCKSP_DLLEXPORT PCFactorSetMatSolverPackage_Factor(PC pc,const MatSolverPackage stype) 169 { 170 PetscErrorCode ierr; 171 PC_Factor *lu = (PC_Factor*)pc->data; 172 173 PetscFunctionBegin; 174 if (lu->fact) { 175 const MatSolverPackage ltype; 176 PetscBool flg; 177 ierr = MatFactorGetSolverPackage(lu->fact,<ype);CHKERRQ(ierr); 178 ierr = PetscStrcmp(stype,ltype,&flg);CHKERRQ(ierr); 179 if (!flg) SETERRQ(((PetscObject)pc)->comm,PETSC_ERR_ARG_WRONGSTATE,"Cannot change solver matrix package after PC has been setup or used"); 180 } else { 181 ierr = MatGetFactor(pc->pmat,stype,lu->factortype,&lu->fact);CHKERRQ(ierr); 182 } 183 ierr = PetscFree(lu->solvertype);CHKERRQ(ierr); 184 ierr = PetscStrallocpy(stype,&lu->solvertype);CHKERRQ(ierr); 185 PetscFunctionReturn(0); 186 } 187 EXTERN_C_END 188 189 EXTERN_C_BEGIN 190 #undef __FUNCT__ 191 #define __FUNCT__ "PCFactorGetMatSolverPackage_Factor" 192 PetscErrorCode PETSCKSP_DLLEXPORT PCFactorGetMatSolverPackage_Factor(PC pc,const MatSolverPackage *stype) 193 { 194 PC_Factor *lu = (PC_Factor*)pc->data; 195 196 PetscFunctionBegin; 197 *stype = lu->solvertype; 198 PetscFunctionReturn(0); 199 } 200 EXTERN_C_END 201 202 EXTERN_C_BEGIN 203 #undef __FUNCT__ 204 #define __FUNCT__ "PCFactorSetColumnPivot_Factor" 205 PetscErrorCode PETSCKSP_DLLEXPORT PCFactorSetColumnPivot_Factor(PC pc,PetscReal dtcol) 206 { 207 PC_Factor *dir = (PC_Factor*)pc->data; 208 209 PetscFunctionBegin; 210 if (dtcol < 0.0 || dtcol > 1.0) SETERRQ1(((PetscObject)pc)->comm,PETSC_ERR_ARG_OUTOFRANGE,"Column pivot tolerance is %G must be between 0 and 1",dtcol); 211 dir->info.dtcol = dtcol; 212 PetscFunctionReturn(0); 213 } 214 EXTERN_C_END 215 216 EXTERN_C_BEGIN 217 #undef __FUNCT__ 218 #define __FUNCT__ "PCSetFromOptions_Factor" 219 PetscErrorCode PETSCKSP_DLLEXPORT PCSetFromOptions_Factor(PC pc) 220 { 221 PC_Factor *factor = (PC_Factor*)pc->data; 222 PetscErrorCode ierr; 223 PetscBool flg = PETSC_FALSE,set; 224 char tname[256], solvertype[64]; 225 PetscFList ordlist; 226 PetscEnum etmp; 227 228 PetscFunctionBegin; 229 if (!MatOrderingRegisterAllCalled) {ierr = MatOrderingRegisterAll(PETSC_NULL);CHKERRQ(ierr);} 230 ierr = PetscOptionsBool("-pc_factor_in_place","Form factored matrix in the same memory as the matrix","PCFactorSetUseInPlace",flg,&flg,PETSC_NULL);CHKERRQ(ierr); 231 if (flg) { 232 ierr = PCFactorSetUseInPlace(pc);CHKERRQ(ierr); 233 } 234 ierr = PetscOptionsReal("-pc_factor_fill","Expected non-zeros in factored matrix","PCFactorSetFill",((PC_Factor*)factor)->info.fill,&((PC_Factor*)factor)->info.fill,0);CHKERRQ(ierr); 235 236 ierr = PetscOptionsEnum("-pc_factor_shift_type","Type of shift to add to diagonal","PCFactorSetShiftType", 237 MatFactorShiftTypes,(PetscEnum)(int)((PC_Factor*)factor)->info.shifttype,&etmp,&flg);CHKERRQ(ierr); 238 if (flg) { 239 ierr = PCFactorSetShiftType(pc,(MatFactorShiftType)etmp);CHKERRQ(ierr); 240 } 241 ierr = PetscOptionsReal("-pc_factor_shift_amount","Shift added to diagonal","PCFactorSetShiftAmount",((PC_Factor*)factor)->info.shiftamount,&((PC_Factor*)factor)->info.shiftamount,0);CHKERRQ(ierr); 242 243 ierr = PetscOptionsReal("-pc_factor_zeropivot","Pivot is considered zero if less than","PCFactorSetZeroPivot",((PC_Factor*)factor)->info.zeropivot,&((PC_Factor*)factor)->info.zeropivot,0);CHKERRQ(ierr); 244 ierr = PetscOptionsReal("-pc_factor_column_pivot","Column pivot tolerance (used only for some factorization)","PCFactorSetColumnPivot",((PC_Factor*)factor)->info.dtcol,&((PC_Factor*)factor)->info.dtcol,&flg);CHKERRQ(ierr); 245 246 flg = ((PC_Factor*)factor)->info.pivotinblocks ? PETSC_TRUE : PETSC_FALSE; 247 ierr = PetscOptionsBool("-pc_factor_pivot_in_blocks","Pivot inside matrix dense blocks for BAIJ and SBAIJ","PCFactorSetPivotInBlocks",flg,&flg,&set);CHKERRQ(ierr); 248 if (set) { 249 ierr = PCFactorSetPivotInBlocks(pc,flg);CHKERRQ(ierr); 250 } 251 252 flg = PETSC_FALSE; 253 ierr = PetscOptionsBool("-pc_factor_reuse_fill","Use fill from previous factorization","PCFactorSetReuseFill",flg,&flg,PETSC_NULL);CHKERRQ(ierr); 254 if (flg) { 255 ierr = PCFactorSetReuseFill(pc,PETSC_TRUE);CHKERRQ(ierr); 256 } 257 flg = PETSC_FALSE; 258 ierr = PetscOptionsBool("-pc_factor_reuse_ordering","Reuse ordering from previous factorization","PCFactorSetReuseOrdering",flg,&flg,PETSC_NULL);CHKERRQ(ierr); 259 if (flg) { 260 ierr = PCFactorSetReuseOrdering(pc,PETSC_TRUE);CHKERRQ(ierr); 261 } 262 263 ierr = MatGetOrderingList(&ordlist);CHKERRQ(ierr); 264 ierr = PetscOptionsList("-pc_factor_mat_ordering_type","Reordering to reduce nonzeros in factored matrix","PCFactorSetMatOrderingType",ordlist,((PC_Factor*)factor)->ordering,tname,256,&flg);CHKERRQ(ierr); 265 if (flg) { 266 ierr = PCFactorSetMatOrderingType(pc,tname);CHKERRQ(ierr); 267 } 268 269 /* maybe should have MatGetSolverTypes(Mat,&list) like the ordering list */ 270 ierr = PetscOptionsString("-pc_factor_mat_solver_package","Specific direct solver to use","MatGetFactor",((PC_Factor*)factor)->solvertype,solvertype,64,&flg);CHKERRQ(ierr); 271 if (flg) { 272 ierr = PCFactorSetMatSolverPackage(pc,solvertype);CHKERRQ(ierr); 273 } 274 PetscFunctionReturn(0); 275 } 276 EXTERN_C_END 277 278 EXTERN_C_BEGIN 279 #undef __FUNCT__ 280 #define __FUNCT__ "PCView_Factor" 281 PetscErrorCode PCView_Factor(PC pc,PetscViewer viewer) 282 { 283 PC_Factor *factor = (PC_Factor*)pc->data; 284 PetscErrorCode ierr; 285 PetscBool isstring,iascii; 286 287 PetscFunctionBegin; 288 ierr = PetscTypeCompare((PetscObject)viewer,PETSCVIEWERSTRING,&isstring);CHKERRQ(ierr); 289 ierr = PetscTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr); 290 if (iascii) { 291 if (factor->factortype == MAT_FACTOR_ILU || factor->factortype == MAT_FACTOR_ICC){ 292 if (factor->info.dt > 0) { 293 ierr = PetscViewerASCIIPrintf(viewer," drop tolerance %G\n",factor->info.dt);CHKERRQ(ierr); 294 ierr = PetscViewerASCIIPrintf(viewer," max nonzeros per row %D\n",factor->info.dtcount);CHKERRQ(ierr); 295 ierr = PetscViewerASCIIPrintf(viewer," column permutation tolerance %G\n",(PetscInt)factor->info.dtcol);CHKERRQ(ierr); 296 } else if (factor->info.levels == 1) { 297 ierr = PetscViewerASCIIPrintf(viewer," %D level of fill\n",(PetscInt)factor->info.levels);CHKERRQ(ierr); 298 } else { 299 ierr = PetscViewerASCIIPrintf(viewer," %D levels of fill\n",(PetscInt)factor->info.levels);CHKERRQ(ierr); 300 } 301 } 302 303 ierr = PetscViewerASCIIPrintf(viewer," tolerance for zero pivot %G\n",factor->info.zeropivot);CHKERRQ(ierr); 304 if (factor->info.shifttype==(PetscReal)MAT_SHIFT_POSITIVE_DEFINITE) { 305 ierr = PetscViewerASCIIPrintf(viewer," using Manteuffel shift\n");CHKERRQ(ierr); 306 } 307 if (factor->info.shifttype==(PetscReal)MAT_SHIFT_NONZERO) { 308 ierr = PetscViewerASCIIPrintf(viewer," using diagonal shift to prevent zero pivot\n");CHKERRQ(ierr); 309 } 310 if (factor->info.shifttype==(PetscReal)MAT_SHIFT_INBLOCKS) { 311 ierr = PetscViewerASCIIPrintf(viewer," using diagonal shift on blocks to prevent zero pivot\n");CHKERRQ(ierr); 312 } 313 314 ierr = PetscViewerASCIIPrintf(viewer," matrix ordering: %s\n",factor->ordering);CHKERRQ(ierr); 315 316 if (factor->fact) { 317 MatInfo info; 318 ierr = MatGetInfo(factor->fact,MAT_LOCAL,&info);CHKERRQ(ierr); 319 ierr = PetscViewerASCIIPrintf(viewer," factor fill ratio given %G, needed %G\n",info.fill_ratio_given,info.fill_ratio_needed);CHKERRQ(ierr); 320 ierr = PetscViewerASCIIPrintf(viewer," Factored matrix follows:\n");CHKERRQ(ierr); 321 ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr); 322 ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr); 323 ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr); 324 ierr = PetscViewerPushFormat(viewer,PETSC_VIEWER_ASCII_INFO);CHKERRQ(ierr); 325 ierr = MatView(factor->fact,viewer);CHKERRQ(ierr); 326 ierr = PetscViewerPopFormat(viewer);CHKERRQ(ierr); 327 ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr); 328 ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr); 329 ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr); 330 } 331 332 } else if (isstring) { 333 MatFactorType t; 334 ierr = MatGetFactorType(factor->fact,&t);CHKERRQ(ierr); 335 if (t == MAT_FACTOR_ILU || t == MAT_FACTOR_ICC){ 336 ierr = PetscViewerStringSPrintf(viewer," lvls=%D,order=%s",(PetscInt)factor->info.levels,factor->ordering);CHKERRQ(ierr);CHKERRQ(ierr); 337 } 338 } else { 339 SETERRQ1(((PetscObject)pc)->comm,PETSC_ERR_SUP,"Viewer type %s not supported for PC_Factor",((PetscObject)viewer)->type_name); 340 } 341 PetscFunctionReturn(0); 342 } 343 EXTERN_C_END 344