Stokhos Package Browser (Single Doxygen Collection) Version of the Day
Loading...
Searching...
No Matches
block_monitor.h
Go to the documentation of this file.
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
21#pragma once
22
23#include <cusp/detail/config.h>
24#include <cusp/array1d.h>
25#include <cusp/blas.h>
26#include <cusp/array2d.h>
27#include <cusp/print.h>
28
29#include <limits>
30#include <iostream>
31#include <iomanip>
32
33// Classes to monitor iterative solver progress, check for convergence, etc.
34// Follows the implementation of Iteration in the ITL:
35// http://www.osl.iu.edu/research/itl/doc/Iteration.html
36
37namespace cusp
38{
53template <typename ValueType>
55{
56public:
57 typedef typename norm_type<ValueType>::type Real;
58
73 template <typename MV>
75 size_t iteration_limit = 500,
78 bool verbose = true) :
79 numRHS(b.num_cols),
84 verbose_(verbose),
85 b_norm(b.num_cols)
86 {
87 for (int i = 0; i < numRHS; i++)
88 b_norm[i] = cusp::blas::nrm2(b.column(i));
89 }
90
93 void operator++(void) { ++iteration_count_; } // prefix increment
94
100 template <typename MV>
101 bool finished(const MV& r)
102 {
103
104 if (converged(r))
105 {
106 if (verbose_) {
107 cusp::array1d<ValueType, cusp::host_memory> resid(numRHS);
108 std::cout << "Successfully converged after " << iteration_count() << " iterations to tolerance " << tolerance(0) << std::endl;
109 std::cout << "with max residual norm ";
110 Real norm_max = 0;
111 for (int i = 0; i < numRHS; i++) {
112 resid[i] = cusp::blas::nrm2(r.column(i));
113 if (resid[i] > norm_max) norm_max = resid[i];
114 }
115 std::cout << norm_max << std::endl;
116
117 //cusp::print(resid);
118 }
119
120 return true;
121 }
122 else if (iteration_count() >= iteration_limit())
123 {
124 if (verbose_) {
125 cusp::array1d<ValueType, cusp::host_memory> resid(numRHS);
126 std::cout << "Failed to converge after " << iteration_count() << " iterations." << std::endl;
127 std::cout << "with max residual norm ";
128 Real norm_max = 0;
129 for (int i = 0; i < numRHS; i++) {
130 resid[i] = cusp::blas::nrm2(r.column(i));
131 if (resid[i] > norm_max) norm_max = resid[i];
132 }
133 std::cout << norm_max << std::endl;
134
135 //cusp::print(resid);
136 }
137
138 return true;
139 }
140 else
141 {
142 return false;
143 }
144
145
146 }
147
150 template <typename MV>
151 bool converged(MV& r) const
152 {
153 for (int i = 0; i < numRHS; i++){
154
155 if (cusp::blas::nrm2(r.column(i)) > tolerance(i)){
156 return false;
157 }
158 }
159
160 return true;
161 }
162
165 size_t iteration_count() const { return iteration_count_; }
166
169 size_t iteration_limit() const { return iteration_limit_; }
170
174
178
184 Real tolerance(int i) const { return absolute_tolerance() + relative_tolerance() * b_norm[i]; }
185
186protected:
187
191 size_t numRHS;
194 cusp::array1d<ValueType, cusp::host_memory> b_norm;
195
196};
197
207/*
208template <typename ValueType>
209class verbose_block_monitor : public default_block_monitor<ValueType>
210{
211 typedef typename norm_type<ValueType>::type Real;
212 typedef cusp::default_block_monitor<ValueType> super;
213
214 public:
215
216 *! Construct a \p verbose_monitor for a given right-hand-side \p b
217 *
218 * The \p verbose_monitor terminates iteration when the residual norm
219 * satisfies the condition
220 * ||b - A x|| <= absolute_tolerance + relative_tolerance * ||b||
221 * or when the iteration limit is reached.
222 *
223 * \param b right-hand-side of the linear system A x = b
224 * \param iteration_limit maximum number of solver iterations to allow
225 * \param relative_tolerance determines convergence criteria
226 * \param absolute_tolerance determines convergence criteria
227 *
228 * \tparam VectorType vector
229 *
230
231
232 template <typename MV>
233 verbose_block_monitor(const MV& b, size_t iteration_limit = 500, Real relative_tolerance = 1e-5, Real absolute_tolerance = 0)
234 : super(b, iteration_limit, relative_tolerance, absolute_tolerance)
235 {
236 std::cout << "Solver will continue until ";
237 std::cout << "residual norm" << super::tolerance() << " or reaching ";
238 std::cout << super::iteration_limit() << " iterations " << std::endl;
239 std::cout << " Iteration Number | Residual Norm" << std::endl;
240 }
241
242 template <typename MV>
243 bool finished(const MV& r)
244 {
245 for (int i = 0; i < r.num_cols; i++)
246 super::r_norm[i] = cusp::blas::nrm2(r.column(i));
247
248 std::cout << " " << std::setw(10) << super::iteration_count();
249 std::cout << " " << std::setw(10) << std::scientific << super::residual_norm_average() << std::endl;
250
251 if (super::converged(r))
252 {
253 std::cout << "Successfully converged after " << super::iteration_count() << " iterations." << std::endl;
254 return true;
255 }
256 else if (super::iteration_count() >= super::iteration_limit())
257 {
258 std::cout << "Failed to converge after " << super::iteration_count() << " iterations." << std::endl;
259 return true;
260 }
261 else
262 {
263 return false;
264 }
265 }
266};
267
268*/
269
273} // end namespace cusp
default_block_monitor(const MV &b, size_t iteration_limit=500, Real absolute_tolerance=1e-6, Real relative_tolerance=1e-6, bool verbose=true)
bool converged(MV &r) const
cusp::array1d< ValueType, cusp::host_memory > b_norm
norm_type< ValueType >::type Real
Real tolerance(int i) const