This file is indexed.

/usr/include/rheolef/puzawa.h is in librheolef-dev 6.7-6.

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
# ifndef _SKIT_PUZAWA_H
# define _SKIT_PUZAWA_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
/// 
/// =========================================================================
#include "rheolef/diststream.h"

namespace rheolef { 
/*D:puzawa
NAME: @code{puzawa} -- Uzawa algorithm.
@findex puzawa
@cindex Uzawa algorithm
@cindex iterative solver
@cindex preconditioner
SYNOPSIS:
  @example
    template <class Matrix, class Vector, class Preconditioner, class Real>
    int puzawa (const Matrix &A, Vector &x, const Vector &b, const Preconditioner &M,
      int &max_iter, Real &tol, const Real& rho, odiststream *p_derr=0);
  @end example

EXAMPLE:
  @noindent
  The simplest call to 'puzawa' has the folling form:
  @example
    size_t max_iter = 100;
    double tol = 1e-7;
    int status = puzawa(A, x, b, EYE, max_iter, tol, 1.0, &derr);
  @end example

DESCRIPTION:       
  @noindent
  @code{puzawa} solves the linear
  system A*x=b using the Uzawa method. The Uzawa method is a
  descent method in the direction opposite to the  gradient,
  with a constant step length 'rho'. The convergence is assured
  when the step length 'rho' is small enough.
  If matrix A is symmetric positive definite, please uses 'pcg' that
  computes automatically the optimal descdnt step length at
  each iteration.

  @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

AUTHOR: 
    Pierre Saramito
    | Pierre.Saramito@imag.fr
    LJK-IMAG, 38041 Grenoble cedex 9, France
DATE: 
    20 april 2009
METHODS: @puzawa
End:
*/
//<puzawa:
template < class Matrix, class Vector, class Preconditioner, class Real, class Size>
int puzawa(const Matrix &A, Vector &x, const Vector &Mb, const Preconditioner &M,
    Size &max_iter, Real &tol, const Real& rho,
    odiststream *p_derr, std::string label)
{
    Vector b = M.solve(Mb);
    Real norm2_b = dot(Mb,b);
    Real norm2_r = norm2_b;
    if (norm2_b == Real(0)) norm2_b = 1;
    if (p_derr) (*p_derr) << "[" << label << "] #iteration residue" << std::endl;
    for (Size n = 0; n <= max_iter; n++) {
	Vector Mr = A*x - Mb;
	Vector r = M.solve(Mr);
        norm2_r = dot(Mr, r);
	if (p_derr) (*p_derr) << "[" << label << "] " << n << " " << sqrt(norm2_r/norm2_b) << std::endl;
        if (norm2_r <= sqr(tol)*norm2_b) {
          tol = sqrt(norm2_r/norm2_b);
          max_iter = n;
          return 0;     
        }
        x  -= rho*r;
    }
    tol = sqrt(norm2_r/norm2_b);
    return 1;
}
//>puzawa:
}// namespace rheolef
# endif // _SKIT_PUZAWA_H