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