1*bf2c3715SXin Li //===================================================== 2*bf2c3715SXin Li // File : action_cholesky.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_CHOLESKY 20*bf2c3715SXin Li #define ACTION_CHOLESKY 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_cholesky { 32*bf2c3715SXin Li 33*bf2c3715SXin Li public : 34*bf2c3715SXin Li 35*bf2c3715SXin Li // Ctor 36*bf2c3715SXin Li Action_cholesky(int size)37*bf2c3715SXin Li Action_cholesky( int size ):_size(size) 38*bf2c3715SXin Li { 39*bf2c3715SXin Li MESSAGE("Action_cholesky Ctor"); 40*bf2c3715SXin Li 41*bf2c3715SXin Li // STL mat/vec initialization 42*bf2c3715SXin Li init_matrix_symm<pseudo_random>(X_stl,_size); 43*bf2c3715SXin Li init_matrix<null_function>(C_stl,_size); 44*bf2c3715SXin Li 45*bf2c3715SXin Li // make sure X is invertible 46*bf2c3715SXin Li for (int i=0; i<_size; ++i) 47*bf2c3715SXin Li X_stl[i][i] = std::abs(X_stl[i][i]) * 1e2 + 100; 48*bf2c3715SXin Li 49*bf2c3715SXin Li // generic matrix and vector initialization 50*bf2c3715SXin Li Interface::matrix_from_stl(X_ref,X_stl); 51*bf2c3715SXin Li Interface::matrix_from_stl(X,X_stl); 52*bf2c3715SXin Li Interface::matrix_from_stl(C,C_stl); 53*bf2c3715SXin Li 54*bf2c3715SXin Li _cost = 0; 55*bf2c3715SXin Li for (int j=0; j<_size; ++j) 56*bf2c3715SXin Li { 57*bf2c3715SXin Li double r = std::max(_size - j -1,0); 58*bf2c3715SXin Li _cost += 2*(r*j+r+j); 59*bf2c3715SXin Li } 60*bf2c3715SXin Li } 61*bf2c3715SXin Li 62*bf2c3715SXin Li // invalidate copy ctor 63*bf2c3715SXin Li Action_cholesky(const Action_cholesky &)64*bf2c3715SXin Li Action_cholesky( const Action_cholesky & ) 65*bf2c3715SXin Li { 66*bf2c3715SXin Li INFOS("illegal call to Action_cholesky Copy Ctor"); 67*bf2c3715SXin Li exit(1); 68*bf2c3715SXin Li } 69*bf2c3715SXin Li 70*bf2c3715SXin Li // Dtor 71*bf2c3715SXin Li ~Action_cholesky(void)72*bf2c3715SXin Li ~Action_cholesky( void ){ 73*bf2c3715SXin Li 74*bf2c3715SXin Li MESSAGE("Action_cholesky Dtor"); 75*bf2c3715SXin Li 76*bf2c3715SXin Li // deallocation 77*bf2c3715SXin Li Interface::free_matrix(X_ref,_size); 78*bf2c3715SXin Li Interface::free_matrix(X,_size); 79*bf2c3715SXin Li Interface::free_matrix(C,_size); 80*bf2c3715SXin Li } 81*bf2c3715SXin Li 82*bf2c3715SXin Li // action name 83*bf2c3715SXin Li name(void)84*bf2c3715SXin Li static inline std::string name( void ) 85*bf2c3715SXin Li { 86*bf2c3715SXin Li return "cholesky_"+Interface::name(); 87*bf2c3715SXin Li } 88*bf2c3715SXin Li nb_op_base(void)89*bf2c3715SXin Li double nb_op_base( void ){ 90*bf2c3715SXin Li return _cost; 91*bf2c3715SXin Li } 92*bf2c3715SXin Li initialize(void)93*bf2c3715SXin Li inline void initialize( void ){ 94*bf2c3715SXin Li Interface::copy_matrix(X_ref,X,_size); 95*bf2c3715SXin Li } 96*bf2c3715SXin Li calculate(void)97*bf2c3715SXin Li inline void calculate( void ) { 98*bf2c3715SXin Li Interface::cholesky(X,C,_size); 99*bf2c3715SXin Li } 100*bf2c3715SXin Li check_result(void)101*bf2c3715SXin Li void check_result( void ){ 102*bf2c3715SXin Li // calculation check 103*bf2c3715SXin Li // STL_interface<typename Interface::real_type>::cholesky(X_stl,C_stl,_size); 104*bf2c3715SXin Li // 105*bf2c3715SXin Li // typename Interface::real_type error= 106*bf2c3715SXin Li // STL_interface<typename Interface::real_type>::norm_diff(C_stl,resu_stl); 107*bf2c3715SXin Li // 108*bf2c3715SXin Li // if (error>1.e-6){ 109*bf2c3715SXin Li // INFOS("WRONG CALCULATION...residual=" << error); 110*bf2c3715SXin Li // exit(0); 111*bf2c3715SXin Li // } 112*bf2c3715SXin Li 113*bf2c3715SXin Li } 114*bf2c3715SXin Li 115*bf2c3715SXin Li private : 116*bf2c3715SXin Li 117*bf2c3715SXin Li typename Interface::stl_matrix X_stl; 118*bf2c3715SXin Li typename Interface::stl_matrix C_stl; 119*bf2c3715SXin Li 120*bf2c3715SXin Li typename Interface::gene_matrix X_ref; 121*bf2c3715SXin Li typename Interface::gene_matrix X; 122*bf2c3715SXin Li typename Interface::gene_matrix C; 123*bf2c3715SXin Li 124*bf2c3715SXin Li int _size; 125*bf2c3715SXin Li double _cost; 126*bf2c3715SXin Li }; 127*bf2c3715SXin Li 128*bf2c3715SXin Li #endif 129