1//=====================================================
2// File   :  action_hessenberg.hh
3// Copyright (C) 2008 Gael Guennebaud <gael.guennebaud@inria.fr>
4//=====================================================
5//
6// This program is free software; you can redistribute it and/or
7// modify it under the terms of the GNU General Public License
8// as published by the Free Software Foundation; either version 2
9// of the License, or (at your option) any later version.
10//
11// This program 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// You should have received a copy of the GNU General Public License
16// along with this program; if not, write to the Free Software
17// Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
18//
19#ifndef ACTION_HESSENBERG
20#define ACTION_HESSENBERG
21#include "utilities.h"
22#include "STL_interface.hh"
23#include <string>
24#include "init/init_function.hh"
25#include "init/init_vector.hh"
26#include "init/init_matrix.hh"
27
28using namespace std;
29
30template<class Interface>
31class Action_hessenberg {
32
33public :
34
35  // Ctor
36
37  Action_hessenberg( int size ):_size(size)
38  {
39    MESSAGE("Action_hessenberg Ctor");
40
41    // STL vector initialization
42    init_matrix<pseudo_random>(X_stl,_size);
43
44    init_matrix<null_function>(C_stl,_size);
45    init_matrix<null_function>(resu_stl,_size);
46
47    // generic matrix and vector initialization
48    Interface::matrix_from_stl(X_ref,X_stl);
49    Interface::matrix_from_stl(X,X_stl);
50    Interface::matrix_from_stl(C,C_stl);
51
52    _cost = 0;
53    for (int j=0; j<_size-2; ++j)
54    {
55      double r = std::max(0,_size-j-1);
56      double b = std::max(0,_size-j-2);
57      _cost += 6 + 3*b + r*r*4 + r*_size*4;
58    }
59  }
60
61  // invalidate copy ctor
62
63  Action_hessenberg( const  Action_hessenberg & )
64  {
65    INFOS("illegal call to Action_hessenberg Copy Ctor");
66    exit(1);
67  }
68
69  // Dtor
70
71  ~Action_hessenberg( void ){
72
73    MESSAGE("Action_hessenberg Dtor");
74
75    // deallocation
76    Interface::free_matrix(X_ref,_size);
77    Interface::free_matrix(X,_size);
78    Interface::free_matrix(C,_size);
79  }
80
81  // action name
82
83  static inline std::string name( void )
84  {
85    return "hessenberg_"+Interface::name();
86  }
87
88  double nb_op_base( void ){
89    return _cost;
90  }
91
92  inline void initialize( void ){
93    Interface::copy_matrix(X_ref,X,_size);
94  }
95
96  inline void calculate( void ) {
97      Interface::hessenberg(X,C,_size);
98  }
99
100  void check_result( void ){
101    // calculation check
102    Interface::matrix_to_stl(C,resu_stl);
103
104//     STL_interface<typename Interface::real_type>::hessenberg(X_stl,C_stl,_size);
105//
106//     typename Interface::real_type error=
107//       STL_interface<typename Interface::real_type>::norm_diff(C_stl,resu_stl);
108//
109//     if (error>1.e-6){
110//       INFOS("WRONG CALCULATION...residual=" << error);
111//       exit(0);
112//     }
113
114  }
115
116private :
117
118  typename Interface::stl_matrix X_stl;
119  typename Interface::stl_matrix C_stl;
120  typename Interface::stl_matrix resu_stl;
121
122  typename Interface::gene_matrix X_ref;
123  typename Interface::gene_matrix X;
124  typename Interface::gene_matrix C;
125
126  int _size;
127  double _cost;
128};
129
130template<class Interface>
131class Action_tridiagonalization {
132
133public :
134
135  // Ctor
136
137  Action_tridiagonalization( int size ):_size(size)
138  {
139    MESSAGE("Action_tridiagonalization Ctor");
140
141    // STL vector initialization
142    init_matrix<pseudo_random>(X_stl,_size);
143
144    for(int i=0; i<_size; ++i)
145    {
146      for(int j=0; j<i; ++j)
147        X_stl[i][j] = X_stl[j][i];
148    }
149
150    init_matrix<null_function>(C_stl,_size);
151    init_matrix<null_function>(resu_stl,_size);
152
153    // generic matrix and vector initialization
154    Interface::matrix_from_stl(X_ref,X_stl);
155    Interface::matrix_from_stl(X,X_stl);
156    Interface::matrix_from_stl(C,C_stl);
157
158    _cost = 0;
159    for (int j=0; j<_size-2; ++j)
160    {
161      double r = std::max(0,_size-j-1);
162      double b = std::max(0,_size-j-2);
163      _cost += 6. + 3.*b + r*r*8.;
164    }
165  }
166
167  // invalidate copy ctor
168
169  Action_tridiagonalization( const  Action_tridiagonalization & )
170  {
171    INFOS("illegal call to Action_tridiagonalization Copy Ctor");
172    exit(1);
173  }
174
175  // Dtor
176
177  ~Action_tridiagonalization( void ){
178
179    MESSAGE("Action_tridiagonalization Dtor");
180
181    // deallocation
182    Interface::free_matrix(X_ref,_size);
183    Interface::free_matrix(X,_size);
184    Interface::free_matrix(C,_size);
185  }
186
187  // action name
188
189  static inline std::string name( void ) { return "tridiagonalization_"+Interface::name(); }
190
191  double nb_op_base( void ){
192    return _cost;
193  }
194
195  inline void initialize( void ){
196    Interface::copy_matrix(X_ref,X,_size);
197  }
198
199  inline void calculate( void ) {
200      Interface::tridiagonalization(X,C,_size);
201  }
202
203  void check_result( void ){
204    // calculation check
205    Interface::matrix_to_stl(C,resu_stl);
206
207//     STL_interface<typename Interface::real_type>::tridiagonalization(X_stl,C_stl,_size);
208//
209//     typename Interface::real_type error=
210//       STL_interface<typename Interface::real_type>::norm_diff(C_stl,resu_stl);
211//
212//     if (error>1.e-6){
213//       INFOS("WRONG CALCULATION...residual=" << error);
214//       exit(0);
215//     }
216
217  }
218
219private :
220
221  typename Interface::stl_matrix X_stl;
222  typename Interface::stl_matrix C_stl;
223  typename Interface::stl_matrix resu_stl;
224
225  typename Interface::gene_matrix X_ref;
226  typename Interface::gene_matrix X;
227  typename Interface::gene_matrix C;
228
229  int _size;
230  double _cost;
231};
232
233#endif
234