← Back to branch summary

~ubuntu-branches/ubuntu/trusty/rheolef/trusty

« back to all changes in this revision

Viewing changes to skit/plib2/puzawa.h

  • Committer: Package Import Robot
  • Author(s): Pierre Saramito
  • Date: 2012-04-06 09:12:21 UTC
  • mfrom: (1.1.5)
  • Revision ID: package-import@ubuntu.com-20120406091221-m58me99p1nxqui49
Tags: 6.0-1
http://bugs.debian.org/667356
http://bugs.debian.org/661689
* New upstream release 6.0 (major changes):
  - massively distributed and parallel support
  - full FEM characteristic method (Lagrange-Gakerkin method) support
  - enhanced users documentation 
  - source code supports g++-4.7 (closes: #667356)
* debian/control: dependencies for MPI distributed solvers added
* debian/rules: build commands simplified
* debian/librheolef-dev.install: man1/* to man9/* added
* debian/changelog: package description rewritted (closes: #661689)

Show diffs side-by-side

added added

removed removed

Lines of Context:
skit/plib2/puzawa.h
 
1
# ifndef _SKIT_PUZAWA_H
 
2
# define _SKIT_PUZAWA_H
 
3
///
 
4
/// This file is part of Rheolef.
 
5
///
 
6
/// Copyright (C) 2000-2009 Pierre Saramito <Pierre.Saramito@imag.fr>
 
7
///
 
8
/// Rheolef is free software; you can redistribute it and/or modify
 
9
/// it under the terms of the GNU General Public License as published by
 
10
/// the Free Software Foundation; either version 2 of the License, or
 
11
/// (at your option) any later version.
 
12
///
 
13
/// Rheolef is distributed in the hope that it will be useful,
 
14
/// but WITHOUT ANY WARRANTY; without even the implied warranty of
 
15
/// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 
16
/// GNU General Public License for more details.
 
17
///
 
18
/// You should have received a copy of the GNU General Public License
 
19
/// along with Rheolef; if not, write to the Free Software
 
20
/// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
 
21
/// 
 
22
/// =========================================================================
 
23
#include "rheolef/diststream.h"
 
24
 
 
25
namespace rheolef { 
 
26
/*D:puzawa
 
27
NAME: @code{puzawa} -- Uzawa algorithm.
 
28
@findex puzawa
 
29
@cindex Uzawa algorithm
 
30
@cindex iterative solver
 
31
@cindex preconditioner
 
32
SYNOPSIS:
 
33
  @example
 
34
    template <class Matrix, class Vector, class Preconditioner, class Real>
 
35
    int puzawa (const Matrix &A, Vector &x, const Vector &b, const Preconditioner &M,
 
36
      int &max_iter, Real &tol, const Real& rho, odiststream *p_derr=0);
 
37
  @end example
 
38
 
 
39
EXAMPLE:
 
40
  @noindent
 
41
  The simplest call to 'puzawa' has the folling form:
 
42
  @example
 
43
    size_t max_iter = 100;
 
44
    double tol = 1e-7;
 
45
    int status = puzawa(A, x, b, EYE, max_iter, tol, 1.0, &derr);
 
46
  @end example
 
47
 
 
48
DESCRIPTION:       
 
49
  @noindent
 
50
  @code{puzawa} solves the linear
 
51
  system A*x=b using the Uzawa method. The Uzawa method is a
 
52
  descent method in the direction opposite to the  gradient,
 
53
  with a constant step length 'rho'. The convergence is assured
 
54
  when the step length 'rho' is small enough.
 
55
  If matrix A is symmetric positive definite, please uses 'pcg' that
 
56
  computes automatically the optimal descdnt step length at
 
57
  each iteration.
 
58
 
 
59
  @noindent
 
60
  The return value indicates convergence within max_iter (input)
 
61
  iterations (0), or no convergence within max_iter iterations (1).
 
62
  Upon successful return, output arguments have the following values:
 
63
  @table @code
 
64
    @itemx x
 
65
        approximate solution to Ax = b
 
66
 
 
67
    @itemx max_iter
 
68
        the number of iterations performed before the tolerance was reached
 
69
 
 
70
    @itemx tol
 
71
        the residual after the final iteration
 
72
  @end table
 
73
 
 
74
AUTHOR: 
 
75
    Pierre Saramito
 
76
    | Pierre.Saramito@imag.fr
 
77
    LJK-IMAG, 38041 Grenoble cedex 9, France
 
78
DATE: 
 
79
    20 april 2009
 
80
METHODS: @puzawa
 
81
End:
 
82
*/
 
83
//<puzawa:
 
84
template < class Matrix, class Vector, class Preconditioner, class Real, class Size>
 
85
int puzawa(const Matrix &A, Vector &x, const Vector &Mb, const Preconditioner &M,
 
86
    Size &max_iter, Real &tol, const Real& rho,
 
87
    odiststream *p_derr, std::string label)
 
88
{
 
89
    Vector b = M.solve(Mb);
 
90
    Real norm2_b = dot(Mb,b);
 
91
    Real norm2_r = norm2_b;
 
92
    if (norm2_b == Real(0)) norm2_b = 1;
 
93
    if (p_derr) (*p_derr) << "[" << label << "] #iteration residue" << std::endl;
 
94
    for (Size n = 0; n <= max_iter; n++) {
 
95
        Vector Mr = A*x - Mb;
 
96
        Vector r = M.solve(Mr);
 
97
        norm2_r = dot(Mr, r);
 
98
        if (p_derr) (*p_derr) << "[" << label << "] " << n << " " << sqrt(norm2_r/norm2_b) << std::endl;
 
99
        if (norm2_r <= sqr(tol)*norm2_b) {
 
100
          tol = sqrt(norm2_r/norm2_b);
 
101
          max_iter = n;
 
102
          return 0;     
 
103
        }
 
104
        x  -= rho*r;
 
105
    }
 
106
    tol = sqrt(norm2_r/norm2_b);
 
107
    return 1;
 
108
}
 
109
//>puzawa:
 
110
}// namespace rheolef
 
111
# endif // _SKIT_PUZAWA_H