1*bf2c3715SXin Li //===================================================== 2*bf2c3715SXin Li // File : action_trisolve.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_TRISOLVE 20*bf2c3715SXin Li #define ACTION_TRISOLVE 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_trisolve { 32*bf2c3715SXin Li 33*bf2c3715SXin Li public : 34*bf2c3715SXin Li 35*bf2c3715SXin Li // Ctor 36*bf2c3715SXin Li Action_trisolve(int size)37*bf2c3715SXin Li Action_trisolve( int size ):_size(size) 38*bf2c3715SXin Li { 39*bf2c3715SXin Li MESSAGE("Action_trisolve Ctor"); 40*bf2c3715SXin Li 41*bf2c3715SXin Li // STL vector initialization 42*bf2c3715SXin Li init_matrix<pseudo_random>(L_stl,_size); 43*bf2c3715SXin Li init_vector<pseudo_random>(B_stl,_size); 44*bf2c3715SXin Li init_vector<null_function>(X_stl,_size); 45*bf2c3715SXin Li for (int j=0; j<_size; ++j) 46*bf2c3715SXin Li { 47*bf2c3715SXin Li for (int i=0; i<j; ++i) 48*bf2c3715SXin Li L_stl[j][i] = 0; 49*bf2c3715SXin Li L_stl[j][j] += 3; 50*bf2c3715SXin Li } 51*bf2c3715SXin Li 52*bf2c3715SXin Li init_vector<null_function>(resu_stl,_size); 53*bf2c3715SXin Li 54*bf2c3715SXin Li // generic matrix and vector initialization 55*bf2c3715SXin Li Interface::matrix_from_stl(L,L_stl); 56*bf2c3715SXin Li Interface::vector_from_stl(X,X_stl); 57*bf2c3715SXin Li Interface::vector_from_stl(B,B_stl); 58*bf2c3715SXin Li 59*bf2c3715SXin Li _cost = 0; 60*bf2c3715SXin Li for (int j=0; j<_size; ++j) 61*bf2c3715SXin Li { 62*bf2c3715SXin Li _cost += 2*j + 1; 63*bf2c3715SXin Li } 64*bf2c3715SXin Li } 65*bf2c3715SXin Li 66*bf2c3715SXin Li // invalidate copy ctor 67*bf2c3715SXin Li Action_trisolve(const Action_trisolve &)68*bf2c3715SXin Li Action_trisolve( const Action_trisolve & ) 69*bf2c3715SXin Li { 70*bf2c3715SXin Li INFOS("illegal call to Action_trisolve Copy Ctor"); 71*bf2c3715SXin Li exit(1); 72*bf2c3715SXin Li } 73*bf2c3715SXin Li 74*bf2c3715SXin Li // Dtor 75*bf2c3715SXin Li ~Action_trisolve(void)76*bf2c3715SXin Li ~Action_trisolve( void ){ 77*bf2c3715SXin Li 78*bf2c3715SXin Li MESSAGE("Action_trisolve Dtor"); 79*bf2c3715SXin Li 80*bf2c3715SXin Li // deallocation 81*bf2c3715SXin Li Interface::free_matrix(L,_size); 82*bf2c3715SXin Li Interface::free_vector(B); 83*bf2c3715SXin Li Interface::free_vector(X); 84*bf2c3715SXin Li } 85*bf2c3715SXin Li 86*bf2c3715SXin Li // action name 87*bf2c3715SXin Li name(void)88*bf2c3715SXin Li static inline std::string name( void ) 89*bf2c3715SXin Li { 90*bf2c3715SXin Li return "trisolve_vector_"+Interface::name(); 91*bf2c3715SXin Li } 92*bf2c3715SXin Li nb_op_base(void)93*bf2c3715SXin Li double nb_op_base( void ){ 94*bf2c3715SXin Li return _cost; 95*bf2c3715SXin Li } 96*bf2c3715SXin Li initialize(void)97*bf2c3715SXin Li inline void initialize( void ){ 98*bf2c3715SXin Li //Interface::copy_vector(X_ref,X,_size); 99*bf2c3715SXin Li } 100*bf2c3715SXin Li calculate(void)101*bf2c3715SXin Li inline void calculate( void ) { 102*bf2c3715SXin Li Interface::trisolve_lower(L,B,X,_size); 103*bf2c3715SXin Li } 104*bf2c3715SXin Li check_result()105*bf2c3715SXin Li void check_result(){ 106*bf2c3715SXin Li if (_size>128) return; 107*bf2c3715SXin Li // calculation check 108*bf2c3715SXin Li Interface::vector_to_stl(X,resu_stl); 109*bf2c3715SXin Li 110*bf2c3715SXin Li STL_interface<typename Interface::real_type>::trisolve_lower(L_stl,B_stl,X_stl,_size); 111*bf2c3715SXin Li 112*bf2c3715SXin Li typename Interface::real_type error= 113*bf2c3715SXin Li STL_interface<typename Interface::real_type>::norm_diff(X_stl,resu_stl); 114*bf2c3715SXin Li 115*bf2c3715SXin Li if (error>1.e-4){ 116*bf2c3715SXin Li INFOS("WRONG CALCULATION...residual=" << error); 117*bf2c3715SXin Li exit(2); 118*bf2c3715SXin Li } //else INFOS("CALCULATION OK...residual=" << error); 119*bf2c3715SXin Li 120*bf2c3715SXin Li } 121*bf2c3715SXin Li 122*bf2c3715SXin Li private : 123*bf2c3715SXin Li 124*bf2c3715SXin Li typename Interface::stl_matrix L_stl; 125*bf2c3715SXin Li typename Interface::stl_vector X_stl; 126*bf2c3715SXin Li typename Interface::stl_vector B_stl; 127*bf2c3715SXin Li typename Interface::stl_vector resu_stl; 128*bf2c3715SXin Li 129*bf2c3715SXin Li typename Interface::gene_matrix L; 130*bf2c3715SXin Li typename Interface::gene_vector X; 131*bf2c3715SXin Li typename Interface::gene_vector B; 132*bf2c3715SXin Li 133*bf2c3715SXin Li int _size; 134*bf2c3715SXin Li double _cost; 135*bf2c3715SXin Li }; 136*bf2c3715SXin Li 137*bf2c3715SXin Li #endif 138