Actual source code: damgsnesad.c

  1: #define PETSCSNES_DLL
  2: 
 3:  #include petscda.h
 4:  #include petscmg.h
 5:  #include petscdmmg.h
 6:  #include src/inline/ilu.h
 7:  #include src/snes/impls/ls/ls.h

 10: EXTERN PetscErrorCode  NLFRelax_DAAD(NLF,MatSORType,PetscInt,Vec);
 11: EXTERN PetscErrorCode  NLFRelax_DAAD4(NLF,MatSORType,PetscInt,Vec);
 12: EXTERN PetscErrorCode  NLFRelax_DAAD9(NLF,MatSORType,PetscInt,Vec);
 13: EXTERN PetscErrorCode  NLFRelax_DAADb(NLF,MatSORType,PetscInt,Vec);
 15: EXTERN PetscErrorCode DMMGFormFunction(SNES,Vec,Vec,void *);
 16: EXTERN PetscErrorCode SNESLSCheckLocalMin_Private(Mat,Vec,Vec,PetscReal,PetscTruth*);

 20: /*
 21:     DMMGComputeJacobianWithAdic - Evaluates the Jacobian via Adic when the user has provided
 22:     a local function evaluation routine.
 23: */
 24: PetscErrorCode DMMGComputeJacobianWithAdic(SNES snes,Vec X,Mat *J,Mat *B,MatStructure *flag,void *ptr)
 25: {
 26:   DMMG           dmmg = (DMMG) ptr;
 28:   Vec            localX;
 29:   DA             da = (DA) dmmg->dm;

 32:   DAGetLocalVector(da,&localX);
 33:   DAGlobalToLocalBegin(da,X,INSERT_VALUES,localX);
 34:   DAGlobalToLocalEnd(da,X,INSERT_VALUES,localX);
 35:   DAComputeJacobian1WithAdic(da,localX,*B,dmmg->user);
 36:   DARestoreLocalVector(da,&localX);
 37:   /* Assemble true Jacobian; if it is different */
 38:   if (*J != *B) {
 39:     MatAssemblyBegin(*J,MAT_FINAL_ASSEMBLY);
 40:     MatAssemblyEnd(*J,MAT_FINAL_ASSEMBLY);
 41:   }
 42:   MatSetOption(*B,MAT_NEW_NONZERO_LOCATION_ERR);
 43:   *flag = SAME_NONZERO_PATTERN;
 44:   return(0);
 45: }

 49: /*@
 50:     SNESDAComputeJacobianWithAdic - This is a universal Jacobian evaluation routine
 51:     that may be used with SNESSetJacobian() as long as the user context has a DA as
 52:     its first record and DASetLocalAdicFunction() has been called.  

 54:    Collective on SNES

 56:    Input Parameters:
 57: +  snes - the SNES context
 58: .  X - input vector
 59: .  J - Jacobian
 60: .  B - Jacobian used in preconditioner (usally same as J)
 61: .  flag - indicates if the matrix changed its structure
 62: -  ptr - optional user-defined context, as set by SNESSetFunction()

 64:    Level: intermediate

 66: .seealso: DASetLocalFunction(), DASetLocalAdicFunction(), SNESSetFunction(), SNESSetJacobian()

 68: @*/
 69: PetscErrorCode  SNESDAComputeJacobianWithAdic(SNES snes,Vec X,Mat *J,Mat *B,MatStructure *flag,void *ptr)
 70: {
 71:   DA             da = *(DA*) ptr;
 73:   Vec            localX;

 76:   DAGetLocalVector(da,&localX);
 77:   DAGlobalToLocalBegin(da,X,INSERT_VALUES,localX);
 78:   DAGlobalToLocalEnd(da,X,INSERT_VALUES,localX);
 79:   DAComputeJacobian1WithAdic(da,localX,*B,ptr);
 80:   DARestoreLocalVector(da,&localX);
 81:   /* Assemble true Jacobian; if it is different */
 82:   if (*J != *B) {
 83:     MatAssemblyBegin(*J,MAT_FINAL_ASSEMBLY);
 84:     MatAssemblyEnd(*J,MAT_FINAL_ASSEMBLY);
 85:   }
 86:   MatSetOption(*B,MAT_NEW_NONZERO_LOCATION_ERR);
 87:   *flag = SAME_NONZERO_PATTERN;
 88:   return(0);
 89: }

 91:  #include src/ksp/pc/impls/mg/mgimpl.h
 92: /*
 93:           This is pre-beta FAS code. It's design should not be taken seriously!

 95:               R is the usual multigrid restriction (e.g. the tranpose of peicewise linear interpolation)
 96:               Q is either a scaled injection or the usual R
 97: */
100: PetscErrorCode DMMGSolveFAS(DMMG *dmmg,PetscInt level)
101: {
103:   PetscInt       i,j,k;
104:   PetscReal      norm;
105:   PC_MG          **mg;
106:   PC             pc;

109:   VecSet(dmmg[level]->r,0.0);
110:   for (j=1; j<=level; j++) {
111:     if (!dmmg[j]->inject) {
112:       DMGetInjection(dmmg[j-1]->dm,dmmg[j]->dm,&dmmg[j]->inject);
113:     }
114:   }

116:   KSPGetPC(dmmg[level]->ksp,&pc);
117:   mg   = ((PC_MG**)pc->data);

119:   for (i=0; i<100; i++) {

121:     for (j=level; j>0; j--) {

123:       /* Relax residual_fine - F(x_fine) = 0 */
124:       for (k=0; k<dmmg[j]->presmooth; k++) {
125:         NLFRelax_DAAD(dmmg[j]->nlf,SOR_SYMMETRIC_SWEEP,1,dmmg[j]->x);
126:       }

128:       /* R*(residual_fine - F(x_fine)) */
129:       DMMGFormFunction(0,dmmg[j]->x,dmmg[j]->w,dmmg[j]);
130:       VecAYPX(dmmg[j]->w,-1.0,dmmg[j]->r);

132:       if (j == level || dmmg[j]->monitorall) {
133:         /* norm( residual_fine - f(x_fine) ) */
134:         VecNorm(dmmg[j]->w,NORM_2,&norm);
135:         if (j == level) {
136:           if (norm < dmmg[level]->abstol) goto theend;
137:           if (i == 0) {
138:             dmmg[level]->rrtol = norm*dmmg[level]->rtol;
139:           } else {
140:             if (norm < dmmg[level]->rrtol) goto theend;
141:           }
142:         }
143:       }

145:       if (dmmg[j]->monitorall) {
146:         for (k=0; k<level-j+1; k++) {PetscPrintf(dmmg[j]->comm,"  ");}
147:         PetscPrintf(dmmg[j]->comm,"FAS function norm %G\n",norm);
148:       }
149:       MatRestrict(mg[j]->restrct,dmmg[j]->w,dmmg[j-1]->r);
150: 
151:       /* F(Q*x_fine) */
152:       VecScatterBegin(dmmg[j]->x,dmmg[j-1]->x,INSERT_VALUES,SCATTER_FORWARD,dmmg[j]->inject);
153:       VecScatterEnd(dmmg[j]->x,dmmg[j-1]->x,INSERT_VALUES,SCATTER_FORWARD,dmmg[j]->inject);
154:       DMMGFormFunction(0,dmmg[j-1]->x,dmmg[j-1]->w,dmmg[j-1]);

156:       /* residual_coarse = F(Q*x_fine) + R*(residual_fine - F(x_fine)) */
157:       VecAYPX(dmmg[j-1]->r,1.0,dmmg[j-1]->w);

159:       /* save Q*x_fine into b (needed when interpolating compute x back up */
160:       VecCopy(dmmg[j-1]->x,dmmg[j-1]->b);
161:     }

163:     for (j=0; j<dmmg[0]->coarsesmooth; j++) {
164:       NLFRelax_DAAD(dmmg[0]->nlf,SOR_SYMMETRIC_SWEEP,1,dmmg[0]->x);
165:     }
166:     if (dmmg[0]->monitorall){
167:       DMMGFormFunction(0,dmmg[0]->x,dmmg[0]->w,dmmg[0]);
168:       VecAXPY(dmmg[0]->w,-1.0,dmmg[0]->r);
169:       VecNorm(dmmg[0]->w,NORM_2,&norm);
170:       for (k=0; k<level+1; k++) {PetscPrintf(dmmg[0]->comm,"  ");}
171:       PetscPrintf(dmmg[0]->comm,"FAS coarse grid function norm %G\n",norm);
172:     }

174:     for (j=1; j<=level; j++) {
175:       /* x_fine = x_fine + R'*(x_coarse - Q*x_fine) */
176:       VecAXPY(dmmg[j-1]->x,-1.0,dmmg[j-1]->b);
177:       MatInterpolateAdd(mg[j]->interpolate,dmmg[j-1]->x,dmmg[j]->x,dmmg[j]->x);

179:       if (dmmg[j]->monitorall) {
180:         /* norm( F(x_fine) - residual_fine ) */
181:         DMMGFormFunction(0,dmmg[j]->x,dmmg[j]->w,dmmg[j]);
182:         VecAXPY(dmmg[j]->w,-1.0,dmmg[j]->r);
183:         VecNorm(dmmg[j]->w,NORM_2,&norm);
184:         for (k=0; k<level-j+1; k++) {PetscPrintf(dmmg[j]->comm,"  ");}
185:         PetscPrintf(dmmg[j]->comm,"FAS function norm %G\n",norm);
186:       }

188:       /* Relax residual_fine - F(x_fine)  = 0 */
189:       for (k=0; k<dmmg[j]->postsmooth; k++) {
190:         NLFRelax_DAAD(dmmg[j]->nlf,SOR_SYMMETRIC_SWEEP,1,dmmg[j]->x);
191:       }

193:       if (dmmg[j]->monitorall) {
194:         /* norm( F(x_fine) - residual_fine ) */
195:         DMMGFormFunction(0,dmmg[j]->x,dmmg[j]->w,dmmg[j]);
196:         VecAXPY(dmmg[j]->w,-1.0,dmmg[j]->r);
197:         VecNorm(dmmg[j]->w,NORM_2,&norm);
198:         for (k=0; k<level-j+1; k++) {PetscPrintf(dmmg[j]->comm,"  ");}
199:         PetscPrintf(dmmg[j]->comm,"FAS function norm %G\n",norm);
200:       }
201:     }

203:     if (dmmg[level]->monitor){
204:       DMMGFormFunction(0,dmmg[level]->x,dmmg[level]->w,dmmg[level]);
205:       VecNorm(dmmg[level]->w,NORM_2,&norm);
206:       PetscPrintf(dmmg[level]->comm,"%D FAS function norm %G\n",i+1,norm);
207:     }
208:   }
209:   theend:
210:   return(0);
211: }

213: /*
214:     This is the point-block version of FAS
215: */
218: PetscErrorCode DMMGSolveFASb(DMMG *dmmg,PetscInt level)
219: {
221:   PetscInt       i,j,k;
222:   PetscReal      norm;
223:   PC_MG          **mg;
224:   PC             pc;

227:   VecSet(dmmg[level]->r,0.0);
228:   for (j=1; j<=level; j++) {
229:     if (!dmmg[j]->inject) {
230:       DMGetInjection(dmmg[j-1]->dm,dmmg[j]->dm,&dmmg[j]->inject);
231:     }
232:   }

234:   KSPGetPC(dmmg[level]->ksp,&pc);
235:   mg   = ((PC_MG**)pc->data);

237:   for (i=0; i<100; i++) {

239:     for (j=level; j>0; j--) {

241:       /* Relax residual_fine - F(x_fine) = 0 */
242:       for (k=0; k<dmmg[j]->presmooth; k++) {
243:         NLFRelax_DAADb(dmmg[j]->nlf,SOR_SYMMETRIC_SWEEP,1,dmmg[j]->x);
244:       }

246:       /* R*(residual_fine - F(x_fine)) */
247:       DMMGFormFunction(0,dmmg[j]->x,dmmg[j]->w,dmmg[j]);
248:       VecAYPX(dmmg[j]->w,-1.0,dmmg[j]->r);

250:       if (j == level || dmmg[j]->monitorall) {
251:         /* norm( residual_fine - f(x_fine) ) */
252:         VecNorm(dmmg[j]->w,NORM_2,&norm);
253:         if (j == level) {
254:           if (norm < dmmg[level]->abstol) goto theend;
255:           if (i == 0) {
256:             dmmg[level]->rrtol = norm*dmmg[level]->rtol;
257:           } else {
258:             if (norm < dmmg[level]->rrtol) goto theend;
259:           }
260:         }
261:       }

263:       if (dmmg[j]->monitorall) {
264:         for (k=0; k<level-j+1; k++) {PetscPrintf(dmmg[j]->comm,"  ");}
265:         PetscPrintf(dmmg[j]->comm,"FAS function norm %G\n",norm);
266:       }
267:       MatRestrict(mg[j]->restrct,dmmg[j]->w,dmmg[j-1]->r);
268: 
269:       /* F(Q*x_fine) */
270:       VecScatterBegin(dmmg[j]->x,dmmg[j-1]->x,INSERT_VALUES,SCATTER_FORWARD,dmmg[j]->inject);
271:       VecScatterEnd(dmmg[j]->x,dmmg[j-1]->x,INSERT_VALUES,SCATTER_FORWARD,dmmg[j]->inject);
272:       DMMGFormFunction(0,dmmg[j-1]->x,dmmg[j-1]->w,dmmg[j-1]);

274:       /* residual_coarse = F(Q*x_fine) + R*(residual_fine - F(x_fine)) */
275:       VecAYPX(dmmg[j-1]->r,1.0,dmmg[j-1]->w);

277:       /* save Q*x_fine into b (needed when interpolating compute x back up */
278:       VecCopy(dmmg[j-1]->x,dmmg[j-1]->b);
279:     }

281:     for (j=0; j<dmmg[0]->coarsesmooth; j++) {
282:       NLFRelax_DAADb(dmmg[0]->nlf,SOR_SYMMETRIC_SWEEP,1,dmmg[0]->x);
283:     }
284:     if (dmmg[0]->monitorall){
285:       DMMGFormFunction(0,dmmg[0]->x,dmmg[0]->w,dmmg[0]);
286:       VecAXPY(dmmg[0]->w,-1.0,dmmg[0]->r);
287:       VecNorm(dmmg[0]->w,NORM_2,&norm);
288:       for (k=0; k<level+1; k++) {PetscPrintf(dmmg[0]->comm,"  ");}
289:       PetscPrintf(dmmg[0]->comm,"FAS coarse grid function norm %G\n",norm);
290:     }

292:     for (j=1; j<=level; j++) {
293:       /* x_fine = x_fine + R'*(x_coarse - Q*x_fine) */
294:       VecAXPY(dmmg[j-1]->x,-1.0,dmmg[j-1]->b);
295:       MatInterpolateAdd(mg[j]->interpolate,dmmg[j-1]->x,dmmg[j]->x,dmmg[j]->x);

297:       if (dmmg[j]->monitorall) {
298:         /* norm( F(x_fine) - residual_fine ) */
299:         DMMGFormFunction(0,dmmg[j]->x,dmmg[j]->w,dmmg[j]);
300:         VecAXPY(dmmg[j]->w,-1.0,dmmg[j]->r);
301:         VecNorm(dmmg[j]->w,NORM_2,&norm);
302:         for (k=0; k<level-j+1; k++) {PetscPrintf(dmmg[j]->comm,"  ");}
303:         PetscPrintf(dmmg[j]->comm,"FAS function norm %G\n",norm);
304:       }

306:       /* Relax residual_fine - F(x_fine)  = 0 */
307:       for (k=0; k<dmmg[j]->postsmooth; k++) {
308:         NLFRelax_DAADb(dmmg[j]->nlf,SOR_SYMMETRIC_SWEEP,1,dmmg[j]->x);
309:       }

311:       if (dmmg[j]->monitorall) {
312:         /* norm( F(x_fine) - residual_fine ) */
313:         DMMGFormFunction(0,dmmg[j]->x,dmmg[j]->w,dmmg[j]);
314:         VecAXPY(dmmg[j]->w,-1.0,dmmg[j]->r);
315:         VecNorm(dmmg[j]->w,NORM_2,&norm);
316:         for (k=0; k<level-j+1; k++) {PetscPrintf(dmmg[j]->comm,"  ");}
317:         PetscPrintf(dmmg[j]->comm,"FAS function norm %G\n",norm);
318:       }
319:     }

321:     if (dmmg[level]->monitor){
322:       DMMGFormFunction(0,dmmg[level]->x,dmmg[level]->w,dmmg[level]);
323:       VecNorm(dmmg[level]->w,NORM_2,&norm);
324:       PetscPrintf(dmmg[level]->comm,"%D FAS function norm %G\n",i+1,norm);
325:     }
326:   }
327:   theend:
328:   return(0);
329: }

332: #include "adic/ad_utils.h"

337: PetscErrorCode PetscADView(PetscInt N,PetscInt nc,double *ptr,PetscViewer viewer)
338: {
339:   PetscInt       i,j,nlen  = PetscADGetDerivTypeSize();
340:   char           *cptr = (char*)ptr;
341:   double         *values;

345:   for (i=0; i<N; i++) {
346:     PetscPrintf(PETSC_COMM_SELF,"Element %D value %G derivatives: ",i,*(double*)cptr);
347:     values = PetscADGetGradArray(cptr);
348:     for (j=0; j<nc; j++) {
349:       PetscPrintf(PETSC_COMM_SELF,"%G ",*values++);
350:     }
351:     PetscPrintf(PETSC_COMM_SELF,"\n");
352:     cptr += nlen;
353:   }

355:   return(0);
356: }

360: PetscErrorCode DMMGSolveFAS4(DMMG *dmmg,PetscInt level)
361: {
363:   PetscInt       i,j,k;
364:   PetscReal      norm;
365:   PetscScalar    zero = 0.0,mone = -1.0,one = 1.0;
366:   PC_MG          **mg;
367:   PC             pc;

370:   VecSet(dmmg[level]->r,zero);
371:   for (j=1; j<=level; j++) {
372:     if (!dmmg[j]->inject) {
373:       DMGetInjection(dmmg[j-1]->dm,dmmg[j]->dm,&dmmg[j]->inject);
374:     }
375:   }

377:   KSPGetPC(dmmg[level]->ksp,&pc);
378:   mg   = ((PC_MG**)pc->data);
379:   for (i=0; i<100; i++) {

381:     for (j=level; j>0; j--) {
382:   PetscPrintf(PETSC_COMM_WORLD,"I am here");
383:       /* Relax residual_fine - F(x_fine) = 0 */
384:       for (k=0; k<dmmg[j]->presmooth; k++) {
385:         NLFRelax_DAAD4(dmmg[j]->nlf,SOR_SYMMETRIC_SWEEP,1,dmmg[j]->x);
386:       }

388:       /* R*(residual_fine - F(x_fine)) */
389:       DMMGFormFunction(0,dmmg[j]->x,dmmg[j]->w,dmmg[j]);
390:       VecAYPX(dmmg[j]->w,mone,dmmg[j]->r);

392:       if (j == level || dmmg[j]->monitorall) {
393:         /* norm( residual_fine - f(x_fine) ) */
394:         VecNorm(dmmg[j]->w,NORM_2,&norm);
395:         if (j == level) {
396:           if (norm < dmmg[level]->abstol) goto theend;
397:           if (i == 0) {
398:             dmmg[level]->rrtol = norm*dmmg[level]->rtol;
399:           } else {
400:             if (norm < dmmg[level]->rrtol) goto theend;
401:           }
402:         }
403:       }  PetscPrintf(PETSC_COMM_WORLD,"I am here");

405:       if (dmmg[j]->monitorall) {
406:         for (k=0; k<level-j+1; k++) {PetscPrintf(dmmg[j]->comm,"  ");}
407:         PetscPrintf(dmmg[j]->comm,"FAS function norm %g\n",norm);
408:       }
409:       MatRestrict(mg[j]->restrct,dmmg[j]->w,dmmg[j-1]->r);
410: 
411:       /* F(R*x_fine) */
412:       VecScatterBegin(dmmg[j]->x,dmmg[j-1]->x,INSERT_VALUES,SCATTER_FORWARD,dmmg[j]->inject);
413:       VecScatterEnd(dmmg[j]->x,dmmg[j-1]->x,INSERT_VALUES,SCATTER_FORWARD,dmmg[j]->inject);
414:       DMMGFormFunction(0,dmmg[j-1]->x,dmmg[j-1]->w,dmmg[j-1]);

416:       /* residual_coarse = F(R*x_fine) + R*(residual_fine - F(x_fine)) */
417:       VecAYPX(dmmg[j-1]->r,one,dmmg[j-1]->w);

419:       /* save R*x_fine into b (needed when interpolating compute x back up */
420:       VecCopy(dmmg[j-1]->x,dmmg[j-1]->b);
421:     }

423:     for (j=0; j<dmmg[0]->presmooth; j++) {
424:       NLFRelax_DAAD4(dmmg[0]->nlf,SOR_SYMMETRIC_SWEEP,1,dmmg[0]->x);
425:     }
426:     if (dmmg[0]->monitorall){
427:       DMMGFormFunction(0,dmmg[0]->x,dmmg[0]->w,dmmg[0]);
428:       VecAXPY(dmmg[0]->w,mone,dmmg[0]->r);
429:       VecNorm(dmmg[0]->w,NORM_2,&norm);
430:       for (k=0; k<level+1; k++) {PetscPrintf(dmmg[0]->comm,"  ");}
431:       PetscPrintf(dmmg[0]->comm,"FAS coarse grid function norm %g\n",norm);
432:     }

434:     for (j=1; j<=level; j++) {
435:       /* x_fine = x_fine + R'*(x_coarse - R*x_fine) */
436:       VecAXPY(dmmg[j-1]->x,mone,dmmg[j-1]->b);
437:       MatInterpolateAdd(mg[j]->interpolate,dmmg[j-1]->x,dmmg[j]->x,dmmg[j]->x);

439:       if (dmmg[j]->monitorall) {
440:         /* norm( F(x_fine) - residual_fine ) */
441:         DMMGFormFunction(0,dmmg[j]->x,dmmg[j]->w,dmmg[j]);
442:         VecAXPY(dmmg[j]->w,mone,dmmg[j]->r);
443:         VecNorm(dmmg[j]->w,NORM_2,&norm);
444:         for (k=0; k<level-j+1; k++) {PetscPrintf(dmmg[j]->comm,"  ");}
445:         PetscPrintf(dmmg[j]->comm,"FAS function norm %g\n",norm);
446:       }

448:       /* Relax residual_fine - F(x_fine)  = 0 */
449:       for (k=0; k<dmmg[j]->postsmooth; k++) {
450:         NLFRelax_DAAD4(dmmg[j]->nlf,SOR_SYMMETRIC_SWEEP,1,dmmg[j]->x);
451:       }

453:       if (dmmg[j]->monitorall) {
454:         /* norm( F(x_fine) - residual_fine ) */
455:         DMMGFormFunction(0,dmmg[j]->x,dmmg[j]->w,dmmg[j]);
456:         VecAXPY(dmmg[j]->w,mone,dmmg[j]->r);
457:         VecNorm(dmmg[j]->w,NORM_2,&norm);
458:         for (k=0; k<level-j+1; k++) {PetscPrintf(dmmg[j]->comm,"  ");}
459:         PetscPrintf(dmmg[j]->comm,"FAS function norm %g\n",norm);
460:       }
461:     }

463:     if (dmmg[level]->monitor){
464:       DMMGFormFunction(0,dmmg[level]->x,dmmg[level]->w,dmmg[level]);
465:       VecNorm(dmmg[level]->w,NORM_2,&norm);
466:       PetscPrintf(dmmg[level]->comm,"%D FAS function norm %g\n",i,norm);
467:     }
468:   }
469:   theend:
470:   return(0);
471: }

473: /*
474:    This function provide several FAS v_cycle iteration 

476:    iter: the number of FAS it run         

478: */
481: PetscErrorCode DMMGSolveFASn(DMMG *dmmg,PetscInt level,PetscInt iter)
482: {
484:   PetscInt       i,j,k;
485:   PetscReal      norm;
486:   PC_MG          **mg;
487:   PC             pc;

490:   VecSet(dmmg[level]->r,0.0);
491:   for (j=1; j<=level; j++) {
492:     if (!dmmg[j]->inject) {
493:       DMGetInjection(dmmg[j-1]->dm,dmmg[j]->dm,&dmmg[j]->inject);
494:     }
495:   }

497:   KSPGetPC(dmmg[level]->ksp,&pc);
498:   mg   = ((PC_MG**)pc->data);

500:   for (i=0; i<iter; i++) {

502:     for (j=level; j>0; j--) {

504:       /* Relax residual_fine - F(x_fine) = 0 */
505:       for (k=0; k<dmmg[j]->presmooth; k++) {
506:         NLFRelax_DAAD(dmmg[j]->nlf,SOR_SYMMETRIC_SWEEP,1,dmmg[j]->x);
507:       }

509:       /* R*(residual_fine - F(x_fine)) */
510:       DMMGFormFunction(0,dmmg[j]->x,dmmg[j]->w,dmmg[j]);
511:       VecAYPX(dmmg[j]->w,-1.0,dmmg[j]->r);

513:       if (j == level || dmmg[j]->monitorall) {
514:         /* norm( residual_fine - f(x_fine) ) */
515:         VecNorm(dmmg[j]->w,NORM_2,&norm);
516:         if (j == level) {
517:           if (norm < dmmg[level]->abstol) goto theend;
518:           if (i == 0) {
519:             dmmg[level]->rrtol = norm*dmmg[level]->rtol;
520:           } else {
521:             if (norm < dmmg[level]->rrtol) goto theend;
522:           }
523:         }
524:       }

526:       if (dmmg[j]->monitorall) {
527:         for (k=0; k<level-j+1; k++) {PetscPrintf(dmmg[j]->comm,"  ");}
528:         PetscPrintf(dmmg[j]->comm,"FAS function norm %g\n",norm);
529:       }
530:       MatRestrict(mg[j]->restrct,dmmg[j]->w,dmmg[j-1]->r);
531: 
532:       /* F(RI*x_fine) */
533:       VecScatterBegin(dmmg[j]->x,dmmg[j-1]->x,INSERT_VALUES,SCATTER_FORWARD,dmmg[j]->inject);
534:       VecScatterEnd(dmmg[j]->x,dmmg[j-1]->x,INSERT_VALUES,SCATTER_FORWARD,dmmg[j]->inject);
535:       DMMGFormFunction(0,dmmg[j-1]->x,dmmg[j-1]->w,dmmg[j-1]);

537:       /* residual_coarse = F(RI*x_fine) + R*(residual_fine - F(x_fine)) */
538:       VecAYPX(dmmg[j-1]->r,1.0,dmmg[j-1]->w);

540:       /* save RI*x_fine into b (needed when interpolating compute x back up */
541:       VecCopy(dmmg[j-1]->x,dmmg[j-1]->b);
542:     }

544:     for (j=0; j<dmmg[0]->coarsesmooth; j++) {
545:       NLFRelax_DAAD(dmmg[0]->nlf,SOR_SYMMETRIC_SWEEP,1,dmmg[0]->x);
546:     }
547:     if (dmmg[0]->monitorall){
548:       DMMGFormFunction(0,dmmg[0]->x,dmmg[0]->w,dmmg[0]);
549:       VecAXPY(dmmg[0]->w,-1.0,dmmg[0]->r);
550:       VecNorm(dmmg[0]->w,NORM_2,&norm);
551:       for (k=0; k<level+1; k++) {PetscPrintf(dmmg[0]->comm,"  ");}
552:       PetscPrintf(dmmg[0]->comm,"FAS coarse grid function norm %g\n",norm);
553:     }

555:     for (j=1; j<=level; j++) {
556:       /* x_fine = x_fine + R'*(x_coarse - RI*x_fine) */
557:       VecAXPY(dmmg[j-1]->x,-1.0,dmmg[j-1]->b);
558:       MatInterpolateAdd(mg[j]->interpolate,dmmg[j-1]->x,dmmg[j]->x,dmmg[j]->x);

560:       if (dmmg[j]->monitorall) {
561:         /* norm( F(x_fine) - residual_fine ) */
562:         DMMGFormFunction(0,dmmg[j]->x,dmmg[j]->w,dmmg[j]);
563:         VecAXPY(dmmg[j]->w,-1.0,dmmg[j]->r);
564:         VecNorm(dmmg[j]->w,NORM_2,&norm);
565:         for (k=0; k<level-j+1; k++) {PetscPrintf(dmmg[j]->comm,"  ");}
566:         PetscPrintf(dmmg[j]->comm,"FAS function norm %g\n",norm);
567:       }

569:       /* Relax residual_fine - F(x_fine)  = 0 */
570:       for (k=0; k<dmmg[j]->postsmooth; k++) {
571:         NLFRelax_DAAD(dmmg[j]->nlf,SOR_SYMMETRIC_SWEEP,1,dmmg[j]->x);
572:       }

574:       if (dmmg[j]->monitorall) {
575:         /* norm( F(x_fine) - residual_fine ) */
576:         DMMGFormFunction(0,dmmg[j]->x,dmmg[j]->w,dmmg[j]);
577:         VecAXPY(dmmg[j]->w,-1.0,dmmg[j]->r);
578:         VecNorm(dmmg[j]->w,NORM_2,&norm);
579:         for (k=0; k<level-j+1; k++) {PetscPrintf(dmmg[j]->comm,"  ");}
580:         PetscPrintf(dmmg[j]->comm,"FAS function norm %g\n",norm);
581:       }
582:     }

584:     if (dmmg[level]->monitor){
585:       DMMGFormFunction(0,dmmg[level]->x,dmmg[level]->w,dmmg[level]);
586:       VecNorm(dmmg[level]->w,NORM_2,&norm);
587:       PetscPrintf(dmmg[level]->comm,"%D FAS function norm %g\n",i+1,norm);
588:     }
589:   }
590:   theend:
591:   return(0);
592: }
593: /*
594:           This is a simple FAS setup function
595: */


600: PetscErrorCode DMMGSolveFASSetUp(DMMG *dmmg,PetscInt level)
601: {
603:   PetscInt       j;//,nlevels=dmmg[0]->nlevels-1;
604:   //PC             pc;

607:   VecSet(dmmg[level]->r,0.0);
608:   for (j=1; j<=level; j++) {
609:     if (!dmmg[j]->inject) {
610:       DMGetInjection(dmmg[j-1]->dm,dmmg[j]->dm,&dmmg[j]->inject);
611:     }
612:   }
613:   VecSet(dmmg[level]->r,0.0);
614:   dmmg[level]->rrtol = 0.0001*dmmg[level]->rtol;//I want to get more precise solution with FAS
615:    return(0);
616: }


619: /*
620:   This is function to implement multiplicative FAS 


623: Options:
624:  
625: -dmmg_fas_cycles 1 : FAS v-cycle
626:                  2 : FAS w-cycle


629: */

633: PetscErrorCode DMMGSolveFASMCycle(DMMG *dmmg,PetscInt level,PetscTruth* converged)
634: {
636:   PetscInt       j,k,cycles=1,nlevels=level;//nlevels=dmmg[0]->nlevels-1;
637:                   // I need to put nlevels=level in order to get grid sequence correctly
638:   PetscReal      norm;
639:   PC_MG          **mg;
640:   PC             pc;
641: 
643: 

645:   PetscOptionsGetInt(PETSC_NULL,"-dmmg_fas_cycles",&cycles,PETSC_NULL);
646: 
647:   KSPGetPC(dmmg[level]->ksp,&pc);
648:   mg   = ((PC_MG**)pc->data);

650:   j=level;

652:   if(j) {/* not the coarsest level */
653:     /* Relax residual_fine - F(x_fine) = 0 */
654:     for (k=0; k<dmmg[j]->presmooth; k++) {
655:       NLFRelax_DAAD(dmmg[j]->nlf,SOR_SYMMETRIC_SWEEP,1,dmmg[j]->x);
656:     }
657: 
658: 

660:     /* R*(residual_fine - F(x_fine)) */
661:     DMMGFormFunction(0,dmmg[j]->x,dmmg[j]->w,dmmg[j]);
662:     VecAYPX(dmmg[j]->w,-1.0,dmmg[j]->r);

664:     if (j == nlevels || dmmg[j]->monitorall) {
665:       /* norm( residual_fine - f(x_fine) ) */
666:       VecNorm(dmmg[j]->w,NORM_2,&norm);
667: 
668:       if (j == nlevels) {
669:         if (norm < dmmg[level]->abstol) {
670:           *converged = PETSC_TRUE;
671:            goto theend;
672:         }

674:           if (norm < dmmg[level]->rrtol){
675:           *converged = PETSC_TRUE;
676:           goto theend;
677: 
678:         }
679:       }
680:     }

682:     if (dmmg[j]->monitorall) {
683:       for (k=0; k<level-j+1; k++) {PetscPrintf(dmmg[j]->comm,"  ");}
684:       PetscPrintf(dmmg[j]->comm,"FAS function norm %g\n",norm);
685:     }
686:     MatRestrict(mg[j]->restrct,dmmg[j]->w,dmmg[j-1]->r);
687: 
688:     /* F(RI*x_fine) */
689:     VecScatterBegin(dmmg[j]->x,dmmg[j-1]->x,INSERT_VALUES,SCATTER_FORWARD,dmmg[j]->inject);
690:     VecScatterEnd(dmmg[j]->x,dmmg[j-1]->x,INSERT_VALUES,SCATTER_FORWARD,dmmg[j]->inject);
691:     DMMGFormFunction(0,dmmg[j-1]->x,dmmg[j-1]->w,dmmg[j-1]);

693:     /* residual_coarse = F(RI*x_fine) + R*(residual_fine - F(x_fine)) */
694:     VecAYPX(dmmg[j-1]->r,1.0,dmmg[j-1]->w);

696:     /* save RI*x_fine into b (needed when interpolating compute x back up */
697:     VecCopy(dmmg[j-1]->x,dmmg[j-1]->b);
698: 
699: 
700:     while (cycles--) {
701:       DMMGSolveFASMCycle(dmmg,level-1,converged);
702:     }
703:   }
704:   else { /* for the coarsest level */
705:     for (k=0; k<dmmg[0]->coarsesmooth; k++) {
706:       NLFRelax_DAAD(dmmg[0]->nlf,SOR_SYMMETRIC_SWEEP,1,dmmg[0]->x);
707:     }
708:     if (dmmg[0]->monitorall){
709:       DMMGFormFunction(0,dmmg[0]->x,dmmg[0]->w,dmmg[0]);
710:       VecAXPY(dmmg[0]->w,-1.0,dmmg[0]->r);
711:       VecNorm(dmmg[0]->w,NORM_2,&norm);
712:       for (k=0; k<level+1; k++) {PetscPrintf(dmmg[0]->comm,"  ");}
713:       PetscPrintf(dmmg[0]->comm,"FAS coarse grid function norm %g\n",norm);
714:     }
715:  if (j == nlevels || dmmg[j]->monitorall) {
716:       /* norm( residual_fine - f(x_fine) ) */
717:       VecNorm(dmmg[j]->w,NORM_2,&norm);
718: 
719:       if (j == nlevels) {
720:         if (norm < dmmg[level]->abstol) {
721:           *converged = PETSC_TRUE;
722:            goto theend;
723:         }
724: 
725:           if (norm < dmmg[level]->rrtol){
726:           *converged = PETSC_TRUE;
727:           goto theend;
728: 
729:         }
730:       }
731:     }

733: 
734:   }
735:   j=level;
736:   if(j) { /* not for the coarsest level */
737:     /* x_fine = x_fine + R'*(x_coarse - RI*x_fine) */
738:     VecAXPY(dmmg[j-1]->x,-1.0,dmmg[j-1]->b);
739:     MatInterpolateAdd(mg[j]->interpolate,dmmg[j-1]->x,dmmg[j]->x,dmmg[j]->x);

741:     if (dmmg[j]->monitorall) {
742:       /* norm( F(x_fine) - residual_fine ) */
743:       DMMGFormFunction(0,dmmg[j]->x,dmmg[j]->w,dmmg[j]);
744:       VecAXPY(dmmg[j]->w,-1.0,dmmg[j]->r);
745:       VecNorm(dmmg[j]->w,NORM_2,&norm);
746:       for (k=0; k<level-j+1; k++) {PetscPrintf(dmmg[j]->comm,"  ");}
747:       PetscPrintf(dmmg[j]->comm,"FAS function norm %g\n",norm);
748:     }

750:     /* Relax residual_fine - F(x_fine)  = 0 */
751:     for (k=0; k<dmmg[j]->postsmooth; k++) {
752:       NLFRelax_DAAD(dmmg[j]->nlf,SOR_SYMMETRIC_SWEEP,1,dmmg[j]->x);
753:     }

755:     if (dmmg[j]->monitorall) {
756:       /* norm( F(x_fine) - residual_fine ) */
757:       DMMGFormFunction(0,dmmg[j]->x,dmmg[j]->w,dmmg[j]);
758:       VecAXPY(dmmg[j]->w,-1.0,dmmg[j]->r);
759:       VecNorm(dmmg[j]->w,NORM_2,&norm);
760:       for (k=0; k<level-j+1; k++) {PetscPrintf(dmmg[j]->comm,"  ");}
761:       PetscPrintf(dmmg[j]->comm,"FAS function norm %g\n",norm);
762:     }
763: 
764: 
765: }
766: theend:
767: return(0);
768: }

770: /*
771:   This is function to implement multiplicative FAS with block smoother


774: Options:
775:  
776: -dmmg_fas_cycles 1 : FAS v-cycle
777:                  2 : FAS w-cycle


780: */


785: PetscErrorCode DMMGSolveFASMCycle9(DMMG *dmmg,PetscInt level,PetscTruth* converged)
786: {
788:   PetscInt       j,k,cycles=1,nlevels=level;//nlevels=dmmg[0]->nlevels-1;
789:                   // I need to put nlevels=level in order to get grid sequence correctly
790:   PetscReal      norm;
791:   PC_MG          **mg;
792:   PC             pc;
793: 
795: 
796: 
797:   PetscOptionsGetInt(PETSC_NULL,"-dmmg_fas_cycles",&cycles,PETSC_NULL);
798: 
799:   KSPGetPC(dmmg[level]->ksp,&pc);
800:   mg   = ((PC_MG**)pc->data);
801:   //   for (j=level; j>0; j--) {
802:   j=level;
803:   //PetscPrintf(dmmg[level]->comm,"j=%d,nlevels=%d",j,nlevels);
804:   if(j) {/* not the coarsest level */
805:     /* Relax residual_fine - F(x_fine) = 0 */
806:     for (k=0; k<dmmg[j]->presmooth; k++) {
807:       NLFRelax_DAAD9(dmmg[j]->nlf,SOR_SYMMETRIC_SWEEP,1,dmmg[j]->x);
808:     }
809: 
810: 

812:     /* R*(residual_fine - F(x_fine)) */
813:     DMMGFormFunction(0,dmmg[j]->x,dmmg[j]->w,dmmg[j]);
814:     VecAYPX(dmmg[j]->w,-1.0,dmmg[j]->r);

816:     if (j == nlevels || dmmg[j]->monitorall) {
817:       /* norm( residual_fine - f(x_fine) ) */
818:       VecNorm(dmmg[j]->w,NORM_2,&norm);
819: 
820:       if (j == nlevels) {
821:         if (norm < dmmg[level]->abstol) {
822:           *converged = PETSC_TRUE;
823:            goto theend;
824:         }
825:         /*        if (i == 0) {
826:           dmmg[level]->rrtol = norm*dmmg[level]->rtol;
827:           } else {*/
828:           if (norm < dmmg[level]->rrtol){
829:           *converged = PETSC_TRUE;
830:           goto theend;
831: 
832:         }
833:       }
834:     }

836:     if (dmmg[j]->monitorall) {
837:       for (k=0; k<level-j+1; k++) {PetscPrintf(dmmg[j]->comm,"  ");}
838:       PetscPrintf(dmmg[j]->comm,"FAS function norm %g\n",norm);
839:     }
840:     MatRestrict(mg[j]->restrct,dmmg[j]->w,dmmg[j-1]->r);
841: 
842:     /* F(RI*x_fine) */
843:     VecScatterBegin(dmmg[j]->x,dmmg[j-1]->x,INSERT_VALUES,SCATTER_FORWARD,dmmg[j]->inject);
844:     VecScatterEnd(dmmg[j]->x,dmmg[j-1]->x,INSERT_VALUES,SCATTER_FORWARD,dmmg[j]->inject);
845:     DMMGFormFunction(0,dmmg[j-1]->x,dmmg[j-1]->w,dmmg[j-1]);

847:     /* residual_coarse = F(RI*x_fine) + R*(residual_fine - F(x_fine)) */
848:     VecAYPX(dmmg[j-1]->r,1.0,dmmg[j-1]->w);

850:     /* save RI*x_fine into b (needed when interpolating compute x back up */
851:     VecCopy(dmmg[j-1]->x,dmmg[j-1]->b);
852: 
853: 
854:     while (cycles--) {
855:       DMMGSolveFASMCycle9(dmmg,level-1,converged);
856:     }
857:   }
858:   else { /* for the coarsest level */
859:     for (k=0; k<dmmg[0]->coarsesmooth; k++) {
860:       NLFRelax_DAAD9(dmmg[0]->nlf,SOR_SYMMETRIC_SWEEP,1,dmmg[0]->x);
861:     }
862:     if (dmmg[0]->monitorall){
863:       DMMGFormFunction(0,dmmg[0]->x,dmmg[0]->w,dmmg[0]);
864:       VecAXPY(dmmg[0]->w,-1.0,dmmg[0]->r);
865:       VecNorm(dmmg[0]->w,NORM_2,&norm);
866:       for (k=0; k<level+1; k++) {PetscPrintf(dmmg[0]->comm,"  ");}
867:       PetscPrintf(dmmg[0]->comm,"FAS coarse grid function norm %g\n",norm);
868:     }
869:  if (j == nlevels || dmmg[j]->monitorall) {
870:       /* norm( residual_fine - f(x_fine) ) */
871:       VecNorm(dmmg[j]->w,NORM_2,&norm);
872: 
873:       if (j == nlevels) {
874:         if (norm < dmmg[level]->abstol) {
875:           *converged = PETSC_TRUE;
876:            goto theend;
877:         }
878:         /*        if (i == 0) {
879:           dmmg[level]->rrtol = norm*dmmg[level]->rtol;
880:           } else {*/
881:           if (norm < dmmg[level]->rrtol){
882:           *converged = PETSC_TRUE;
883:           goto theend;
884: 
885:         }
886:       }
887:     }

889: 
890:   }
891:   j=level;
892:   if(j) { /* not for the coarsest level */
893:     /* x_fine = x_fine + R'*(x_coarse - RI*x_fine) */
894:     VecAXPY(dmmg[j-1]->x,-1.0,dmmg[j-1]->b);
895:     MatInterpolateAdd(mg[j]->interpolate,dmmg[j-1]->x,dmmg[j]->x,dmmg[j]->x);

897:     if (dmmg[j]->monitorall) {
898:       /* norm( F(x_fine) - residual_fine ) */
899:       DMMGFormFunction(0,dmmg[j]->x,dmmg[j]->w,dmmg[j]);
900:       VecAXPY(dmmg[j]->w,-1.0,dmmg[j]->r);
901:       VecNorm(dmmg[j]->w,NORM_2,&norm);
902:       for (k=0; k<level-j+1; k++) {PetscPrintf(dmmg[j]->comm,"  ");}
903:       PetscPrintf(dmmg[j]->comm,"FAS function norm %g\n",norm);
904:     }

906:     /* Relax residual_fine - F(x_fine)  = 0 */
907:     for (k=0; k<dmmg[j]->postsmooth; k++) {
908:       NLFRelax_DAAD9(dmmg[j]->nlf,SOR_SYMMETRIC_SWEEP,1,dmmg[j]->x);
909:     }

911:     if (dmmg[j]->monitorall) {
912:       /* norm( F(x_fine) - residual_fine ) */
913:       DMMGFormFunction(0,dmmg[j]->x,dmmg[j]->w,dmmg[j]);
914:       VecAXPY(dmmg[j]->w,-1.0,dmmg[j]->r);
915:       VecNorm(dmmg[j]->w,NORM_2,&norm);
916:       for (k=0; k<level-j+1; k++) {PetscPrintf(dmmg[j]->comm,"  ");}
917:       PetscPrintf(dmmg[j]->comm,"FAS function norm %g\n",norm);
918:     }
919: 
920:     /* if(j==nlevels){
921:      if (dmmg[level]->monitor){
922:     DMMGFormFunction(0,dmmg[level]->x,dmmg[level]->w,dmmg[level]);
923:     VecNorm(dmmg[level]->w,NORM_2,&norm);
924:       
925:     PetscPrintf(dmmg[level]->comm," FAS function norm %g\n",norm);
926:     
927:     }
928:     }*/
929: }
930: theend:
931: return(0);
932: }

934: /*
935:   This is function to implement full FAS with block smoother(9 points together)


938: Options:
939:  
940: -dmmg_fas_cycles 1 : FAS v-cycle
941:                  2 : FAS w-cycle


944: */


949: PetscErrorCode DMMGSolveFASFCycle(DMMG *dmmg,PetscInt l,PetscTruth* converged)
950: {
952:   PetscInt       j;//l = dmmg[0]->nlevels-1;
953:   PC_MG          **mg;
954:   PC             pc;

957:   KSPGetPC(dmmg[l]->ksp,&pc);
958:   mg   = ((PC_MG**)pc->data);
959:   // restriction all the way down to the coarse level
960:   if(l>0) {
961:     for(j=l;j>0;j--) {
962: 
963:       /* R*(residual_fine - F(x_fine)) */
964:       DMMGFormFunction(0,dmmg[j]->x,dmmg[j]->w,dmmg[j]);
965:       VecAYPX(dmmg[j]->w,-1.0,dmmg[j]->r);

967:       MatRestrict(mg[j]->restrct,dmmg[j]->w,dmmg[j-1]->r);
968: 
969:       /* F(RI*x_fine) */
970:       VecScatterBegin(dmmg[j]->x,dmmg[j-1]->x,INSERT_VALUES,SCATTER_FORWARD,dmmg[j]->inject);
971:       VecScatterEnd(dmmg[j]->x,dmmg[j-1]->x,INSERT_VALUES,SCATTER_FORWARD,dmmg[j]->inject);
972:       DMMGFormFunction(0,dmmg[j-1]->x,dmmg[j-1]->w,dmmg[j-1]);

974:       /* residual_coarse = F(RI*x_fine) + R*(residual_fine - F(x_fine)) */
975:       VecAYPX(dmmg[j-1]->r,1.0,dmmg[j-1]->w);

977:       /* save RI*x_fine into b (needed when interpolating compute x back up */
978:       VecCopy(dmmg[j-1]->x,dmmg[j-1]->b);

980:     }

982:     // all the way up to the finest level
983:     for (j=0; j<l; j++) {
984:       DMMGSolveFASMCycle(dmmg,j,PETSC_NULL);
985:       /* x_fine = x_fine + R'*(x_coarse - RI*x_fine) */
986:       VecAXPY(dmmg[j]->x,-1.0,dmmg[j]->b);
987:       MatInterpolateAdd(mg[j+1]->interpolate,dmmg[j]->x,dmmg[j+1]->x,dmmg[j+1]->x);

989:     }
990:   }
991:   DMMGSolveFASMCycle(dmmg,l,converged);
992:   return(0);
993: }



997: /*
998:   This is function to implement full FAS  with block smoother ( 9 points together)


1001: Options:
1002:  
1003: -dmmg_fas_cycles 1 : FAS v-cycle
1004:                  2 : FAS w-cycle


1007: */
1010: PetscErrorCode DMMGSolveFASFCycle9(DMMG *dmmg,PetscInt l,PetscTruth* converged)
1011: {
1013:   PetscInt       j;//l = dmmg[0]->nlevels-1;
1014:   PC_MG          **mg;
1015:   PC             pc;

1018:   KSPGetPC(dmmg[l]->ksp,&pc);
1019:   mg   = ((PC_MG**)pc->data);
1020:   // restriction all the way down to the coarse level
1021:   if(l>0) {
1022:     for(j=l;j>0;j--) {
1023: 
1024:       /* R*(residual_fine - F(x_fine)) */
1025:       DMMGFormFunction(0,dmmg[j]->x,dmmg[j]->w,dmmg[j]);
1026:       VecAYPX(dmmg[j]->w,-1.0,dmmg[j]->r);

1028:       MatRestrict(mg[j]->restrct,dmmg[j]->w,dmmg[j-1]->r);
1029: 
1030:       /* F(RI*x_fine) */
1031:       VecScatterBegin(dmmg[j]->x,dmmg[j-1]->x,INSERT_VALUES,SCATTER_FORWARD,dmmg[j]->inject);
1032:       VecScatterEnd(dmmg[j]->x,dmmg[j-1]->x,INSERT_VALUES,SCATTER_FORWARD,dmmg[j]->inject);
1033:       DMMGFormFunction(0,dmmg[j-1]->x,dmmg[j-1]->w,dmmg[j-1]);

1035:       /* residual_coarse = F(RI*x_fine) + R*(residual_fine - F(x_fine)) */
1036:       VecAYPX(dmmg[j-1]->r,1.0,dmmg[j-1]->w);

1038:       /* save RI*x_fine into b (needed when interpolating compute x back up */
1039:       VecCopy(dmmg[j-1]->x,dmmg[j-1]->b);

1041:     }

1043:     // all the way up to the finest level
1044:     for (j=0; j<l; j++) {
1045:       DMMGSolveFASMCycle9(dmmg,j,PETSC_NULL);
1046:       /* x_fine = x_fine + R'*(x_coarse - RI*x_fine) */
1047:       VecAXPY(dmmg[j]->x,-1.0,dmmg[j]->b);
1048:       MatInterpolateAdd(mg[j+1]->interpolate,dmmg[j]->x,dmmg[j+1]->x,dmmg[j+1]->x);

1050:     }
1051:   }
1052:   DMMGSolveFASMCycle9(dmmg,l,converged);
1053:   return(0);
1054: }

1056: /*
1057:           This is function is to solve nonlinear system with FAS

1059: Options:

1061: -dmmg_fas_9:    using block smoother 
1062:  
1063: -dmmg_fas_full: using full FAS


1066: */
1069: PetscErrorCode DMMGSolveFASCycle(DMMG *dmmg,PetscInt level)
1070: {
1072:   PetscInt       i;
1073:   PetscTruth     converged = PETSC_FALSE, flg,flgb;
1074:   PetscReal      norm;

1077:    DMMGSolveFASSetUp(dmmg,level);
1078:   PetscOptionsHasName(PETSC_NULL,"-dmmg_fas_9",&flgb);
1079: 
1080:   if(flgb){

1082:     PetscOptionsHasName(PETSC_NULL,"-dmmg_fas_full",&flg);
1083:     if (flg) {
1084:       for(i=0;i<1000;i++){
1085:         PetscPrintf(dmmg[level]->comm,"%D ",i+1);
1086:         DMMGSolveFASFCycle9(dmmg,level,&converged);
1087:         if (dmmg[level]->monitor){
1088:           DMMGFormFunction(0,dmmg[level]->x,dmmg[level]->w,dmmg[level]);
1089:           VecNorm(dmmg[level]->w,NORM_2,&norm);
1090:           PetscPrintf(dmmg[level]->comm," FAS function norm %g\n",norm);
1091:         }
1092:         if (converged) return(0);
1093:       }
1094:     }
1095:     else{
1096:       for(i=0;i<1000;i++){
1097:         PetscPrintf(dmmg[level]->comm,"%D ",i+1);
1098:         DMMGSolveFASMCycle9(dmmg,level,&converged);
1099:         if (dmmg[level]->monitor){
1100:           DMMGFormFunction(0,dmmg[level]->x,dmmg[level]->w,dmmg[level]);
1101:           VecNorm(dmmg[level]->w,NORM_2,&norm);
1102:           PetscPrintf(dmmg[level]->comm," FAS function norm %g\n",norm);
1103:         }

1105:         if (converged) return(0);
1106:       }
1107:     }
1108:   }
1109:   else {
1110:     PetscOptionsHasName(PETSC_NULL,"-dmmg_fas_full",&flg);
1111:     if (flg) {
1112:       for(i=0;i<1000;i++){
1113:         PetscPrintf(dmmg[level]->comm,"%D ",i+1);
1114:         DMMGSolveFASFCycle(dmmg,level,&converged);
1115:         if (dmmg[level]->monitor){
1116:           DMMGFormFunction(0,dmmg[level]->x,dmmg[level]->w,dmmg[level]);
1117:           VecNorm(dmmg[level]->w,NORM_2,&norm);
1118:           PetscPrintf(dmmg[level]->comm," FAS function norm %g\n",norm);
1119:         }
1120:         if (converged) return(0);
1121:       }
1122:     }
1123:     else{
1124:       for(i=0;i<1000;i++){
1125:         PetscPrintf(dmmg[level]->comm,"%D ",i+1);
1126:         DMMGSolveFASMCycle(dmmg,level,&converged);
1127:         if (dmmg[level]->monitor){
1128:           DMMGFormFunction(0,dmmg[level]->x,dmmg[level]->w,dmmg[level]);
1129:           VecNorm(dmmg[level]->w,NORM_2,&norm);
1130:           PetscPrintf(dmmg[level]->comm," FAS function norm %g\n",norm);
1131:         }

1133:         if (converged) return(0);
1134:       }

1136:     }
1137:   }
1138:   return(0);
1139: }

1141: /*
1142:           This is function is to implement one  FAS iteration 

1144: Options:

1146: -dmmg_fas_9:    using block smoother 
1147:  
1148: -dmmg_fas_full: using full FAS

1150: */
1153: PetscErrorCode DMMGSolveFASCyclen(DMMG *dmmg,PetscInt level)
1154: {
1156:   PetscTruth     converged = PETSC_FALSE, flg,flgb;
1157:   PetscReal      norm;
1158:   // PC_MG          **mg;
1159:   //PC             pc;

1162:    DMMGSolveFASSetUp(dmmg,level);
1163:      PetscOptionsHasName(PETSC_NULL,"-dmmg_fas_9",&flgb);
1164: 
1165:   if(flgb){

1167:     PetscOptionsHasName(PETSC_NULL,"-dmmg_fas_full",&flg);
1168:     if (flg) {
1169: 
1170:         DMMGSolveFASFCycle9(dmmg,level,&converged);
1171:         if (dmmg[level]->monitor){
1172:           DMMGFormFunction(0,dmmg[level]->x,dmmg[level]->w,dmmg[level]);
1173:           VecNorm(dmmg[level]->w,NORM_2,&norm);
1174:           PetscPrintf(dmmg[level]->comm," FAS function norm %g\n",norm);
1175:         }

1177:     }
1178:     else{
1179: 
1180:         DMMGSolveFASMCycle9(dmmg,level,&converged);
1181:         if (dmmg[level]->monitor){
1182:           DMMGFormFunction(0,dmmg[level]->x,dmmg[level]->w,dmmg[level]);
1183:           VecNorm(dmmg[level]->w,NORM_2,&norm);
1184:           PetscPrintf(dmmg[level]->comm," FAS function norm %g\n",norm);
1185:         }


1188:     }
1189:   }
1190:   else {
1191:     PetscOptionsHasName(PETSC_NULL,"-dmmg_fas_full",&flg);
1192:     if (flg) {
1193: 
1194:         DMMGSolveFASFCycle(dmmg,level,&converged);
1195:         if (dmmg[level]->monitor){
1196:           DMMGFormFunction(0,dmmg[level]->x,dmmg[level]->w,dmmg[level]);
1197:           VecNorm(dmmg[level]->w,NORM_2,&norm);
1198:           PetscPrintf(dmmg[level]->comm," FAS function norm %g\n",norm);
1199:         }

1201:     }
1202:     else{
1203: 
1204:         DMMGSolveFASMCycle(dmmg,level,&converged);
1205:         if (dmmg[level]->monitor){
1206:           DMMGFormFunction(0,dmmg[level]->x,dmmg[level]->w,dmmg[level]);
1207:           VecNorm(dmmg[level]->w,NORM_2,&norm);
1208:           PetscPrintf(dmmg[level]->comm," FAS function norm %g\n",norm);
1209:         }



1213:     }
1214:   }
1215: 

1217:   return(0);
1218: }


1221: /*

1223: This is function to implement Nonlinear CG to accelerate FAS

1225: In order to use this acceleration, the option is

1227: -dmmg_fas_NCG



1231: */


1236: PetscErrorCode DMMGSolveFAS_NCG(DMMG *dmmg, PetscInt level)
1237: {
1238:   SNES           snes = dmmg[level]->snes;
1239:   SNES_LS        *neP = (SNES_LS*)snes->data;
1241:   PetscInt       maxits,i,lits;
1242:   PetscTruth     lssucceed;
1243:   // MatStructure   flg = DIFFERENT_NONZERO_PATTERN;
1244:   PetscReal      fnorm,gnorm,xnorm,ynorm,betaFR,betaPR,beta,betaHS,betaDY;
1245:   Vec            Y,X,F,G,W,Gradold,Sk;
1246:   //KSP            ksp;



1251:  if (!snes->solve) SETERRQ(PETSC_ERR_ORDER,"SNESSetType() or SNESSetFromOptions() must be called before SNESSolve()");

1253:   VecDuplicate(dmmg[level]->x,&Sk);
1254:   snes->vec_sol = snes->vec_sol_always = dmmg[level]->x;
1255:   if (!snes->setupcalled) {
1256:     SNESSetUp(snes);
1257:   }
1258:   if (snes->conv_hist_reset) snes->conv_hist_len = 0;
1260:   snes->nfuncs = 0; snes->linear_its = 0; snes->numFailures = 0;


1263:   snes->reason  = SNES_CONVERGED_ITERATING;

1265:   maxits        = snes->max_its;        /* maximum number of iterations */
1266:   X                = snes->vec_sol;        /* solution vector */
1267:   F                = snes->vec_func;        /* residual vector */
1268:   Y                = snes->work[0];        /* work vectors */
1269:   G                = snes->work[1];
1270:   W                = snes->work[2];

1272:   PetscObjectTakeAccess(snes);
1273:   snes->iter = 0;
1274:   PetscObjectGrantAccess(snes);
1275: 
1276:   X    = dmmg[level]->x;
1277:   VecCopy(X,Y);
1278:   VecCopy(X,G);

1280:   // to get the residual for the F
1281:   SNESComputeFunction(snes,X,F);
1282: 
1283:   VecNorm(F,NORM_2,&fnorm);        /* fnorm <- ||F||  */
1284:   if (fnorm != fnorm) SETERRQ(PETSC_ERR_FP,"User provided compute function generated a Not-a-Number");
1285:   PetscObjectTakeAccess(snes);
1286:   snes->norm = fnorm;
1287:   PetscObjectGrantAccess(snes);
1288:   SNESLogConvHistory(snes,fnorm,0);
1289:   SNESMonitor(snes,0,fnorm);

1291:   if (fnorm < snes->abstol) {snes->reason = SNES_CONVERGED_FNORM_ABS; return(0);}

1293:   /* set parameter for default relative tolerance convergence test */
1294:   snes->ttol = fnorm*snes->rtol;
1295:   // dmmg[level]->rrtol= snes->ttol;

1297:   // set this to store the old grad
1298:   Gradold=snes->vec_sol_update_always;
1299: 
1300:   // compute the search direction Y
1301:   // I need to put Q(x)=x-FAS(x) here
1302:   DMMGSolveFASCyclen(dmmg,level);
1303:   //  F    = X - dmmg[level]->x; this is the gradient direction
1304:   VecAXPY(Y,-1.0,X);
1305:   // copy the gradient to the old
1306:   VecCopy(Y,Gradold);
1307:   // copy X back
1308:   VecCopy(G,X);
1309:   VecWAXPY(Sk,-1.0,X,X);
1310:   // so far I put X=X_c, F= F(x_c),  Gradold= Y=grad(x_c)

1312:   //  for (i=0; i<maxits; i++) {

1314:    for (i=0; i<10000; i++) {


1317:     PetscPrintf(PETSC_COMM_WORLD,"iter=%d",i+1);
1318: 
1319: 
1320:     // X=x_c, F=F(x_c),Y search direction; G=F(x_new), W=x_new,
1321:     VecNorm(X,NORM_2,&xnorm);
1322:     (*neP->LineSearch)(snes,neP->lsP,X,F,G,Y,W,fnorm,&ynorm,&gnorm,&lssucceed);
1323:     PetscInfo4(snes,"SNESSolve_LS: fnorm=%18.16e, gnorm=%18.16e, ynorm=%18.16e, lssucceed=%d\n",fnorm,gnorm,ynorm,(int)lssucceed);
1324:     if (snes->reason == SNES_DIVERGED_FUNCTION_COUNT) break;
1325:     PetscPrintf(PETSC_COMM_WORLD,"step=%g,oldnorm=%g,norm=%g ",ynorm,fnorm,gnorm);
1326: 
1327:     fnorm=gnorm; //copy the new function norm; this will effect the line_search
1328:      VecWAXPY(Sk,-1.0,X,W);
1329:     //update the new solution
1330:     ierr=VecCopy(W,X);
1331:     ierr=VecCopy(G,F);

1333: 
1334:   // compute the new search direction G
1335:   // I need to put Q(x)=x-FAS(x) here

1337:   DMMGSolveFASCyclen(dmmg,level);
1338:   //G    = X - dmmg[level]->x; G is the new gradient, Y is old gradient
1339: 
1340:   VecWAXPY(G,-1.0,X,W);
1341:   // copy W back to X
1342:   VecCopy(W,X);
1343:   VecNorm(G,NORM_2,&gnorm);
1344:   VecNorm(Gradold,NORM_2,&ynorm);
1345:   betaFR = gnorm*gnorm/(ynorm*ynorm); //FR_beta
1346: 
1347:   VecWAXPY(W,-1.0,Gradold,G);
1348:   VecDot(W,G,&gnorm);
1349:   //  VecNorm(G,NORM_2,&gnorm);
1350:   VecNorm(Gradold,NORM_2,&ynorm);
1351:   betaPR = gnorm/(ynorm*ynorm); //PR_beta
1352: 
1353:   if ( betaPR<-betaFR)
1354:     {
1355:       beta =- betaFR;
1356:     }
1357:   else {
1358:     if (betaPR>betaFR)
1359:       {beta=betaFR;}
1360:     else{

1362:       beta=betaPR;
1363:     }
1364:   }
1365:   //  beta=betaFR;
1366:   //beta=betaPR;

1368:   // try another beta
1369: 
1370:      VecWAXPY(W,-1.0,Gradold,G);
1371:      VecDot(W,G,&betaHS);
1372:      VecDot(W,Y,&gnorm);
1373:      betaHS=-betaHS/gnorm;
1374:      VecDot(G,G,&betaDY);
1375:      betaDY=-betaDY/gnorm;
1376:      if(betaHS<betaDY)
1377:        beta=betaHS;
1378:      else
1379:        beta=betaDY;
1380:      if(beta<0)
1381:        beta=0;
1382: 
1383:      PetscPrintf(PETSC_COMM_WORLD,"betaHS=%g,betaDY=%g\n",betaHS,betaDY);
1384: 

1386: 
1387:     // compute the c_2
1388:     VecDot(G,Y,&gnorm);
1389:     VecDot(Gradold,Y,&ynorm);
1390:     PetscPrintf(PETSC_COMM_WORLD,"beta=%g,c_2=%g\n",beta,fabs(gnorm/ynorm));
1391:      VecNorm(G,NORM_2,&gnorm);
1392:     VecNorm(Y,NORM_2,&ynorm);
1393:     PetscPrintf(PETSC_COMM_WORLD,"size=%g\n",fabs(gnorm/(beta*ynorm)));

1395:     // update the direction: Y= G + beta * Y
1396:     VecAXPBY(Y,1.0,beta,G);
1397:     VecCopy(G,Gradold);
1398:     //ierr =VecCopy(G,Y);
1399:     snes->iter = i+1;
1400:     snes->norm = fnorm;
1401:     PetscObjectGrantAccess(snes);
1402:     SNESLogConvHistory(snes,fnorm,lits);
1403:     SNESMonitor(snes,i+1,fnorm);
1404: 
1405:      if (!lssucceed) {
1406:       PetscTruth ismin;
1407:       beta=0;
1408:       if (++snes->numFailures >= snes->maxFailures) {
1409:       snes->reason = SNES_DIVERGED_LS_FAILURE;
1410:         SNESLSCheckLocalMin_Private(snes->jacobian,F,W,fnorm,&ismin);
1411:         if (ismin) snes->reason = SNES_DIVERGED_LOCAL_MIN;
1412:         break;
1413:       }
1414:       }
1415: 

1417:     /* Test for convergence */
1418:     if (snes->converged) {
1419:       VecNorm(X,NORM_2,&xnorm);        /* xnorm = || X || */
1420:       (*snes->converged)(snes,snes->iter,xnorm,1.0,fnorm,&snes->reason,snes->cnvP);
1421:       if (snes->reason) {
1422:         break;
1423:       }
1424:     }
1425:   }
1426:   if (X != snes->vec_sol) {
1427:     VecCopy(X,snes->vec_sol);
1428:   }
1429:   if (F != snes->vec_func) {
1430:     VecCopy(F,snes->vec_func);
1431:   }
1432:   snes->vec_sol_always  = snes->vec_sol;
1433:   snes->vec_func_always = snes->vec_func;
1434:   if (i == maxits) {
1435:     PetscInfo1(snes,"SNESSolve_LS: Maximum number of iterations has been reached: %D\n",maxits);
1436:     snes->reason = SNES_DIVERGED_MAX_IT;
1437:   }
1438:   PetscPrintf(PETSC_COMM_WORLD,"reason=%d\n",snes->reason);
1439:   // VecView(X,PETSC_VIEWER_STDOUT_WORLD);
1440:   return(0);
1441: }



1445: /*

1447: This is function to implement NonGMRES  to accelerate FAS

1449: In order to use this acceleration, the option is

1451: -dmmg_fas_NGMRES


1454: Options:

1456: -dmmg_fas_ngmres_m : the number of previous solutions to keep

1458: */

1462: PetscErrorCode DMMGSolveFAS_NGMRES(DMMG *dmmg, PetscInt level)
1463: {
1464:   SNES           snes = dmmg[level]->snes;
1466:   PetscInt       maxits=10000,i,k,l,j,subm=3,iter;
1467:  PetscOptionsGetInt(PETSC_NULL,"-dmmg_fas_ngmres_m",&subm,PETSC_NULL);

1469:   PetscTruth     restart=PETSC_FALSE, selectA=PETSC_FALSE;
1470:   PetscReal      fnorm,gnorm,dnorm,dnormtemp,dminnorm,fminnorm,tol=1.e-12,gammaA=2,epsilonB=0.1,deltaB=0.9,gammaC;
1471:   Vec            X,F,G,W,D,u[subm],res[subm];
1472:    PetscScalar    H[subm][subm],q[subm][subm],beta[subm],xi[subm],alpha[subm],alphasum,det,Hinv[16];


1476:  gammaC=2; if (gammaA>gammaC) gammaC=gammaA;
1477:  VecDuplicate(dmmg[level]->x,&X);
1478:  VecDuplicate(dmmg[level]->x,&F);
1479:  VecDuplicate(dmmg[level]->x,&W);
1480:  VecDuplicate(dmmg[level]->x,&G);
1481:  VecDuplicate(dmmg[level]->x,&D);

1483:  for(i=0;i<subm;i++) {/* get the space for the solution */
1484:    VecDuplicate(dmmg[level]->x,&u[i]);
1485:    VecDuplicate(dmmg[level]->x,&res[i]);
1486:  }

1488:   X    = dmmg[level]->x;
1489:   VecCopy(X,u[0]);
1490: 
1491:   // to get the residual for the F
1492:   SNESComputeFunction(snes,X,F);
1493:   VecCopy(F,res[0]);
1494:   VecNorm(F,NORM_2,&fnorm);        /* fnorm <- ||F||  */
1495:   fnorm=fnorm*fnorm;
1496:   iter=1;
1497:   restartline:
1498: 
1499:   q[0][0] = fnorm; fminnorm=fnorm;


1502:    for (k=1; k<maxits; k++) {

1504: 
1505:      PetscPrintf(PETSC_COMM_WORLD,"\n k=%d,iter=%d fmin=%g ",k,iter++,sqrt(fminnorm));
1506: 
1507:     /* compute the X=u^M , F=r, fnorm =||F||*/

1509:     DMMGSolveFASCyclen(dmmg,level);
1510:     SNESComputeFunction(snes,X,F);
1511:     VecNorm(F,NORM_2,&fnorm);
1512:     if (fnorm < tol) { return(0);}
1513:     fnorm =fnorm*fnorm;
1514:     if (fnorm<fminnorm) fminnorm=fnorm;
1515:     //PetscPrintf(PETSC_COMM_WORLD,"fmin=%g",sqrt(fminnorm));
1516:     /* compute u^A */
1517:     //l=min(subm,k)
1518:     l=subm;
1519:     if (k<l) l=k;
1520: 
1521:     /* compute the matrix H and RHS */
1522: 
1523: 
1524:     for(i=0;i<l;i++){
1525:       VecDot(F,res[i],&xi[i]);
1526:       beta[i]=fnorm-xi[i];
1527:     }
1528: 
1529:     for(i=0;i<l;i++){
1530:       for(j=0;j<l;j++){
1531:        H[i][j]=q[i][j]-xi[i]-xi[j]+fnorm;
1532:       }
1533:     }
1534:     /* Here is special for subm=2 */
1535:     if(l==1){
1536:       //   H[0][0] = q[0][0]-xi[0]-xi[0]+fnorm;
1537:       alpha[0]= beta[0]/H[0][0];
1538:     }
1539:     if(l==2){
1540: 
1541:       alpha[0]= ( beta[0]*H[1][1] -beta[1]*H[0][1] )/( H[0][0]*H[1][1]- H[0][1]*H[1][0]);
1542:       alpha[1]= ( beta[1]*H[0][0] -beta[0]*H[1][0] )/( H[0][0]*H[1][1]- H[0][1]*H[1][0]);
1543: 
1544:     }
1545:     if(l==3) {
1546: 
1547:   det = H[0][0]*(H[1][1]*H[2][2]-H[1][2]*H[2][1])-H[0][1]*(H[1][0]*H[2][2]-H[1][2]*H[2][0])+H[0][2]*(H[1][0]*H[2][1]-H[1][1]*H[2][0]);
1548:       alpha[0]= (beta[0]*(H[1][1]*H[2][2]-H[1][2]*H[2][1])-H[0][1]*(beta[1]*H[2][2]-H[1][2]*beta[2])+H[0][2]*(beta[1]*H[2][1]-H[1][1]*beta[2]))/det;
1549:       alpha[1]=(H[0][0]*(beta[1]*H[2][2]-H[1][2]*beta[2])-beta[0]*(H[1][0]*H[2][2]-H[1][2]*H[2][0])+H[0][2]*(H[1][0]*beta[2]-beta[1]*H[2][0]))/det;
1550:       alpha[2]=(H[0][0]*(H[1][1]*beta[2]-beta[1]*H[2][1])-H[0][1]*(H[1][0]*beta[2]-beta[1]*H[2][0])+beta[0]*(H[1][0]*H[2][1]-H[1][1]*H[2][0]))/det;


1553:     }
1554: 
1555:     if(l==4){
1556:         Hinv[0]=H[0][0];Hinv[1]=H[0][1];Hinv[2]=H[0][2];Hinv[3]=H[0][3];
1557:       Hinv[4]=H[1][0];Hinv[5]=H[1][1];Hinv[6]=H[1][2];Hinv[7]=H[1][3];
1558:       Hinv[8]=H[2][0];Hinv[9]=H[2][1];Hinv[10]=H[2][2];Hinv[11]=H[2][3];
1559:       Hinv[12]=H[3][0];Hinv[13]=H[3][1];Hinv[14]=H[3][2];Hinv[15]=H[3][3];
1560:       Kernel_A_gets_inverse_A_4(Hinv);
1561:       for(i=0;i<l;i++)
1562:         alpha[i]=Hinv[4*i]*beta[0]+Hinv[4*i+1]*beta[1]+Hinv[4*i+2]*beta[2]+Hinv[4*i+3]*beta[3];
1563: 
1564:     }
1565:     alphasum=0;
1566:     for (i=0;i<l;i++)
1567:       alphasum=alphasum+alpha[i];
1568: 
1569:     /* W= u^A */
1570:     VecCopy(X,W);
1571:     VecAXPBY(W,0.0,1-alphasum,X);
1572: 
1573:     for(i=0;i<l;i++)
1574:          VecAXPY(W,alpha[i],u[i]);
1575: 
1576:     /* W= F(G) */
1577:     SNESComputeFunction(snes,W,G);
1578:     VecNorm(G,NORM_2,&gnorm);
1579:     gnorm=gnorm*gnorm;


1582: 
1583:     /* select the uA or uM */
1584:     // Criterion A
1585:     if(sqrt(gnorm)<gammaA*sqrt(fminnorm)){
1586:       //PetscPrintf(PETSC_COMM_WORLD,"Crite A\n");
1587:        selectA=PETSC_TRUE;
1588:     }
1589:     // Criterion B
1590: 
1591:     ierr=VecCopy(W,D);
1592:     ierr=VecAXPY(D,-1,X);
1593:     ierr=VecNorm(D,NORM_2,&dnorm);
1594:     dminnorm=10000000;
1595:     for(i=0;i<l;i++) {
1596:       ierr=VecCopy(W,D);
1597:       ierr=VecAXPY(D,-1,u[i]);
1598:       ierr=VecNorm(D,NORM_2,&dnormtemp);
1599:       if(dnormtemp<dminnorm) dminnorm=dnormtemp;
1600:     }
1601:      if( epsilonB*dnorm<dminnorm || sqrt(gnorm)<deltaB*sqrt(fminnorm))
1602:       selectA =PETSC_TRUE;
1603:     else
1604:       selectA=PETSC_FALSE;
1605: 
1606:     if(selectA){  /* uA selected */
1607:       selectA=PETSC_FALSE;
1608:       //   PetscPrintf(PETSC_COMM_WORLD,"Select A\n");
1609:       VecCopy(W,X);
1610:       VecCopy(G,F);
1611:       fnorm=gnorm;
1612:     }
1613:      if (fnorm<fminnorm) fminnorm=fnorm;
1614: 

1616:     /* Test for convergence */
1617:     if (sqrt(fnorm)<tol) {
1618:       return(0);
1619:     }
1620: 
1621:     /* Test for restart */
1622:     if(sqrt(gnorm)>gammaC*sqrt(fminnorm)) {
1623:       PetscPrintf(PETSC_COMM_WORLD,"restart for C ");
1624:       restart=PETSC_TRUE;
1625:     }
1626:     if(epsilonB*dnorm>dminnorm && sqrt(gnorm)>deltaB*sqrt(fminnorm)) {
1627:       PetscPrintf(PETSC_COMM_WORLD,"restart for D ");
1628:       restart=PETSC_TRUE;
1629:       }
1630:     /* Prepation for the next iteration */
1631: 
1632:      //turn off restart
1633:     //restart=PETSC_FALSE;
1634:     if(restart){
1635:       restart=PETSC_FALSE;
1636:       goto restartline;
1637:     }
1638:     else {
1639:       j=k%subm;
1640:       VecCopy(F,res[j]);
1641:       VecCopy(X,u[j]);
1642:       for(i=0;i<l;i++){
1643:           ierr= VecDot(F,res[i],&q[j][i]);
1644:         q[i][j]=q[j][i];
1645:       }
1646:       if(l<subm)
1647:       q[j][j]=fnorm;
1648:     }
1649:    }

1651:   return(0);
1652: }

1654: /*
1655: This is a function to provide the function value:

1657: x-FAS(x)

1659: */



1665: PetscErrorCode DMMGFASFunction(SNES snes,Vec x,Vec f,void *ptr)
1666: {
1667:       DMMG*          dmmg=(DMMG*)ptr;
1669:       Vec            temp;
1670:       PetscInt       level=dmmg[0]->nlevels-1;
1672: 
1673:       VecDuplicate(dmmg[level]->x,&temp);
1674: 
1675: 
1676:       VecCopy(dmmg[level]->x,temp);
1677:       VecCopy(x,dmmg[level]->x);
1678:       VecCopy(x,f);
1679: 
1680:       // I need to put -F(x)=x-FAS(x) here
1681:       DMMGSolveFASCyclen(dmmg,level);
1682:       VecAXPY(f,-1.0,dmmg[level]->x);
1683:       // y = alpha x + y.
1684:       ierr=VecCopy(temp,dmmg[level]->x);
1685:       //copy W back to X
1686: 
1687:  return(0);
1688: }



1692: /* this function is to implement Quasi-Newton method with implicit Broyden updating methods(limit memory version)

1694: In order to use this method, the option is -dmmg_fas_QNewton

1696: Options:

1698: -dmmg_fas_QNewton_m: the number of the vectors to keep for inverse of Jacobian

1700: -dmmg_fas_initialJacobian: will use matrix-free GMRES to solve the initial Jacobian 

1702:                         with  options -snes_mf -snes_max_it 1 -ksp_max_it n


1705: In this function, does not have line search and nonlinear gmres acceleration

1707: */

1711: PetscErrorCode DMMGSolveFAS_QNewton(DMMG *dmmg, PetscInt level)
1712: {

1714: 

1716:   SNES           snes = dmmg[level]->snes, snes0;
1718:   PetscInt       maxits=10000,i,k,l,subm=3,subm01;
1719:   PetscOptionsGetInt(PETSC_NULL,"-dmmg_fas_QNewton_m",&subm,PETSC_NULL);
1720:   subm01=subm-1;
1721:    PetscTruth   flg;
1722:    PetscReal      fnorm,gnorm,tol=1.e-12;
1723:   Vec            X,F,G,W,D,Y,v[subm],w[subm],s0,s1,F0,F1;
1724: 


1728:  PetscOptionsHasName(PETSC_NULL,"-dmmg_fas_initialJacobian",&flg);
1729:  VecDuplicate(dmmg[level]->x,&X);
1730:  VecDuplicate(dmmg[level]->x,&F);
1731:  VecDuplicate(dmmg[level]->x,&W);
1732:  VecDuplicate(dmmg[level]->x,&G);
1733:  VecDuplicate(dmmg[level]->x,&D);
1734:  VecDuplicate(dmmg[level]->x,&Y);
1735:  VecDuplicate(dmmg[level]->x,&s0);
1736:  VecDuplicate(dmmg[level]->x,&s1);
1737:  VecDuplicate(dmmg[level]->x,&F0);
1738:  VecDuplicate(dmmg[level]->x,&F1);
1739: 
1740:  // creat a snes for solve the initial Jacobian
1741:  SNESCreate(dmmg[level]->comm,&snes0);
1742:  SNESSetFunction(snes0,F,DMMGFASFunction,dmmg);
1743:  SNESSetFromOptions(snes0);

1745:  for(i=0;i<subm;i++) {/* get the space for the solution */
1746:    VecDuplicate(dmmg[level]->x,&v[i]);
1747:    VecDuplicate(dmmg[level]->x,&w[i]);
1748:  }

1750:  //We first try B0==I
1751:    X    = dmmg[level]->x;
1752: 
1753:    if(flg){
1754:      ierr= VecAXPBY(Y,0.0,0.0,X);
1755:      ierr= VecCopy(X,s0);
1756:      ierr= SNESSolve(snes0,Y,s0);
1757:      ierr= VecAXPY(s0,-1.0,X);
1758:    }
1759:    else{
1760:      ierr=VecCopy(X,W);
1761:      ierr=VecCopy(X,Y);
1762: 
1763:      // I need to put -F(x)=x-FAS(x) here

1765:      DMMGSolveFASCyclen(dmmg,level);
1766:      VecAXPY(Y,-1.0,X);
1767:      // y = alpha x + y.
1768:      ierr=VecCopy(W,X);
1769:      //copy W back to X
1770: 
1771:      // Y stores the -F(x)
1772:      ierr= VecAXPBY(Y,0.0,-1.0,X);
1773:      ierr= VecCopy(Y,s0);

1775:    }
1776: 
1777:    VecAXPY(X,1.0,s0);
1778: 
1779: for(k=0; k<maxits; k++){
1780: 
1781:  /* Test for convergence */
1782:    SNESComputeFunction(snes,X,F);
1783:    VecNorm(F,NORM_2,&fnorm);        /* fnorm <- ||F||  */
1784:    PetscPrintf(PETSC_COMM_WORLD,"k=%d, fnorm=%g\n",
1785:                        k,fnorm);
1786:     if (fnorm<tol) {
1787:       return(0);
1788:     }
1789: 

1791:     if(flg){
1792:       //     ierr= SNESSolve(snes0,Y,s1);
1793:       ierr= VecAXPBY(Y,0.0,0.0,X);
1794:       ierr= VecCopy(dmmg[level]->x,s1);
1795:       ierr= SNESSolve(snes0,Y,s1);
1796:       ierr= VecAXPY(s1,-1.0,X);
1797: 
1798:     }
1799:     else{
1800:       ierr=VecCopy(X,W);
1801:       ierr=VecCopy(X,Y);
1802: 
1803:       // I need to put -F(x)=x-FAS(x) here

1805:       DMMGSolveFASCyclen(dmmg,level);
1806:       VecAXPY(Y,-1.0,X);
1807:       // y = alpha x + y.
1808:       ierr=VecCopy(W,X);
1809:       //copy W back to X
1810: 
1811:       //So far, I got X=x_k, Y=-F(x_k)
1812:       // I should solve the G=-B_0^{-1}F(x_k) first, but I choose B_0^{-1}=I,
1813:       ierr=VecCopy(Y,F1);
1814:       ierr= VecAXPBY(Y,0.0,-1.0,X);
1815:       ierr=VecCopy(Y,s1);

1817:     }
1818: 
1819:    l=subm;
1820:    if (k<l) l=k;
1821: 
1822: 
1823:    for (i=0;i<l;i++){
1824:      // compute [I+v(i)w(i)^T]*s(k)
1825:      ierr= VecDot(w[i],s1,&gnorm);
1826:      VecAXPY(s1,gnorm,v[i]);
1827: 
1828:    }
1829:    if(l==subm) {
1830:      for(i=0;i<subm01;i++){
1831:        ierr= VecCopy(w[i+1],w[i]);
1832:        ierr= VecCopy(v[i+1],v[i]);
1833:      }
1834:      l--;
1835:    }
1836:      ierr= VecCopy(s0,w[l]);
1837:      ierr= VecCopy(s0,Y);
1838:      ierr= VecCopy(s1,v[l]);
1839:      ierr= VecAXPY(Y,-1.0,s1);
1840:      ierr= VecDot(w[l],Y,&gnorm);
1841:      ierr= VecAXPBY(v[l],0.0,1.0/gnorm,w[l]);

1843:      ierr= VecDot(s1,w[l],&gnorm);
1844:      ierr= VecAXPY(s1,gnorm,v[l]);
1845:      ierr= VecCopy(s1,s0);
1846:      ierr=VecAXPY(X,1.0,s1);
1847:  }
1848:   return(0);
1849: }