← Back to branch summary

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

« back to all changes in this revision

Viewing changes to skit/tst/solve_iterative_tst.cc

  • 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/tst/solve_iterative_tst.cc
 
1
///
 
2
/// This file is part of Rheolef.
 
3
///
 
4
/// Copyright (C) 2000-2009 Pierre Saramito <Pierre.Saramito@imag.fr>
 
5
///
 
6
/// Rheolef is free software; you can redistribute it and/or modify
 
7
/// it under the terms of the GNU General Public License as published by
 
8
/// the Free Software Foundation; either version 2 of the License, or
 
9
/// (at your option) any later version.
 
10
///
 
11
/// Rheolef is distributed in the hope that it will be useful,
 
12
/// but WITHOUT ANY WARRANTY; without even the implied warranty of
 
13
/// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 
14
/// GNU General Public License for more details.
 
15
///
 
16
/// You should have received a copy of the GNU General Public License
 
17
/// along with Rheolef; if not, write to the Free Software
 
18
/// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
 
19
///
 
20
/// =========================================================================
 
21
//
 
22
// This code presents various solvers
 
23
//
 
24
// also print matlab code for non-regression test purpose, i.e.
 
25
//
 
26
// usage:
 
27
//    makefish 4 | solve_iterative_tst | octave -q
 
28
//  machine-precision independent computation:
 
29
//    makefish 4 | solve_iterative_tst -ndigits 12 | octave -q
 
30
//
 
31
# include "rheolef/skit.h"
 
32
using namespace std;
 
33
 
 
34
extern int solve_bicg     (const csr<Float>&, const vec<Float>&);
 
35
extern int solve_bicgstab (const csr<Float>&, const vec<Float>&);
 
36
extern int solve_cg       (const csr<Float>&, const vec<Float>&);
 
37
extern int solve_cgs      (const csr<Float>&, const vec<Float>&);
 
38
extern int solve_cheby    (const csr<Float>&, const vec<Float>&);
 
39
extern int solve_gmres    (const csr<Float>&, const vec<Float>&);
 
40
extern int solve_ir       (const csr<Float>&, const vec<Float>&);
 
41
extern int solve_qmr      (const csr<Float>&, const vec<Float>&);
 
42
 
 
43
bool do_trunc;
 
44
 
 
45
Float eps1 = 1e-7;
 
46
 
 
47
template <class T>
 
48
inline T my_round (const T& x) { return eps1*T(int(x/eps1+0.5)); }
 
49
 
 
50
template <class Iterator>
 
51
void
 
52
round_value (Iterator first, Iterator last)
 
53
{
 
54
    while (first != last) {
 
55
        *first = my_round(*first);
 
56
        first++;
 
57
    }
 
58
}
 
59
void check (const char* method, int max_iter, const csr<Float>& a, 
 
60
            const vec<Float>& x, const vec<Float>& b)
 
61
{
 
62
    if (do_trunc) {
 
63
        vec<Float> x1 = x;
 
64
        round_value(x1.begin(), x1.end());
 
65
        cout << "x=" <<  x1 << ";"    << endl;
 
66
    } else {
 
67
        cout << "x=" <<  x << ";"    << endl;
 
68
    }
 
69
    cout << "error_" << method << "_" << max_iter;
 
70
    if (!do_trunc) {
 
71
         cout << "=norm(a*x-b)\n\n";
 
72
    } else {
 
73
         cout << "=eps1*round(norm(a*x-b)/eps1)\n\n";
 
74
    }
 
75
}
 
76
int main (int argc, char** argv)
 
77
{
 
78
    do_trunc = false;
 
79
    int digits10 = numeric_limits<Float>::digits10;
 
80
    // cout.setf(ios::scientific, !(ios::fixed | ios::floatfield) );
 
81
 
 
82
    if (argc == 3) {
 
83
        // non-regression mode:
 
84
        // avoid machine-dependent output digits in non-regression mode
 
85
 
 
86
        // avoid variation of output formats 
 
87
        // when double, doubledouble or bigfloat
 
88
        numeric_flags<Float>::output_as_double(true);
 
89
 
 
90
        do_trunc = true;
 
91
        // cout.setf(ios::scientific);
 
92
        cout.setf(ios::fixed);
 
93
        digits10 = atoi(argv[2]);
 
94
        cout << "eps1=sqrt(10^(-" << digits10 << "));\n";
 
95
        eps1=pow(Float(10), Float(-digits10/2));
 
96
 
 
97
    }
 
98
    cout << setprecision(digits10);
 
99
 
 
100
    csr<Float> a;
 
101
    cin >> a;
 
102
    vec<Float> x(a.ncol());
 
103
    x = 1;
 
104
    vec<Float> b = a*x;
 
105
 
 
106
    cout << ml << "a=" << a << ";" << endl;
 
107
    cout << "b=" << b << ";" << endl;
 
108
 
 
109
    solve_bicg(a, b);
 
110
    solve_bicgstab(a, b);
 
111
    solve_cg(a, b);
 
112
    solve_cgs(a, b);
 
113
    solve_cheby(a, b);
 
114
    solve_gmres(a, b);
 
115
    solve_qmr(a, b);
 
116
    
 
117
    return 0;
 
118
}
 
119