fgmres.hpp 5.85 KB
Newer Older
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
/******************************************************************************
 *
 * AMDiS - Adaptive multidimensional simulations
 *
 * Copyright (C) 2013 Dresden University of Technology. All Rights Reserved.
 * Web: https://fusionforge.zih.tu-dresden.de/projects/amdis
 *
 * Authors: 
 * Simon Vey, Thomas Witkowski, Andreas Naumann, Simon Praetorius, et al.
 *
 * This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE
 * WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.
 *
 *
 * This file is part of AMDiS
 *
 * See also license.opensource.txt in the distribution.
 * 
 ******************************************************************************/

// Written by Simon Praetorius


#ifndef ITL_FGMRES_INCLUDE
#define ITL_FGMRES_INCLUDE

#include <algorithm>
#include <boost/numeric/mtl/concept/collection.hpp>
#include <boost/numeric/mtl/vector/dense_vector.hpp>
#include <boost/numeric/mtl/matrix/dense2D.hpp>
#include <boost/numeric/mtl/matrix/multi_vector.hpp>
#include <boost/numeric/mtl/operation/givens.hpp>
#include <boost/numeric/mtl/operation/two_norm.hpp>
#include <boost/numeric/mtl/utility/exception.hpp>
#include <boost/numeric/mtl/utility/irange.hpp>

#include "details.hpp"

namespace itl {
  
/// Flexible Generalized Minimal Residual method (without restart)
/// Cite: Youcef Saad, A Flexible Inner-Outer Preconditioned GMRES Algorithm, 1993
/** It computes at most kmax_in iterations (or size(x) depending on what is smaller) 
    regardless on whether the termination criterion is reached or not.   **/
template < typename Matrix, typename Vector, typename Preconditioner, typename Iteration >
int fgmres_full(const Matrix &A, Vector &x, const Vector &b,
               Preconditioner &P, Iteration& iter)
{
    using mtl::irange; using mtl::iall; using std::abs; using std::sqrt; using math::reciprocal;
    typedef typename mtl::Collection<Vector>::value_type Scalar;
    typedef typename mtl::Collection<Vector>::size_type  Size;

    if (size(b) == 0) throw mtl::logic_error("empty rhs vector");
    
55
    const Scalar                zero= math::zero(Scalar()), breakdown_tol= 1.e-16, kappa = 10.0;
56
57
    Scalar                      rho, bnrm2, temp, hr;
    Size                        k, kmax(std::min(size(x), Size(iter.max_iterations() - iter.iterations())));
58
    Vector                      w(resource(x)), r(b - A*x);
59
60
61
62
63
64
    mtl::matrix::multi_vector<Vector>   V(Vector(resource(x), zero), kmax+1); 
    mtl::matrix::multi_vector<Vector>   Z(Vector(resource(x), zero), kmax+1); 
    mtl::matrix::dense2D<Scalar>        H(kmax+1, kmax);
    
    mtl::vector::dense_vector<Scalar>   sn(kmax, zero), cs(kmax, zero), s(kmax+1, zero), y(kmax, zero);  // replicated in distributed solvers 

65
66
67
68
69
70
    bnrm2 = two_norm(b);
    if (bnrm2 == zero) {
      set_to_zero(x);
      // b == 0 => solution = 0
      return iter.terminate(bnrm2);
    }
71
    
72
73
    rho = two_norm(r);				// norm of preconditioned residual
    if (iter.finished(rho))			// initial guess is good enough solution
74
75
      return iter;
    
76
    V.vector(0) = r * reciprocal(rho);
77
    H = zero;
78
    s[0] = rho;
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

    // FGMRES iteration
    for (k= 0; k < kmax && !iter.finished(rho) ; ++k, ++iter) {
      Z.vector(k) = solve(P, V.vector(k));
      w = A * Z.vector(k);
      temp = two_norm(w);
      
      for (Size j= 0; j < k+1; j++) {
	H[j][k] = dot(V.vector(j), w);
	w -= H[j][k] * V.vector(j);
      }
      H[k+1][k]= two_norm(w);
            
      // reorthogonalization, only if "heuristic" condition is fulfilled
      if (H[k+1][k] < temp * reciprocal(kappa)) {
	for (Size i= 0; i < k+1; i++) {
	    hr = dot(w, V.vector(i));
            H[i][k] += hr;
            w -= hr * V.vector(i);
	}
        temp = two_norm(w);
	
	if (temp < H[k+1][k] * reciprocal(kappa)) {
	  set_to_zero(w);
	  H[k+1][k] = 0.0;
	} else {
	  H[k+1][k] = temp;
	}
      }
      
109
      if (H[k+1][k] < breakdown_tol)
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
	return iter.fail(2, "FGMRES: Singular matrix - nearly hard breakdown");
      
      V.vector(k+1) = w * reciprocal(H[k+1][k]);

      // k Given's rotations
      for(Size i= 0; i < k; i++) {
	temp =        cs[i]*H[i][k] + sn[i]*H[i+1][k];
	H[i+1][k] = - sn[i]*H[i][k] + cs[i]*H[i+1][k];
	H[i][k] = temp;
      }
      
      details::rotmat(H[k][k], H[k+1][k], cs[k], sn[k]);
      
      s[k+1] = -sn[k]*s[k];
      s[k]   =  cs[k]*s[k];
      H[k][k] = cs[k]*H[k][k] + sn[k]*H[k+1][k];
      H[k+1][k] = 0.0;
      
128
      rho = std::abs(s[k+1]);
129
130
131
    }
    
    // reduce k, to get regular matrix
132
//     while (k > 0 && abs(s[k-1]) <= iter.atol()) k--;
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150

    // iteration is finished -> compute x: solve H*y=g as far as rank of H allows
    irange range(k);
    for (; !range.empty(); --range) {
      try {
	  y[range] = upper_trisolve(H[range][range], s[range]); 
      } catch (mtl::matrix_singular) { continue; } // if singular then try with sub-matrix
      break;
    }

    if (range.finish() < k)
  	std::cerr << "GMRES orhogonalized with " << k << " vectors but matrix singular, can only use " 
		  << range.finish() << " vectors!\n";
    if (range.empty())
        return iter.fail(3, "GMRES did not find any direction to correct x");
    
    x += Z.vector(range) * y[range];
    
151
152
//     r = b - A*x;
    return iter.terminate(rho);
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
}

/// Flexible Generalized Minimal Residual method with restart
template < typename Matrix, typename Vector, typename LeftPreconditioner,
           typename RightPreconditioner, typename Iteration >
int fgmres(const Matrix &A, Vector &x, const Vector &b,
          LeftPreconditioner &L, RightPreconditioner &R,
	  Iteration& iter, typename mtl::Collection<Vector>::size_type restart)
{   
    do {
	Iteration inner(iter);
	inner.set_max_iterations(std::min(int(iter.iterations()+restart), iter.max_iterations()));
	inner.suppress_resume(true);
	fgmres_full(A, x, b, R, inner);
	iter.update_progress(inner);
    } while (!iter.finished());

    return iter;
}



} // namespace itl

#endif // ITL_FGMRES_INCLUDE