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