IterativeSolverT.cc 5.65 KB
Newer Older
David Bommes's avatar
David Bommes committed
1
2
3
4
5
6
7
8
9
10
11
12
13
14
//=============================================================================
//
//  CLASS IterativeSolverT - IMPLEMENTATION
//
//=============================================================================

#define COMISO_ITERATIVESOLVERT_C

//== INCLUDES =================================================================

#include "IterativeSolverT.hh"

//== NAMESPACES ===============================================================

15
16
namespace COMISO
{
David Bommes's avatar
David Bommes committed
17
18
19
20

//== IMPLEMENTATION ==========================================================

template <class RealT>
21
22
23
bool IterativeSolverT<RealT>::gauss_seidel_local(const Matrix& _A, Vector& _x,
    const Vector& _rhs, const std::vector<unsigned int>& _idxs,
    const int _max_iter, const Real& _tolerance)
David Bommes's avatar
David Bommes committed
24
{
25
26
  if (_max_iter == 0)
    return false;
David Bommes's avatar
David Bommes committed
27

28
  typedef typename gmm::linalg_traits<Matrix>::const_sub_col_type ColT;
David Bommes's avatar
David Bommes committed
29
30
31
32
33
34
  typedef typename gmm::linalg_traits<ColT>::const_iterator CIter;

  // clear old data
  i_temp.clear();
  q.clear();

35
36
  for (unsigned int i = 0; i < _idxs.size(); ++i)
    q.push_back(_idxs[i]);
David Bommes's avatar
David Bommes committed
37
38
39

  int it_count = 0;

40
  while (!q.empty() && it_count < _max_iter)
David Bommes's avatar
David Bommes committed
41
42
  {
    ++it_count;
43
    const auto i = q.front();
David Bommes's avatar
David Bommes committed
44
45
46
    q.pop_front();
    i_temp.clear();

47
48
49
    double res_i = -_rhs[i];
    double x_i_new = _rhs[i];
    double diag = 1.0;
David Bommes's avatar
David Bommes committed
50

51
52
53
    const ColT col = mat_const_col(_A, i);
    for (auto it = gmm::vect_const_begin(col), ite = gmm::vect_const_end(col);
         it != ite; ++it)
David Bommes's avatar
David Bommes committed
54
    {
55
56
57
58
59
      const auto j = static_cast<unsigned>(it.index());
      res_i += (*it) * _x[j];
      x_i_new -= (*it) * _x[j];
      if (j != i)
        i_temp.push_back(j);
David Bommes's avatar
David Bommes committed
60
      else
61
        diag = *it;
David Bommes's avatar
David Bommes committed
62
63
64
    }

    // take inverse of diag
65
    diag = 1.0 / diag;
David Bommes's avatar
David Bommes committed
66
67

    // compare relative residuum normalized by diagonal entry
68
    if (fabs(res_i * diag) > _tolerance)
David Bommes's avatar
David Bommes committed
69
    {
70
      _x[i] += x_i_new * diag;
David Bommes's avatar
David Bommes committed
71

72
73
      for (unsigned int j = 0; j < i_temp.size(); ++j)
        q.push_back(i_temp[j]);
David Bommes's avatar
David Bommes committed
74
75
76
    }
  }

77
  return q.empty(); // converged?
David Bommes's avatar
David Bommes committed
78
79
80
81
82
}

//-----------------------------------------------------------------------------

template <class RealT>
83
84
85
bool IterativeSolverT<RealT>::gauss_seidel_local2(const Matrix& _A, Vector& _x,
    const Vector& _rhs, const std::vector<unsigned int>& _idxs,
    const int _max_iter, const Real& _tolerance)
David Bommes's avatar
David Bommes committed
86
{
87
  typedef typename gmm::linalg_traits<Matrix>::const_sub_col_type ColT;
David Bommes's avatar
David Bommes committed
88
89
  typedef typename gmm::linalg_traits<ColT>::const_iterator CIter;

90
  double t2 = _tolerance * _tolerance;
David Bommes's avatar
David Bommes committed
91
92
93
94
95

  // clear old data
  i_temp.clear();
  s.clear();

96
97
  for (unsigned int i = 0; i < _idxs.size(); ++i)
    s.insert(_idxs[i]);
David Bommes's avatar
David Bommes committed
98
99
100
101
102

  int it_count = 0;

  bool finished = false;

103
  while (!finished && it_count < _max_iter)
David Bommes's avatar
David Bommes committed
104
105
106
  {
    finished = true;
    std::set<int>::iterator s_it = s.begin();
107
    for (; s_it != s.end(); ++s_it)
David Bommes's avatar
David Bommes committed
108
109
110
111
112
    {
      ++it_count;
      unsigned int cur_i = *s_it;
      i_temp.clear();

113
      ColT col = mat_const_col(_A, cur_i);
David Bommes's avatar
David Bommes committed
114

115
116
117
118
      CIter it = gmm::vect_const_begin(col);
      CIter ite = gmm::vect_const_end(col);

      double res_i = -_rhs[cur_i];
David Bommes's avatar
David Bommes committed
119
      double x_i_new = _rhs[cur_i];
120
121
      double diag = 1.0;
      for (; it != ite; ++it)
David Bommes's avatar
David Bommes committed
122
      {
123
124
125
126
127
128
        res_i += (*it) * _x[it.index()];
        x_i_new -= (*it) * _x[it.index()];
        if (it.index() != cur_i)
          i_temp.push_back(static_cast<int>(it.index()));
        else
          diag = *it;
David Bommes's avatar
David Bommes committed
129
130
131
      }

      // compare relative residuum normalized by diagonal entry
132
      if (res_i * res_i / diag > t2)
David Bommes's avatar
David Bommes committed
133
      {
134
135
136
137
138
139
140
        _x[cur_i] += x_i_new / _A(cur_i, cur_i);

        for (unsigned int j = 0; j < i_temp.size(); ++j)
        {
          finished = false;
          s.insert(i_temp[j]);
        }
David Bommes's avatar
David Bommes committed
141
142
143
144
145
146
147
      }
    }
  }

  // converged?
  return finished;
}
148

David Bommes's avatar
David Bommes committed
149
150
151
//-----------------------------------------------------------------------------

template <class RealT>
152
153
bool IterativeSolverT<RealT>::conjugate_gradient(const Matrix& _A, Vector& _x,
    const Vector& _rhs, int& _max_iter, Real& _tolerance)
David Bommes's avatar
David Bommes committed
154
155
156
157
158
159
160
161
{
  Real rho, rho_1(0), a;

  // initialize vectors
  p_.resize(_x.size());
  q_.resize(_x.size());
  r_.resize(_x.size());
  d_.resize(_x.size());
162
  gmm::copy(_x, p_);
David Bommes's avatar
David Bommes committed
163
164

  // initialize diagonal (for relative norm)
165
166
  for (unsigned int i = 0; i < _x.size(); ++i)
    d_[i] = 1.0 / _A(i, i);
David Bommes's avatar
David Bommes committed
167
168
169
170
171

  // start with iteration 0
  int cur_iter(0);

  gmm::mult(_A, gmm::scaled(_x, Real(-1)), _rhs, r_);
172
  rho = gmm::vect_sp(r_, r_);
David Bommes's avatar
David Bommes committed
173
  gmm::copy(r_, p_);
174

David Bommes's avatar
David Bommes committed
175
176
177
178
  bool not_converged = true;
  Real res_norm(0);

  // while not converged
179
180
  while ((not_converged = ((res_norm = vect_norm_rel(r_, d_)) > _tolerance)) &&
         cur_iter < _max_iter)
David Bommes's avatar
David Bommes committed
181
182
183
  {
    //    std::cerr << "iter " << cur_iter << "  res " << res_norm << std::endl;

184
185
186
    if (cur_iter != 0)
    {
      rho = gmm::vect_sp(r_, r_);
David Bommes's avatar
David Bommes committed
187
188
189
190
191
      gmm::add(r_, gmm::scaled(p_, rho / rho_1), p_);
    }

    gmm::mult(_A, p_, q_);

192
    a = rho / gmm::vect_sp(q_, p_);
David Bommes's avatar
David Bommes committed
193
194
195
196
197
198
    gmm::add(gmm::scaled(p_, a), _x);
    gmm::add(gmm::scaled(q_, -a), r_);
    rho_1 = rho;

    ++cur_iter;
  }
199
200

  _max_iter = cur_iter;
David Bommes's avatar
David Bommes committed
201
202
203
204
205
206
207
208
  _tolerance = res_norm;

  return (!not_converged);
}

//-----------------------------------------------------------------------------

template <class RealT>
209
210
typename IterativeSolverT<RealT>::Real IterativeSolverT<RealT>::vect_norm_rel(
    const Vector& _v, const Vector& _diag) const
David Bommes's avatar
David Bommes committed
211
212
213
{
  Real res = 0.0;

214
  for (unsigned int i = 0; i < _v.size(); ++i)
David Bommes's avatar
David Bommes committed
215
  {
216
    res = std::max(fabs(_v[i] * _diag[i]), res);
David Bommes's avatar
David Bommes committed
217

218
219
220
    //     Real cur = fabs(_v[i]*_diag[i]);
    //     if(cur > res)
    //       res = cur;
David Bommes's avatar
David Bommes committed
221
222
223
224
225
226
227
228
229
  }
  return res;
}

//-----------------------------------------------------------------------------

//=============================================================================
} // namespace COMISO
//=============================================================================