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