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: }