OOPS
DRPCGMinimizer.h
Go to the documentation of this file.
1 /*
2  * (C) Copyright 2009-2016 ECMWF.
3  *
4  * This software is licensed under the terms of the Apache Licence Version 2.0
5  * which can be obtained at http://www.apache.org/licenses/LICENSE-2.0.
6  * In applying this licence, ECMWF does not waive the privileges and immunities
7  * granted to it by virtue of its status as an intergovernmental organisation nor
8  * does it submit to any jurisdiction.
9  */
10 
11 #ifndef OOPS_ASSIMILATION_DRPCGMINIMIZER_H_
12 #define OOPS_ASSIMILATION_DRPCGMINIMIZER_H_
13 
14 #include <cmath>
15 #include <string>
16 #include <vector>
17 
18 #include "eckit/config/Configuration.h"
19 #include "oops/assimilation/BMatrix.h"
20 #include "oops/assimilation/ControlIncrement.h"
21 #include "oops/assimilation/CostFunction.h"
22 #include "oops/assimilation/DRMinimizer.h"
23 #include "oops/assimilation/HtRinvHMatrix.h"
24 #include "oops/assimilation/QNewtonLMP.h"
25 #include "oops/util/dot_product.h"
26 #include "oops/util/formats.h"
27 #include "oops/util/Logger.h"
28 
29 namespace oops {
30 
31 /// DRPCG Minimizer
32 /*!
33  * \brief Derber-Rosati Preconditioned Conjugate Gradients solver.
34  *
35  * This solver is based on the standard Preconditioned Conjugate
36  * Gradients solver for Ax=b (G. H. Golub and C. F. Van Loan, Matrix
37  * Computations), and on the Derber and Rosati double
38  * PCG algorithm (J. Derber and A. Rosati, 1989, J. Phys. Oceanog. 1333-1347).
39  * For details see S. Gurol, PhD Manuscript, 2013.
40  * It solves \f$ Ax=b\f$ for the particular case \f$ A=B^{-1}+C\f$,
41  * without requiring the application of \f$ B^{-1}\f$. This algorithm
42  * is similar to DRIPCG except it includes standard PCG instead IPCG
43  * and stopping criteria is based on the preconditioner norm.
44  *
45  * A must be square, symmetric, positive definite.
46  *
47  * A preconditioner must be supplied that, given a vector q, returns an
48  * approximation to \f$ (AB)^{-1} q\f$. Possible preconditioning
49  * is detailed in S. Gurol, PhD Manuscript, 2013.
50  * Note that the traditional \f$ B\f$-preconditioning corresponds to
51  * precond=\f$I\f$.
52  *
53  * On entry:
54  * - dx = starting point, \f$ dx_{0} \f$.
55  * - dxh = \f$ B^{-1} dx_{0} \f$.
56  * - rr = \f$ (sum B^-1 dx^{b}_{i} + ) H^T R^{-1} d \f$
57  * - B = \f$ B \f$.
58  * - C = \f$ C \f$.
59  * - precond = preconditioner \f$ F_k \approx (AB)^{-1} \f$.
60  *
61  * On exit, dx will contain the solution \f$ dx \f$ and dxh will contain
62  * \f$ B^{-1} dx\f$.
63  * The return value is the achieved reduction in preconditioned residual norm.
64  *
65  * Iteration will stop if the maximum iteration limit "maxiter" is reached
66  * or if the residual norm reduces by a factor of "tolerance".
67  *
68  * Each matrix must implement a method:
69  * - void multiply(const VECTOR&, VECTOR&) const
70  *
71  * which applies the matrix to the first argument, and returns the
72  * matrix-vector product in the second. (Note: the const is optional, but
73  * recommended.)
74  */
75 
76 // -----------------------------------------------------------------------------
77 
78 template<typename MODEL, typename OBS> class DRPCGMinimizer : public DRMinimizer<MODEL, OBS> {
79  typedef BMatrix<MODEL, OBS> Bmat_;
80  typedef CostFunction<MODEL, OBS> CostFct_;
81  typedef ControlIncrement<MODEL, OBS> CtrlInc_;
82  typedef HtRinvHMatrix<MODEL, OBS> HtRinvH_;
83 
84  public:
85  const std::string classname() const override {return "DRPCGMinimizer";}
86  DRPCGMinimizer(const eckit::Configuration &, const CostFct_ &);
87  ~DRPCGMinimizer() {}
88 
89  private:
90  double solve(CtrlInc_ &, CtrlInc_ &, CtrlInc_ &, const Bmat_ &, const HtRinvH_ &,
91  const double, const double, const int, const double) override;
92 
93  QNewtonLMP<CtrlInc_, Bmat_> lmp_;
94 };
95 
96 // =============================================================================
97 
98 template<typename MODEL, typename OBS>
99 DRPCGMinimizer<MODEL, OBS>::DRPCGMinimizer(const eckit::Configuration & conf, const CostFct_ & J)
100  : DRMinimizer<MODEL, OBS>(J), lmp_(conf)
101 {}
102 
103 // -----------------------------------------------------------------------------
104 
105 template<typename MODEL, typename OBS>
106 double DRPCGMinimizer<MODEL, OBS>::solve(CtrlInc_ & dx, CtrlInc_ & dxh, CtrlInc_ & rr,
107  const Bmat_ & B, const HtRinvH_ & HtRinvH,
108  const double costJ0Jb, const double costJ0JoJc,
109  const int maxiter, const double tolerance) {
110  // dx increment
111  // dxh B^{-1} dx
112  // rr (sum B^{-1} dx_i^{b} +) G^T H^{-1} d
113 
114  CtrlInc_ qq(dxh);
115  CtrlInc_ pp(dxh);
116  CtrlInc_ hh(dxh);
117  CtrlInc_ zz(dxh);
118  CtrlInc_ pr(dxh);
119  CtrlInc_ r0(dxh);
120  CtrlInc_ ww(dxh);
121 
122  // vectors for re-orthogonalization
123  std::vector<CtrlInc_> vvecs;
124  std::vector<CtrlInc_> zvecs;
125  std::vector<double> scals;
126  // reserve space in vectors to avoid extra copies
127  vvecs.reserve(maxiter+1);
128  zvecs.reserve(maxiter+1);
129  scals.reserve(maxiter+1);
130 
131  // J0
132  const double costJ0 = costJ0Jb + costJ0JoJc;
133 
134  // r_{0}
135  r0 = rr;
136 
137  // r_{0}^T r_{0}
138  Log::info() << "normr0 " << dot_product(rr, rr) << std::endl;
139 
140  // z_{0} = B LMP r_{0}
141  lmp_.multiply(rr, pr);
142  B.multiply(pr, zz);
143  // p_{0} = z_{0}
144  pp = zz;
145  // h_{0} = LMP r_{0}
146  hh = pr;
147 
148  // r_{0}^T z_{0}
149  double dotRr0 = dot_product(rr, zz);
150  double normReduction = 1.0;
151  double rdots = dotRr0;
152  double rdots_old = 0.0;
153 
154  vvecs.push_back(rr);
155  zvecs.push_back(zz);
156  scals.push_back(1.0/dotRr0);
157 
158  Log::info() << std::endl;
159  for (int jiter = 0; jiter < maxiter; ++jiter) {
160  Log::info() << " DRPCG Starting Iteration " << jiter+1 << std::endl;
161 
162  if (jiter > 0) {
163  // beta_{i} = r_{i+1}^T z_{i+1} / r_{i}^T z_{i}
164  double beta = rdots/rdots_old;
165 
166  // p_{i+1} = z_{i+1} + beta*p_{i}
167  pp *= beta;
168  pp += zz;
169 
170  // h_{i+1} = LMP r_{i+1} + beta*h_{i}
171  hh *= beta;
172  hh += pr;
173  }
174 
175  // q_{i} = h_{i} + H^T R^{-1} H p_{i}
176  HtRinvH.multiply(pp, qq);
177  qq += hh;
178 
179  // alpha_{i} = r_{i}^T z_{i} / q_{i}^T p_{i}
180  double rho = dot_product(pp, qq);
181  double alpha = rdots/rho;
182 
183  // dx_{i+1} = dx_{i} + alpha * p_{i}
184  dx.axpy(alpha, pp);
185  // dxh_{i+1} = dxh_{i} + alpha * h_{i} ! for diagnosing Jb
186  dxh.axpy(alpha, hh);
187  // r_{i+1} = r_{i} - alpha * q_{i}
188  rr.axpy(-alpha, qq);
189 
190  // Compute the quadratic cost function
191  // J[dx_{i}] = J[0] - 0.5 dx_{i}^T r_{0}
192  double costJ = costJ0 - 0.5 * dot_product(dx, r0);
193  // Jb[dx_{i}] = 0.5 dx_{i}^T f_{i}
194  double costJb = costJ0Jb + 0.5 * dot_product(dx, dxh);
195  // Jo[dx_{i}] + Jc[dx_{i}] = J[dx_{i}] - Jb[dx_{i}]
196  double costJoJc = costJ - costJb;
197 
198  // Re-orthogonalization
199  for (int jj = 0; jj < jiter; ++jj) {
200  double proj = scals[jj] * dot_product(rr, zvecs[jj]);
201  rr.axpy(-proj, vvecs[jj]);
202  }
203 
204  // z_{i+1} = B LMP r_{i+1}
205  lmp_.multiply(rr, pr);
206  B.multiply(pr, zz);
207 
208  // r_{i}^T z_{i}
209  rdots_old = rdots;
210  // r_{i+1}^T z_{i+1}
211  rdots = dot_product(rr, zz);
212 
213  Log::info() << "rdots " << rdots << " iteration " << jiter << std::endl;
214 
215  // r_{i+1}^T z_{i+1} / r_{0}^T z_{0}
216  normReduction = sqrt(rdots/dotRr0);
217 
218  Log::info() << "DRPCG end of iteration " << jiter+1 << std::endl
219  << " Norm reduction (" << std::setw(2) << jiter+1 << ") = "
220  << util::full_precision(normReduction) << std::endl
221  << " Quadratic cost function: J (" << std::setw(2) << jiter+1 << ") = "
222  << util::full_precision(costJ) << std::endl
223  << " Quadratic cost function: Jb (" << std::setw(2) << jiter+1 << ") = "
224  << util::full_precision(costJb) << std::endl
225  << " Quadratic cost function: JoJc(" << std::setw(2) << jiter+1 << ") = "
226  << util::full_precision(costJoJc) << std::endl << std::endl;
227 
228  // Save the pairs for preconditioning
229  lmp_.push(pp, hh, qq, rho);
230 
231  if (normReduction < tolerance) {
232  Log::info() << "DRPCG: Achieved required reduction in residual norm." << std::endl;
233  break;
234  }
235 
236  vvecs.push_back(rr);
237  zvecs.push_back(zz);
238  scals.push_back(1.0/rdots);
239  }
240 
241 // Generate the (second-level) Limited Memory Preconditioner
242  lmp_.update(B);
243 
244  return normReduction;
245 }
246 
247 // -----------------------------------------------------------------------------
248 
249 } // namespace oops
250 
251 #endif // OOPS_ASSIMILATION_DRPCGMINIMIZER_H_
DRMinimizer.h
oops
The namespace for the main oops code.
Definition: ErrorCovarianceL95.cc:22
HtRinvHMatrix.h
oops::BMatrix::multiply
void multiply(const CtrlInc_ &x, CtrlInc_ &bx) const
Definition: BMatrix.h:33
oops::DRPCGMinimizer::~DRPCGMinimizer
~DRPCGMinimizer()
Definition: DRPCGMinimizer.h:87
oops::DRPCGMinimizer::CtrlInc_
ControlIncrement< MODEL, OBS > CtrlInc_
Definition: DRPCGMinimizer.h:81
CostFunction.h
oops::DRPCGMinimizer::classname
const std::string classname() const override
Definition: DRPCGMinimizer.h:85
oops::QNewtonLMP
Definition: QNewtonLMP.h:40
oops::ControlIncrement
Definition: ControlIncrement.h:46
oops::DRPCGMinimizer::DRPCGMinimizer
DRPCGMinimizer(const eckit::Configuration &, const CostFct_ &)
Definition: DRPCGMinimizer.h:99
oops::DRPCGMinimizer::Bmat_
BMatrix< MODEL, OBS > Bmat_
Definition: DRPCGMinimizer.h:79
QNewtonLMP.h
oops::DRPCGMinimizer::solve
double solve(CtrlInc_ &, CtrlInc_ &, CtrlInc_ &, const Bmat_ &, const HtRinvH_ &, const double, const double, const int, const double) override
Definition: DRPCGMinimizer.h:106
oops::HtRinvHMatrix
The matrix.
Definition: HtRinvHMatrix.h:35
BMatrix.h
oops::DRMinimizer
DR (Derber and Rosati) Minimizers.
Definition: DRMinimizer.h:44
oops::DRPCGMinimizer::lmp_
QNewtonLMP< CtrlInc_, Bmat_ > lmp_
Definition: DRPCGMinimizer.h:93
oops::DRPCGMinimizer::HtRinvH_
HtRinvHMatrix< MODEL, OBS > HtRinvH_
Definition: DRPCGMinimizer.h:82
ControlIncrement.h
oops::DRPCGMinimizer
DRPCG Minimizer.
Definition: DRPCGMinimizer.h:78
oops::BMatrix
The matrix.
Definition: BMatrix.h:27
oops::CostFunction
Cost Function.
Definition: CostFunction.h:53
oops::HtRinvHMatrix::multiply
void multiply(const CtrlInc_ &dx, CtrlInc_ &dz) const
Definition: HtRinvHMatrix.h:60
oops::ControlIncrement::axpy
void axpy(const double, const ControlIncrement &)
Definition: ControlIncrement.h:196
oops::DRPCGMinimizer::CostFct_
CostFunction< MODEL, OBS > CostFct_
Definition: DRPCGMinimizer.h:80