1 /* 2 Defines a Eisenstat trick SSOR preconditioner. This uses about 3 %50 of the usual amount of floating point ops used for SSOR + Krylov 4 method. But it requires actually solving the preconditioned problem 5 with both left and right preconditioning. 6 */ 7 #include "src/ksp/pc/pcimpl.h" /*I "petscpc.h" I*/ 8 9 typedef struct { 10 Mat shell,A; 11 Vec b,diag; /* temporary storage for true right hand side */ 12 PetscReal omega; 13 PetscTruth usediag; /* indicates preconditioner should include diagonal scaling*/ 14 } PC_Eisenstat; 15 16 17 #undef __FUNCT__ 18 #define __FUNCT__ "PCMult_Eisenstat" 19 static int PCMult_Eisenstat(Mat mat,Vec b,Vec x) 20 { 21 int ierr; 22 PC pc; 23 PC_Eisenstat *eis; 24 25 PetscFunctionBegin; 26 ierr = MatShellGetContext(mat,(void **)&pc);CHKERRQ(ierr); 27 eis = (PC_Eisenstat*)pc->data; 28 ierr = MatRelax(eis->A,b,eis->omega,SOR_EISENSTAT,0.0,1,1,x);CHKERRQ(ierr); 29 PetscFunctionReturn(0); 30 } 31 32 #undef __FUNCT__ 33 #define __FUNCT__ "PCApply_Eisenstat" 34 static int PCApply_Eisenstat(PC pc,Vec x,Vec y) 35 { 36 PC_Eisenstat *eis = (PC_Eisenstat*)pc->data; 37 int ierr; 38 39 PetscFunctionBegin; 40 if (eis->usediag) {ierr = VecPointwiseMult(x,eis->diag,y);CHKERRQ(ierr);} 41 else {ierr = VecCopy(x,y);CHKERRQ(ierr);} 42 PetscFunctionReturn(0); 43 } 44 45 #undef __FUNCT__ 46 #define __FUNCT__ "PCPre_Eisenstat" 47 static int PCPre_Eisenstat(PC pc,KSP ksp,Vec x,Vec b) 48 { 49 PC_Eisenstat *eis = (PC_Eisenstat*)pc->data; 50 PetscTruth nonzero; 51 int ierr; 52 53 PetscFunctionBegin; 54 if (pc->mat != pc->pmat) SETERRQ(PETSC_ERR_SUP,"Cannot have different mat and pmat"); 55 56 /* swap shell matrix and true matrix */ 57 eis->A = pc->mat; 58 pc->mat = eis->shell; 59 60 if (!eis->b) { 61 ierr = VecDuplicate(b,&eis->b);CHKERRQ(ierr); 62 PetscLogObjectParent(pc,eis->b); 63 } 64 65 /* save true b, other option is to swap pointers */ 66 ierr = VecCopy(b,eis->b);CHKERRQ(ierr); 67 68 /* if nonzero initial guess, modify x */ 69 ierr = KSPGetInitialGuessNonzero(ksp,&nonzero);CHKERRQ(ierr); 70 if (nonzero) { 71 ierr = MatRelax(eis->A,x,eis->omega,SOR_APPLY_UPPER,0.0,1,1,x);CHKERRQ(ierr); 72 } 73 74 /* modify b by (L + D)^{-1} */ 75 ierr = MatRelax(eis->A,b,eis->omega,(MatSORType)(SOR_ZERO_INITIAL_GUESS | 76 SOR_FORWARD_SWEEP),0.0,1,1,b);CHKERRQ(ierr); 77 PetscFunctionReturn(0); 78 } 79 80 #undef __FUNCT__ 81 #define __FUNCT__ "PCPost_Eisenstat" 82 static int PCPost_Eisenstat(PC pc,KSP ksp,Vec x,Vec b) 83 { 84 PC_Eisenstat *eis = (PC_Eisenstat*)pc->data; 85 int ierr; 86 87 PetscFunctionBegin; 88 ierr = MatRelax(eis->A,x,eis->omega,(MatSORType)(SOR_ZERO_INITIAL_GUESS | 89 SOR_BACKWARD_SWEEP),0.0,1,1,x);CHKERRQ(ierr); 90 pc->mat = eis->A; 91 /* get back true b */ 92 ierr = VecCopy(eis->b,b);CHKERRQ(ierr); 93 PetscFunctionReturn(0); 94 } 95 96 #undef __FUNCT__ 97 #define __FUNCT__ "PCDestroy_Eisenstat" 98 static int PCDestroy_Eisenstat(PC pc) 99 { 100 PC_Eisenstat *eis = (PC_Eisenstat *)pc->data; 101 int ierr; 102 103 PetscFunctionBegin; 104 if (eis->b) {ierr = VecDestroy(eis->b);CHKERRQ(ierr);} 105 if (eis->shell) {ierr = MatDestroy(eis->shell);CHKERRQ(ierr);} 106 if (eis->diag) {ierr = VecDestroy(eis->diag);CHKERRQ(ierr);} 107 ierr = PetscFree(eis);CHKERRQ(ierr); 108 PetscFunctionReturn(0); 109 } 110 111 #undef __FUNCT__ 112 #define __FUNCT__ "PCSetFromOptions_Eisenstat" 113 static int PCSetFromOptions_Eisenstat(PC pc) 114 { 115 PC_Eisenstat *eis = (PC_Eisenstat*)pc->data; 116 int ierr; 117 PetscTruth flg; 118 119 PetscFunctionBegin; 120 ierr = PetscOptionsHead("Eisenstat SSOR options");CHKERRQ(ierr); 121 ierr = PetscOptionsReal("-pc_eisenstat_omega","Relaxation factor 0 < omega < 2","PCEisenstatSetOmega",eis->omega,&eis->omega,0);CHKERRQ(ierr); 122 ierr = PetscOptionsName("-pc_eisenstat_no_diagonal_scaling","Do not use standard diagonal scaling","PCEisenstatNoDiagonalScaling",&flg);CHKERRQ(ierr); 123 if (flg) { 124 ierr = PCEisenstatNoDiagonalScaling(pc);CHKERRQ(ierr); 125 } 126 ierr = PetscOptionsTail();CHKERRQ(ierr); 127 PetscFunctionReturn(0); 128 } 129 130 #undef __FUNCT__ 131 #define __FUNCT__ "PCView_Eisenstat" 132 static int PCView_Eisenstat(PC pc,PetscViewer viewer) 133 { 134 PC_Eisenstat *eis = (PC_Eisenstat*)pc->data; 135 int ierr; 136 PetscTruth iascii; 137 138 PetscFunctionBegin; 139 ierr = PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_ASCII,&iascii);CHKERRQ(ierr); 140 if (iascii) { 141 ierr = PetscViewerASCIIPrintf(viewer,"Eisenstat: omega = %g\n",eis->omega);CHKERRQ(ierr); 142 if (eis->usediag) { 143 ierr = PetscViewerASCIIPrintf(viewer,"Eisenstat: Using diagonal scaling (default)\n");CHKERRQ(ierr); 144 } else { 145 ierr = PetscViewerASCIIPrintf(viewer,"Eisenstat: Not using diagonal scaling\n");CHKERRQ(ierr); 146 } 147 } else { 148 SETERRQ1(1,"Viewer type %s not supported for Eisenstat PC",((PetscObject)viewer)->type_name); 149 } 150 PetscFunctionReturn(0); 151 } 152 153 #undef __FUNCT__ 154 #define __FUNCT__ "PCSetUp_Eisenstat" 155 static int PCSetUp_Eisenstat(PC pc) 156 { 157 int ierr,M,N,m,n; 158 PC_Eisenstat *eis = (PC_Eisenstat*)pc->data; 159 160 PetscFunctionBegin; 161 if (!pc->setupcalled) { 162 ierr = MatGetSize(pc->mat,&M,&N);CHKERRQ(ierr); 163 ierr = MatGetLocalSize(pc->mat,&m,&n);CHKERRQ(ierr); 164 ierr = MatCreate(pc->comm,m,N,M,N,&eis->shell);CHKERRQ(ierr); 165 ierr = MatSetType(eis->shell,MATSHELL);CHKERRQ(ierr); 166 ierr = MatShellSetContext(eis->shell,(void*)pc);CHKERRQ(ierr); 167 PetscLogObjectParent(pc,eis->shell); 168 ierr = MatShellSetOperation(eis->shell,MATOP_MULT,(void(*)(void))PCMult_Eisenstat);CHKERRQ(ierr); 169 } 170 if (!eis->usediag) PetscFunctionReturn(0); 171 if (!pc->setupcalled) { 172 ierr = MatGetVecs(pc->pmat,&eis->diag,0);CHKERRQ(ierr); 173 PetscLogObjectParent(pc,eis->diag); 174 } 175 ierr = MatGetDiagonal(pc->pmat,eis->diag);CHKERRQ(ierr); 176 PetscFunctionReturn(0); 177 } 178 179 /* --------------------------------------------------------------------*/ 180 181 EXTERN_C_BEGIN 182 #undef __FUNCT__ 183 #define __FUNCT__ "PCEisenstatSetOmega_Eisenstat" 184 int PCEisenstatSetOmega_Eisenstat(PC pc,PetscReal omega) 185 { 186 PC_Eisenstat *eis; 187 188 PetscFunctionBegin; 189 if (omega >= 2.0 || omega <= 0.0) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,"Relaxation out of range"); 190 eis = (PC_Eisenstat*)pc->data; 191 eis->omega = omega; 192 PetscFunctionReturn(0); 193 } 194 EXTERN_C_END 195 196 EXTERN_C_BEGIN 197 #undef __FUNCT__ 198 #define __FUNCT__ "PCEisenstatNoDiagonalScaling_Eisenstat" 199 int PCEisenstatNoDiagonalScaling_Eisenstat(PC pc) 200 { 201 PC_Eisenstat *eis; 202 203 PetscFunctionBegin; 204 eis = (PC_Eisenstat*)pc->data; 205 eis->usediag = PETSC_FALSE; 206 PetscFunctionReturn(0); 207 } 208 EXTERN_C_END 209 210 #undef __FUNCT__ 211 #define __FUNCT__ "PCEisenstatSetOmega" 212 /*@ 213 PCEisenstatSetOmega - Sets the SSOR relaxation coefficient, omega, 214 to use with Eisenstat's trick (where omega = 1.0 by default). 215 216 Collective on PC 217 218 Input Parameters: 219 + pc - the preconditioner context 220 - omega - relaxation coefficient (0 < omega < 2) 221 222 Options Database Key: 223 . -pc_eisenstat_omega <omega> - Sets omega 224 225 Notes: 226 The Eisenstat trick implementation of SSOR requires about 50% of the 227 usual amount of floating point operations used for SSOR + Krylov method; 228 however, the preconditioned problem must be solved with both left 229 and right preconditioning. 230 231 To use SSOR without the Eisenstat trick, employ the PCSOR preconditioner, 232 which can be chosen with the database options 233 $ -pc_type sor -pc_sor_symmetric 234 235 Level: intermediate 236 237 .keywords: PC, Eisenstat, set, SOR, SSOR, relaxation, omega 238 239 .seealso: PCSORSetOmega() 240 @*/ 241 int PCEisenstatSetOmega(PC pc,PetscReal omega) 242 { 243 int ierr,(*f)(PC,PetscReal); 244 245 PetscFunctionBegin; 246 PetscValidHeaderSpecific(pc,PC_COOKIE,1); 247 ierr = PetscObjectQueryFunction((PetscObject)pc,"PCEisenstatSetOmega_C",(void (**)(void))&f);CHKERRQ(ierr); 248 if (f) { 249 ierr = (*f)(pc,omega);CHKERRQ(ierr); 250 } 251 PetscFunctionReturn(0); 252 } 253 254 #undef __FUNCT__ 255 #define __FUNCT__ "PCEisenstatNoDiagonalScaling" 256 /*@ 257 PCEisenstatNoDiagonalScaling - Causes the Eisenstat preconditioner 258 not to do additional diagonal preconditioning. For matrices with a constant 259 along the diagonal, this may save a small amount of work. 260 261 Collective on PC 262 263 Input Parameter: 264 . pc - the preconditioner context 265 266 Options Database Key: 267 . -pc_eisenstat_no_diagonal_scaling - Activates PCEisenstatNoDiagonalScaling() 268 269 Level: intermediate 270 271 Note: 272 If you use the KPSSetDiagonalScaling() or -ksp_diagonal_scale option then you will 273 likley want to use this routine since it will save you some unneeded flops. 274 275 .keywords: PC, Eisenstat, use, diagonal, scaling, SSOR 276 277 .seealso: PCEisenstatSetOmega() 278 @*/ 279 int PCEisenstatNoDiagonalScaling(PC pc) 280 { 281 int ierr,(*f)(PC); 282 283 PetscFunctionBegin; 284 PetscValidHeaderSpecific(pc,PC_COOKIE,1); 285 ierr = PetscObjectQueryFunction((PetscObject)pc,"PCEisenstatNoDiagonalScaling_C",(void (**)(void))&f);CHKERRQ(ierr); 286 if (f) { 287 ierr = (*f)(pc);CHKERRQ(ierr); 288 } 289 PetscFunctionReturn(0); 290 } 291 292 /* --------------------------------------------------------------------*/ 293 294 /*MC 295 PCEISENSTAT - An implementation of SSOR (symmetric successive over relaxation, symmetric Gauss-Seidel) 296 preconditioning that incorporates Eisenstat's trick to reduce the amount of computation needed. 297 298 Options Database Keys: 299 + -pc_eisenstat_omega <omega> - Sets omega 300 - -pc_eisenstat_no_diagonal_scaling - Activates PCEisenstatNoDiagonalScaling() 301 302 Level: beginner 303 304 Concepts: SOR, preconditioners, Gauss-Seidel, Eisenstat's trick 305 306 Notes: Only implemented for the SeqAIJ matrix format. 307 Not a true parallel SOR, in parallel this implementation corresponds to block 308 Jacobi with SOR on each block. 309 310 .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PC, 311 PCEisenstatNoDiagonalScaling(), PCEisenstatSetOmega(), PCSOR 312 M*/ 313 314 EXTERN_C_BEGIN 315 #undef __FUNCT__ 316 #define __FUNCT__ "PCCreate_Eisenstat" 317 int PCCreate_Eisenstat(PC pc) 318 { 319 int ierr; 320 PC_Eisenstat *eis; 321 322 PetscFunctionBegin; 323 ierr = PetscNew(PC_Eisenstat,&eis);CHKERRQ(ierr); 324 PetscLogObjectMemory(pc,sizeof(PC_Eisenstat)); 325 326 pc->ops->apply = PCApply_Eisenstat; 327 pc->ops->presolve = PCPre_Eisenstat; 328 pc->ops->postsolve = PCPost_Eisenstat; 329 pc->ops->applyrichardson = 0; 330 pc->ops->setfromoptions = PCSetFromOptions_Eisenstat; 331 pc->ops->destroy = PCDestroy_Eisenstat; 332 pc->ops->view = PCView_Eisenstat; 333 pc->ops->setup = PCSetUp_Eisenstat; 334 335 pc->data = (void*)eis; 336 eis->omega = 1.0; 337 eis->b = 0; 338 eis->diag = 0; 339 eis->usediag = PETSC_TRUE; 340 341 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCEisenstatSetOmega_C","PCEisenstatSetOmega_Eisenstat", 342 PCEisenstatSetOmega_Eisenstat);CHKERRQ(ierr); 343 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCEisenstatNoDiagonalScaling_C", 344 "PCEisenstatNoDiagonalScaling_Eisenstat", 345 PCEisenstatNoDiagonalScaling_Eisenstat);CHKERRQ(ierr); 346 PetscFunctionReturn(0); 347 } 348 EXTERN_C_END 349