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