gauss_seidel.hpp 4.74 KB
Newer Older
1
// Software License for MTL
2
3
//
// Copyright (c) 2007 The Trustees of Indiana University.
4
//               2008 Dresden University of Technology and the Trustees of Indiana University.
5
//               2010 SimuNova UG (haftungsbeschränkt), www.simunova.com.
6
7
// All rights reserved.
// Authors: Peter Gottschling and Andrew Lumsdaine
8
//
9
// This file is part of the Matrix Template Library
10
//
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
// See also license.mtl.txt in the distribution.

#ifndef ITL_GAUSS_SEIDEL_INCLUDE
#define ITL_GAUSS_SEIDEL_INCLUDE

#include <boost/assert.hpp>
#include <boost/numeric/mtl/concept/collection.hpp>
#include <boost/numeric/mtl/utility/exception.hpp>
#include <boost/numeric/mtl/utility/is_row_major.hpp>
#include <boost/numeric/mtl/utility/property_map.hpp>
#include <boost/numeric/mtl/utility/range_generator.hpp>
#include <boost/numeric/mtl/utility/tag.hpp>
#include <boost/numeric/mtl/matrix/compressed2D.hpp>
#include <boost/numeric/mtl/interface/vpt.hpp>

namespace itl {

/// Gauss-Seidel smoother
/** Constructor takes references to a matrix and a right-hand side vector.
30
    operator() is applied on a vector and changes it in place.
31
    Matrix must be square, stored row-major and free of zero entries in the diagonal.
32
    Vectors b and x must have the same number of rows as A.
33
34
35
36
37
38
39
40
**/
template <typename Matrix>
class gauss_seidel
{
    typedef typename mtl::Collection<Matrix>::value_type Scalar;
    typedef typename mtl::Collection<Matrix>::size_type  size_type;
  public:
    /// Construct with constant references to matrix and RHS vector
41
    gauss_seidel(const Matrix& A) : A(A), dia_inv(num_rows(A))
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
    {
	BOOST_STATIC_ASSERT((mtl::traits::is_row_major<Matrix>::value)); // No CCS
	assert(num_rows(A) == num_cols(A)); // Matrix must be square
	for (size_type i= 0; i < num_rows(A); ++i) {
	    Scalar a= A[i][i];
	    MTL_THROW_IF(a == 0, mtl::missing_diagonal());
	    dia_inv[i]= 1.0 / a;
	}
    }

    /// Apply Gauss-Seidel on vector \p x, i.e. \p x is changed
    template <typename Vector, typename RHSVector>
    Vector& operator()(Vector& x, const RHSVector& b) const
    {
	mtl::vampir_trace<8551> tracer;
	namespace tag= mtl::tag; using namespace mtl::traits;
58
	using mtl::begin; using mtl::end;
59

60
61
62
63
        typedef typename range_generator<tag::row, Matrix>::type       a_cur_type;
        typedef typename range_generator<tag::nz, a_cur_type>::type    a_icur_type;
	typename col<Matrix>::type                   col_a(A);
	typename const_value<Matrix>::type           value_a(A);
64
65
66
67
68
69

	typedef typename mtl::Collection<Vector>::value_type           value_type;

	a_cur_type ac= begin<tag::row>(A), aend= end<tag::row>(A);
	for (unsigned i= 0; ac != aend; ++ac, ++i) {
	    value_type tmp= b[i];
70
	    for (a_icur_type aic= begin<tag::nz>(ac), aiend= end<tag::nz>(ac); aic != aiend; ++aic)
71
		if (col_a(*aic) != i)
72
		    tmp-= value_a(*aic) * x[col_a(*aic)];
73
74
75
76
77
78
79
80
81
82
	    x[i]= dia_inv[i] * tmp;
	}
 	return x;
    }

   private:
    const Matrix&    A;
    mtl::vector::dense_vector<Scalar>  dia_inv;
};

83

84
85
86
87
88
89
90
91
template <typename Value, typename Parameters>
class gauss_seidel<mtl::matrix::compressed2D<Value, Parameters> >
{
    typedef mtl::matrix::compressed2D<Value, Parameters> Matrix;
    typedef typename mtl::Collection<Matrix>::value_type Scalar;
    typedef typename mtl::Collection<Matrix>::size_type  size_type;
  public:
    /// Construct with constant references to matrix and RHS vector
92
    gauss_seidel(const Matrix& A)
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
      : A(A), dia_inv(num_rows(A)), dia_pos(num_rows(A))
    {
	BOOST_STATIC_ASSERT((mtl::traits::is_row_major<Matrix>::value)); // No CCS
	assert(num_rows(A) == num_cols(A)); // Matrix must be square
	for (size_type i= 0; i < num_rows(A); ++i) {
	    mtl::utilities::maybe<size_type> pos = A.indexer(A, i, i);
	    MTL_THROW_IF(!pos, mtl::missing_diagonal());
	    dia_inv[i]= 1.0 / A.value_from_offset(pos);
	    dia_pos[i]= pos;
	}
    }

    /// Apply Gauss-Seidel on vector \p x, i.e. \p x is changed
    template <typename Vector, typename RHSVector>
    Vector& operator()(Vector& x, const RHSVector& b) const
    {
	mtl::vampir_trace<8551> tracer;
	typedef typename mtl::Collection<Vector>::value_type           value_type;
111
	typedef typename mtl::Collection<Matrix>::size_type            size_type;
112
113
114
115
116
117
118
119
120
121
	const size_type nr= num_rows(A);
	size_type cj1= A.ref_major()[0];
	for (size_type i= 0; i < nr; ++i) {
	    value_type tmp= b[i];
	    size_type cj0= cj1, cjm= dia_pos[i];
	    cj1= A.ref_major()[i+1];
	    for (; cj0 < cjm; cj0++)
		tmp-= A.data[cj0] * x[A.ref_minor()[cj0]];
	    for (size_type j= cjm+1; j < cj1; j++)
		tmp-= A.data[j] * x[A.ref_minor()[j]];
122
123
	    x[i]= dia_inv[i] * tmp;
	}
124
125
126
127
128
129
130
131
132
133
134
135
136
137
 	return x;
    }


  private:
    const Matrix&    A;
    mtl::vector::dense_vector<Scalar>     dia_inv;
    mtl::vector::dense_vector<size_type>  dia_pos;
};


} // namespace itl

#endif // ITL_GAUSS_SEIDEL_INCLUDE