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