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 PetscTruth 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,PetscTruth 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 EXTERN_C_BEGIN 154 #undef __FUNCT__ 155 #define __FUNCT__ "PCFactorGetMatrix_Factor" 156 PetscErrorCode PETSCKSP_DLLEXPORT PCFactorGetMatrix_Factor(PC pc,Mat *mat) 157 { 158 PC_Factor *ilu = (PC_Factor*)pc->data; 159 160 PetscFunctionBegin; 161 if (!ilu->fact) SETERRQ(((PetscObject)pc)->comm,PETSC_ERR_ORDER,"Matrix not yet factored; call after KSPSetUp() or PCSetUp()"); 162 *mat = ilu->fact; 163 PetscFunctionReturn(0); 164 } 165 EXTERN_C_END 166 167 EXTERN_C_BEGIN 168 #undef __FUNCT__ 169 #define __FUNCT__ "PCFactorSetMatSolverPackage_Factor" 170 PetscErrorCode PETSCKSP_DLLEXPORT PCFactorSetMatSolverPackage_Factor(PC pc,const MatSolverPackage stype) 171 { 172 PetscErrorCode ierr; 173 PC_Factor *lu = (PC_Factor*)pc->data; 174 175 PetscFunctionBegin; 176 if (lu->fact) { 177 const MatSolverPackage ltype; 178 PetscTruth flg; 179 ierr = MatFactorGetSolverPackage(lu->fact,<ype);CHKERRQ(ierr); 180 ierr = PetscStrcmp(stype,ltype,&flg);CHKERRQ(ierr); 181 if (!flg) SETERRQ(((PetscObject)pc)->comm,PETSC_ERR_ARG_WRONGSTATE,"Cannot change solver matrix package after PC has been setup or used"); 182 } else { 183 ierr = MatGetFactor(pc->pmat,stype,lu->factortype,&lu->fact);CHKERRQ(ierr); 184 } 185 ierr = PetscFree(lu->solvertype);CHKERRQ(ierr); 186 ierr = PetscStrallocpy(stype,&lu->solvertype);CHKERRQ(ierr); 187 PetscFunctionReturn(0); 188 } 189 EXTERN_C_END 190 191 EXTERN_C_BEGIN 192 #undef __FUNCT__ 193 #define __FUNCT__ "PCFactorGetMatSolverPackage_Factor" 194 PetscErrorCode PETSCKSP_DLLEXPORT PCFactorGetMatSolverPackage_Factor(PC pc,const MatSolverPackage *stype) 195 { 196 PC_Factor *lu = (PC_Factor*)pc->data; 197 198 PetscFunctionBegin; 199 *stype = lu->solvertype; 200 PetscFunctionReturn(0); 201 } 202 EXTERN_C_END 203 204 EXTERN_C_BEGIN 205 #undef __FUNCT__ 206 #define __FUNCT__ "PCFactorSetColumnPivot_Factor" 207 PetscErrorCode PETSCKSP_DLLEXPORT PCFactorSetColumnPivot_Factor(PC pc,PetscReal dtcol) 208 { 209 PC_Factor *dir = (PC_Factor*)pc->data; 210 211 PetscFunctionBegin; 212 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); 213 dir->info.dtcol = dtcol; 214 PetscFunctionReturn(0); 215 } 216 EXTERN_C_END 217 218 EXTERN_C_BEGIN 219 #undef __FUNCT__ 220 #define __FUNCT__ "PCSetFromOptions_Factor" 221 PetscErrorCode PETSCKSP_DLLEXPORT PCSetFromOptions_Factor(PC pc) 222 { 223 PC_Factor *factor = (PC_Factor*)pc->data; 224 PetscErrorCode ierr; 225 PetscTruth flg = PETSC_FALSE,set; 226 char tname[256], solvertype[64]; 227 PetscFList ordlist; 228 PetscEnum etmp; 229 230 PetscFunctionBegin; 231 if (!MatOrderingRegisterAllCalled) {ierr = MatOrderingRegisterAll(PETSC_NULL);CHKERRQ(ierr);} 232 ierr = PetscOptionsTruth("-pc_factor_in_place","Form factored matrix in the same memory as the matrix","PCFactorSetUseInPlace",flg,&flg,PETSC_NULL);CHKERRQ(ierr); 233 if (flg) { 234 ierr = PCFactorSetUseInPlace(pc);CHKERRQ(ierr); 235 } 236 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); 237 238 ierr = PetscOptionsEnum("-pc_factor_shift_type","Type of shift to add to diagonal","PCFactorSetShiftType", 239 MatFactorShiftTypes,(PetscEnum)(int)((PC_Factor*)factor)->info.shifttype,&etmp,&flg);CHKERRQ(ierr); 240 if (flg) { 241 ierr = PCFactorSetShiftType(pc,(MatFactorShiftType)etmp);CHKERRQ(ierr); 242 } 243 ierr = PetscOptionsReal("-pc_factor_shift_amount","Shift added to diagonal","PCFactorSetShiftAmount",((PC_Factor*)factor)->info.shiftamount,&((PC_Factor*)factor)->info.shiftamount,0);CHKERRQ(ierr); 244 245 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); 246 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); 247 248 flg = ((PC_Factor*)factor)->info.pivotinblocks ? PETSC_TRUE : PETSC_FALSE; 249 ierr = PetscOptionsTruth("-pc_factor_pivot_in_blocks","Pivot inside matrix dense blocks for BAIJ and SBAIJ","PCFactorSetPivotInBlocks",flg,&flg,&set);CHKERRQ(ierr); 250 if (set) { 251 ierr = PCFactorSetPivotInBlocks(pc,flg);CHKERRQ(ierr); 252 } 253 254 flg = PETSC_FALSE; 255 ierr = PetscOptionsTruth("-pc_factor_reuse_fill","Use fill from previous factorization","PCFactorSetReuseFill",flg,&flg,PETSC_NULL);CHKERRQ(ierr); 256 if (flg) { 257 ierr = PCFactorSetReuseFill(pc,PETSC_TRUE);CHKERRQ(ierr); 258 } 259 flg = PETSC_FALSE; 260 ierr = PetscOptionsTruth("-pc_factor_reuse_ordering","Reuse ordering from previous factorization","PCFactorSetReuseOrdering",flg,&flg,PETSC_NULL);CHKERRQ(ierr); 261 if (flg) { 262 ierr = PCFactorSetReuseOrdering(pc,PETSC_TRUE);CHKERRQ(ierr); 263 } 264 265 ierr = MatGetOrderingList(&ordlist);CHKERRQ(ierr); 266 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); 267 if (flg) { 268 ierr = PCFactorSetMatOrderingType(pc,tname);CHKERRQ(ierr); 269 } 270 271 /* maybe should have MatGetSolverTypes(Mat,&list) like the ordering list */ 272 ierr = PetscOptionsString("-pc_factor_mat_solver_package","Specific direct solver to use","MatGetFactor",((PC_Factor*)factor)->solvertype,solvertype,64,&flg);CHKERRQ(ierr); 273 if (flg) { 274 ierr = PCFactorSetMatSolverPackage(pc,solvertype);CHKERRQ(ierr); 275 } 276 PetscFunctionReturn(0); 277 } 278 EXTERN_C_END 279 280 EXTERN_C_BEGIN 281 #undef __FUNCT__ 282 #define __FUNCT__ "PCView_Factor" 283 PetscErrorCode PCView_Factor(PC pc,PetscViewer viewer) 284 { 285 PC_Factor *factor = (PC_Factor*)pc->data; 286 PetscErrorCode ierr; 287 PetscTruth isstring,iascii; 288 289 PetscFunctionBegin; 290 ierr = PetscTypeCompare((PetscObject)viewer,PETSCVIEWERSTRING,&isstring);CHKERRQ(ierr); 291 ierr = PetscTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr); 292 if (iascii) { 293 if (factor->factortype == MAT_FACTOR_ILU || factor->factortype == MAT_FACTOR_ICC){ 294 if (factor->info.dt > 0) { 295 ierr = PetscViewerASCIIPrintf(viewer," drop tolerance %G\n",factor->info.dt);CHKERRQ(ierr); 296 ierr = PetscViewerASCIIPrintf(viewer," max nonzeros per row %D\n",factor->info.dtcount);CHKERRQ(ierr); 297 ierr = PetscViewerASCIIPrintf(viewer," column permutation tolerance %G\n",(PetscInt)factor->info.dtcol);CHKERRQ(ierr); 298 } else if (factor->info.levels == 1) { 299 ierr = PetscViewerASCIIPrintf(viewer," %D level of fill\n",(PetscInt)factor->info.levels);CHKERRQ(ierr); 300 } else { 301 ierr = PetscViewerASCIIPrintf(viewer," %D levels of fill\n",(PetscInt)factor->info.levels);CHKERRQ(ierr); 302 } 303 } 304 305 ierr = PetscViewerASCIIPrintf(viewer," tolerance for zero pivot %G\n",factor->info.zeropivot);CHKERRQ(ierr); 306 if (factor->info.shifttype==(PetscReal)MAT_SHIFT_POSITIVE_DEFINITE) { 307 ierr = PetscViewerASCIIPrintf(viewer," using Manteuffel shift\n");CHKERRQ(ierr); 308 } 309 if (factor->info.shifttype==(PetscReal)MAT_SHIFT_NONZERO) { 310 ierr = PetscViewerASCIIPrintf(viewer," using diagonal shift to prevent zero pivot\n");CHKERRQ(ierr); 311 } 312 if (factor->info.shifttype==(PetscReal)MAT_SHIFT_INBLOCKS) { 313 ierr = PetscViewerASCIIPrintf(viewer," using diagonal shift on blocks to prevent zero pivot\n");CHKERRQ(ierr); 314 } 315 316 ierr = PetscViewerASCIIPrintf(viewer," matrix ordering: %s\n",factor->ordering);CHKERRQ(ierr); 317 318 if (factor->fact) { 319 MatInfo info; 320 ierr = MatGetInfo(factor->fact,MAT_LOCAL,&info);CHKERRQ(ierr); 321 ierr = PetscViewerASCIIPrintf(viewer," factor fill ratio given %G, needed %G\n",info.fill_ratio_given,info.fill_ratio_needed);CHKERRQ(ierr); 322 ierr = PetscViewerASCIIPrintf(viewer," Factored matrix follows:\n");CHKERRQ(ierr); 323 ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr); 324 ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr); 325 ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr); 326 ierr = PetscViewerPushFormat(viewer,PETSC_VIEWER_ASCII_INFO);CHKERRQ(ierr); 327 ierr = MatView(factor->fact,viewer);CHKERRQ(ierr); 328 ierr = PetscViewerPopFormat(viewer);CHKERRQ(ierr); 329 ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr); 330 ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr); 331 ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr); 332 } 333 334 } else if (isstring) { 335 MatFactorType t; 336 ierr = MatGetFactorType(factor->fact,&t);CHKERRQ(ierr); 337 if (t == MAT_FACTOR_ILU || t == MAT_FACTOR_ICC){ 338 ierr = PetscViewerStringSPrintf(viewer," lvls=%D,order=%s",(PetscInt)factor->info.levels,factor->ordering);CHKERRQ(ierr);CHKERRQ(ierr); 339 } 340 } else { 341 SETERRQ1(((PetscObject)pc)->comm,PETSC_ERR_SUP,"Viewer type %s not supported for PC_Factor",((PetscObject)viewer)->type_name); 342 } 343 PetscFunctionReturn(0); 344 } 345 EXTERN_C_END 346