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