action_lu_solve.hh
Go to the documentation of this file.
00001 //=====================================================
00002 // File   :  action_lu_solve.hh
00003 // Author :  L. Plagne <laurent.plagne@edf.fr)>        
00004 // Copyright (C) EDF R&D,  lun sep 30 14:23:19 CEST 2002
00005 //=====================================================
00006 // 
00007 // This program is free software; you can redistribute it and/or
00008 // modify it under the terms of the GNU General Public License
00009 // as published by the Free Software Foundation; either version 2
00010 // of the License, or (at your option) any later version.
00011 // 
00012 // This program is distributed in the hope that it will be useful,
00013 // but WITHOUT ANY WARRANTY; without even the implied warranty of
00014 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
00015 // GNU General Public License for more details.
00016 // You should have received a copy of the GNU General Public License
00017 // along with this program; if not, write to the Free Software
00018 // Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
00019 // 
00020 #ifndef ACTION_LU_SOLVE
00021 #define ACTION_LU_SOLVE
00022 #include "utilities.h"
00023 #include "STL_interface.hh"
00024 #include <string>
00025 #include "init/init_function.hh"
00026 #include "init/init_vector.hh"
00027 #include "init/init_matrix.hh"
00028 
00029 using namespace std;
00030 
00031 template<class Interface>
00032 class Action_lu_solve 
00033 {
00034 
00035 public :
00036 
00037   static inline std::string name( void )
00038   {
00039     return "lu_solve_"+Interface::name();
00040   }
00041   
00042   static double nb_op_base(int size){
00043     return 2.0*size*size*size/3.0;  // questionable but not really important
00044   }
00045 
00046 
00047   static double calculate( int nb_calc, int size ) {
00048 
00049     // STL matrix and vector initialization
00050     
00051     typename Interface::stl_matrix A_stl;
00052     typename Interface::stl_vector B_stl;
00053     typename Interface::stl_vector X_stl;
00054 
00055     init_matrix<pseudo_random>(A_stl,size);
00056     init_vector<pseudo_random>(B_stl,size);
00057     init_vector<null_function>(X_stl,size);
00058 
00059     // generic matrix and vector initialization
00060 
00061     typename Interface::gene_matrix A;
00062     typename Interface::gene_vector B;
00063     typename Interface::gene_vector X;
00064 
00065     typename Interface::gene_matrix LU; 
00066 
00067     Interface::matrix_from_stl(A,A_stl);
00068     Interface::vector_from_stl(B,B_stl);
00069     Interface::vector_from_stl(X,X_stl);
00070     Interface::matrix_from_stl(LU,A_stl);
00071   
00072     // local variable :
00073 
00074     typename Interface::Pivot_Vector pivot; // pivot vector
00075     Interface::new_Pivot_Vector(pivot,size);
00076     
00077     // timer utilities
00078 
00079     Portable_Timer chronos;
00080 
00081     // time measurement
00082 
00083     chronos.start();
00084     
00085     for (int ii=0;ii<nb_calc;ii++){
00086 
00087       // LU factorization
00088       Interface::copy_matrix(A,LU,size);
00089       Interface::LU_factor(LU,pivot,size);
00090       
00091       // LU solve
00092 
00093       Interface::LU_solve(LU,pivot,B,X,size);
00094 
00095     }
00096 
00097     // Time stop
00098 
00099     chronos.stop();
00100 
00101     double time=chronos.user_time();
00102   
00103     // check result :
00104 
00105     typename Interface::stl_vector B_new_stl(size);
00106     Interface::vector_to_stl(X,X_stl);
00107 
00108     STL_interface<typename Interface::real_type>::matrix_vector_product(A_stl,X_stl,B_new_stl,size); 
00109   
00110     typename Interface::real_type error=
00111       STL_interface<typename Interface::real_type>::norm_diff(B_stl,B_new_stl);
00112     
00113     if (error>1.e-5){
00114       INFOS("WRONG CALCULATION...residual=" << error);
00115       STL_interface<typename Interface::real_type>::display_vector(B_stl);
00116       STL_interface<typename Interface::real_type>::display_vector(B_new_stl);
00117       exit(0);
00118     }
00119     
00120     // deallocation and return time
00121     
00122     Interface::free_matrix(A,size);
00123     Interface::free_vector(B);
00124     Interface::free_vector(X);
00125     Interface::free_Pivot_Vector(pivot);
00126 
00127     return time;
00128   }
00129 
00130 };
00131   
00132 
00133 #endif
00134 
00135 
00136 


libicr
Author(s): Robert Krug
autogenerated on Mon Jan 6 2014 11:32:28