← Back to branch summary

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

« back to all changes in this revision

Viewing changes to skit/lib/cg.h

  • Committer: Bazaar Package Importer
  • Author(s): Christophe Prud'homme
  • Date: 2010-06-12 09:08:59 UTC
  • Revision ID: james.westby@ubuntu.com-20100612090859-8gpm2gc7j3ab43et
Tags: upstream-5.89
Import upstream version 5.89

Show diffs side-by-side

added added

removed removed

Lines of Context:
skit/lib/cg.h
 
1
# ifndef _SKIT_CG_H
 
2
# define _SKIT_CG_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
/*D:cg
 
24
NAME: @code{cg} -- conjugate gradient algorithm. 
 
25
@findex cg
 
26
@cindex conjugate gradient algorithm
 
27
@cindex iterative solver
 
28
SYNOPSIS:
 
29
    @example
 
30
        template <class Matrix, class Vector, class Preconditioner, class Real>
 
31
        int cg (const Matrix &A, Vector &x, const Vector &b,
 
32
                const Preconditioner &M, int &max_iter, Real &tol);
 
33
    @end example
 
34
 
 
35
EXAMPLE:
 
36
  @noindent
 
37
  The simplest call to 'cg' has the folling form:
 
38
    @example
 
39
        int status = cg(a, x, b, EYE, 100, 1e-7);
 
40
    @end example
 
41
 
 
42
DESCRIPTION:       
 
43
  @noindent
 
44
  @code{cg} solves the symmetric positive definite linear
 
45
  system Ax=b using the Conjugate Gradient method.
 
46
 
 
47
  @noindent
 
48
  The return value indicates convergence within max_iter (input)
 
49
  iterations (0), or no convergence within max_iter iterations (1).
 
50
  Upon successful return, output arguments have the following values:
 
51
  @table @code
 
52
    @itemx x
 
53
        approximate solution to Ax = b
 
54
 
 
55
    @itemx max_iter
 
56
        the number of iterations performed before the tolerance was reached
 
57
 
 
58
    @itemx tol
 
59
        the residual after the final iteration
 
60
  @end table
 
61
 
 
62
NOTE: 
 
63
 
 
64
  @noindent
 
65
  @code{cg} is an iterative template routine.
 
66
 
 
67
  @noindent
 
68
  @code{cg} follows the algorithm described on p. 15 in
 
69
 
 
70
  @quotation
 
71
        Templates for the Solution of Linear Systems: Building Blocks for Iterative Methods, 
 
72
        2nd Edition, 
 
73
        R. Barrett, M. Berry, T. F. Chan, J. Demmel, J. Donato, J. Dongarra, V. Eijkhout,
 
74
        R. Pozo, C. Romine, H. Van der Vorst,
 
75
        SIAM, 1994, 
 
76
        @url{ftp.netlib.org/templates/templates.ps}.
 
77
  @end quotation
 
78
 
 
79
  @noindent
 
80
  The present implementation is inspired from 
 
81
  @code{IML++ 1.2} iterative method library,
 
82
  @url{http://math.nist.gov/iml++}.
 
83
 
 
84
AUTHOR: 
 
85
 
 
86
    Pierre Saramito
 
87
    | Pierre.Saramito@imag.fr
 
88
    LMC-IMAG, 38041 Grenoble cedex 9, France
 
89
 
 
90
DATE: 
 
91
    
 
92
    28 february 1997
 
93
 
 
94
METHODS: @cg
 
95
End:
 
96
*/
 
97
// ========== [ method implementation ] ====================================
 
98
 
 
99
template < class Matrix, class Vector, class Preconditioner, class Real>
 
100
int cg(const Matrix &A, Vector &x, const Vector &b,
 
101
       const Preconditioner &M, int &max_iter, Real &tol)
 
102
{
 
103
    Real resid;
 
104
    Vector p, q;
 
105
    Real  alpha, beta, rho, rho_1 = 0;
 
106
 
 
107
    Real normb = norm(b);
 
108
    Vector r = b - A*x;
 
109
 
 
110
    if (normb == Real(0)) 
 
111
        normb = 1;
 
112
  
 
113
    if ((resid = norm(r) / normb) <= tol) {
 
114
        tol = resid;
 
115
        max_iter = 0;
 
116
        return 0;
 
117
    }
 
118
    for (int i = 1; i <= max_iter; i++) {
 
119
    
 
120
        if (i == 1) {
 
121
            Vector z = M.solve(r);
 
122
            rho = dot(r, z);
 
123
            p = z;
 
124
        } else {
 
125
            Vector z = M.solve(r);
 
126
            rho = dot(r, z);
 
127
            beta = rho / rho_1;
 
128
            p = z + beta * p;
 
129
        }
 
130
        q = A*p;
 
131
        alpha = rho / dot(p, q);
 
132
    
 
133
        x += alpha * p;
 
134
        r -= alpha * q;
 
135
    
 
136
        if ((resid = norm(r) / normb) <= tol) {
 
137
          tol = resid;
 
138
          max_iter = i;
 
139
          return 0;     
 
140
        }
 
141
        rho_1 = rho;
 
142
    }
 
143
    tol = resid;
 
144
    return 1;
 
145
}
 
146
# endif // _SKIT_CG_H