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