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 #if defined(PETSC_MISSING_LAPACK_GELSS) 65 SETERRQ(PetscObjectComm((PetscObject)snes),PETSC_ERR_SUP,"NGMRES with LS requires the LAPACK GELSS routine."); 66 #else 67 ierr = PetscBLASIntCast(l,&ngmres->m);CHKERRQ(ierr); 68 ierr = PetscBLASIntCast(l,&ngmres->n);CHKERRQ(ierr); 69 ngmres->info = 0; 70 ngmres->rcond = -1.; 71 ierr = PetscFPTrapPush(PETSC_FP_TRAP_OFF);CHKERRQ(ierr); 72 #if defined(PETSC_USE_COMPLEX) 73 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)); 74 #else 75 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)); 76 #endif 77 ierr = PetscFPTrapPop();CHKERRQ(ierr); 78 if (ngmres->info < 0) SETERRQ(PetscObjectComm((PetscObject)snes),PETSC_ERR_LIB,"Bad argument to GELSS"); 79 if (ngmres->info > 0) SETERRQ(PetscObjectComm((PetscObject)snes),PETSC_ERR_LIB,"SVD failed to converge"); 80 #endif 81 } 82 for (i=0; i<l; i++) { 83 if (PetscIsInfOrNanScalar(beta[i])) SETERRQ(PetscObjectComm((PetscObject)snes),PETSC_ERR_LIB,"SVD generated inconsistent output"); 84 } 85 alph_total = 0.; 86 for (i = 0; i < l; i++) alph_total += beta[i]; 87 88 ierr = VecCopy(XM,XA);CHKERRQ(ierr); 89 ierr = VecScale(XA,1.-alph_total);CHKERRQ(ierr); 90 ierr = VecMAXPY(XA,l,beta,Xdot);CHKERRQ(ierr); 91 /* check the validity of the step */ 92 ierr = VecCopy(XA,Y);CHKERRQ(ierr); 93 ierr = VecAXPY(Y,-1.0,X);CHKERRQ(ierr); 94 ierr = SNESLineSearchPostCheck(snes->linesearch,X,Y,XA,&changed_y,&changed_w);CHKERRQ(ierr); 95 if (!ngmres->approxfunc) { 96 if (snes->npc && snes->npcside== PC_LEFT) { 97 ierr = SNESApplyNPC(snes,XA,NULL,FA);CHKERRQ(ierr); 98 } else { 99 ierr =SNESComputeFunction(snes,XA,FA);CHKERRQ(ierr); 100 } 101 } else { 102 ierr = VecCopy(FM,FA);CHKERRQ(ierr); 103 ierr = VecScale(FA,1.-alph_total);CHKERRQ(ierr); 104 ierr = VecMAXPY(FA,l,beta,Fdot);CHKERRQ(ierr); 105 } 106 PetscFunctionReturn(0); 107 } 108 109 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) 110 { 111 PetscErrorCode ierr; 112 SNES_NGMRES *ngmres = (SNES_NGMRES*) snes->data; 113 PetscReal dcurnorm,dmin = -1.0; 114 Vec *Xdot = ngmres->Xdot; 115 PetscInt i; 116 117 PetscFunctionBegin; 118 if (xMnorm) { 119 ierr = VecNormBegin(XM,NORM_2,xMnorm);CHKERRQ(ierr); 120 } 121 if (fMnorm) { 122 ierr = VecNormBegin(FM,NORM_2,fMnorm);CHKERRQ(ierr); 123 } 124 if (yMnorm) { 125 ierr = VecCopy(X,D);CHKERRQ(ierr); 126 ierr = VecAXPY(D,-1.0,XM);CHKERRQ(ierr); 127 ierr = VecNormBegin(D,NORM_2,yMnorm);CHKERRQ(ierr); 128 } 129 if (xAnorm) { 130 ierr = VecNormBegin(XA,NORM_2,xAnorm);CHKERRQ(ierr); 131 } 132 if (fAnorm) { 133 ierr = VecNormBegin(FA,NORM_2,fAnorm);CHKERRQ(ierr); 134 } 135 if (yAnorm) { 136 ierr = VecCopy(X,D);CHKERRQ(ierr); 137 ierr = VecAXPY(D,-1.0,XA);CHKERRQ(ierr); 138 ierr = VecNormBegin(D,NORM_2,yAnorm);CHKERRQ(ierr); 139 } 140 if (dnorm) { 141 ierr = VecCopy(XA,D);CHKERRQ(ierr); 142 ierr = VecAXPY(D,-1.0,XM);CHKERRQ(ierr); 143 ierr = VecNormBegin(D,NORM_2,dnorm);CHKERRQ(ierr); 144 } 145 if (dminnorm) { 146 for (i=0; i<l; i++) { 147 ierr = VecCopy(Xdot[i],D);CHKERRQ(ierr); 148 ierr = VecAXPY(D,-1.0,XA);CHKERRQ(ierr); 149 ierr = VecNormBegin(D,NORM_2,&ngmres->xnorms[i]);CHKERRQ(ierr); 150 } 151 } 152 if (xMnorm) {ierr = VecNormEnd(XM,NORM_2,xMnorm);CHKERRQ(ierr);} 153 if (fMnorm) {ierr = VecNormEnd(FM,NORM_2,fMnorm);CHKERRQ(ierr);} 154 if (yMnorm) {ierr = VecNormEnd(D,NORM_2,yMnorm);CHKERRQ(ierr);} 155 if (xAnorm) {ierr = VecNormEnd(XA,NORM_2,xAnorm);CHKERRQ(ierr);} 156 if (fAnorm) {ierr = VecNormEnd(FA,NORM_2,fAnorm);CHKERRQ(ierr);} 157 if (yAnorm) {ierr = VecNormEnd(D,NORM_2,yAnorm);CHKERRQ(ierr);} 158 if (dnorm) {ierr = VecNormEnd(D,NORM_2,dnorm);CHKERRQ(ierr);} 159 if (dminnorm) { 160 for (i=0; i<l; i++) { 161 ierr = VecNormEnd(D,NORM_2,&ngmres->xnorms[i]);CHKERRQ(ierr); 162 dcurnorm = ngmres->xnorms[i]; 163 if ((dcurnorm < dmin) || (dmin < 0.0)) dmin = dcurnorm; 164 } 165 *dminnorm = dmin; 166 } 167 PetscFunctionReturn(0); 168 } 169 170 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) 171 { 172 SNES_NGMRES *ngmres = (SNES_NGMRES*) snes->data; 173 PetscErrorCode ierr; 174 SNESLineSearchReason lssucceed; 175 PetscBool selectA; 176 177 PetscFunctionBegin; 178 if (ngmres->select_type == SNES_NGMRES_SELECT_LINESEARCH) { 179 /* X = X + \lambda(XA - X) */ 180 if (ngmres->monitor) { 181 ierr = PetscViewerASCIIPrintf(ngmres->monitor,"||F_A||_2 = %e, ||F_M||_2 = %e\n",fAnorm,fMnorm);CHKERRQ(ierr); 182 } 183 ierr = VecCopy(FM,F);CHKERRQ(ierr); 184 ierr = VecCopy(XM,X);CHKERRQ(ierr); 185 ierr = VecCopy(XA,Y);CHKERRQ(ierr); 186 ierr = VecAYPX(Y,-1.0,X);CHKERRQ(ierr); 187 *fnorm = fMnorm; 188 ierr = SNESLineSearchApply(ngmres->additive_linesearch,X,F,fnorm,Y);CHKERRQ(ierr); 189 ierr = SNESLineSearchGetReason(ngmres->additive_linesearch,&lssucceed);CHKERRQ(ierr); 190 ierr = SNESLineSearchGetNorms(ngmres->additive_linesearch,xnorm,fnorm,ynorm);CHKERRQ(ierr); 191 if (lssucceed) { 192 if (++snes->numFailures >= snes->maxFailures) { 193 snes->reason = SNES_DIVERGED_LINE_SEARCH; 194 PetscFunctionReturn(0); 195 } 196 } 197 if (ngmres->monitor) { 198 ierr = PetscViewerASCIIPrintf(ngmres->monitor,"Additive solution: ||F||_2 = %e\n",*fnorm);CHKERRQ(ierr); 199 } 200 } else if (ngmres->select_type == SNES_NGMRES_SELECT_DIFFERENCE) { 201 selectA = PETSC_TRUE; 202 /* Conditions for choosing the accelerated answer */ 203 /* Criterion A -- the norm of the function isn't increased above the minimum by too much */ 204 if (fAnorm >= ngmres->gammaA*fminnorm) selectA = PETSC_FALSE; 205 206 /* Criterion B -- the choice of x^A isn't too close to some other choice */ 207 if (ngmres->epsilonB*dnorm<dminnorm || PetscSqrtReal(*fnorm)<ngmres->deltaB*PetscSqrtReal(fminnorm)) { 208 } else selectA=PETSC_FALSE; 209 210 if (selectA) { 211 if (ngmres->monitor) { 212 ierr = PetscViewerASCIIPrintf(ngmres->monitor,"picked X_A, ||F_A||_2 = %e, ||F_M||_2 = %e\n",fAnorm,fMnorm);CHKERRQ(ierr); 213 } 214 /* copy it over */ 215 *xnorm = xAnorm; 216 *fnorm = fAnorm; 217 *ynorm = yAnorm; 218 ierr = VecCopy(FA,F);CHKERRQ(ierr); 219 ierr = VecCopy(XA,X);CHKERRQ(ierr); 220 } else { 221 if (ngmres->monitor) { 222 ierr = PetscViewerASCIIPrintf(ngmres->monitor,"picked X_M, ||F_A||_2 = %e, ||F_M||_2 = %e\n",fAnorm,fMnorm);CHKERRQ(ierr); 223 } 224 *xnorm = xMnorm; 225 *fnorm = fMnorm; 226 *ynorm = yMnorm; 227 ierr = VecCopy(XM,Y);CHKERRQ(ierr); 228 ierr = VecAXPY(Y,-1.0,X);CHKERRQ(ierr); 229 ierr = VecCopy(FM,F);CHKERRQ(ierr); 230 ierr = VecCopy(XM,X);CHKERRQ(ierr); 231 } 232 } else { /* none */ 233 *xnorm = xAnorm; 234 *fnorm = fAnorm; 235 *ynorm = yAnorm; 236 ierr = VecCopy(FA,F);CHKERRQ(ierr); 237 ierr = VecCopy(XA,X);CHKERRQ(ierr); 238 } 239 PetscFunctionReturn(0); 240 } 241 242 PetscErrorCode SNESNGMRESSelectRestart_Private(SNES snes,PetscInt l,PetscReal fMnorm, PetscReal fAnorm,PetscReal dnorm,PetscReal fminnorm,PetscReal dminnorm,PetscBool *selectRestart) 243 { 244 SNES_NGMRES *ngmres = (SNES_NGMRES*)snes->data; 245 PetscErrorCode ierr; 246 247 PetscFunctionBegin; 248 *selectRestart = PETSC_FALSE; 249 /* difference stagnation restart */ 250 if ((ngmres->epsilonB*dnorm > dminnorm) && (PetscSqrtReal(fAnorm) > ngmres->deltaB*PetscSqrtReal(fminnorm)) && l > 0) { 251 if (ngmres->monitor) { 252 ierr = PetscViewerASCIIPrintf(ngmres->monitor,"difference restart: %e > %e\n",ngmres->epsilonB*dnorm,dminnorm);CHKERRQ(ierr); 253 } 254 *selectRestart = PETSC_TRUE; 255 } 256 /* residual stagnation restart */ 257 if (PetscSqrtReal(fAnorm) > ngmres->gammaC*PetscSqrtReal(fminnorm)) { 258 if (ngmres->monitor) { 259 ierr = PetscViewerASCIIPrintf(ngmres->monitor,"residual restart: %e > %e\n",PetscSqrtReal(fAnorm),ngmres->gammaC*PetscSqrtReal(fminnorm));CHKERRQ(ierr); 260 } 261 *selectRestart = PETSC_TRUE; 262 } 263 264 /* F_M stagnation restart */ 265 if (ngmres->restart_fm_rise && fMnorm > snes->norm) { 266 if (ngmres->monitor) { 267 ierr = PetscViewerASCIIPrintf(ngmres->monitor,"F_M rise restart: %e > %e\n",fMnorm,snes->norm);CHKERRQ(ierr); 268 } 269 *selectRestart = PETSC_TRUE; 270 } 271 272 PetscFunctionReturn(0); 273 } 274