1 //=====================================================
2 // File   :  action_cholesky.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_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 
Action_cholesky(int size)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 
Action_cholesky(const Action_cholesky &)64   Action_cholesky( const  Action_cholesky & )
65   {
66     INFOS("illegal call to Action_cholesky Copy Ctor");
67     exit(1);
68   }
69 
70   // Dtor
71 
~Action_cholesky(void)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 
name(void)84   static inline std::string name( void )
85   {
86     return "cholesky_"+Interface::name();
87   }
88 
nb_op_base(void)89   double nb_op_base( void ){
90     return _cost;
91   }
92 
initialize(void)93   inline void initialize( void ){
94     Interface::copy_matrix(X_ref,X,_size);
95   }
96 
calculate(void)97   inline void calculate( void ) {
98       Interface::cholesky(X,C,_size);
99   }
100 
check_result(void)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