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