This file is indexed.

/usr/include/rheolef/pgmres.h is in librheolef-dev 6.5-1+b1.

This file is owned by root:root, with mode 0o644.

The actual contents of the file can be viewed below.

  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
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
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
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
# ifndef _SKIT_GMRES_H
# define _SKIT_GMRES_H
///
/// This file is part of Rheolef.
///
/// Copyright (C) 2000-2009 Pierre Saramito <Pierre.Saramito@imag.fr>
///
/// Rheolef is free software; you can redistribute it and/or modify
/// it under the terms of the GNU General Public License as published by
/// the Free Software Foundation; either version 2 of the License, or
/// (at your option) any later version.
///
/// Rheolef is distributed in the hope that it will be useful,
/// but WITHOUT ANY WARRANTY; without even the implied warranty of
/// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
/// GNU General Public License for more details.
///
/// You should have received a copy of the GNU General Public License
/// along with Rheolef; if not, write to the Free Software
/// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
/// 
/// =========================================================================

#ifdef HAVE_CONFIG_H
#include "rheolef/compiler.h"
#endif
#include <cmath> 

/*Class:pgmres
NAME: @code{pgmres} -- generalized minimum residual method 
@findex qmr
@cindex generalized minimum residual method
@cindex iterative solver
SYNOPSIS:
  @example
   template <class Matrix, class Vector, class Preconditioner,
     class SmallMatrix, class SmallVector, class Real, class Size>
   int pgmres (const Matrix &A, Vector &x, const Vector &b, const Preconditioner &M,
     SmallMatrix &H, const SmallVector& dummy,
     Size m, Size &max_iter, Real &tol);
  @end example
EXAMPLE:
  @noindent
  The simplest call to @code{pgmres} has the folling form:
  @example
        double tol = 1e-7;
        size_t max_iter = 100;
  	size_t m = 6;
  	boost::numeric::ublas::matrix<double> H(m+1,m+1);
  	vec<double,sequential> dummy;
  	int status = pgmres (a, x, b, ic0(a), H, dummy, m, max_iter, tol);
  @end example
DESCRIPTION:       
  @noindent
  @code{pgmres} solves the unsymmetric linear system Ax = b 
  using the generalized minimum residual method.

  @noindent
  The return value indicates convergence within max_iter (input)
  iterations (0), or no convergence within max_iter iterations (1).
  Upon successful return, output arguments have the following values:
  @table @code
    @item x
	approximate solution to Ax = b

    @item max_iter
	the number of iterations performed before the tolerance was reached

    @item tol
	the residual after the final iteration
  @end table
  @noindent
  In addition, M specifies a preconditioner, H specifies a matrix
  to hold the coefficients of the upper Hessenberg matrix constructed
  by the @code{pgmres} iterations, @code{m} specifies the number of iterations
  for each restart.

  @noindent
  @code{pgmres} requires two matrices as input, A and H.
  The matrix A, which will typically be a sparse matrix) corresponds
  to the matrix in the linear system Ax=b.
  The matrix H, which will be typically a dense matrix, corresponds
  to the upper Hessenberg matrix H that is constructed during the 
  @code{pgmres} iterations. Within @code{pgmres}, H is used in a different way
  than A, so its class must supply different functionality.
  That is, A is only accessed though its matrix-vector and 
  transpose-matrix-vector multiplication functions.
  On the other hand, @code{pgmres} solves a dense upper triangular linear 
  system of equations on H. Therefore, the class
  to which H belongs must provide H(i,j) operator for element acess.

NOTE: 
  @noindent
  It is important to remember that we use the convention that indices
  are 0-based. That is H(0,0) is the first component of the
  matrix H. Also, the type of the matrix must be compatible with the
  type of single vector entry. That is, operations such as
  H(i,j)*x(j) must be able to be carried out.
  
  @noindent
  @code{pgmres} is an iterative template routine.

  @noindent
  @code{pgmres} follows the algorithm described on p. 20 in
  @emph{Templates for the solution of linear systems: building blocks for iterative methods}, 
	2nd Edition, 
        R. Barrett, M. Berry, T. F. Chan, J. Demmel, J. Donato, J. Dongarra, V. Eijkhout,
	R. Pozo, C. Romine, H. Van der Vorst,
        SIAM, 1994, 
	@url{ftp.netlib.org/templates/templates.ps}.

  @noindent
  The present implementation is inspired from @code{IML++ 1.2} iterative method library,
  @url{http://math.nist.gov/iml++}.

AUTHOR: 

    Pierre Saramito
    | Pierre.Saramito@imag.fr
    LMC-IMAG, 38041 Grenoble cedex 9, France

DATE: 
    
    12 march 1997

METHODS: @pgmres
End:
*/
// ========== [ method implementation ] ====================================
namespace rheolef { 

//<pgmres:
template <class SmallMatrix, class Vector, class SmallVector, class Size>
void 
Update(Vector &x, Size k, SmallMatrix &h, SmallVector &s, Vector v[])
{
  SmallVector y = s;
  // Back solve:  
  for (int i = k; i >= 0; i--) {
    y(i) /= h(i,i);
    for (int j = i - 1; j >= 0; j--)
      y(j) -= h(j,i) * y(i);
  }
  for (Size j = 0; j <= k; j++) {
    x += v[j] * y(j);
  }
}
template<class Real> 
void GeneratePlaneRotation(Real &dx, Real &dy, Real &cs, Real &sn)
{
  if (dy == Real(0)) {
    cs = 1.0;
    sn = 0.0;
  } else if (abs(dy) > abs(dx)) {
    Real temp = dx / dy;
    sn = 1.0 / ::sqrt( 1.0 + temp*temp );
    cs = temp * sn;
  } else {
    Real temp = dy / dx;
    cs = 1.0 / ::sqrt( 1.0 + temp*temp );
    sn = temp * cs;
  }
}
template<class Real> 
void ApplyPlaneRotation(Real &dx, Real &dy, Real &cs, Real &sn)
{
  Real temp  =  cs * dx + sn * dy;
  dy = -sn * dx + cs * dy;
  dx = temp;
}
template <class Matrix, class Vector, class Preconditioner,
          class SmallMatrix, class SmallVector, class Real, class Size>
int
pgmres (const Matrix &A, Vector &x, const Vector &b, const Preconditioner &M,
      SmallMatrix &H, const SmallVector&, const Size &m, Size &max_iter, Real &tol)
{
  odiststream* p_derr = &derr;
  std::string label = "pgmres";
  if (p_derr) (*p_derr) << "[" << label << "] #iteration residue" << std::endl;
  Vector w;
  SmallVector s(m+1), cs(m+1), sn(m+1);
  Size i;
  Size j = 1;
  Size k;
  Real resid;
  Real normb = norm(M.solve(b));
  Vector r = M.solve(b - A * x);
  Real beta = norm(r);
  if (normb == Real(0)) {
    normb = 1;
  } 
  resid = norm(r) / normb;
  if (resid  <= tol) {
    tol = resid;
    max_iter = 0;
    return 0;
  }
  Vector *v = new Vector[m+1];
  while (j <= max_iter) {
    v[0] = r * (1.0 / beta);    // ??? r / beta
    s = 0.0;
    s(0) = beta;
    for (i = 0; i < m && j <= max_iter; i++, j++) {
      w = M.solve(A * v[i]);
      for (k = 0; k <= i; k++) {
        H(k, i) = dot(w, v[k]);
        w -= H(k, i) * v[k];
      }
      H(i+1, i) = norm(w);
      v[i+1] = w * (1.0 / H(i+1, i)); // ??? w / H(i+1, i)
      for (k = 0; k < i; k++) {
        ApplyPlaneRotation(H(k,i), H(k+1,i), cs(k), sn(k));
      }
      GeneratePlaneRotation(H(i,i), H(i+1,i), cs(i), sn(i));
      ApplyPlaneRotation(H(i,i), H(i+1,i), cs(i), sn(i));
      ApplyPlaneRotation(s(i), s(i+1), cs(i), sn(i));
      resid = abs(s(i+1)) / normb;
      if (p_derr) (*p_derr) << "[" << label << "] " << j << " " << resid << std::endl;
      if (resid  < tol) {
        Update(x, i, H, s, v);
        tol = resid;
        max_iter = j;
        delete [] v;
        return 0;
      }
    }
    Update(x, m - 1, H, s, v);
    r = M.solve(b - A * x);
    beta = norm(r);
    resid = beta / normb;
    if (p_derr) (*p_derr) << "[" << label << "] " << j << " " << resid << std::endl;
    if (resid < tol) {
      tol = resid;
      max_iter = j;
      delete [] v;
      return 0;
    }
  }
  tol = resid;
  delete [] v;
  return 1;
}
//>pgmres:
}// namespace rheolef
# endif // _SKIT_GMRES_H