1 2 #include <private/snesimpl.h> 3 4 typedef struct { 5 PetscReal lambda; /* The default step length for the update */ 6 Vec * dX; /* The change in X */ 7 Vec * dF; /* The change in F */ 8 PetscInt m; /* the number of kept previous steps */ 9 } QNContext; 10 11 #undef __FUNCT__ 12 #define __FUNCT__ "LBGFSApplyJinv_Private" 13 PetscErrorCode LBGFSApplyJinv_Private(SNES snes, PetscInt it, Vec g, Vec z) { 14 15 PetscErrorCode ierr; 16 17 QNContext * qn = (QNContext *)snes->data; 18 19 Vec * dX = qn->dX; 20 Vec * dF = qn->dF; 21 22 PetscInt k, i; 23 PetscInt m = qn->m; 24 PetscScalar * alpha = PETSC_NULL; 25 PetscScalar * beta = PETSC_NULL; 26 PetscScalar * rho = PETSC_NULL; 27 PetscScalar t; 28 PetscInt l = m; 29 30 PetscFunctionBegin; 31 32 if (it < m) l = it+1; 33 ierr = PetscMalloc3(m, PetscScalar, &alpha, m, PetscScalar, &beta, m, PetscScalar, &rho);CHKERRQ(ierr); 34 35 /* precalculate alpha, beta, rho corresponding to the normal indices*/ 36 for (i = 0; i < l; i++) { 37 ierr = VecDot(dX[i], dF[i], &t);CHKERRQ(ierr); 38 rho[i] = 1. / t; 39 beta[i] = 0; 40 alpha[i] = 0; 41 } 42 43 ierr = VecCopy(g, z);CHKERRQ(ierr); 44 45 /* outward recursion starting at iteration k's update and working back */ 46 for (i = 0; i < l; i++) { 47 k = (it - i) % m; 48 ierr = VecDot(dX[k], z, &t);CHKERRQ(ierr); 49 alpha[k] = t / rho[k]; 50 ierr = VecAXPY(z, -alpha[k], dF[k]);CHKERRQ(ierr); 51 } 52 53 /* inner application of the initial inverse jacobian approximation */ 54 /* right now it's just the identity. Nothing needs to go here. */ 55 56 /* inward recursion starting at the first update and working forward*/ 57 for (i = l - 1; i >= 0; i--) { 58 k = (it - i) % m; 59 ierr = VecDot(dF[k], z, &t);CHKERRQ(ierr); 60 beta[k] = rho[k]*t; 61 ierr = VecAXPY(z, (alpha[k] - beta[k]), dX[k]); 62 } 63 ierr = PetscFree3(alpha, beta, rho);CHKERRQ(ierr); 64 65 PetscFunctionReturn(0); 66 } 67 68 #undef __FUNCT__ 69 #define __FUNCT__ "SNESSolve_QN" 70 static PetscErrorCode SNESSolve_QN(SNES snes) 71 { 72 73 PetscErrorCode ierr; 74 QNContext * qn = (QNContext*) snes->data; 75 76 Vec x; 77 Vec f; 78 Vec p, pold; 79 80 PetscInt i, j, k, l; 81 82 PetscReal fnorm; 83 PetscScalar gdot; 84 PetscInt m = qn->m; 85 86 Vec * W = qn->dX; 87 Vec * V = qn->dF; 88 89 /* basically just a regular newton's method except for the application of the jacobian */ 90 PetscFunctionBegin; 91 92 x = snes->vec_sol; 93 f = snes->vec_func; 94 p = snes->vec_sol_update; 95 pold = snes->work[0]; 96 97 snes->reason = SNES_CONVERGED_ITERATING; 98 99 ierr = PetscObjectTakeAccess(snes);CHKERRQ(ierr); 100 snes->iter = 0; 101 snes->norm = 0.; 102 ierr = PetscObjectGrantAccess(snes);CHKERRQ(ierr); 103 ierr = SNESComputeFunction(snes,x,f);CHKERRQ(ierr); 104 if (snes->domainerror) { 105 snes->reason = SNES_DIVERGED_FUNCTION_DOMAIN; 106 PetscFunctionReturn(0); 107 } 108 ierr = VecNorm(f, NORM_2, &fnorm);CHKERRQ(ierr); /* fnorm <- ||F|| */ 109 if (PetscIsInfOrNanReal(fnorm)) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FP,"Infinite or not-a-number generated in norm"); 110 ierr = PetscObjectTakeAccess(snes);CHKERRQ(ierr); 111 snes->norm = fnorm; 112 ierr = PetscObjectGrantAccess(snes);CHKERRQ(ierr); 113 SNESLogConvHistory(snes,fnorm,0); 114 ierr = SNESMonitor(snes,0,fnorm);CHKERRQ(ierr); 115 116 /* set parameter for default relative tolerance convergence test */ 117 snes->ttol = fnorm*snes->rtol; 118 /* test convergence */ 119 ierr = (*snes->ops->converged)(snes,0,0.0,0.0,fnorm,&snes->reason,snes->cnvP);CHKERRQ(ierr); 120 if (snes->reason) PetscFunctionReturn(0); 121 ierr = VecCopy(f, pold);CHKERRQ(ierr); 122 ierr = VecAXPY(x, -1.0, pold);CHKERRQ(ierr); 123 for(i = 0; i < snes->max_its; i++) { 124 125 ierr = SNESComputeFunction(snes, x, f);CHKERRQ(ierr); 126 ierr = VecNorm(f, NORM_2, &fnorm);CHKERRQ(ierr); 127 if (PetscIsInfOrNanReal(fnorm)) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FP,"Infinite or not-a-number generated in norm"); 128 ierr = PetscObjectTakeAccess(snes);CHKERRQ(ierr); 129 snes->norm = fnorm; 130 ierr = PetscObjectGrantAccess(snes);CHKERRQ(ierr); 131 SNESLogConvHistory(snes,fnorm,i+1); 132 ierr = SNESMonitor(snes,i+1,fnorm);CHKERRQ(ierr); 133 134 /* set parameter for default relative tolerance convergence test */ 135 snes->ttol = fnorm*snes->rtol; 136 /* test convergence */ 137 ierr = (*snes->ops->converged)(snes,i+1,0.0,0.0,fnorm,&snes->reason,snes->cnvP);CHKERRQ(ierr); 138 if (snes->reason) PetscFunctionReturn(0); 139 k = (i) % m; 140 l = m; 141 if (i < l) l = i; 142 ierr = VecCopy(f, p);CHKERRQ(ierr); 143 for (j=0; j<k; j++) { /* p = product_{j<i} [I+v(j)w(j)^T]*p */ 144 ierr = VecDot(W[j],p,&gdot);CHKERRQ(ierr); 145 ierr = VecAXPY(p,gdot,V[j]);CHKERRQ(ierr); 146 } 147 ierr = VecCopy(pold,W[k]);CHKERRQ(ierr); /* w[i] = pold */ 148 ierr = VecAXPY(pold,-1.0,p);CHKERRQ(ierr); /* v[i] = p */ 149 ierr = VecDot(W[k],pold,&gdot);CHKERRQ(ierr); /* ----------------- */ 150 ierr = VecCopy(p,V[k]);CHKERRQ(ierr); /* w[i]'*(Pold - p) */ 151 ierr = VecScale(V[k],1.0/gdot);CHKERRQ(ierr); 152 153 ierr = VecDot(W[k],p,&gdot);CHKERRQ(ierr); /* p = (I + v[i]*w[i]')*p */ 154 ierr = VecAXPY(p,gdot,V[k]);CHKERRQ(ierr); 155 ierr = VecCopy(p,pold);CHKERRQ(ierr); 156 157 ierr = VecAXPY(x, -1.0, p);CHKERRQ(ierr); 158 } 159 if (i == snes->max_its) { 160 ierr = PetscInfo1(snes, "Maximum number of iterations has been reached: %D\n", snes->max_its);CHKERRQ(ierr); 161 if (!snes->reason) snes->reason = SNES_DIVERGED_MAX_IT; 162 } 163 PetscFunctionReturn(0); 164 } 165 166 167 #undef __FUNCT__ 168 #define __FUNCT__ "SNESSetUp_QN" 169 static PetscErrorCode SNESSetUp_QN(SNES snes) 170 { 171 QNContext * qn = (QNContext *)snes->data; 172 PetscErrorCode ierr; 173 PetscFunctionBegin; 174 ierr = VecDuplicateVecs(snes->vec_sol, qn->m, &qn->dX);CHKERRQ(ierr); 175 ierr = VecDuplicateVecs(snes->vec_sol, qn->m, &qn->dF);CHKERRQ(ierr); 176 ierr = SNESDefaultGetWork(snes,1);CHKERRQ(ierr); 177 PetscFunctionReturn(0); 178 } 179 180 #undef __FUNCT__ 181 #define __FUNCT__ "SNESReset_QN" 182 static PetscErrorCode SNESReset_QN(SNES snes) 183 { 184 PetscErrorCode ierr; 185 QNContext * qn; 186 PetscFunctionBegin; 187 if (snes->data) { 188 qn = (QNContext *)snes->data; 189 if (qn->dX) { 190 ierr = VecDestroyVecs(qn->m, &qn->dX);CHKERRQ(ierr); 191 } 192 if (qn->dF) { 193 ierr = VecDestroyVecs(qn->m, &qn->dF);CHKERRQ(ierr); 194 } 195 } 196 if (snes->work) {ierr = VecDestroyVecs(snes->nwork,&snes->work);CHKERRQ(ierr);} 197 PetscFunctionReturn(0); 198 } 199 200 #undef __FUNCT__ 201 #define __FUNCT__ "SNESDestroy_QN" 202 static PetscErrorCode SNESDestroy_QN(SNES snes) 203 { 204 PetscErrorCode ierr; 205 PetscFunctionBegin; 206 ierr = SNESReset_QN(snes);CHKERRQ(ierr); 207 ierr = PetscFree(snes->data);CHKERRQ(ierr); 208 PetscFunctionReturn(0); 209 } 210 211 #undef __FUNCT__ 212 #define __FUNCT__ "SNESSetFromOptions_QN" 213 static PetscErrorCode SNESSetFromOptions_QN(SNES snes) 214 { 215 216 PetscErrorCode ierr; 217 QNContext * qn; 218 219 PetscFunctionBegin; 220 221 qn = (QNContext *)snes->data; 222 223 ierr = PetscOptionsHead("SNES QN options");CHKERRQ(ierr); 224 ierr = PetscOptionsReal("-snes_qn_lambda", "SOR mixing parameter", "SNES", qn->lambda, &qn->lambda, PETSC_NULL);CHKERRQ(ierr); 225 ierr = PetscOptionsInt("-snes_qn_m", "Number of past states saved for L-Broyden methods", "SNES", qn->m, &qn->m, PETSC_NULL);CHKERRQ(ierr); 226 ierr = PetscOptionsTail();CHKERRQ(ierr); 227 PetscFunctionReturn(0); 228 } 229 230 /* -------------------------------------------------------------------------- */ 231 /*MC 232 SNESQN - Limited-Memory Quasi-Newton methods for the solution of nonlinear systems. 233 234 Implements a limited-memory "good" Broyden update method. 235 236 Level: beginner 237 238 .seealso: SNESCreate(), SNES, SNESSetType(), SNESLS, SNESTR 239 M*/ 240 EXTERN_C_BEGIN 241 #undef __FUNCT__ 242 #define __FUNCT__ "SNESCreate_QN" 243 PetscErrorCode SNESCreate_QN(SNES snes) 244 { 245 246 PetscErrorCode ierr; 247 QNContext * qn; 248 249 PetscFunctionBegin; 250 snes->ops->setup = SNESSetUp_QN; 251 snes->ops->solve = SNESSolve_QN; 252 snes->ops->destroy = SNESDestroy_QN; 253 snes->ops->setfromoptions = SNESSetFromOptions_QN; 254 snes->ops->view = 0; 255 snes->ops->reset = SNESReset_QN; 256 257 ierr = PetscNewLog(snes, QNContext, &qn);CHKERRQ(ierr); 258 snes->data = (void *) qn; 259 qn->m = 10; 260 qn->lambda = 1.; 261 qn->dX = PETSC_NULL; 262 qn->dF = PETSC_NULL; 263 PetscFunctionReturn(0); 264 } 265 EXTERN_C_END 266