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