jor.hpp 4.84 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
30
// See also license.mtl.txt in the distribution.

#ifndef ITL_JOR_INCLUDE
#define ITL_JOR_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 "./relaxation_parameter.hpp"

namespace itl {

/// Jacobi smoother with relaxation
/** Constructor takes references to a matrix and a right-hand side vector.
31
    operator() is applied on a vector and changes it in place.
32
    Matrix must be square, stored row-major and free of zero entries in the diagonal.
33
    Vectors b and x must have the same number of rows as A.
34
35
36
37
38
39
40
41
**/
template <typename Matrix, typename Omega = default_omega>
class jor
{
    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
42
    jor(const Matrix& A) : A(A), dia_inv(num_rows(A))
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
    {
	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 JOR on vector \p x, i.e. \p x is changed
    template <typename Vector, typename RHSVector>
    Vector& operator()(Vector& x, const RHSVector& b) const
    {
	namespace tag= mtl::tag; using namespace mtl::traits;
	using mtl::begin; using mtl::end; using std::swap;

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

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

67
	static Vector x0(resource(x));
68
69
70
	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];
71
	    for (a_icur_type aic= begin<tag::nz>(ac), aiend= end<tag::nz>(ac); aic != aiend; ++aic)
72
		if (col_a(*aic) != i)
73
		    tmp-= value_a(*aic) * x[col_a(*aic)];
74
75
76
77
78
79
80
81
82
83
84
	    x0[i]= Omega::value * (dia_inv[i] * tmp) + (1. - Omega::value)*x0[i];
	}
	swap(x0, x);
 	return x;
    }

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

85

86
87
88
89
90
91
92
93
template <typename Value, typename Parameters, typename Omega>
class jor<mtl::matrix::compressed2D<Value, Parameters> , Omega>
{
    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
94
    jor(const Matrix& A)
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
      : 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 JOR on vector \p x, i.e. \p x is changed
    template <typename Vector, typename RHSVector>
    Vector& operator()(Vector& x, const RHSVector& b) const
    {
	typedef typename mtl::Collection<Vector>::value_type           value_type;
112
	typedef typename mtl::Collection<Matrix>::size_type            size_type;
113
	const size_type nr= num_rows(A);
114
	static Vector x0(resource(x));
115
116
117
118
119
120
121
122
123
	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]];
124
125
	    x0[i] = Omega::value * (dia_inv[i] * tmp) + (1. - Omega::value)*x0[i];
	}
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
	swap(x0, x);
 	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_JOR_INCLUDE