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