Am Montag, 13. Mai 2022, finden Wartungsarbeiten am Gitlab-Server (Update auf neue Version statt). Der Dienst wird daher am Montag für einige Zeit nicht verfügbar sein.
On Monday, May 13th 2022, the Gitlab server will be updated. The service will therefore not be accessible for some time on Monday.

gauss_seidel.hpp 4.78 KB
 Praetorius, Simon committed Oct 13, 2016 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 ``````// Software License for MTL // // Copyright (c) 2007 The Trustees of Indiana University. // 2008 Dresden University of Technology and the Trustees of Indiana University. // 2010 SimuNova UG (haftungsbeschränkt), www.simunova.com. // All rights reserved. // Authors: Peter Gottschling and Andrew Lumsdaine // // This file is part of the Matrix Template Library // // See also license.mtl.txt in the distribution. #ifndef ITL_GAUSS_SEIDEL_INCLUDE #define ITL_GAUSS_SEIDEL_INCLUDE #include #include #include #include #include #include #include #include #include namespace itl { /// Gauss-Seidel smoother /** Constructor takes references to a matrix and a right-hand side vector. operator() is applied on a vector and changes it in place. Matrix must be square, stored row-major and free of zero entries in the diagonal. Vectors b and x must have the same number of rows as A. **/ template class gauss_seidel { typedef typename mtl::Collection::value_type Scalar; typedef typename mtl::Collection::size_type size_type; public: /// Construct with constant references to matrix and RHS vector gauss_seidel(const Matrix& A) : A(A), dia_inv(num_rows(A)) { BOOST_STATIC_ASSERT((mtl::traits::is_row_major::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 Vector& operator()(Vector& x, const RHSVector& b) const { mtl::vampir_trace<8551> tracer; namespace tag= mtl::tag; using namespace mtl::traits; using mtl::begin; using mtl::end; typedef typename range_generator::type a_cur_type; typedef typename range_generator::type a_icur_type; typename col::type col_a(A); typename const_value::type value_a(A); typedef typename mtl::Collection::value_type value_type; a_cur_type ac= begin(A), aend= end(A); for (unsigned i= 0; ac != aend; ++ac, ++i) { value_type tmp= b[i]; for (a_icur_type aic= begin(ac), aiend= end(ac); aic != aiend; ++aic) if (col_a(*aic) != i) tmp-= value_a(*aic) * x[col_a(*aic)]; x[i]= dia_inv[i] * tmp; } return x; } private: const Matrix& A; mtl::vector::dense_vector dia_inv; }; template class gauss_seidel > { typedef mtl::matrix::compressed2D Matrix; typedef typename mtl::Collection::value_type Scalar; typedef typename mtl::Collection::size_type size_type; public: /// Construct with constant references to matrix and RHS vector gauss_seidel(const Matrix& A) : A(A), dia_inv(num_rows(A)), dia_pos(num_rows(A)) { BOOST_STATIC_ASSERT((mtl::traits::is_row_major::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 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 Vector& operator()(Vector& x, const RHSVector& b) const { mtl::vampir_trace<8551> tracer; typedef typename mtl::Collection::value_type value_type; typedef typename mtl::Collection::size_type size_type; 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]]; x[i]= dia_inv[i] * tmp; } return x; } private: const Matrix& A; mtl::vector::dense_vector dia_inv; mtl::vector::dense_vector dia_pos; }; } // namespace itl #endif // ITL_GAUSS_SEIDEL_INCLUDE``````