Actual source code: ex43-44.h
1: /*
2: Accelerates Newton's method by solving a small problem defined by those elements with large residual plus one level of overlap
4: This is a toy code for playing around
6: Counts residual entries as small if they are less then .2 times the maximum
7: Decides to solve a reduced problem if the number of large entries is less than 20 percent of all entries (and this has been true for criteria_reduce iterations)
8: */
9: #include petscsnes.h
15: typedef struct {
16: Vec xwork,fwork;
17: VecScatter scatter;
18: SNES snes;
19: IS is;
20: } SubCtx;
24: PetscErrorCode FormFunctionSub(SNES snes,Vec x,Vec f,void *ictx)
25: {
27: SubCtx *ctx = (SubCtx*) ictx;
30: VecScatterBegin(ctx->scatter,x,ctx->xwork,INSERT_VALUES,SCATTER_REVERSE);
31: VecScatterEnd(ctx->scatter,x,ctx->xwork,INSERT_VALUES,SCATTER_REVERSE);
32: SNESComputeFunction(ctx->snes,ctx->xwork,ctx->fwork);
33: VecScatterBegin(ctx->scatter,ctx->fwork,f,INSERT_VALUES,SCATTER_FORWARD);
34: VecScatterEnd(ctx->scatter,ctx->fwork,f,INSERT_VALUES,SCATTER_FORWARD);
35: return(0);
36: }
40: PetscErrorCode FormJacobianSub(SNES snes,Vec x,Mat *A, Mat *B, MatStructure *str,void *ictx)
41: {
43: SubCtx *ctx = (SubCtx*) ictx;
44: Mat As,Bs;
47: VecScatterBegin(ctx->scatter,x,ctx->xwork,INSERT_VALUES,SCATTER_REVERSE);
48: VecScatterEnd(ctx->scatter,x,ctx->xwork,INSERT_VALUES,SCATTER_REVERSE);
49: SNESGetJacobian(ctx->snes,&As,&Bs,PETSC_NULL,PETSC_NULL);
50: SNESComputeJacobian(ctx->snes,ctx->xwork,&As,&Bs,str);
51: if (*B) {
52: MatGetSubMatrix(Bs,ctx->is,ctx->is,MAT_REUSE_MATRIX,B);
53: } else {
54: MatGetSubMatrix(Bs,ctx->is,ctx->is,MAT_INITIAL_MATRIX,B);
55: }
56: if (!*A) {
57: *A = *B;
58: PetscObjectReference((PetscObject)*A);
59: }
60: if (*A != *B) {
61: MatAssemblyBegin(*A,MAT_FINAL_ASSEMBLY);
62: MatAssemblyEnd(*A,MAT_FINAL_ASSEMBLY);
63: }
64: return(0);
65: }
69: PetscErrorCode SolveSubproblem(SNES snes)
70: {
72: Vec residual,solution;
73: PetscReal rmax;
74: PetscInt i,n,cnt,*indices;
75: PetscScalar *r;
76: SNES snessub;
77: Vec x,f;
78: SubCtx ctx;
79: Mat mat;
82: ctx.snes = snes;
83: SNESGetSolution(snes,&solution);
84: SNESGetFunction(snes,&residual,0,0);
85: VecNorm(residual,NORM_INFINITY,&rmax);
86: VecGetLocalSize(residual,&n);
87: VecGetArray(residual,&r);
88: cnt = 0;
89: for (i=0; i<n; i++) {
90: if (PetscAbsScalar(r[i]) > .20*rmax ) cnt++;
91: }
92: PetscMalloc(cnt*sizeof(PetscInt),&indices);
93: cnt = 0;
94: for (i=0; i<n; i++) {
95: if (PetscAbsScalar(r[i]) > .20*rmax ) indices[cnt++] = i;
96: }
97: if (cnt > .2*n) return(0);
99: printf("number in subproblem %d\n",cnt);
100: VecRestoreArray(residual,&r);
101: ISCreateGeneral(PETSC_COMM_WORLD,cnt,indices,&ctx.is);
102: PetscFree(indices);
104: SNESGetJacobian(snes,0,&mat,0,0);
105: MatIncreaseOverlap(mat,1,&ctx.is,2);
106: ISSort(ctx.is);
107: ISGetLocalSize(ctx.is,&cnt);
108: printf("number in subproblem %d\n",cnt);
109: VecCreate(PETSC_COMM_WORLD,&x);
110: VecSetSizes(x,cnt,cnt);
111: VecSetType(x,VECSEQ);
112: VecDuplicate(x,&f);
113: VecScatterCreate(solution,ctx.is,x,PETSC_NULL,&ctx.scatter);
115: VecDuplicate(solution,&ctx.xwork);
116: VecCopy(solution,ctx.xwork);
117: VecDuplicate(residual,&ctx.fwork);
118: VecScatterBegin(ctx.scatter,solution,x,INSERT_VALUES,SCATTER_FORWARD);
119: VecScatterEnd(ctx.scatter,solution,x,INSERT_VALUES,SCATTER_FORWARD);
120: SNESCreate(PETSC_COMM_WORLD,&snessub);
121: SNESSetOptionsPrefix(snessub,"sub_");
122: SNESSetFunction(snessub,f,FormFunctionSub,(void*)&ctx);
123: SNESSetJacobian(snessub,0,0,FormJacobianSub,(void*)&ctx);
124: SNESSetFromOptions(snessub);
125: SNESSolve(snessub,PETSC_NULL,x);
126: VecScatterBegin(ctx.scatter,x,solution,INSERT_VALUES,SCATTER_REVERSE);
127: VecScatterEnd(ctx.scatter,x,solution,INSERT_VALUES,SCATTER_REVERSE);
129: ISDestroy(ctx.is);
130: VecDestroy(ctx.xwork);
131: VecDestroy(ctx.fwork);
132: VecScatterDestroy(ctx.scatter);
133: VecDestroy(x);
134: VecDestroy(f);
135: SNESDestroy(snessub);
136: return(0);
137: }
141: static PetscInt CountGood = 0;
145: PetscErrorCode MonitorRange(SNES snes,PetscInt it,PetscReal rnorm,void *dummy)
146: {
147: PetscErrorCode ierr;
148: PetscReal perc;
150: SNESMonitorRange_Private(snes,it,&perc);
151: if (perc < .20) CountGood++;
152: else CountGood = 0;
153: return(0);
154: }