/[escript]/branches/diaplayground/cusplibrary/cusp/monitor.h
ViewVC logotype

Annotation of /branches/diaplayground/cusplibrary/cusp/monitor.h

Parent Directory Parent Directory | Revision Log Revision Log


Revision 4955 - (hide annotations)
Tue May 20 04:33:15 2014 UTC (6 years, 1 month ago) by caltinay
File MIME type: text/plain
File size: 11348 byte(s)
added pristine copy of cusplibrary (apache license) to be used by ripley.

1 caltinay 4955 /*
2     * Copyright 2008-2009 NVIDIA Corporation
3     *
4     * Licensed under the Apache License, Version 2.0 (the "License");
5     * you may not use this file except in compliance with the License.
6     * You may obtain a copy of the License at
7     *
8     * http://www.apache.org/licenses/LICENSE-2.0
9     *
10     * Unless required by applicable law or agreed to in writing, software
11     * distributed under the License is distributed on an "AS IS" BASIS,
12     * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
13     * See the License for the specific language governing permissions and
14     * limitations under the License.
15     */
16    
17     /*! \file monitor.h
18     * \brief Monitor iterative solver convergence
19     */
20    
21     #pragma once
22    
23     #include <cusp/detail/config.h>
24    
25     #include <cusp/blas.h>
26    
27     #include <limits>
28     #include <iostream>
29     #include <iomanip>
30    
31     // Classes to monitor iterative solver progress, check for convergence, etc.
32     // Follows the implementation of Iteration in the ITL:
33     // http://www.osl.iu.edu/research/itl/doc/Iteration.html
34    
35     namespace cusp
36     {
37     /*! \addtogroup iterative_solvers Iterative Solvers
38     * \addtogroup monitors Monitors
39     * \ingroup iterative_solvers
40     * \{
41     */
42    
43     /*! \p default_monitor : Implements standard convergence criteria
44     * and reporting for iterative solvers.
45     *
46     * \tparam ValueType scalar type used in the solver (e.g. \c float or \c cusp::complex<double>).
47     *
48     * The following code snippet demonstrates how to configure
49     * the \p default_monitor and use it with an iterative solver.
50     *
51     * \code
52     * #include <cusp/csr_matrix.h>
53     * #include <cusp/monitor.h>
54     * #include <cusp/krylov/cg.h>
55     * #include <cusp/gallery/poisson.h>
56     *
57     * int main(void)
58     * {
59     * // create an empty sparse matrix structure (CSR format)
60     * cusp::csr_matrix<int, float, cusp::device_memory> A;
61     *
62     * // initialize matrix
63     * cusp::gallery::poisson5pt(A, 10, 10);
64     *
65     * // allocate storage for solution (x) and right hand side (b)
66     * cusp::array1d<float, cusp::device_memory> x(A.num_rows, 0);
67     * cusp::array1d<float, cusp::device_memory> b(A.num_rows, 1);
68     *
69     * // set stopping criteria:
70     * // iteration_limit = 100
71     * // relative_tolerance = 1e-6
72     * cusp::default_monitor<float> monitor(b, 100, 1e-6);
73     *
74     * // solve the linear system A x = b
75     * cusp::krylov::cg(A, x, b, monitor);
76     *
77     * // report solver results
78     * if (monitor.converged())
79     * {
80     * std::cout << "Solver converged to " << monitor.relative_tolerance() << " relative tolerance";
81     * std::cout << " after " << monitor.iteration_count() << " iterations" << std::endl;
82     * }
83     * else
84     * {
85     * std::cout << "Solver reached iteration limit " << monitor.iteration_limit() << " before converging";
86     * std::cout << " to " << monitor.relative_tolerance() << " relative tolerance " << std::endl;
87     * }
88     *
89     * return 0;
90     * }
91     * \endcode
92     *
93     * \see \p verbose_monitor
94     *
95     */
96     template <typename ValueType>
97     class default_monitor
98     {
99     public:
100     typedef typename norm_type<ValueType>::type Real;
101    
102     /*! Construct a \p default_monitor for a given right-hand-side \p b
103     *
104     * The \p default_monitor terminates iteration when the residual norm
105     * satisfies the condition
106     * ||b - A x|| <= absolute_tolerance + relative_tolerance * ||b||
107     * or when the iteration limit is reached.
108     *
109     * \param b right-hand-side of the linear system A x = b
110     * \param iteration_limit maximum number of solver iterations to allow
111     * \param relative_tolerance determines convergence criteria
112     * \param absolute_tolerance determines convergence criteria
113     *
114     * \tparam VectorType vector
115     */
116     template <typename Vector>
117     default_monitor(const Vector& b, size_t iteration_limit = 500, Real relative_tolerance = 1e-5, Real absolute_tolerance = 0)
118     : b_norm(cusp::blas::nrm2(b)),
119     r_norm(std::numeric_limits<Real>::max()),
120     iteration_limit_(iteration_limit),
121     iteration_count_(0),
122     relative_tolerance_(relative_tolerance),
123     absolute_tolerance_(absolute_tolerance)
124     {}
125    
126     /*! increment the iteration count
127     */
128     void operator++(void) { ++iteration_count_; } // prefix increment
129    
130     /*! applies convergence criteria to determine whether iteration is finished
131     *
132     * \param r residual vector of the linear system (r = b - A x)
133     * \tparam Vector vector
134     */
135     template <typename Vector>
136     bool finished(const Vector& r)
137     {
138     r_norm = cusp::blas::nrm2(r);
139    
140     return converged() || iteration_count() >= iteration_limit();
141     }
142    
143     /*! whether the last tested residual satifies the convergence tolerance
144     */
145     bool converged() const
146     {
147     return residual_norm() <= tolerance();
148     }
149    
150     /*! Euclidean norm of last residual
151     */
152     Real residual_norm() const { return r_norm; }
153    
154     /*! number of iterations
155     */
156     size_t iteration_count() const { return iteration_count_; }
157    
158     /*! maximum number of iterations
159     */
160     size_t iteration_limit() const { return iteration_limit_; }
161    
162     /*! relative tolerance
163     */
164     Real relative_tolerance() const { return relative_tolerance_; }
165    
166     /*! absolute tolerance
167     */
168     Real absolute_tolerance() const { return absolute_tolerance_; }
169    
170     /*! tolerance
171     *
172     * Equal to absolute_tolerance() + relative_tolerance() * ||b||
173     *
174     */
175     Real tolerance() const { return absolute_tolerance() + relative_tolerance() * b_norm; }
176    
177     protected:
178    
179     Real r_norm;
180     Real b_norm;
181     Real relative_tolerance_;
182     Real absolute_tolerance_;
183    
184     size_t iteration_limit_;
185     size_t iteration_count_;
186     };
187    
188     /*! \p verbose_monitor is similar to \p default monitor except that
189     * it displays the solver status during iteration and reports a
190     * summary after iteration has stopped.
191     *
192     * \tparam ValueType scalar type used in the solver (e.g. \c float or \c cusp::complex<double>).
193     *
194     * \see \p default_monitor
195     */
196     template <typename ValueType>
197     class verbose_monitor : public default_monitor<ValueType>
198     {
199     typedef typename norm_type<ValueType>::type Real;
200     typedef cusp::default_monitor<ValueType> super;
201    
202     public:
203     /*! Construct a \p verbose_monitor for a given right-hand-side \p b
204     *
205     * The \p verbose_monitor terminates iteration when the residual norm
206     * satisfies the condition
207     * ||b - A x|| <= absolute_tolerance + relative_tolerance * ||b||
208     * or when the iteration limit is reached.
209     *
210     * \param b right-hand-side of the linear system A x = b
211     * \param iteration_limit maximum number of solver iterations to allow
212     * \param relative_tolerance determines convergence criteria
213     * \param absolute_tolerance determines convergence criteria
214     *
215     * \tparam VectorType vector
216     */
217     template <typename Vector>
218     verbose_monitor(const Vector& b, size_t iteration_limit = 500, Real relative_tolerance = 1e-5, Real absolute_tolerance = 0)
219     : super(b, iteration_limit, relative_tolerance, absolute_tolerance)
220     {
221     std::cout << "Solver will continue until ";
222     std::cout << "residual norm " << super::tolerance() << " or reaching ";
223     std::cout << super::iteration_limit() << " iterations " << std::endl;
224     std::cout << " Iteration Number | Residual Norm" << std::endl;
225     }
226    
227     template <typename Vector>
228     bool finished(const Vector& r)
229     {
230     super::r_norm = cusp::blas::nrm2(r);
231    
232     std::cout << " " << std::setw(10) << super::iteration_count();
233     std::cout << " " << std::setw(10) << std::scientific << super::residual_norm() << std::endl;
234    
235     if (super::converged())
236     {
237     std::cout << "Successfully converged after " << super::iteration_count() << " iterations." << std::endl;
238     return true;
239     }
240     else if (super::iteration_count() >= super::iteration_limit())
241     {
242     std::cout << "Failed to converge after " << super::iteration_count() << " iterations." << std::endl;
243     return true;
244     }
245     else
246     {
247     return false;
248     }
249     }
250     };
251     /*! \}
252     */
253    
254    
255     /*! \p convergence_monitor is similar to \p default monitor except that
256     * it displays the solver status during iteration and reports a
257     * summary after iteration has stopped.
258     *
259     * \tparam ValueType scalar type used in the solver (e.g. \c float or \c cusp::complex<double>).
260     *
261     * \see \p default_monitor
262     */
263     template <typename ValueType>
264     class convergence_monitor : public default_monitor<ValueType>
265     {
266     typedef typename norm_type<ValueType>::type Real;
267     typedef cusp::default_monitor<ValueType> super;
268    
269     public:
270     /*! Construct a \p convergence_monitor for a given right-hand-side \p b
271     *
272     * The \p convergence_monitor terminates iteration when the residual norm
273     * satisfies the condition
274     * ||b - A x|| <= absolute_tolerance + relative_tolerance * ||b||
275     * or when the iteration limit is reached.
276     *
277     * \param b right-hand-side of the linear system A x = b
278     * \param iteration_limit maximum number of solver iterations to allow
279     * \param relative_tolerance determines convergence criteria
280     * \param absolute_tolerance determines convergence criteria
281     *
282     * \tparam VectorType vector
283     */
284    
285     cusp::array1d<Real,cusp::host_memory> residuals;
286    
287     template <typename Vector>
288     convergence_monitor(const Vector& b, size_t iteration_limit = 500, Real relative_tolerance = 1e-5, Real absolute_tolerance = 0)
289     : super(b, iteration_limit, relative_tolerance, absolute_tolerance)
290     {
291     residuals.reserve(iteration_limit);
292     }
293    
294     template <typename Vector>
295     bool finished(const Vector& r)
296     {
297     super::r_norm = cusp::blas::nrm2(r);
298     residuals.push_back(super::r_norm);
299    
300     return super::converged() || super::iteration_count() >= super::iteration_limit();
301     }
302    
303     void print(void)
304     {
305     std::cout << "Solver will continue until ";
306     std::cout << "residual norm " << super::tolerance() << " or reaching ";
307     std::cout << super::iteration_limit() << " iterations " << std::endl;
308    
309     std::cout << "Ran " << super::iteration_count();
310     std::cout << " iterations with a final residual of ";
311     std::cout << super::r_norm << std::endl;
312    
313     std::cout << "geometric convergence factor : " << geometric_rate() << std::endl;
314     std::cout << "immediate convergence factor : " << immediate_rate() << std::endl;
315     std::cout << "average convergence factor : " << average_rate() << std::endl;
316     }
317    
318     Real immediate_rate(void)
319     {
320     size_t num = residuals.size();
321     return residuals[num-1] / residuals[num-2];
322     }
323    
324     Real geometric_rate(void)
325     {
326     size_t num = residuals.size();
327     return std::pow(residuals[num-1] / residuals[0], Real(1.0)/num);
328     }
329    
330     Real average_rate(void)
331     {
332     size_t num = residuals.size();
333     cusp::array1d<Real,cusp::host_memory> avg_vec(num-1);
334     thrust::transform(residuals.begin() + 1, residuals.end(), residuals.begin(), avg_vec.begin(), thrust::divides<Real>());
335     Real sum = thrust::reduce(avg_vec.begin(), avg_vec.end(), Real(0), thrust::plus<Real>());
336     return sum / Real(avg_vec.size());
337     }
338     };
339     /*! \}
340     */
341    
342     } // end namespace cusp
343    

  ViewVC Help
Powered by ViewVC 1.1.26