1 //=====================================================
2 // File   :  action_lu_decomp.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_PARTIAL_LU
20 #define ACTION_PARTIAL_LU
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_partial_lu {
32 
33 public :
34 
35   // Ctor
36 
Action_partial_lu(int size)37   Action_partial_lu( int size ):_size(size)
38   {
39     MESSAGE("Action_partial_lu Ctor");
40 
41     // STL vector initialization
42     init_matrix<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] = X_stl[i][i] * 1e2 + 1;
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 = 2.0*size*size*size/3.0 + size*size;
55   }
56 
57   // invalidate copy ctor
58 
Action_partial_lu(const Action_partial_lu &)59   Action_partial_lu( const  Action_partial_lu & )
60   {
61     INFOS("illegal call to Action_partial_lu Copy Ctor");
62     exit(1);
63   }
64 
65   // Dtor
66 
~Action_partial_lu(void)67   ~Action_partial_lu( void ){
68 
69     MESSAGE("Action_partial_lu Dtor");
70 
71     // deallocation
72     Interface::free_matrix(X_ref,_size);
73     Interface::free_matrix(X,_size);
74     Interface::free_matrix(C,_size);
75   }
76 
77   // action name
78 
name(void)79   static inline std::string name( void )
80   {
81     return "partial_lu_decomp_"+Interface::name();
82   }
83 
nb_op_base(void)84   double nb_op_base( void ){
85     return _cost;
86   }
87 
initialize(void)88   inline void initialize( void ){
89     Interface::copy_matrix(X_ref,X,_size);
90   }
91 
calculate(void)92   inline void calculate( void ) {
93       Interface::partial_lu_decomp(X,C,_size);
94   }
95 
check_result(void)96   void check_result( void ){
97     // calculation check
98 //     Interface::matrix_to_stl(C,resu_stl);
99 
100 //     STL_interface<typename Interface::real_type>::lu_decomp(X_stl,C_stl,_size);
101 //
102 //     typename Interface::real_type error=
103 //       STL_interface<typename Interface::real_type>::norm_diff(C_stl,resu_stl);
104 //
105 //     if (error>1.e-6){
106 //       INFOS("WRONG CALCULATION...residual=" << error);
107 //       exit(0);
108 //     }
109 
110   }
111 
112 private :
113 
114   typename Interface::stl_matrix X_stl;
115   typename Interface::stl_matrix C_stl;
116 
117   typename Interface::gene_matrix X_ref;
118   typename Interface::gene_matrix X;
119   typename Interface::gene_matrix C;
120 
121   int _size;
122   double _cost;
123 };
124 
125 #endif
126