MADNESS  0.10.1
gmres.h
Go to the documentation of this file.
1 /*
2  This file is part of MADNESS.
3 
4  Copyright (C) 2007,2010 Oak Ridge National Laboratory
5 
6  This program is free software; you can redistribute it and/or modify
7  it under the terms of the GNU General Public License as published by
8  the Free Software Foundation; either version 2 of the License, or
9  (at your option) any later version.
10 
11  This program is distributed in the hope that it will be useful,
12  but WITHOUT ANY WARRANTY; without even the implied warranty of
13  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14  GNU General Public License for more details.
15 
16  You should have received a copy of the GNU General Public License
17  along with this program; if not, write to the Free Software
18  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
19 
20  For more information please contact:
21 
22  Robert J. Harrison
23  Oak Ridge National Laboratory
24  One Bethel Valley Road
25  P.O. Box 2008, MS-6367
26 
27  email: harrisonrj@ornl.gov
28  tel: 865-241-3937
29  fax: 865-572-0680
30 
31 
32  $Id$
33 */
34 
35 #ifndef MADNESS_LINALG_GMRES_H__INCLUDED
36 #define MADNESS_LINALG_GMRES_H__INCLUDED
37 
38 #include <madness/tensor/tensor.h>
39 #include <madness/world/print.h>
40 #include <iostream>
42 
43 /** \file gmres.h
44  \brief Defines a general operator interface and a templated GMRES solver
45  for solving linear equations.
46 
47  \ingroup solvers
48 
49  An AbstractVectorSpace class is also defined to guarantee a uniform way
50  for linear algebra routines to access common operations like norm, inner
51  product etc. Implementations for MADNESS Vectors and Functions are
52  provided in the VectorSpace<floating_point, NDIM> and
53  FunctionSpace<floating_point, NDIM> classes. */
54 
55 namespace madness {
56 
57  /** \brief A generic operator: takes in one \c T and produces another \c T.
58 
59  Override the protected action() function to implement an Operator.
60  */
61  template <typename T>
62  class Operator {
63  protected:
64  /** \brief The action of the operator
65 
66  \param[in] in The input vector
67  \param[out] out The action of the operator on the input vector
68  */
69  virtual void action(const T &in, T &out) const = 0;
70 
71  public:
72  /** \brief Public access to the operator's action, returns out for
73  convenience.
74 
75  \param[in] in The input vector
76  \param[out] out The action of the operator on the input vector
77  */
78  T & applyOp(const T &in, T &out) const {
79  action(in, out);
80  return out;
81  }
82 
83  virtual ~Operator() {}
84  };
85 
86  /** \brief A generic vector space which provides common operations needed
87  by linear algebra routines (norm, inner product, etc.)
88 
89  Most of these routines will presumedly be defined, but the names may
90  not be the same (i.e. madness:Function uses norm2, madness::Tensor
91  uses normf)
92 
93  - \c T is the vector type
94  - \c real_type the real type used for norms, etc.
95  - \c scalar_type is the type used for inner products and gaxpy
96  coefficients, etc.
97 
98  When implementing a child class, real_type and scalar_type will
99  probably be obtained from TensorTypeData<T> (see the VectorSpace
100  and FunctionSpace classes below). */
101  template <typename T, typename real_type, typename scalar_type>
103  private:
104  /** \brief Bury the default constructor to prevent its use. */
106 
107  public:
108  /// \brief The world
110 
111  /** \brief Make a vector space.
112 
113  The World is needed to limit output, and may be needed
114  for spaces working with MADNESS functions.
115 
116  \param[in] world The world.
117  */
119 
120  virtual ~AbstractVectorSpace() {}
121 
122  /// \brief The norm of a vector
123  virtual real_type norm(const T &) const = 0;
124 
125  /** \brief Scales the vector (in-place) by a constant
126  \f[ \vec{x} \leftarrow c \vec{x} \f]
127 
128  \return The scaled vector */
129  virtual T & scale(T &, const scalar_type &) const = 0;
130 
131  /** \brief Standard bilinear gaxpy
132  \f[ \vec{x} \leftarrow a \vec{x} + b \vec{y} \f]
133 
134  \return The new vector, \f$\vec{x}\f$ */
135  virtual T & gaxpy(T &x, const scalar_type &a, const T &y,
136  const scalar_type &b) const = 0;
137 
138  /// \brief The inner product between two vectors
139  virtual scalar_type inner(const T &, const T &) const = 0;
140 
141  /** \brief Any special instructions to be executed when a vector
142  is no longer needed.
143 
144  Unless otherwise specified, do nothing. */
145  virtual void destroy(T &) const {}
146  };
147 
148  // Real part of a real number
149  // This is necessary to make templating nightmares disappear in
150  // VectorSpace::norm
151  static inline double real(double x) { return x; }
152 
153  // Imaginary part of a real number
154  // This is necessary to make templating nightmares disappear in
155  // VectorSpace::norm
156  static inline double imag(double x) { return 0.0; }
157 
158  /// \brief A vector space using MADNESS Vectors
159  template <typename T, int NDIM>
160  class VectorSpace : public AbstractVectorSpace<Vector<T, NDIM>,
161  typename TensorTypeData<T>::float_scalar_type, T> {
162 
163  protected:
165 
166  public:
168  typedef T scalar_type;
169 
172 
173  virtual ~VectorSpace() {}
174 
176  const {
177 
178  real_type ret, mag;
179  int i;
180 
181  mag = real(vec[0]);
182  ret = mag*mag;
183  if(iscplx) {
184  mag = imag(vec[0]);
185  ret += mag*mag;
186  }
187 
188  for(i = 1; i < NDIM; ++i) {
189  mag = real(vec[i]);
190  ret += mag*mag;
191  if(iscplx) {
192  mag = imag(vec[i]);
193  ret += mag*mag;
194  }
195  }
196 
197  return sqrt(ret);
198  }
199 
201  Vector<scalar_type, NDIM> &vec, const scalar_type &c) const {
202 
203  vec *= c;
204  return vec;
205  }
206 
209  const Vector<scalar_type, NDIM> &y, const scalar_type &b)
210  const {
211 
212  for(int i = 0; i < NDIM; ++i)
213  x[i] = a * x[i] + b * y[i];
214  return x;
215  }
216 
218  const Vector<scalar_type, NDIM> &r) const {
219 
220  scalar_type ret;
221 
222  if(iscplx) {
223  ret = conj(l[0]) * r[0];
224  for(int i = 1; i < NDIM; ++i)
225  ret += conj(l[i]) * r[i];
226  }
227  else {
228  ret = l[0] * r[0];
229  for(int i = 1; i < NDIM; ++i)
230  ret += l[i] * r[i];
231  }
232 
233  return ret;
234  }
235  };
236 
237  /// \brief A vector space using MADNESS Functions
238  template <typename T, int NDIM>
239  class FunctionSpace : public AbstractVectorSpace<Function<T, NDIM>,
240  typename TensorTypeData<T>::float_scalar_type, T> {
241 
242  public:
244  typedef T scalar_type;
245 
248 
249  virtual ~FunctionSpace() {}
250 
252  const {
253 
254  return vec.norm2();
255  }
256 
258  Function<scalar_type, NDIM> &vec, const scalar_type &c) const {
259 
260  vec.scale(c);
261  return vec;
262  }
263 
266  const Function<scalar_type, NDIM> &y, const scalar_type &b)
267  const {
268 
269  x.gaxpy(a, y, b);
270  return x;
271  }
272 
274  const Function<scalar_type, NDIM> &r) const {
275 
276  return l.inner(r);
277  }
278 
279  virtual void destroy(Function<scalar_type, NDIM> &f) const {
280  f.clear();
281  }
282  };
283 
284  /// \brief A vector space using MADNESS Vectors of MADNESS Functions
285  template <typename T, int VDIM, int FDIM>
287  Vector<Function<T, FDIM>, VDIM>,
288  typename TensorTypeData<T>::float_scalar_type, T> {
289 
290  public:
292  typedef T scalar_type;
293 
295  <Vector<Function<T, FDIM>, VDIM>, real_type, scalar_type>
296  (world) {}
297 
299 
300  virtual real_type norm(
301  const Vector<Function<scalar_type, FDIM>, VDIM> &vec)
302  const {
303 
304  real_type temp, ret = vec[0].norm2();
305  ret *= ret;
306  for(int i = 1; i < VDIM; ++i) {
307  temp = vec[i].norm2();
308  ret += temp*temp;
309  }
310  return sqrt(ret);
311  }
312 
315  const scalar_type &c) const {
316 
317  for(int i = 0; i < VDIM; ++i)
318  vec[i].scale(c);
319  return vec;
320  }
321 
324  const scalar_type &a,
325  const Vector<Function<scalar_type, FDIM>, VDIM> &y,
326  const scalar_type &b) const {
327 
328  for(int i = 0; i < VDIM; ++i)
329  x[i].gaxpy(a, y[i], b);
330  return x;
331  }
332 
334  const Vector<Function<scalar_type, FDIM>, VDIM> &l,
335  const Vector<Function<scalar_type, FDIM>, VDIM> &r) const {
336 
337  scalar_type ret = l[0].inner(r[0]);
338  for(int i = 0; i < VDIM; ++i)
339  ret += l[i].inner(r[i]);
340  return ret;
341  }
342 
344  const {
345  for(int i = 0; i < VDIM; ++i)
346  f[i].clear();
347  }
348  };
349 
350  /** \brief A GMRES solver routine for linear systems,
351  \f$ \mathbf{A} \vec{x} = \vec{b} \f$.
352 
353  Requires the vector space, the operator, A, the inhomogeneity b, an
354  initial guess x, the maximum number of iterations, the convergence
355  threshold for the solution residual, the convergence threshold for
356  the norm of the update vector, and a flag for producing output.
357 
358  The vector space object provides a standard way to compute the needed
359  linear algebra operations (norm, inner product, etc.).
360 
361  The printed output, if desired, shows iteration number, the residual,
362  the norm of the change in solution vectors from one iteration to the
363  next, and the effective rank of the GMRES matrix at that iteration.
364 
365  On output, x has the computed solution, resid_thresh has its residual,
366  update_thresh has the latest updatenorm, and maxiters has the number
367  of iterations needed.
368 
369  \param[in] space The AbstractVectorSpace interface for the type
370  \param[in] op The Operator \f$\mathbf{A}\f$
371  \param[in] b The right-hand-side vector \f$\vec{b}\f$
372  \param[in,out] x Input: The initial guess for \f$\vec{x}\f$. Output:
373  The computed solution for \f$\vec{x}\f$.
374  \param[in,out] maxiters Input: Maximum number of iterations to perform.
375  Output: Actual iterations performed.
376  \param[in,out] resid_thresh Input: Convergence threshold for the
377  residual \f$||\mathbf{A} \vec{x} - \vec{b}||\f$.
378  Output: The residual after the final iteration.
379  \param[in,out] update_thresh Input: Convergence threshold for the
380  update vector \f$|| \vec{x}_{iter} -
381  \vec{x}_{iter-1}||\f$. Output: The value after the
382  final iteration.
383  \param[in] outp True if output to stdout is desired, false otherwise,
384  defaults to false if unspecified.
385  */
386  template <typename T, typename real_type, typename scalar_type>
388  const Operator<T> &op, const T &b, T &x, int &maxiters,
389  real_type &resid_thresh, real_type &update_thresh,
390  const bool outp = false) {
391 
392  int iter, i;
393  long rank;
394  std::vector<T> V;
395  T r;
396  Tensor<scalar_type> H(maxiters+1, maxiters);
397  Tensor<scalar_type> betae(maxiters+1);
398  Tensor<scalar_type> y, yold;
400  Tensor<real_type> sumsq;
401  real_type resid, norm, updatenorm;
402  World &world = space.world;
403 
404  // initialize
405  H = 0.0;
406  betae = 0.0;
407 
408  // construct the first subspace basis vector
409  iter = 0;
410  op.applyOp(x, r);
411  space.gaxpy(r, -1.0, b, 1.0);
412  betae[0] = resid = space.norm(r);
413  if(outp && world.rank() == 0)
414  printf("itr rnk update_norm resid\n%.3d N/A N/A %.6e\n",
415  iter, resid);
416  if(resid < resid_thresh) {
417  maxiters = 1;
418  resid_thresh = resid;
419  update_thresh = 0.0;
420  return;
421  }
422  space.scale(r, 1.0 / resid);
423  V.push_back(r);
424  ++iter;
425 
426  do {
427  // compute the new vector
428  op.applyOp(V[iter - 1], r);
429 
430  // orthogonalize the new vector
431  for(i = 0; i < iter; ++i) {
432  H(i, iter-1) = space.inner(V[i], r);
433  space.gaxpy(r, 1.0, V[i], -H(i, iter-1));
434  }
435  H(iter, iter-1) = norm = space.norm(r);
436 
437  // normalize the vector and save it
438  space.scale(r, 1.0 / norm);
439  V.push_back(r);
440 
441  // solve Hy == betae for y
442  gelss(H(Slice(0, iter), Slice(0, iter-1)), betae(Slice(0, iter)),
443  1.0e-12, y, s, rank, sumsq);
444 
445  // residual from the least-squares fit
446  resid = sumsq[0];
447 
448  // compute the update norm,
449  // || x_n - x_{n-1} ||
450  // = || (x_0 + V_n y_n) - (x_0 + V_{n-1} y_{n-1}) ||
451  // = || V_n (y_n - y_{n-1}) ||, assuming y_{n-1}[n] = 0
452  // = || y_n - y_{n-1} ||
453  if(iter == 1)
454  updatenorm = y.normf();
455  else {
456  scalar_type temp = y[0] - yold[0];
457  updatenorm = real(temp)*real(temp) + imag(temp)*imag(temp);
458  for(i = 1; i < iter-1; ++i) {
459  temp = y[i] - yold[i];
460  updatenorm += real(temp)*real(temp) +
461  imag(temp)*imag(temp);
462  }
463  updatenorm += real(y[iter-1]) * real(y[iter-1]) +
464  imag(y[iter-1]) * imag(y[iter-1]);
465  updatenorm = sqrt(updatenorm);
466  }
467  yold = copy(y);
468 
469  if(norm < 1.0e-10) {
470  // we just got the zero vector -> no more progress
471 
472  // if this is the first iteration, compute the residual
473  if(iter == 1) {
474  space.destroy(r);
475  space.gaxpy(x, 1.0, V[0], betae[0]);
476  op.applyOp(x, r);
477  space.gaxpy(r, -1.0, b, 1.0);
478  resid_thresh = space.norm(r);
479  update_thresh = updatenorm;
480  space.destroy(r);
481  space.destroy(V[0]);
482 
483  maxiters = 1;
484  if(outp && world.rank() == 0)
485  printf("%.3d N/A %.6e %.6e ** Zero Vector Encount" \
486  "ered **\n", iter, updatenorm, resid_thresh);
487  return;
488  }
489 
490  if(outp && world.rank() == 0)
491  printf("%.3d %.3ld %.6e %.6e ** Zero Vector Encount" \
492  "ered **\n", iter, rank, updatenorm, resid);
493  break;
494  }
495 
496  if(outp && world.rank() == 0) {
497  printf("%.3d %.3ld %.6e %.6e", iter, rank, updatenorm, resid);
498  if(iter != rank)
499  printf(" ** Questionable Progress **");
500  printf("\n");
501  }
502 
503  ++iter;
504  } while(iter <= maxiters && resid > resid_thresh &&
505  updatenorm > update_thresh);
506 
507  // build the solution vector and destroy the basis vectors
508  for(i = 0; i < iter; ++i) {
509  space.gaxpy(x, 1.0, V[i], y[i]);
510  space.destroy(V[i]);
511  }
512 
513  resid_thresh = resid;
514  update_thresh = updatenorm;
515  maxiters = iter;
516  }
517 
518 } // end of madness namespace
519 
520 #endif // MADNESS_LINALG_GMRES_H__INCLUDED
A generic vector space which provides common operations needed by linear algebra routines (norm,...
Definition: gmres.h:102
AbstractVectorSpace()
Bury the default constructor to prevent its use.
AbstractVectorSpace(World &world)
Make a vector space.
Definition: gmres.h:118
virtual T & gaxpy(T &x, const scalar_type &a, const T &y, const scalar_type &b) const =0
Standard bilinear gaxpy.
virtual T & scale(T &, const scalar_type &) const =0
Scales the vector (in-place) by a constant.
virtual ~AbstractVectorSpace()
Definition: gmres.h:120
virtual scalar_type inner(const T &, const T &) const =0
The inner product between two vectors.
World & world
The world.
Definition: gmres.h:109
virtual void destroy(T &) const
Any special instructions to be executed when a vector is no longer needed.
Definition: gmres.h:145
virtual real_type norm(const T &) const =0
The norm of a vector.
A vector space using MADNESS Functions.
Definition: gmres.h:240
virtual scalar_type inner(const Function< scalar_type, NDIM > &l, const Function< scalar_type, NDIM > &r) const
The inner product between two vectors.
Definition: gmres.h:273
virtual real_type norm(const Function< scalar_type, NDIM > &vec) const
The norm of a vector.
Definition: gmres.h:251
virtual Function< scalar_type, NDIM > & gaxpy(Function< scalar_type, NDIM > &x, const scalar_type &a, const Function< scalar_type, NDIM > &y, const scalar_type &b) const
Standard bilinear gaxpy.
Definition: gmres.h:264
virtual Function< scalar_type, NDIM > & scale(Function< scalar_type, NDIM > &vec, const scalar_type &c) const
Scales the vector (in-place) by a constant.
Definition: gmres.h:257
virtual ~FunctionSpace()
Definition: gmres.h:249
FunctionSpace(World &world)
Definition: gmres.h:246
virtual void destroy(Function< scalar_type, NDIM > &f) const
Any special instructions to be executed when a vector is no longer needed.
Definition: gmres.h:279
T scalar_type
Definition: gmres.h:244
TensorTypeData< T >::float_scalar_type real_type
Definition: gmres.h:243
A multiresolution adaptive numerical function.
Definition: mra.h:122
Function< T, NDIM > & gaxpy(const T &alpha, const Function< Q, NDIM > &other, const R &beta, bool fence=true)
Inplace, general bi-linear operation in wavelet basis. No communication except for optional fence.
Definition: mra.h:981
A generic operator: takes in one T and produces another T.
Definition: gmres.h:62
virtual ~Operator()
Definition: gmres.h:83
virtual void action(const T &in, T &out) const =0
The action of the operator.
T & applyOp(const T &in, T &out) const
Public access to the operator's action, returns out for convenience.
Definition: gmres.h:78
A slice defines a sub-range or patch of a dimension.
Definition: slice.h:103
Traits class to specify support of numeric types.
Definition: type_data.h:56
A tensor is a multidimension array.
Definition: tensor.h:317
float_scalar_type normf() const
Returns the Frobenius norm of the tensor.
Definition: tensor.h:1726
A vector space using MADNESS Vectors of MADNESS Functions.
Definition: gmres.h:288
virtual Vector< Function< scalar_type, FDIM >, VDIM > & scale(Vector< Function< scalar_type, FDIM >, VDIM > &vec, const scalar_type &c) const
Scales the vector (in-place) by a constant.
Definition: gmres.h:313
virtual scalar_type inner(const Vector< Function< scalar_type, FDIM >, VDIM > &l, const Vector< Function< scalar_type, FDIM >, VDIM > &r) const
The inner product between two vectors.
Definition: gmres.h:333
TensorTypeData< T >::float_scalar_type real_type
Definition: gmres.h:291
VectorOfFunctionsSpace(World &world)
Definition: gmres.h:294
virtual ~VectorOfFunctionsSpace()
Definition: gmres.h:298
virtual Vector< Function< scalar_type, FDIM >, VDIM > & gaxpy(Vector< Function< scalar_type, FDIM >, VDIM > &x, const scalar_type &a, const Vector< Function< scalar_type, FDIM >, VDIM > &y, const scalar_type &b) const
Standard bilinear gaxpy.
Definition: gmres.h:322
virtual void destroy(Vector< Function< scalar_type, FDIM >, VDIM > &f) const
Any special instructions to be executed when a vector is no longer needed.
Definition: gmres.h:343
T scalar_type
Definition: gmres.h:292
virtual real_type norm(const Vector< Function< scalar_type, FDIM >, VDIM > &vec) const
The norm of a vector.
Definition: gmres.h:300
A vector space using MADNESS Vectors.
Definition: gmres.h:161
virtual Vector< scalar_type, NDIM > & scale(Vector< scalar_type, NDIM > &vec, const scalar_type &c) const
Scales the vector (in-place) by a constant.
Definition: gmres.h:200
TensorTypeData< T >::float_scalar_type real_type
Definition: gmres.h:167
virtual Vector< scalar_type, NDIM > & gaxpy(Vector< scalar_type, NDIM > &x, const scalar_type &a, const Vector< scalar_type, NDIM > &y, const scalar_type &b) const
Standard bilinear gaxpy.
Definition: gmres.h:207
virtual real_type norm(const Vector< scalar_type, NDIM > &vec) const
The norm of a vector.
Definition: gmres.h:175
VectorSpace(World &world)
Definition: gmres.h:170
virtual scalar_type inner(const Vector< scalar_type, NDIM > &l, const Vector< scalar_type, NDIM > &r) const
The inner product between two vectors.
Definition: gmres.h:217
static const bool iscplx
Definition: gmres.h:164
T scalar_type
Definition: gmres.h:168
virtual ~VectorSpace()
Definition: gmres.h:173
A simple, fixed dimension vector.
Definition: vector.h:64
A parallel world class.
Definition: world.h:132
ProcessID rank() const
Returns the process rank in this World (same as MPI_Comm_rank()).
Definition: world.h:318
auto T(World &world, response_space &f) -> response_space
Definition: global_functions.cc:34
Tensor< double > op(const Tensor< double > &x)
Definition: kain.cc:508
double norm(const T &t)
Definition: adquad.h:42
File holds all helper structures necessary for the CC_Operator and CC2 class.
Definition: DFParameters.h:10
Function< T, NDIM > conj(const Function< T, NDIM > &f, bool fence=true)
Return the complex conjugate of the input function with the same distribution and optional fence.
Definition: mra.h:2046
Function< T, NDIM > copy(const Function< T, NDIM > &f, const std::shared_ptr< WorldDCPmapInterface< Key< NDIM > > > &pmap, bool fence=true)
Create a new copy of the function with different distribution and optional fence.
Definition: mra.h:2002
void gelss(const Tensor< T > &a, const Tensor< T > &b, double rcond, Tensor< T > &x, Tensor< typename Tensor< T >::scalar_type > &s, long &rank, Tensor< typename Tensor< T >::scalar_type > &sumsq)
Solve Ax = b for general A using the LAPACK *gelss routines.
Definition: lapack.cc:884
NDIM & f
Definition: mra.h:2416
double imag(double x)
Definition: complexfun.h:56
double real(double x)
Definition: complexfun.h:52
void GMRES(const AbstractVectorSpace< T, real_type, scalar_type > &space, const Operator< T > &op, const T &b, T &x, int &maxiters, real_type &resid_thresh, real_type &update_thresh, const bool outp=false)
A GMRES solver routine for linear systems, .
Definition: gmres.h:387
Vector< T, sizeof...(Ts)+1 > vec(T t, Ts... ts)
Factory function for creating a madness::Vector.
Definition: vector.h:711
static const double b
Definition: nonlinschro.cc:119
static const double a
Definition: nonlinschro.cc:118
Defines simple templates for printing to std::cout "a la Python".
static const double c
Definition: relops.cc:10
static double V(const coordT &r)
Definition: tdse.cc:288
Defines and implements most of Tensor.
Prototypes for a partial interface from Tensor to LAPACK.
static const std::size_t NDIM
Definition: testpdiff.cc:42