1 #include <../src/snes/impls/gs/gsimpl.h> 2 3 PETSC_EXTERN PetscErrorCode SNESComputeNGSDefaultSecant(SNES snes,Vec X,Vec B,void *ctx) 4 { 5 PetscErrorCode ierr; 6 SNES_NGS *gs = (SNES_NGS*)snes->data; 7 PetscInt i,j,k,ncolors; 8 DM dm; 9 PetscBool flg; 10 ISColoring coloring = gs->coloring; 11 MatColoring mc; 12 Vec W,G,F; 13 PetscScalar h=gs->h; 14 IS *coloris; 15 PetscScalar f,g,x,w,d; 16 PetscReal dxt,xt,ft,ft1=0; 17 const PetscInt *idx; 18 PetscInt size,s; 19 PetscReal atol,rtol,stol; 20 PetscInt its; 21 PetscErrorCode (*func)(SNES,Vec,Vec,void*); 22 void *fctx; 23 PetscBool mat = gs->secant_mat,equal,isdone,alldone; 24 PetscScalar *xa,*fa,*wa,*ga; 25 26 PetscFunctionBegin; 27 if (snes->nwork < 3) { 28 ierr = SNESSetWorkVecs(snes,3);CHKERRQ(ierr); 29 } 30 W = snes->work[0]; 31 G = snes->work[1]; 32 F = snes->work[2]; 33 ierr = VecGetOwnershipRange(X,&s,NULL);CHKERRQ(ierr); 34 ierr = SNESNGSGetTolerances(snes,&atol,&rtol,&stol,&its);CHKERRQ(ierr); 35 ierr = SNESGetDM(snes,&dm);CHKERRQ(ierr); 36 ierr = SNESGetFunction(snes,NULL,&func,&fctx);CHKERRQ(ierr); 37 if (!coloring) { 38 /* create the coloring */ 39 ierr = DMHasColoring(dm,&flg);CHKERRQ(ierr); 40 if (flg && !mat) { 41 ierr = DMCreateColoring(dm,IS_COLORING_GLOBAL,&coloring);CHKERRQ(ierr); 42 } else { 43 if (!snes->jacobian) {ierr = SNESSetUpMatrices(snes);CHKERRQ(ierr);} 44 ierr = MatColoringCreate(snes->jacobian,&mc);CHKERRQ(ierr); 45 ierr = MatColoringSetDistance(mc,1);CHKERRQ(ierr); 46 ierr = MatColoringSetFromOptions(mc);CHKERRQ(ierr); 47 ierr = MatColoringApply(mc,&coloring);CHKERRQ(ierr); 48 ierr = MatColoringDestroy(&mc);CHKERRQ(ierr); 49 } 50 gs->coloring = coloring; 51 } 52 ierr = ISColoringGetIS(coloring,&ncolors,&coloris);CHKERRQ(ierr); 53 ierr = VecEqual(X,snes->vec_sol,&equal);CHKERRQ(ierr); 54 if (equal && snes->normschedule == SNES_NORM_ALWAYS) { 55 /* assume that the function is already computed */ 56 ierr = VecCopy(snes->vec_func,F);CHKERRQ(ierr); 57 } else { 58 ierr = PetscLogEventBegin(SNES_NGSFuncEval,snes,X,B,0);CHKERRQ(ierr); 59 ierr = (*func)(snes,X,F,fctx);CHKERRQ(ierr); 60 ierr = PetscLogEventEnd(SNES_NGSFuncEval,snes,X,B,0);CHKERRQ(ierr); 61 if (B) {ierr = VecAXPY(F,-1.0,B);CHKERRQ(ierr);} 62 } 63 ierr = VecGetArray(X,&xa);CHKERRQ(ierr); 64 ierr = VecGetArray(F,&fa);CHKERRQ(ierr); 65 ierr = VecGetArray(G,&ga);CHKERRQ(ierr); 66 ierr = VecGetArray(W,&wa);CHKERRQ(ierr); 67 for (i=0;i<ncolors;i++) { 68 ierr = ISGetIndices(coloris[i],&idx);CHKERRQ(ierr); 69 ierr = ISGetLocalSize(coloris[i],&size);CHKERRQ(ierr); 70 ierr = VecCopy(X,W);CHKERRQ(ierr); 71 for (j=0;j<size;j++) { 72 wa[idx[j]-s] += h; 73 } 74 ierr = PetscLogEventBegin(SNES_NGSFuncEval,snes,X,B,0);CHKERRQ(ierr); 75 ierr = (*func)(snes,W,G,fctx);CHKERRQ(ierr); 76 ierr = PetscLogEventEnd(SNES_NGSFuncEval,snes,X,B,0);CHKERRQ(ierr); 77 if (B) {ierr = VecAXPY(G,-1.0,B);CHKERRQ(ierr);} 78 for (k=0;k<its;k++) { 79 dxt = 0.; 80 xt = 0.; 81 ft = 0.; 82 for (j=0;j<size;j++) { 83 f = fa[idx[j]-s]; 84 x = xa[idx[j]-s]; 85 g = ga[idx[j]-s]; 86 w = wa[idx[j]-s]; 87 if (PetscAbsScalar(g-f) > atol) { 88 /* This is equivalent to d = x - (h*f) / PetscRealPart(g-f) */ 89 d = (x*g-w*f) / PetscRealPart(g-f); 90 } else { 91 d = x; 92 } 93 dxt += PetscRealPart(PetscSqr(d-x)); 94 xt += PetscRealPart(PetscSqr(x)); 95 ft += PetscRealPart(PetscSqr(f)); 96 xa[idx[j]-s] = d; 97 } 98 99 if (k == 0) ft1 = PetscSqrtReal(ft); 100 if (k<its-1) { 101 isdone = PETSC_FALSE; 102 if (stol*PetscSqrtReal(xt) > PetscSqrtReal(dxt)) isdone = PETSC_TRUE; 103 if (PetscSqrtReal(ft) < atol) isdone = PETSC_TRUE; 104 if (rtol*ft1 > PetscSqrtReal(ft)) isdone = PETSC_TRUE; 105 ierr = MPIU_Allreduce(&isdone,&alldone,1,MPIU_BOOL,MPI_BAND,PetscObjectComm((PetscObject)snes));CHKERRQ(ierr); 106 if (alldone) break; 107 } 108 if (i < ncolors-1 || k < its-1) { 109 ierr = PetscLogEventBegin(SNES_NGSFuncEval,snes,X,B,0);CHKERRQ(ierr); 110 ierr = (*func)(snes,X,F,fctx);CHKERRQ(ierr); 111 ierr = PetscLogEventEnd(SNES_NGSFuncEval,snes,X,B,0);CHKERRQ(ierr); 112 if (B) {ierr = VecAXPY(F,-1.0,B);CHKERRQ(ierr);} 113 } 114 if (k<its-1) { 115 ierr = VecSwap(X,W);CHKERRQ(ierr); 116 ierr = VecSwap(F,G);CHKERRQ(ierr); 117 } 118 } 119 } 120 ierr = VecRestoreArray(X,&xa);CHKERRQ(ierr); 121 ierr = VecRestoreArray(F,&fa);CHKERRQ(ierr); 122 ierr = VecRestoreArray(G,&ga);CHKERRQ(ierr); 123 ierr = VecRestoreArray(W,&wa);CHKERRQ(ierr); 124 ierr = ISColoringRestoreIS(coloring,&coloris);CHKERRQ(ierr); 125 PetscFunctionReturn(0); 126 } 127