OOPS
RPCGMinimizer.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_RPCGMINIMIZER_H_
12 #define OOPS_ASSIMILATION_RPCGMINIMIZER_H_
13 
14 #include <string>
15 #include <vector>
16 
17 #include "oops/assimilation/CostFunction.h"
18 #include "oops/assimilation/DualMinimizer.h"
19 #include "oops/assimilation/DualVector.h"
20 #include "oops/assimilation/HBHtMatrix.h"
21 #include "oops/assimilation/RinvMatrix.h"
22 #include "oops/base/IdentityMatrix.h"
23 #include "oops/util/dot_product.h"
24 #include "oops/util/formats.h"
25 #include "oops/util/Logger.h"
26 
27 namespace oops {
28 
29 /// RPCG Minimizer
30 /*!
31  * \brief Augmented Restricted Preconditioned Conjugate Gradients.
32  *
33  * This solver is based on the algorithm proposed in Gratton and Tshimanga,
34  * QJRMS, 135: 1573-1585 (2009). It performs minimization in observation space.
35  *
36  * It solves the linear system \f$ (I + R^{-1}HBH^T) lambda = H^T R^{-1}d \f$
37  * with \f$ HBH^T \f$ inner-product in the augmented observation space.
38  *
39  * Note that the traditional \f$ B\f$-preconditioning in model space
40  * corresponds to \f$I\f$ for this algorithm.
41  *
42  * A second-level preconditioner, \f$ G \f$, must be symmetric and
43  * positive definite with respect to \f$ HBH^T \f$ inner-product.
44  * Possible preconditioning is detailed in S. Gurol, PhD Manuscript, 2013.
45  *
46  * On entry:
47  * - vv = starting point, \f$ v_0 \f$.
48  * - vvp = starting point augmented part, \f$ vp_0 \f$
49  * - rr = right hand side.
50  * - dy = \f$ H (xb - xk) \f$
51  * - sigma = \f$ (xb - xk)^T B^{-1} (xb - xk)\f$.
52  * - HBHt = \f$ HBH^T \f$.
53  * - Rinv = \f$ R^{-1} \f$.
54  * - G = preconditioner \f$ G \f$.
55  * - Gt = preconditioner transpose \f$ G^T \f$.
56  *
57  * On exit, vv and vvp will contain the solution \f$ lambda = [vv; vp] \f$
58  *
59  * The return value is the achieved reduction in preconditioned residual norm.
60  *
61  * Iteration will stop if the maximum iteration limit "maxiter" is reached
62  * or if the preconditioned residual norm reduces by a factor of "tolerance".
63  */
64 
65 // -----------------------------------------------------------------------------
66 
67 template<typename MODEL, typename OBS> class RPCGMinimizer : public DualMinimizer<MODEL, OBS> {
68  typedef CostFunction<MODEL, OBS> CostFct_;
69  typedef DualVector<MODEL, OBS> Dual_;
70  typedef HBHtMatrix<MODEL, OBS> HBHt_;
71  typedef RinvMatrix<MODEL, OBS> Rinv_;
72 
73  public:
74  const std::string classname() const override {return "RPCGMinimizer";}
75  RPCGMinimizer(const eckit::Configuration &, const CostFct_ & J): DualMinimizer<MODEL, OBS>(J) {}
76  ~RPCGMinimizer() {}
77 
78  private:
79  double solve(Dual_ &, double &, Dual_ &, const HBHt_ &, const Rinv_ &,
80  const int &, const double &, Dual_ &, const double &) override;
81 };
82 
83 // =============================================================================
84 
85 template<typename MODEL, typename OBS>
86 double RPCGMinimizer<MODEL, OBS>::solve(Dual_ & vv, double & vvp, Dual_ & rr,
87  const HBHt_ & HBHt, const Rinv_ & Rinv,
88  const int & maxiter, const double & tolerance,
89  Dual_ & dy, const double & sigma) {
90  IdentityMatrix<Dual_> precond;
91  IdentityMatrix<Dual_> precondt;
92 
93  Dual_ zz(vv);
94  Dual_ ww(vv);
95  Dual_ pp(vv);
96  Dual_ tt(vv);
97  Dual_ ll(vv);
98  Dual_ qq(vv);
99  Dual_ w(vv);
100  Dual_ v(vv);
101 
102  // augmented part of the vectors
103  double rrp = 1.0;
104  double llp;
105  double wwp;
106  double zzp;
107  double ppp = 0.0;
108  double ttp = 0.0;
109  double qqp;
110  double vp = 1.0;
111  double wp;
112 
113  std::vector<Dual_> vVEC; // required for re-orthogonalization
114  std::vector<Dual_> wVEC; // required for re-orthogonalization
115  // reserve space to avoid extra copies
116  vVEC.reserve(maxiter+1);
117  wVEC.reserve(maxiter+1);
118 
119  std::vector<double> vpVEC;
120  std::vector<double> wpVEC;
121 
122  // llaug = HBHtaug rraug
123  HBHt.multiply(rr, ll);
124  ll.axpy(rrp, dy);
125  llp = dot_product(dy, rr) + sigma*rrp;
126 
127  // wwaug = Gtaug llaug
128  precondt.multiply(ll, ww); // UPDATE WHEN G IS NOT IDENTITY
129  wwp = llp;
130 
131  // zzaug = Gaug rraug
132  precond.multiply(rr, zz); // UPDATE WHEN G IS NOT IDENTITY
133  zzp = rrp;
134 
135  double dotwr = dot_product(rr, ww) + rrp*wwp;
136  double normReduction = 1.0;
137  double dotw0r0 = dotwr;
138  double dotwr_old = 0.0;
139 
140  v = rr;
141  v *= 1/sqrt(dotwr);
142  vp = rrp;
143  vp *= 1/sqrt(dotwr);
144  w = ww;
145  w *= 1/sqrt(dotwr);
146  wp = wwp;
147  wp *= 1/sqrt(dotwr);
148 
149  vVEC.clear();
150  vpVEC.clear();
151  wVEC.clear();
152  wpVEC.clear();
153 
154  vVEC.push_back(v);
155  vpVEC.push_back(vp);
156  wVEC.push_back(w);
157  wpVEC.push_back(wp);
158 
159  Log::info() << std::endl;
160  for (int jiter = 0; jiter < maxiter; ++jiter) {
161  Log::info() << "RPCG Starting Iteration " << jiter+1 << std::endl;
162 
163  if (jiter > 0) {
164  double beta = dotwr/dotwr_old;
165  Log::info() << "RPCG beta = " << beta << std::endl;
166 
167  pp *= beta;
168  ppp *= beta;
169 
170  tt *= beta;
171  ttp *= beta;
172  }
173 
174  pp += zz;
175  ppp += zzp; // ppaug = zzaug + beta*ppaug
176 
177  tt += ww;
178  ttp += wwp; // ttaug = wwaug + beta*ttaug
179 
180  // (RinvHBHt + I) pp
181  Rinv.multiply(tt, qq);
182  qq += pp;
183  qqp = ppp; // qqaug = ppaug + Rinv_aug ttaug
184 
185  double alpha = dotwr/(dot_product(qq, tt) + qqp * ttp);
186  Log::info() << "RPCG alpha = " << alpha << std::endl;
187 
188  vv.axpy(alpha, pp);
189  vvp = vvp + alpha*ppp; // vvaug = vvaug + alpha*ppaug
190 
191  rr.axpy(-alpha, qq);
192  rrp = rrp - alpha*qqp; // rraug = rraug - alpha*qqaug
193 
194  // Re-orthogonalization
195  for (int iiter = 0; iiter < jiter; ++iiter) {
196  double proj = dot_product(rr, wVEC[iiter]) + rrp * wpVEC[iiter];
197  rr.axpy(-proj, vVEC[iiter]);
198  rrp -= proj*vpVEC[iiter];
199  }
200 
201  // llaug = HBHt_aug rraug
202  HBHt.multiply(rr, ll);
203  ll.axpy(rrp, dy);
204  llp = dot_product(dy, rr) + sigma * rrp;
205 
206  // wwaug = Gt_aug llaug
207  precondt.multiply(ll, ww); // UPDATE WHEN G IS NOT IDENTITY
208  wwp = llp;
209 
210  // zzaug = Gaug rraug
211  precond.multiply(rr, zz); // UPDATE WHEN G IS NOT IDENTITY
212  zzp = rrp;
213 
214  dotwr_old = dotwr;
215  dotwr = dot_product(ww, rr) + rrp * wwp;
216 
217  v = rr;
218  v *= 1/sqrt(dotwr);
219  vp = rrp;
220  vp *= 1/sqrt(dotwr);
221  w = ww;
222  w *= 1/sqrt(dotwr);
223  wp = wwp;
224  wp *= 1/sqrt(dotwr);
225 
226  vVEC.push_back(v);
227  vpVEC.push_back(vp);
228  wVEC.push_back(w);
229  wpVEC.push_back(wp);
230 
231  normReduction = sqrt(dotwr/dotw0r0);
232 
233  Log::info() << "RPCG end of iteration " << jiter+1 << ". Norm reduction= "
234  << util::full_precision(normReduction) << std::endl << std::endl;
235 
236  if (normReduction < tolerance) {
237  Log::info() << "RPCG: Achieved required reduction in residual norm." << std::endl;
238  break;
239  }
240  }
241  return normReduction;
242 }
243 
244 // -----------------------------------------------------------------------------
245 
246 } // namespace oops
247 
248 #endif // OOPS_ASSIMILATION_RPCGMINIMIZER_H_
oops
The namespace for the main oops code.
Definition: ErrorCovarianceL95.cc:22
oops::RinvMatrix
The matrix.
Definition: RinvMatrix.h:29
IdentityMatrix.h
oops::RPCGMinimizer::Rinv_
RinvMatrix< MODEL, OBS > Rinv_
Definition: RPCGMinimizer.h:71
oops::DualVector::axpy
void axpy(const double, const DualVector &)
Definition: DualVector.h:212
oops::RinvMatrix::multiply
void multiply(const Dual_ &dx, Dual_ &dy) const
Definition: RinvMatrix.h:36
DualMinimizer.h
CostFunction.h
oops::DualVector
Container of dual space vectors for all terms of the cost function.
Definition: DualVector.h:34
oops::RPCGMinimizer
RPCG Minimizer.
Definition: RPCGMinimizer.h:67
RinvMatrix.h
oops::DualMinimizer
Dual Minimizer.
Definition: DualMinimizer.h:40
oops::RPCGMinimizer::classname
const std::string classname() const override
Definition: RPCGMinimizer.h:74
oops::RPCGMinimizer::solve
double solve(Dual_ &, double &, Dual_ &, const HBHt_ &, const Rinv_ &, const int &, const double &, Dual_ &, const double &) override
Definition: RPCGMinimizer.h:86
oops::HBHtMatrix
The matrix.
Definition: HBHtMatrix.h:32
oops::IdentityMatrix::multiply
void multiply(const VECTOR &a, VECTOR &b) const
Definition: IdentityMatrix.h:20
oops::RPCGMinimizer::RPCGMinimizer
RPCGMinimizer(const eckit::Configuration &, const CostFct_ &J)
Definition: RPCGMinimizer.h:75
HBHtMatrix.h
oops::RPCGMinimizer::CostFct_
CostFunction< MODEL, OBS > CostFct_
Definition: RPCGMinimizer.h:68
oops::IdentityMatrix
Identity matrix.
Definition: IdentityMatrix.h:18
oops::RPCGMinimizer::HBHt_
HBHtMatrix< MODEL, OBS > HBHt_
Definition: RPCGMinimizer.h:70
oops::RPCGMinimizer::Dual_
DualVector< MODEL, OBS > Dual_
Definition: RPCGMinimizer.h:69
oops::RPCGMinimizer::~RPCGMinimizer
~RPCGMinimizer()
Definition: RPCGMinimizer.h:76
oops::CostFunction
Cost Function.
Definition: CostFunction.h:53
oops::HBHtMatrix::multiply
void multiply(const Dual_ &dy, Dual_ &dz) const
Definition: HBHtMatrix.h:60
DualVector.h