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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 4956 - (show annotations)
Tue May 20 12:34:41 2014 UTC (5 years ago) by caltinay
File MIME type: text/plain
File size: 11348 byte(s)
minor changes to cusp to allow compilation with -pedantic and -Werror.
One bug fix in spmv_dia.

1 /*
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 relative_tolerance_(relative_tolerance),
121 absolute_tolerance_(absolute_tolerance),
122 iteration_limit_(iteration_limit),
123 iteration_count_(0)
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 b_norm;
180 Real r_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