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