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 = MatCreateShell(pc->comm,m,N,M,N,(void*)pc,&eis->shell);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 = VecDuplicate(pc->vec,&eis->diag);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