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