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