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 PetscCheck(ivec <= l,PetscObjectComm((PetscObject)snes),PETSC_ERR_ARG_WRONGSTATE,"Cannot update vector %" PetscInt_FMT " with space size %" PetscInt_FMT "!",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 PetscCallBLAS("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 PetscCallBLAS("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 PetscCheck(ngmres->info >= 0,PetscObjectComm((PetscObject)snes),PETSC_ERR_LIB,"Bad argument to GELSS"); 74 PetscCheck(ngmres->info <= 0,PetscObjectComm((PetscObject)snes),PETSC_ERR_LIB,"SVD failed to converge"); 75 } 76 for (i=0; i<l; i++) { 77 PetscCheck(!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) PetscCall(VecNormBegin(XM,NORM_2,xMnorm)); 112 if (fMnorm) PetscCall(VecNormBegin(FM,NORM_2,fMnorm)); 113 if (yMnorm) { 114 PetscCall(VecCopy(X,D)); 115 PetscCall(VecAXPY(D,-1.0,XM)); 116 PetscCall(VecNormBegin(D,NORM_2,yMnorm)); 117 } 118 if (xAnorm) PetscCall(VecNormBegin(XA,NORM_2,xAnorm)); 119 if (fAnorm) PetscCall(VecNormBegin(FA,NORM_2,fAnorm)); 120 if (yAnorm) { 121 PetscCall(VecCopy(X,D)); 122 PetscCall(VecAXPY(D,-1.0,XA)); 123 PetscCall(VecNormBegin(D,NORM_2,yAnorm)); 124 } 125 if (dnorm) { 126 PetscCall(VecCopy(XA,D)); 127 PetscCall(VecAXPY(D,-1.0,XM)); 128 PetscCall(VecNormBegin(D,NORM_2,dnorm)); 129 } 130 if (dminnorm) { 131 for (i=0; i<l; i++) { 132 PetscCall(VecCopy(Xdot[i],D)); 133 PetscCall(VecAXPY(D,-1.0,XA)); 134 PetscCall(VecNormBegin(D,NORM_2,&ngmres->xnorms[i])); 135 } 136 } 137 if (xMnorm) PetscCall(VecNormEnd(XM,NORM_2,xMnorm)); 138 if (fMnorm) PetscCall(VecNormEnd(FM,NORM_2,fMnorm)); 139 if (yMnorm) PetscCall(VecNormEnd(D,NORM_2,yMnorm)); 140 if (xAnorm) PetscCall(VecNormEnd(XA,NORM_2,xAnorm)); 141 if (fAnorm) PetscCall(VecNormEnd(FA,NORM_2,fAnorm)); 142 if (yAnorm) PetscCall(VecNormEnd(D,NORM_2,yAnorm)); 143 if (dnorm) PetscCall(VecNormEnd(D,NORM_2,dnorm)); 144 if (dminnorm) { 145 for (i=0; i<l; i++) { 146 PetscCall(VecNormEnd(D,NORM_2,&ngmres->xnorms[i])); 147 dcurnorm = ngmres->xnorms[i]; 148 if ((dcurnorm < dmin) || (dmin < 0.0)) dmin = dcurnorm; 149 } 150 *dminnorm = dmin; 151 } 152 PetscFunctionReturn(0); 153 } 154 155 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) 156 { 157 SNES_NGMRES *ngmres = (SNES_NGMRES*) snes->data; 158 SNESLineSearchReason lssucceed; 159 PetscBool selectA; 160 161 PetscFunctionBegin; 162 if (ngmres->select_type == SNES_NGMRES_SELECT_LINESEARCH) { 163 /* X = X + \lambda(XA - X) */ 164 if (ngmres->monitor) { 165 PetscCall(PetscViewerASCIIPrintf(ngmres->monitor,"||F_A||_2 = %e, ||F_M||_2 = %e\n",(double)fAnorm,(double)fMnorm)); 166 } 167 PetscCall(VecCopy(FM,F)); 168 PetscCall(VecCopy(XM,X)); 169 PetscCall(VecCopy(XA,Y)); 170 PetscCall(VecAYPX(Y,-1.0,X)); 171 *fnorm = fMnorm; 172 PetscCall(SNESLineSearchApply(ngmres->additive_linesearch,X,F,fnorm,Y)); 173 PetscCall(SNESLineSearchGetReason(ngmres->additive_linesearch,&lssucceed)); 174 PetscCall(SNESLineSearchGetNorms(ngmres->additive_linesearch,xnorm,fnorm,ynorm)); 175 if (lssucceed) { 176 if (++snes->numFailures >= snes->maxFailures) { 177 snes->reason = SNES_DIVERGED_LINE_SEARCH; 178 PetscFunctionReturn(0); 179 } 180 } 181 if (ngmres->monitor) { 182 PetscCall(PetscViewerASCIIPrintf(ngmres->monitor,"Additive solution: ||F||_2 = %e\n",(double)*fnorm)); 183 } 184 } else if (ngmres->select_type == SNES_NGMRES_SELECT_DIFFERENCE) { 185 selectA = PETSC_TRUE; 186 /* Conditions for choosing the accelerated answer */ 187 /* Criterion A -- the norm of the function isn't increased above the minimum by too much */ 188 if (fAnorm >= ngmres->gammaA*fminnorm) selectA = PETSC_FALSE; 189 190 /* Criterion B -- the choice of x^A isn't too close to some other choice */ 191 if (ngmres->epsilonB*dnorm<dminnorm || PetscSqrtReal(*fnorm)<ngmres->deltaB*PetscSqrtReal(fminnorm)) { 192 } else selectA=PETSC_FALSE; 193 194 if (selectA) { 195 if (ngmres->monitor) { 196 PetscCall(PetscViewerASCIIPrintf(ngmres->monitor,"picked X_A, ||F_A||_2 = %e, ||F_M||_2 = %e\n",(double)fAnorm,(double)fMnorm)); 197 } 198 /* copy it over */ 199 *xnorm = xAnorm; 200 *fnorm = fAnorm; 201 *ynorm = yAnorm; 202 PetscCall(VecCopy(FA,F)); 203 PetscCall(VecCopy(XA,X)); 204 } else { 205 if (ngmres->monitor) { 206 PetscCall(PetscViewerASCIIPrintf(ngmres->monitor,"picked X_M, ||F_A||_2 = %e, ||F_M||_2 = %e\n",(double)fAnorm,(double)fMnorm)); 207 } 208 *xnorm = xMnorm; 209 *fnorm = fMnorm; 210 *ynorm = yMnorm; 211 PetscCall(VecCopy(XM,Y)); 212 PetscCall(VecAXPY(Y,-1.0,X)); 213 PetscCall(VecCopy(FM,F)); 214 PetscCall(VecCopy(XM,X)); 215 } 216 } else { /* none */ 217 *xnorm = xAnorm; 218 *fnorm = fAnorm; 219 *ynorm = yAnorm; 220 PetscCall(VecCopy(FA,F)); 221 PetscCall(VecCopy(XA,X)); 222 } 223 PetscFunctionReturn(0); 224 } 225 226 PetscErrorCode SNESNGMRESSelectRestart_Private(SNES snes,PetscInt l,PetscReal fMnorm, PetscReal fAnorm,PetscReal dnorm,PetscReal fminnorm,PetscReal dminnorm,PetscBool *selectRestart) 227 { 228 SNES_NGMRES *ngmres = (SNES_NGMRES*)snes->data; 229 230 PetscFunctionBegin; 231 *selectRestart = PETSC_FALSE; 232 /* difference stagnation restart */ 233 if ((ngmres->epsilonB*dnorm > dminnorm) && (PetscSqrtReal(fAnorm) > ngmres->deltaB*PetscSqrtReal(fminnorm)) && l > 0) { 234 if (ngmres->monitor) { 235 PetscCall(PetscViewerASCIIPrintf(ngmres->monitor,"difference restart: %e > %e\n",(double)(ngmres->epsilonB*dnorm),(double)dminnorm)); 236 } 237 *selectRestart = PETSC_TRUE; 238 } 239 /* residual stagnation restart */ 240 if (PetscSqrtReal(fAnorm) > ngmres->gammaC*PetscSqrtReal(fminnorm)) { 241 if (ngmres->monitor) { 242 PetscCall(PetscViewerASCIIPrintf(ngmres->monitor,"residual restart: %e > %e\n",(double)PetscSqrtReal(fAnorm),(double)(ngmres->gammaC*PetscSqrtReal(fminnorm)))); 243 } 244 *selectRestart = PETSC_TRUE; 245 } 246 247 /* F_M stagnation restart */ 248 if (ngmres->restart_fm_rise && fMnorm > snes->norm) { 249 if (ngmres->monitor) { 250 PetscCall(PetscViewerASCIIPrintf(ngmres->monitor,"F_M rise restart: %e > %e\n",(double)fMnorm,(double)snes->norm)); 251 } 252 *selectRestart = PETSC_TRUE; 253 } 254 255 PetscFunctionReturn(0); 256 } 257