1 #include <../src/snes/impls/ngmres/snesngmres.h> /*I "petscsnes.h" I*/ 2 #include <petscblaslapack.h> 3 4 PetscErrorCode SNESNGMRESUpdateSubspace_Private(SNES snes,PetscInt ivec,PetscInt l,Vec F,PetscReal fnorm,Vec X) 5 { 6 SNES_NGMRES *ngmres = (SNES_NGMRES*) snes->data; 7 Vec *Fdot = ngmres->Fdot; 8 Vec *Xdot = ngmres->Xdot; 9 PetscErrorCode ierr; 10 11 PetscFunctionBegin; 12 if (ivec > l) SETERRQ2(PetscObjectComm((PetscObject)snes),PETSC_ERR_ARG_WRONGSTATE,"Cannot update vector %D with space size %D!",ivec,l); 13 ierr = VecCopy(F,Fdot[ivec]);CHKERRQ(ierr); 14 ierr = VecCopy(X,Xdot[ivec]);CHKERRQ(ierr); 15 16 ngmres->fnorms[ivec] = fnorm; 17 PetscFunctionReturn(0); 18 } 19 20 PetscErrorCode SNESNGMRESFormCombinedSolution_Private(SNES snes,PetscInt ivec,PetscInt l,Vec XM,Vec FM,PetscReal fMnorm,Vec X,Vec XA,Vec FA) 21 { 22 SNES_NGMRES *ngmres = (SNES_NGMRES*) snes->data; 23 PetscInt i,j; 24 Vec *Fdot = ngmres->Fdot; 25 Vec *Xdot = ngmres->Xdot; 26 PetscScalar *beta = ngmres->beta; 27 PetscScalar *xi = ngmres->xi; 28 PetscScalar alph_total = 0.; 29 PetscErrorCode ierr; 30 PetscReal nu; 31 Vec Y = snes->work[2]; 32 PetscBool changed_y,changed_w; 33 34 PetscFunctionBegin; 35 nu = fMnorm*fMnorm; 36 37 /* construct the right hand side and xi factors */ 38 if (l > 0) { 39 ierr = VecMDotBegin(FM,l,Fdot,xi);CHKERRQ(ierr); 40 ierr = VecMDotBegin(Fdot[ivec],l,Fdot,beta);CHKERRQ(ierr); 41 ierr = VecMDotEnd(FM,l,Fdot,xi);CHKERRQ(ierr); 42 ierr = VecMDotEnd(Fdot[ivec],l,Fdot,beta);CHKERRQ(ierr); 43 for (i = 0; i < l; i++) { 44 Q(i,ivec) = beta[i]; 45 Q(ivec,i) = beta[i]; 46 } 47 } else { 48 Q(0,0) = ngmres->fnorms[ivec]*ngmres->fnorms[ivec]; 49 } 50 51 for (i = 0; i < l; i++) beta[i] = nu - xi[i]; 52 53 /* construct h */ 54 for (j = 0; j < l; j++) { 55 for (i = 0; i < l; i++) { 56 H(i,j) = Q(i,j)-xi[i]-xi[j]+nu; 57 } 58 } 59 if (l == 1) { 60 /* simply set alpha[0] = beta[0] / H[0, 0] */ 61 if (H(0,0) != 0.) beta[0] = beta[0]/H(0,0); 62 else beta[0] = 0.; 63 } else { 64 ierr = PetscBLASIntCast(l,&ngmres->m);CHKERRQ(ierr); 65 ierr = PetscBLASIntCast(l,&ngmres->n);CHKERRQ(ierr); 66 ngmres->info = 0; 67 ngmres->rcond = -1.; 68 ierr = PetscFPTrapPush(PETSC_FP_TRAP_OFF);CHKERRQ(ierr); 69 #if defined(PETSC_USE_COMPLEX) 70 PetscStackCallBLAS("LAPACKgelss",LAPACKgelss_(&ngmres->m,&ngmres->n,&ngmres->nrhs,ngmres->h,&ngmres->lda,ngmres->beta,&ngmres->ldb,ngmres->s,&ngmres->rcond,&ngmres->rank,ngmres->work,&ngmres->lwork,ngmres->rwork,&ngmres->info)); 71 #else 72 PetscStackCallBLAS("LAPACKgelss",LAPACKgelss_(&ngmres->m,&ngmres->n,&ngmres->nrhs,ngmres->h,&ngmres->lda,ngmres->beta,&ngmres->ldb,ngmres->s,&ngmres->rcond,&ngmres->rank,ngmres->work,&ngmres->lwork,&ngmres->info)); 73 #endif 74 ierr = PetscFPTrapPop();CHKERRQ(ierr); 75 if (ngmres->info < 0) SETERRQ(PetscObjectComm((PetscObject)snes),PETSC_ERR_LIB,"Bad argument to GELSS"); 76 if (ngmres->info > 0) SETERRQ(PetscObjectComm((PetscObject)snes),PETSC_ERR_LIB,"SVD failed to converge"); 77 } 78 for (i=0; i<l; i++) { 79 if (PetscIsInfOrNanScalar(beta[i])) SETERRQ(PetscObjectComm((PetscObject)snes),PETSC_ERR_LIB,"SVD generated inconsistent output"); 80 } 81 alph_total = 0.; 82 for (i = 0; i < l; i++) alph_total += beta[i]; 83 84 ierr = VecCopy(XM,XA);CHKERRQ(ierr); 85 ierr = VecScale(XA,1.-alph_total);CHKERRQ(ierr); 86 ierr = VecMAXPY(XA,l,beta,Xdot);CHKERRQ(ierr); 87 /* check the validity of the step */ 88 ierr = VecCopy(XA,Y);CHKERRQ(ierr); 89 ierr = VecAXPY(Y,-1.0,X);CHKERRQ(ierr); 90 ierr = SNESLineSearchPostCheck(snes->linesearch,X,Y,XA,&changed_y,&changed_w);CHKERRQ(ierr); 91 if (!ngmres->approxfunc) { 92 if (snes->npc && snes->npcside== PC_LEFT) { 93 ierr = SNESApplyNPC(snes,XA,NULL,FA);CHKERRQ(ierr); 94 } else { 95 ierr = SNESComputeFunction(snes,XA,FA);CHKERRQ(ierr); 96 } 97 } else { 98 ierr = VecCopy(FM,FA);CHKERRQ(ierr); 99 ierr = VecScale(FA,1.-alph_total);CHKERRQ(ierr); 100 ierr = VecMAXPY(FA,l,beta,Fdot);CHKERRQ(ierr); 101 } 102 PetscFunctionReturn(0); 103 } 104 105 PetscErrorCode SNESNGMRESNorms_Private(SNES snes,PetscInt l,Vec X,Vec F,Vec XM,Vec FM,Vec XA,Vec FA,Vec D,PetscReal *dnorm,PetscReal *dminnorm,PetscReal *xMnorm,PetscReal *fMnorm,PetscReal *yMnorm, PetscReal *xAnorm,PetscReal *fAnorm,PetscReal *yAnorm) 106 { 107 PetscErrorCode ierr; 108 SNES_NGMRES *ngmres = (SNES_NGMRES*) snes->data; 109 PetscReal dcurnorm,dmin = -1.0; 110 Vec *Xdot = ngmres->Xdot; 111 PetscInt i; 112 113 PetscFunctionBegin; 114 if (xMnorm) { 115 ierr = VecNormBegin(XM,NORM_2,xMnorm);CHKERRQ(ierr); 116 } 117 if (fMnorm) { 118 ierr = VecNormBegin(FM,NORM_2,fMnorm);CHKERRQ(ierr); 119 } 120 if (yMnorm) { 121 ierr = VecCopy(X,D);CHKERRQ(ierr); 122 ierr = VecAXPY(D,-1.0,XM);CHKERRQ(ierr); 123 ierr = VecNormBegin(D,NORM_2,yMnorm);CHKERRQ(ierr); 124 } 125 if (xAnorm) { 126 ierr = VecNormBegin(XA,NORM_2,xAnorm);CHKERRQ(ierr); 127 } 128 if (fAnorm) { 129 ierr = VecNormBegin(FA,NORM_2,fAnorm);CHKERRQ(ierr); 130 } 131 if (yAnorm) { 132 ierr = VecCopy(X,D);CHKERRQ(ierr); 133 ierr = VecAXPY(D,-1.0,XA);CHKERRQ(ierr); 134 ierr = VecNormBegin(D,NORM_2,yAnorm);CHKERRQ(ierr); 135 } 136 if (dnorm) { 137 ierr = VecCopy(XA,D);CHKERRQ(ierr); 138 ierr = VecAXPY(D,-1.0,XM);CHKERRQ(ierr); 139 ierr = VecNormBegin(D,NORM_2,dnorm);CHKERRQ(ierr); 140 } 141 if (dminnorm) { 142 for (i=0; i<l; i++) { 143 ierr = VecCopy(Xdot[i],D);CHKERRQ(ierr); 144 ierr = VecAXPY(D,-1.0,XA);CHKERRQ(ierr); 145 ierr = VecNormBegin(D,NORM_2,&ngmres->xnorms[i]);CHKERRQ(ierr); 146 } 147 } 148 if (xMnorm) {ierr = VecNormEnd(XM,NORM_2,xMnorm);CHKERRQ(ierr);} 149 if (fMnorm) {ierr = VecNormEnd(FM,NORM_2,fMnorm);CHKERRQ(ierr);} 150 if (yMnorm) {ierr = VecNormEnd(D,NORM_2,yMnorm);CHKERRQ(ierr);} 151 if (xAnorm) {ierr = VecNormEnd(XA,NORM_2,xAnorm);CHKERRQ(ierr);} 152 if (fAnorm) {ierr = VecNormEnd(FA,NORM_2,fAnorm);CHKERRQ(ierr);} 153 if (yAnorm) {ierr = VecNormEnd(D,NORM_2,yAnorm);CHKERRQ(ierr);} 154 if (dnorm) {ierr = VecNormEnd(D,NORM_2,dnorm);CHKERRQ(ierr);} 155 if (dminnorm) { 156 for (i=0; i<l; i++) { 157 ierr = VecNormEnd(D,NORM_2,&ngmres->xnorms[i]);CHKERRQ(ierr); 158 dcurnorm = ngmres->xnorms[i]; 159 if ((dcurnorm < dmin) || (dmin < 0.0)) dmin = dcurnorm; 160 } 161 *dminnorm = dmin; 162 } 163 PetscFunctionReturn(0); 164 } 165 166 PetscErrorCode SNESNGMRESSelect_Private(SNES snes,PetscInt k_restart,Vec XM,Vec FM,PetscReal xMnorm,PetscReal fMnorm,PetscReal yMnorm,Vec XA,Vec FA,PetscReal xAnorm,PetscReal fAnorm,PetscReal yAnorm,PetscReal dnorm,PetscReal fminnorm,PetscReal dminnorm,Vec X,Vec F,Vec Y,PetscReal *xnorm,PetscReal *fnorm,PetscReal *ynorm) 167 { 168 SNES_NGMRES *ngmres = (SNES_NGMRES*) snes->data; 169 PetscErrorCode ierr; 170 SNESLineSearchReason lssucceed; 171 PetscBool selectA; 172 173 PetscFunctionBegin; 174 if (ngmres->select_type == SNES_NGMRES_SELECT_LINESEARCH) { 175 /* X = X + \lambda(XA - X) */ 176 if (ngmres->monitor) { 177 ierr = PetscViewerASCIIPrintf(ngmres->monitor,"||F_A||_2 = %e, ||F_M||_2 = %e\n",fAnorm,fMnorm);CHKERRQ(ierr); 178 } 179 ierr = VecCopy(FM,F);CHKERRQ(ierr); 180 ierr = VecCopy(XM,X);CHKERRQ(ierr); 181 ierr = VecCopy(XA,Y);CHKERRQ(ierr); 182 ierr = VecAYPX(Y,-1.0,X);CHKERRQ(ierr); 183 *fnorm = fMnorm; 184 ierr = SNESLineSearchApply(ngmres->additive_linesearch,X,F,fnorm,Y);CHKERRQ(ierr); 185 ierr = SNESLineSearchGetReason(ngmres->additive_linesearch,&lssucceed);CHKERRQ(ierr); 186 ierr = SNESLineSearchGetNorms(ngmres->additive_linesearch,xnorm,fnorm,ynorm);CHKERRQ(ierr); 187 if (lssucceed) { 188 if (++snes->numFailures >= snes->maxFailures) { 189 snes->reason = SNES_DIVERGED_LINE_SEARCH; 190 PetscFunctionReturn(0); 191 } 192 } 193 if (ngmres->monitor) { 194 ierr = PetscViewerASCIIPrintf(ngmres->monitor,"Additive solution: ||F||_2 = %e\n",*fnorm);CHKERRQ(ierr); 195 } 196 } else if (ngmres->select_type == SNES_NGMRES_SELECT_DIFFERENCE) { 197 selectA = PETSC_TRUE; 198 /* Conditions for choosing the accelerated answer */ 199 /* Criterion A -- the norm of the function isn't increased above the minimum by too much */ 200 if (fAnorm >= ngmres->gammaA*fminnorm) selectA = PETSC_FALSE; 201 202 /* Criterion B -- the choice of x^A isn't too close to some other choice */ 203 if (ngmres->epsilonB*dnorm<dminnorm || PetscSqrtReal(*fnorm)<ngmres->deltaB*PetscSqrtReal(fminnorm)) { 204 } else selectA=PETSC_FALSE; 205 206 if (selectA) { 207 if (ngmres->monitor) { 208 ierr = PetscViewerASCIIPrintf(ngmres->monitor,"picked X_A, ||F_A||_2 = %e, ||F_M||_2 = %e\n",fAnorm,fMnorm);CHKERRQ(ierr); 209 } 210 /* copy it over */ 211 *xnorm = xAnorm; 212 *fnorm = fAnorm; 213 *ynorm = yAnorm; 214 ierr = VecCopy(FA,F);CHKERRQ(ierr); 215 ierr = VecCopy(XA,X);CHKERRQ(ierr); 216 } else { 217 if (ngmres->monitor) { 218 ierr = PetscViewerASCIIPrintf(ngmres->monitor,"picked X_M, ||F_A||_2 = %e, ||F_M||_2 = %e\n",fAnorm,fMnorm);CHKERRQ(ierr); 219 } 220 *xnorm = xMnorm; 221 *fnorm = fMnorm; 222 *ynorm = yMnorm; 223 ierr = VecCopy(XM,Y);CHKERRQ(ierr); 224 ierr = VecAXPY(Y,-1.0,X);CHKERRQ(ierr); 225 ierr = VecCopy(FM,F);CHKERRQ(ierr); 226 ierr = VecCopy(XM,X);CHKERRQ(ierr); 227 } 228 } else { /* none */ 229 *xnorm = xAnorm; 230 *fnorm = fAnorm; 231 *ynorm = yAnorm; 232 ierr = VecCopy(FA,F);CHKERRQ(ierr); 233 ierr = VecCopy(XA,X);CHKERRQ(ierr); 234 } 235 PetscFunctionReturn(0); 236 } 237 238 PetscErrorCode SNESNGMRESSelectRestart_Private(SNES snes,PetscInt l,PetscReal fMnorm, PetscReal fAnorm,PetscReal dnorm,PetscReal fminnorm,PetscReal dminnorm,PetscBool *selectRestart) 239 { 240 SNES_NGMRES *ngmres = (SNES_NGMRES*)snes->data; 241 PetscErrorCode ierr; 242 243 PetscFunctionBegin; 244 *selectRestart = PETSC_FALSE; 245 /* difference stagnation restart */ 246 if ((ngmres->epsilonB*dnorm > dminnorm) && (PetscSqrtReal(fAnorm) > ngmres->deltaB*PetscSqrtReal(fminnorm)) && l > 0) { 247 if (ngmres->monitor) { 248 ierr = PetscViewerASCIIPrintf(ngmres->monitor,"difference restart: %e > %e\n",ngmres->epsilonB*dnorm,dminnorm);CHKERRQ(ierr); 249 } 250 *selectRestart = PETSC_TRUE; 251 } 252 /* residual stagnation restart */ 253 if (PetscSqrtReal(fAnorm) > ngmres->gammaC*PetscSqrtReal(fminnorm)) { 254 if (ngmres->monitor) { 255 ierr = PetscViewerASCIIPrintf(ngmres->monitor,"residual restart: %e > %e\n",PetscSqrtReal(fAnorm),ngmres->gammaC*PetscSqrtReal(fminnorm));CHKERRQ(ierr); 256 } 257 *selectRestart = PETSC_TRUE; 258 } 259 260 /* F_M stagnation restart */ 261 if (ngmres->restart_fm_rise && fMnorm > snes->norm) { 262 if (ngmres->monitor) { 263 ierr = PetscViewerASCIIPrintf(ngmres->monitor,"F_M rise restart: %e > %e\n",fMnorm,snes->norm);CHKERRQ(ierr); 264 } 265 *selectRestart = PETSC_TRUE; 266 } 267 268 PetscFunctionReturn(0); 269 } 270