/usr/include/rheolef/puzawa.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 | # 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
|