1 //=====================================================
2 // File   :  action_lu_solve.hh
3 // Author :  L. Plagne <laurent.plagne@edf.fr)>
4 // Copyright (C) EDF R&D,  lun sep 30 14:23:19 CEST 2002
5 //=====================================================
6 //
7 // This program is free software; you can redistribute it and/or
8 // modify it under the terms of the GNU General Public License
9 // as published by the Free Software Foundation; either version 2
10 // of the License, or (at your option) any later version.
11 //
12 // This program is distributed in the hope that it will be useful,
13 // but WITHOUT ANY WARRANTY; without even the implied warranty of
14 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
15 // GNU General Public License for more details.
16 // You should have received a copy of the GNU General Public License
17 // along with this program; if not, write to the Free Software
18 // Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
19 //
20 #ifndef ACTION_LU_SOLVE
21 #define ACTION_LU_SOLVE
22 #include "utilities.h"
23 #include "STL_interface.hh"
24 #include <string>
25 #include "init/init_function.hh"
26 #include "init/init_vector.hh"
27 #include "init/init_matrix.hh"
28 
29 using namespace std;
30 
31 template<class Interface>
32 class Action_lu_solve
33 {
34 
35 public :
36 
name(void)37   static inline std::string name( void )
38   {
39     return "lu_solve_"+Interface::name();
40   }
41 
nb_op_base(int size)42   static double nb_op_base(int size){
43     return 2.0*size*size*size/3.0;  // questionable but not really important
44   }
45 
46 
calculate(int nb_calc,int size)47   static double calculate( int nb_calc, int size ) {
48 
49     // STL matrix and vector initialization
50 
51     typename Interface::stl_matrix A_stl;
52     typename Interface::stl_vector B_stl;
53     typename Interface::stl_vector X_stl;
54 
55     init_matrix<pseudo_random>(A_stl,size);
56     init_vector<pseudo_random>(B_stl,size);
57     init_vector<null_function>(X_stl,size);
58 
59     // generic matrix and vector initialization
60 
61     typename Interface::gene_matrix A;
62     typename Interface::gene_vector B;
63     typename Interface::gene_vector X;
64 
65     typename Interface::gene_matrix LU;
66 
67     Interface::matrix_from_stl(A,A_stl);
68     Interface::vector_from_stl(B,B_stl);
69     Interface::vector_from_stl(X,X_stl);
70     Interface::matrix_from_stl(LU,A_stl);
71 
72     // local variable :
73 
74     typename Interface::Pivot_Vector pivot; // pivot vector
75     Interface::new_Pivot_Vector(pivot,size);
76 
77     // timer utilities
78 
79     Portable_Timer chronos;
80 
81     // time measurement
82 
83     chronos.start();
84 
85     for (int ii=0;ii<nb_calc;ii++){
86 
87       // LU factorization
88       Interface::copy_matrix(A,LU,size);
89       Interface::LU_factor(LU,pivot,size);
90 
91       // LU solve
92 
93       Interface::LU_solve(LU,pivot,B,X,size);
94 
95     }
96 
97     // Time stop
98 
99     chronos.stop();
100 
101     double time=chronos.user_time();
102 
103     // check result :
104 
105     typename Interface::stl_vector B_new_stl(size);
106     Interface::vector_to_stl(X,X_stl);
107 
108     STL_interface<typename Interface::real_type>::matrix_vector_product(A_stl,X_stl,B_new_stl,size);
109 
110     typename Interface::real_type error=
111       STL_interface<typename Interface::real_type>::norm_diff(B_stl,B_new_stl);
112 
113     if (error>1.e-5){
114       INFOS("WRONG CALCULATION...residual=" << error);
115       STL_interface<typename Interface::real_type>::display_vector(B_stl);
116       STL_interface<typename Interface::real_type>::display_vector(B_new_stl);
117       exit(0);
118     }
119 
120     // deallocation and return time
121 
122     Interface::free_matrix(A,size);
123     Interface::free_vector(B);
124     Interface::free_vector(X);
125     Interface::free_Pivot_Vector(pivot);
126 
127     return time;
128   }
129 
130 };
131 
132 
133 #endif
134 
135 
136 
137