1: #include <petsc-private/snesimpl.h>             /*I   "petscsnes.h"   I*/

  3: typedef struct {
  4:   PetscInt  sweeps;
  5: } SNES_GS;

  9: PetscErrorCode SNESReset_GS(SNES snes)
 10: {
 12:   return(0);
 13: }

 17: PetscErrorCode SNESDestroy_GS(SNES snes)
 18: {

 22:   SNESReset_GS(snes);
 23:   PetscFree(snes->data);
 24:   return(0);
 25: }

 29: PetscErrorCode SNESSetUp_GS(SNES snes)
 30: {
 32:   return(0);
 33: }

 37: PetscErrorCode SNESSetFromOptions_GS(SNES snes)
 38: {

 42:   PetscOptionsHead("SNES GS options");
 43:   PetscOptionsTail();
 44:   return(0);
 45: }

 49: PetscErrorCode SNESView_GS(SNES snes, PetscViewer viewer)
 50: {
 52:   return(0);
 53: }

 57: PetscErrorCode SNESSolve_GS(SNES snes)
 58: {
 59:   Vec            F;
 60:   Vec            X;
 61:   Vec            B;
 62:   PetscInt       i;
 63:   PetscReal      fnorm;
 65:   SNESNormType   normtype;


 69:   X = snes->vec_sol;
 70:   F = snes->vec_func;
 71:   B = snes->vec_rhs;

 73:   PetscObjectTakeAccess(snes);
 74:   snes->iter = 0;
 75:   snes->norm = 0.;
 76:   PetscObjectGrantAccess(snes);
 77:   snes->reason = SNES_CONVERGED_ITERATING;

 79:   SNESGetNormType(snes, &normtype);
 80:   if (normtype == SNES_NORM_FUNCTION || normtype == SNES_NORM_INITIAL_ONLY || normtype == SNES_NORM_INITIAL_FINAL_ONLY) {
 81:     /* compute the initial function and preconditioned update delX */
 82:     if (!snes->vec_func_init_set) {
 83:       SNESComputeFunction(snes,X,F);
 84:       if (snes->domainerror) {
 85:         snes->reason = SNES_DIVERGED_FUNCTION_DOMAIN;
 86:         return(0);
 87:       }
 88:     } else {
 89:       snes->vec_func_init_set = PETSC_FALSE;
 90:     }

 92:     /* convergence test */
 93:     if (!snes->norm_init_set) {
 94:       VecNorm(F, NORM_2, &fnorm); /* fnorm <- ||F||  */
 95:       if (PetscIsInfOrNanReal(fnorm)) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FP,"Infinite or not-a-number generated in norm");
 96:     } else {
 97:       fnorm = snes->norm_init;
 98:       snes->norm_init_set = PETSC_FALSE;
 99:     }
100:     PetscObjectTakeAccess(snes);
101:     snes->iter = 0;
102:     snes->norm = fnorm;
103:     PetscObjectGrantAccess(snes);
104:     SNESLogConvHistory(snes,snes->norm,0);
105:     SNESMonitor(snes,0,snes->norm);

107:     /* set parameter for default relative tolerance convergence test */
108:     snes->ttol = fnorm*snes->rtol;

110:     /* test convergence */
111:     (*snes->ops->converged)(snes,0,0.0,0.0,fnorm,&snes->reason,snes->cnvP);
112:     if (snes->reason) return(0);
113:   } else {
114:     PetscObjectGrantAccess(snes);
115:     SNESLogConvHistory(snes,snes->norm,0);
116:     SNESMonitor(snes,0,snes->norm);
117:   }


120:   /* Call general purpose update function */
121:   if (snes->ops->update) {
122:     (*snes->ops->update)(snes, snes->iter);
123:   }

125:   for (i = 0; i < snes->max_its; i++) {
126:     SNESComputeGS(snes, B, X);
127:     /* only compute norms if requested or about to exit due to maximum iterations */
128:     if (normtype == SNES_NORM_FUNCTION || ((i == snes->max_its - 1) && (normtype == SNES_NORM_INITIAL_FINAL_ONLY || normtype == SNES_NORM_FINAL_ONLY))) {
129:       SNESComputeFunction(snes,X,F);
130:       if (snes->domainerror) {
131:         snes->reason = SNES_DIVERGED_FUNCTION_DOMAIN;
132:         return(0);
133:       }
134:       VecNorm(F, NORM_2, &fnorm); /* fnorm <- ||F||  */
135:       if (PetscIsInfOrNanReal(fnorm)) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FP,"Infinite or not-a-number generated in norm");
136:     }
137:     /* Monitor convergence */
138:     PetscObjectTakeAccess(snes);
139:     snes->iter = i+1;
140:     snes->norm = fnorm;
141:     PetscObjectGrantAccess(snes);
142:     SNESLogConvHistory(snes,snes->norm,0);
143:     SNESMonitor(snes,snes->iter,snes->norm);
144:     /* Test for convergence */
145:     if (normtype == SNES_NORM_FUNCTION) (*snes->ops->converged)(snes,snes->iter,0.0,0.0,fnorm,&snes->reason,snes->cnvP);
146:     if (snes->reason) return(0);
147:     /* Call general purpose update function */
148:     if (snes->ops->update) {
149:       (*snes->ops->update)(snes, snes->iter);
150:     }
151:   }
152:   if (normtype == SNES_NORM_FUNCTION) {
153:     if (i == snes->max_its) {
154:       PetscInfo1(snes,"Maximum number of iterations has been reached: %D\n",snes->max_its);
155:       if(!snes->reason) snes->reason = SNES_DIVERGED_MAX_IT;
156:     }
157:   } else {
158:     if (!snes->reason) snes->reason = SNES_CONVERGED_ITS; /* GS is meant to be used as a preconditioner */
159:   }
160:   return(0);
161: }

163: /*MC
164:   SNESGS - Just calls the user-provided solution routine provided with SNESSetGS()

166:    Level: advanced

168:   Notes:
169:   the Gauss-Seidel smoother is inherited through composition.  If a solver has been created with SNESGetPC(), it will have
170:   its parent's Gauss-Seidel routine associated with it.

172: .seealso: SNESCreate(), SNES, SNESSetType(), SNESSetGS(), SNESType (for list of available types)
173: M*/

175: EXTERN_C_BEGIN
178: PetscErrorCode SNESCreate_GS(SNES snes)
179: {
180:   SNES_GS        *gs;

184:   snes->ops->destroy        = SNESDestroy_GS;
185:   snes->ops->setup          = SNESSetUp_GS;
186:   snes->ops->setfromoptions = SNESSetFromOptions_GS;
187:   snes->ops->view           = SNESView_GS;
188:   snes->ops->solve          = SNESSolve_GS;
189:   snes->ops->reset          = SNESReset_GS;

191:   snes->usesksp             = PETSC_FALSE;
192:   snes->usespc              = PETSC_FALSE;

194:   if (!snes->tolerancesset) {
195:     snes->max_its             = 10000;
196:     snes->max_funcs           = 10000;
197:   }

199:   PetscNewLog(snes, SNES_GS, &gs);
200:   snes->data = (void*) gs;
201:   return(0);
202: }
203: EXTERN_C_END