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