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 
28 using namespace std;
29 
30 template<class Interface>
31 class Action_hessenberg {
32 
33 public :
34 
35   // Ctor
36 
Action_hessenberg(int size)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 
Action_hessenberg(const Action_hessenberg &)63   Action_hessenberg( const  Action_hessenberg & )
64   {
65     INFOS("illegal call to Action_hessenberg Copy Ctor");
66     exit(1);
67   }
68 
69   // Dtor
70 
~Action_hessenberg(void)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 
name(void)83   static inline std::string name( void )
84   {
85     return "hessenberg_"+Interface::name();
86   }
87 
nb_op_base(void)88   double nb_op_base( void ){
89     return _cost;
90   }
91 
initialize(void)92   inline void initialize( void ){
93     Interface::copy_matrix(X_ref,X,_size);
94   }
95 
calculate(void)96   inline void calculate( void ) {
97       Interface::hessenberg(X,C,_size);
98   }
99 
check_result(void)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 
116 private :
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 
130 template<class Interface>
131 class Action_tridiagonalization {
132 
133 public :
134 
135   // Ctor
136 
Action_tridiagonalization(int size)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 
Action_tridiagonalization(const Action_tridiagonalization &)169   Action_tridiagonalization( const  Action_tridiagonalization & )
170   {
171     INFOS("illegal call to Action_tridiagonalization Copy Ctor");
172     exit(1);
173   }
174 
175   // Dtor
176 
~Action_tridiagonalization(void)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 
name(void)189   static inline std::string name( void ) { return "tridiagonalization_"+Interface::name(); }
190 
nb_op_base(void)191   double nb_op_base( void ){
192     return _cost;
193   }
194 
initialize(void)195   inline void initialize( void ){
196     Interface::copy_matrix(X_ref,X,_size);
197   }
198 
calculate(void)199   inline void calculate( void ) {
200       Interface::tridiagonalization(X,C,_size);
201   }
202 
check_result(void)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 
219 private :
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