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 PetscErrorCode PCMult_Eisenstat(Mat mat,Vec b,Vec x) 20 { 21 PetscErrorCode 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 PetscErrorCode PCApply_Eisenstat(PC pc,Vec x,Vec y) 35 { 36 PC_Eisenstat *eis = (PC_Eisenstat*)pc->data; 37 PetscErrorCode 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 PetscErrorCode PCPre_Eisenstat(PC pc,KSP ksp,Vec x,Vec b) 48 { 49 PC_Eisenstat *eis = (PC_Eisenstat*)pc->data; 50 PetscTruth nonzero; 51 PetscErrorCode 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 PetscErrorCode PCPost_Eisenstat(PC pc,KSP ksp,Vec x,Vec b) 83 { 84 PC_Eisenstat *eis = (PC_Eisenstat*)pc->data; 85 PetscErrorCode 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 PetscErrorCode PCDestroy_Eisenstat(PC pc) 99 { 100 PC_Eisenstat *eis = (PC_Eisenstat *)pc->data; 101 PetscErrorCode 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 PetscErrorCode PCSetFromOptions_Eisenstat(PC pc) 114 { 115 PC_Eisenstat *eis = (PC_Eisenstat*)pc->data; 116 PetscErrorCode 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 PetscErrorCode PCView_Eisenstat(PC pc,PetscViewer viewer) 133 { 134 PC_Eisenstat *eis = (PC_Eisenstat*)pc->data; 135 PetscErrorCode 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 PetscErrorCode PCSetUp_Eisenstat(PC pc) 156 { 157 PetscErrorCode ierr; 158 int M,N,m,n; 159 PC_Eisenstat *eis = (PC_Eisenstat*)pc->data; 160 161 PetscFunctionBegin; 162 if (!pc->setupcalled) { 163 ierr = MatGetSize(pc->mat,&M,&N);CHKERRQ(ierr); 164 ierr = MatGetLocalSize(pc->mat,&m,&n);CHKERRQ(ierr); 165 ierr = MatCreate(pc->comm,m,N,M,N,&eis->shell);CHKERRQ(ierr); 166 ierr = MatSetType(eis->shell,MATSHELL);CHKERRQ(ierr); 167 ierr = MatShellSetContext(eis->shell,(void*)pc);CHKERRQ(ierr); 168 PetscLogObjectParent(pc,eis->shell); 169 ierr = MatShellSetOperation(eis->shell,MATOP_MULT,(void(*)(void))PCMult_Eisenstat);CHKERRQ(ierr); 170 } 171 if (!eis->usediag) PetscFunctionReturn(0); 172 if (!pc->setupcalled) { 173 ierr = MatGetVecs(pc->pmat,&eis->diag,0);CHKERRQ(ierr); 174 PetscLogObjectParent(pc,eis->diag); 175 } 176 ierr = MatGetDiagonal(pc->pmat,eis->diag);CHKERRQ(ierr); 177 PetscFunctionReturn(0); 178 } 179 180 /* --------------------------------------------------------------------*/ 181 182 EXTERN_C_BEGIN 183 #undef __FUNCT__ 184 #define __FUNCT__ "PCEisenstatSetOmega_Eisenstat" 185 PetscErrorCode PCEisenstatSetOmega_Eisenstat(PC pc,PetscReal omega) 186 { 187 PC_Eisenstat *eis; 188 189 PetscFunctionBegin; 190 if (omega >= 2.0 || omega <= 0.0) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,"Relaxation out of range"); 191 eis = (PC_Eisenstat*)pc->data; 192 eis->omega = omega; 193 PetscFunctionReturn(0); 194 } 195 EXTERN_C_END 196 197 EXTERN_C_BEGIN 198 #undef __FUNCT__ 199 #define __FUNCT__ "PCEisenstatNoDiagonalScaling_Eisenstat" 200 PetscErrorCode PCEisenstatNoDiagonalScaling_Eisenstat(PC pc) 201 { 202 PC_Eisenstat *eis; 203 204 PetscFunctionBegin; 205 eis = (PC_Eisenstat*)pc->data; 206 eis->usediag = PETSC_FALSE; 207 PetscFunctionReturn(0); 208 } 209 EXTERN_C_END 210 211 #undef __FUNCT__ 212 #define __FUNCT__ "PCEisenstatSetOmega" 213 /*@ 214 PCEisenstatSetOmega - Sets the SSOR relaxation coefficient, omega, 215 to use with Eisenstat's trick (where omega = 1.0 by default). 216 217 Collective on PC 218 219 Input Parameters: 220 + pc - the preconditioner context 221 - omega - relaxation coefficient (0 < omega < 2) 222 223 Options Database Key: 224 . -pc_eisenstat_omega <omega> - Sets omega 225 226 Notes: 227 The Eisenstat trick implementation of SSOR requires about 50% of the 228 usual amount of floating point operations used for SSOR + Krylov method; 229 however, the preconditioned problem must be solved with both left 230 and right preconditioning. 231 232 To use SSOR without the Eisenstat trick, employ the PCSOR preconditioner, 233 which can be chosen with the database options 234 $ -pc_type sor -pc_sor_symmetric 235 236 Level: intermediate 237 238 .keywords: PC, Eisenstat, set, SOR, SSOR, relaxation, omega 239 240 .seealso: PCSORSetOmega() 241 @*/ 242 PetscErrorCode PCEisenstatSetOmega(PC pc,PetscReal omega) 243 { 244 PetscErrorCode ierr,(*f)(PC,PetscReal); 245 246 PetscFunctionBegin; 247 PetscValidHeaderSpecific(pc,PC_COOKIE,1); 248 ierr = PetscObjectQueryFunction((PetscObject)pc,"PCEisenstatSetOmega_C",(void (**)(void))&f);CHKERRQ(ierr); 249 if (f) { 250 ierr = (*f)(pc,omega);CHKERRQ(ierr); 251 } 252 PetscFunctionReturn(0); 253 } 254 255 #undef __FUNCT__ 256 #define __FUNCT__ "PCEisenstatNoDiagonalScaling" 257 /*@ 258 PCEisenstatNoDiagonalScaling - Causes the Eisenstat preconditioner 259 not to do additional diagonal preconditioning. For matrices with a constant 260 along the diagonal, this may save a small amount of work. 261 262 Collective on PC 263 264 Input Parameter: 265 . pc - the preconditioner context 266 267 Options Database Key: 268 . -pc_eisenstat_no_diagonal_scaling - Activates PCEisenstatNoDiagonalScaling() 269 270 Level: intermediate 271 272 Note: 273 If you use the KPSSetDiagonalScaling() or -ksp_diagonal_scale option then you will 274 likley want to use this routine since it will save you some unneeded flops. 275 276 .keywords: PC, Eisenstat, use, diagonal, scaling, SSOR 277 278 .seealso: PCEisenstatSetOmega() 279 @*/ 280 PetscErrorCode PCEisenstatNoDiagonalScaling(PC pc) 281 { 282 PetscErrorCode ierr,(*f)(PC); 283 284 PetscFunctionBegin; 285 PetscValidHeaderSpecific(pc,PC_COOKIE,1); 286 ierr = PetscObjectQueryFunction((PetscObject)pc,"PCEisenstatNoDiagonalScaling_C",(void (**)(void))&f);CHKERRQ(ierr); 287 if (f) { 288 ierr = (*f)(pc);CHKERRQ(ierr); 289 } 290 PetscFunctionReturn(0); 291 } 292 293 /* --------------------------------------------------------------------*/ 294 295 /*MC 296 PCEISENSTAT - An implementation of SSOR (symmetric successive over relaxation, symmetric Gauss-Seidel) 297 preconditioning that incorporates Eisenstat's trick to reduce the amount of computation needed. 298 299 Options Database Keys: 300 + -pc_eisenstat_omega <omega> - Sets omega 301 - -pc_eisenstat_no_diagonal_scaling - Activates PCEisenstatNoDiagonalScaling() 302 303 Level: beginner 304 305 Concepts: SOR, preconditioners, Gauss-Seidel, Eisenstat's trick 306 307 Notes: Only implemented for the SeqAIJ matrix format. 308 Not a true parallel SOR, in parallel this implementation corresponds to block 309 Jacobi with SOR on each block. 310 311 .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PC, 312 PCEisenstatNoDiagonalScaling(), PCEisenstatSetOmega(), PCSOR 313 M*/ 314 315 EXTERN_C_BEGIN 316 #undef __FUNCT__ 317 #define __FUNCT__ "PCCreate_Eisenstat" 318 PetscErrorCode PCCreate_Eisenstat(PC pc) 319 { 320 PetscErrorCode ierr; 321 PC_Eisenstat *eis; 322 323 PetscFunctionBegin; 324 ierr = PetscNew(PC_Eisenstat,&eis);CHKERRQ(ierr); 325 PetscLogObjectMemory(pc,sizeof(PC_Eisenstat)); 326 327 pc->ops->apply = PCApply_Eisenstat; 328 pc->ops->presolve = PCPre_Eisenstat; 329 pc->ops->postsolve = PCPost_Eisenstat; 330 pc->ops->applyrichardson = 0; 331 pc->ops->setfromoptions = PCSetFromOptions_Eisenstat; 332 pc->ops->destroy = PCDestroy_Eisenstat; 333 pc->ops->view = PCView_Eisenstat; 334 pc->ops->setup = PCSetUp_Eisenstat; 335 336 pc->data = (void*)eis; 337 eis->omega = 1.0; 338 eis->b = 0; 339 eis->diag = 0; 340 eis->usediag = PETSC_TRUE; 341 342 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCEisenstatSetOmega_C","PCEisenstatSetOmega_Eisenstat", 343 PCEisenstatSetOmega_Eisenstat);CHKERRQ(ierr); 344 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCEisenstatNoDiagonalScaling_C", 345 "PCEisenstatNoDiagonalScaling_Eisenstat", 346 PCEisenstatNoDiagonalScaling_Eisenstat);CHKERRQ(ierr); 347 PetscFunctionReturn(0); 348 } 349 EXTERN_C_END 350