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