MPQC  3.0.0-alpha
conjgrad.h
1 //
2 // conjgrad.h
3 //
4 // Copyright (C) 2012 Edward Valeev
5 //
6 // Author: Edward Valeev <evaleev@vt.edu>
7 // Maintainer: EV
8 //
9 // This file is part of the SC Toolkit.
10 //
11 // The SC Toolkit is free software; you can redistribute it and/or modify
12 // it under the terms of the GNU Library General Public License as published by
13 // the Free Software Foundation; either version 2, or (at your option)
14 // any later version.
15 //
16 // The SC Toolkit is distributed in the hope that it will be useful,
17 // but WITHOUT ANY WARRANTY; without even the implied warranty of
18 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
19 // GNU Library General Public License for more details.
20 //
21 // You should have received a copy of the GNU Library General Public License
22 // along with the SC Toolkit; see the file COPYING.LIB. If not, write to
23 // the Free Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA.
24 //
25 // The U.S. Government is granted a limited license as per AL 91-7.
26 //
27 
28 #ifdef __GNUG__
29 #pragma interface
30 #endif
31 
32 #ifndef _mpqc_src_lib_math_optimize_conjgrad_h
33 #define _mpqc_src_lib_math_optimize_conjgrad_h
34 
35 #include <assert.h>
36 #include <math/scmat/abstract.h>
37 #include <limits>
38 
39 namespace sc {
40 
52  template <typename F>
54  const RefSCMatrix& b,
55  RefSCMatrix& x,
56  const RefSCMatrix& preconditioner,
57  double convergence_target = -1.0) {
58 
59  const size_t n = x.nrow() * x.ncol();
60  MPQC_ASSERT(n == (preconditioner.nrow() * preconditioner.ncol()));
61 
62  // solution vector
63  RefSCMatrix XX_i;
64  // residual vector
65  RefSCMatrix RR_i = b.clone();
66  // preconditioned residual vector
67  RefSCMatrix ZZ_i;
68  // direction vector
69  RefSCMatrix PP_i;
70  RefSCMatrix APP_i = b.clone();
71 
72  // approximate the condition number as the ratio of the min and max elements of the preconditioner
73  // assuming that preconditioner is the approximate inverse of A in Ax - b =0
74  typedef SCElementFindExtremum<std::greater<double>, SCMatrixIterationRanges::AllElements> MinOp;
75  typedef SCElementFindExtremum<std::less<double>, SCMatrixIterationRanges::AllElements> MaxOp;
76  Ref<MinOp> findmin_op = new MinOp(1e10);
77  Ref<MaxOp> findmax_op = new MaxOp;
78  preconditioner.element_op(findmin_op);
79  preconditioner.element_op(findmax_op);
80  const double precond_min = findmin_op->result().at(0).value;
81  const double precond_max = findmax_op->result().at(0).value;
82  const double cond_number = precond_max / precond_min;
83  // if convergence target is given, estimate of how tightly the system can be converged
84  if (convergence_target < 0.0) {
85  convergence_target = 1e-15 * cond_number;
86  }
87  else { // else warn if the given system is not sufficiently well conditioned
88  if (convergence_target < 1e-15 * cond_number)
89  ExEnv::out0() << indent << "WARNING: linsolv_conjugate_gradient conv. target (" << convergence_target
90  << ") may be too low for 64-bit precision" << std::endl;
91  }
92 
93  bool converged = false;
94  const unsigned int max_niter = 500;
95  double rnorm2 = 0.0;
96  const unsigned int rhs_size = b.nrow() * b.ncol();
97 
98  // starting guess: x_0 = D^-1 . b
100  XX_i = b.copy();
101  XX_i.element_op(multiply_op, preconditioner);
102 
103  // r_0 = b - a(x)
104  a(XX_i, RR_i); // RR_i = a(XX_i)
105  RR_i.scale(-1.0);
106  RR_i.accumulate(b); // RR_i = b - a(XX_i)
107 
108  // z_0 = D^-1 . r_0
109  ZZ_i = RR_i.copy();
110  ZZ_i.element_op(multiply_op, preconditioner);
111 
112  // p_0 = z_0
113  PP_i = ZZ_i.copy();
114 
116 
117  unsigned int iter = 0;
118  while (not converged) {
119 
120  // alpha_i = (r_i . z_i) / (p_i . A . p_i)
121  scalarprod_op->init();
122  RR_i.element_op(scalarprod_op, ZZ_i);
123  const double rz_norm2 = scalarprod_op->result();
124  a(PP_i,APP_i);
125 
126  scalarprod_op->init();
127  PP_i.element_op(scalarprod_op, APP_i);
128  const double pAp_i = scalarprod_op->result();
129  const double alpha_i = rz_norm2 / pAp_i;
130 
131  // x_i += alpha_i p_i
132  Ref<SCElementDAXPY> daxpy_op = new SCElementDAXPY(alpha_i);
133  XX_i.element_op(daxpy_op, PP_i);
134 
135  // r_i -= alpha_i Ap_i
136  Ref<SCElementDAXPY> daxpy_op2 = new SCElementDAXPY(-alpha_i);
137  RR_i.element_op(daxpy_op2, APP_i);
138 
139  Ref<SCElementKNorm> norm2_op = new SCElementKNorm(2);
140  RR_i.element_op(norm2_op);
141  const double r_ip1_norm = norm2_op->result() / rhs_size;
142  if (r_ip1_norm < convergence_target) {
143  converged = true;
144  rnorm2 = r_ip1_norm;
145  }
146 
147  // z_i = D^-1 . r_i
148  ZZ_i.assign(RR_i);
149  ZZ_i.element_op(multiply_op, preconditioner);
150  scalarprod_op->init();
151  ZZ_i.element_op(scalarprod_op, RR_i);
152  const double rz_ip1_norm2 = scalarprod_op->result();
153 
154  const double beta_i = rz_ip1_norm2 / rz_norm2;
155 
156  // p_i = z_i+1 + beta_i p_i
157  // 1) scale p_i by beta_i
158  // 2) add z_i+1 (i.e. current contents of z_i)
159  PP_i.scale(beta_i);
160  Ref<SCElementDAXPY> daxpy_op3 = new SCElementDAXPY( 1.0 );
161  PP_i.element_op(daxpy_op3, ZZ_i);
162 
163  ++iter;
164  //std::cout << "iter=" << iter << " dnorm=" << r_ip1_norm << std::endl;
165 
166  if (iter >= max_niter) {
167  x.assign(XX_i);
168  throw MaxIterExceeded("linsolv_conjugate_gradient", __FILE__, __LINE__, max_niter);
169  }
170  } // solver loop
171 
172  x.assign(XX_i);
173 
174  return rnorm2;
175  }
176 
177  inline size_t size(const RefSCMatrix& m) {
178  return m.nrow() * m.ncol();
179  }
180 
181  inline RefSCMatrix clone(const RefSCMatrix& m) {
182  return m.clone();
183  }
184 
185  inline RefSCMatrix copy(const RefSCMatrix& m) {
186  return m.copy();
187  }
188 
189  inline double minabs_value(const RefSCMatrix& m) {
190  typedef SCElementFindExtremum<sc::abs_greater<double>, SCMatrixIterationRanges::AllElements> MinOp;
191  Ref<MinOp> findmin_op = new MinOp(std::numeric_limits<double>::max());
192  m.element_op(findmin_op);
193  return findmin_op->result().at(0).value;
194  }
195 
196  inline double maxabs_value(const RefSCMatrix& m) {
197  typedef SCElementFindExtremum<sc::abs_less<double>, SCMatrixIterationRanges::AllElements> MaxOp;
198  Ref<MaxOp> findmax_op = new MaxOp;
199  m.element_op(findmax_op);
200  return findmax_op->result().at(0).value;
201  }
202 
203  inline void vec_multiply(RefSCMatrix& m1, const RefSCMatrix& m2) {
204  Ref<SCElementOp2> multiply_op = new SCElementDestructiveProduct;
205  m1.element_op(multiply_op, m2);
206  }
207 
208  inline double dot_product(const RefSCMatrix& m1, const RefSCMatrix& m2) {
209  Ref<SCElementScalarProduct> scalarprod_op = new SCElementScalarProduct;
210  scalarprod_op->init();
211  m1.element_op(scalarprod_op, m2);
212  return scalarprod_op->result();
213  }
214 
215  inline void scale(RefSCMatrix& m, double scaling_factor) {
216  m.scale(scaling_factor);
217  }
218 
219  inline void daxpy(RefSCMatrix& y, double a, const RefSCMatrix& x) {
220  Ref<SCElementDAXPY> daxpy_op = new SCElementDAXPY(a);
221  y.element_op(daxpy_op, x);
222  }
223 
224  inline void assign(RefSCMatrix& m1, const RefSCMatrix& m2) {
225  m1.assign(m2);
226  }
227 
228  inline double norm2(const RefSCMatrix& m) {
229  Ref<SCElementKNorm> norm2_op = new SCElementKNorm(2);
230  m.element_op(norm2_op);
231  return norm2_op->result();
232  }
233 
234  inline void print(const RefSCMatrix& m, const char* label) {
235  m.print(label);
236  }
237 
263  template <typename D, typename F>
265  typedef typename D::element_type value_type;
266 
267  value_type operator()(F& a,
268  const D& b,
269  D& x,
270  const D& preconditioner,
271  value_type convergence_target = -1.0) {
272 
273  const size_t n = size(x);
274  MPQC_ASSERT(n == size(preconditioner));
275 
276  // solution vector
277  D XX_i;
278  // residual vector
279  D RR_i = clone(b);
280  // preconditioned residual vector
281  D ZZ_i;
282  // direction vector
283  D PP_i;
284  D APP_i = clone(b);
285 
286  // approximate the condition number as the ratio of the min and max elements of the preconditioner
287  // assuming that preconditioner is the approximate inverse of A in Ax - b =0
288  const value_type precond_min = minabs_value(preconditioner);
289  const value_type precond_max = maxabs_value(preconditioner);
290  const value_type cond_number = precond_max / precond_min;
291  //std::cout << "condition number = " << precond_max << " / " << precond_min << " = " << cond_number << std::endl;
292  // if convergence target is given, estimate of how tightly the system can be converged
293  if (convergence_target < 0.0) {
294  convergence_target = 1e-15 * cond_number;
295  }
296  else { // else warn if the given system is not sufficiently well conditioned
297  if (convergence_target < 1e-15 * cond_number)
298  std::cout << "WARNING: ConjugateGradient convergence target (" << convergence_target
299  << ") may be too low for 64-bit precision" << std::endl;
300  }
301 
302  bool converged = false;
303  const unsigned int max_niter = 500;
304  value_type rnorm2 = 0.0;
305  const size_t rhs_size = size(b);
306 
307  // starting guess: x_0 = D^-1 . b
309  XX_i = copy(b);
310  vec_multiply(XX_i, preconditioner);
311 
312  // r_0 = b - a(x)
313  a(XX_i, RR_i); // RR_i = a(XX_i)
314  scale(RR_i, -1.0);
315  daxpy(RR_i, 1.0, b); // RR_i = b - a(XX_i)
316 
317  // z_0 = D^-1 . r_0
318  ZZ_i = copy(RR_i);
319  vec_multiply(ZZ_i, preconditioner);
320 
321  // p_0 = z_0
322  PP_i = copy(ZZ_i);
323 
324  unsigned int iter = 0;
325  while (not converged) {
326 
327  // alpha_i = (r_i . z_i) / (p_i . A . p_i)
328  value_type rz_norm2 = dot_product(RR_i, ZZ_i);
329  a(PP_i,APP_i);
330 
331  const value_type pAp_i = dot_product(PP_i, APP_i);
332  const value_type alpha_i = rz_norm2 / pAp_i;
333 
334  // x_i += alpha_i p_i
335  daxpy(XX_i, alpha_i, PP_i);
336 
337  // r_i -= alpha_i Ap_i
338  daxpy(RR_i, -alpha_i, APP_i);
339 
340  const value_type r_ip1_norm = norm2(RR_i) / rhs_size;
341  if (r_ip1_norm < convergence_target) {
342  converged = true;
343  rnorm2 = r_ip1_norm;
344  }
345 
346  // z_i = D^-1 . r_i
347  ZZ_i = copy(RR_i);
348  vec_multiply(ZZ_i, preconditioner);
349 
350  const value_type rz_ip1_norm2 = dot_product(ZZ_i, RR_i);
351 
352  const value_type beta_i = rz_ip1_norm2 / rz_norm2;
353 
354  // p_i = z_i+1 + beta_i p_i
355  // 1) scale p_i by beta_i
356  // 2) add z_i+1 (i.e. current contents of z_i)
357  scale(PP_i, beta_i);
358  daxpy(PP_i, 1.0, ZZ_i);
359 
360  ++iter;
361  //std::cout << "iter=" << iter << " dnorm=" << r_ip1_norm << std::endl;
362 
363  if (iter >= max_niter) {
364  assign(x, XX_i);
365  throw std::domain_error("ConjugateGradient: max # of iterations exceeded");
366  }
367  } // solver loop
368 
369  assign(x, XX_i);
370 
371  return rnorm2;
372  }
373  };
374 
375 } // end of namespace sc
376 
377 #endif // end of header guard
378 
379 
380 // Local Variables:
381 // mode: c++
382 // c-file-style: "CLJ-CONDENSED"
383 // End:
sc::linsolv_conjugate_gradient
double linsolv_conjugate_gradient(F &a, const RefSCMatrix &b, RefSCMatrix &x, const RefSCMatrix &preconditioner, double convergence_target=-1.0)
Solves linear system a(x) = b using conjugate gradient solver where a is a linear function of x.
Definition: conjgrad.h:53
sc::RefSCMatrix
The RefSCMatrix class is a smart pointer to an SCMatrix specialization.
Definition: matrix.h:135
sc::Ref
A template class that maintains references counts.
Definition: ref.h:361
sc::RefSCMatrix::clone
RefSCMatrix clone() const
These call the SCMatrix members of the same name after checking for references to 0.
sc::SCElementDestructiveProduct
does
Definition: elemop.h:234
sc::ConjugateGradientSolver
Solves linear system a(x) = b using conjugate gradient solver where a is a linear function of x.
Definition: conjgrad.h:264
sc::SCElementKNorm
Computes k-norm of matrix.
Definition: elemop.h:499
sc::MaxIterExceeded
This is thrown when an iterative algorithm attempts to use more iterations than allowed.
Definition: scexception.h:371
sc::SCElementDAXPY
does
Definition: elemop.h:247
sc::SCElementScalarProduct
evaluates
Definition: elemop.h:213
sc::SCElementFindExtremum
Searches each range in IterationRanges for element i so that there is no element j in that Range for ...
Definition: elemop.h:358
sc::ExEnv::out0
static std::ostream & out0()
Return an ostream that writes from node 0.
sc
Contains all MPQC code up to version 3.
Definition: mpqcin.h:14

Generated at Sun Jan 26 2020 23:24:00 for MPQC 3.0.0-alpha using the documentation package Doxygen 1.8.16.