Home | History | Annotate | Download | only in actions
      1 //=====================================================
      2 // File   :  action_matrix_matrix_product.hh
      3 // Author :  L. Plagne <laurent.plagne (at) 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_TRISOLVE_MATRIX_PRODUCT
     21 #define ACTION_TRISOLVE_MATRIX_PRODUCT
     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_trisolve_matrix {
     33 
     34 public :
     35 
     36   // Ctor
     37 
     38   Action_trisolve_matrix( int size ):_size(size)
     39   {
     40     MESSAGE("Action_trisolve_matrix Ctor");
     41 
     42     // STL matrix and vector initialization
     43 
     44     init_matrix<pseudo_random>(A_stl,_size);
     45     init_matrix<pseudo_random>(B_stl,_size);
     46     init_matrix<null_function>(X_stl,_size);
     47     init_matrix<null_function>(resu_stl,_size);
     48 
     49     for (int j=0; j<_size; ++j)
     50     {
     51       for (int i=0; i<j; ++i)
     52         A_stl[j][i] = 0;
     53       A_stl[j][j] += 3;
     54     }
     55 
     56     // generic matrix and vector initialization
     57 
     58     Interface::matrix_from_stl(A_ref,A_stl);
     59     Interface::matrix_from_stl(B_ref,B_stl);
     60     Interface::matrix_from_stl(X_ref,X_stl);
     61 
     62     Interface::matrix_from_stl(A,A_stl);
     63     Interface::matrix_from_stl(B,B_stl);
     64     Interface::matrix_from_stl(X,X_stl);
     65 
     66     _cost = 0;
     67     for (int j=0; j<_size; ++j)
     68     {
     69       _cost += 2*j + 1;
     70     }
     71     _cost *= _size;
     72   }
     73 
     74   // invalidate copy ctor
     75 
     76   Action_trisolve_matrix( const  Action_trisolve_matrix & )
     77   {
     78     INFOS("illegal call to Action_trisolve_matrix Copy Ctor");
     79     exit(0);
     80   }
     81 
     82   // Dtor
     83 
     84   ~Action_trisolve_matrix( void ){
     85 
     86     MESSAGE("Action_trisolve_matrix Dtor");
     87 
     88     // deallocation
     89 
     90     Interface::free_matrix(A,_size);
     91     Interface::free_matrix(B,_size);
     92     Interface::free_matrix(X,_size);
     93 
     94     Interface::free_matrix(A_ref,_size);
     95     Interface::free_matrix(B_ref,_size);
     96     Interface::free_matrix(X_ref,_size);
     97 
     98   }
     99 
    100   // action name
    101 
    102   static inline std::string name( void )
    103   {
    104     return "trisolve_matrix_"+Interface::name();
    105   }
    106 
    107   double nb_op_base( void ){
    108     return _cost;
    109   }
    110 
    111   inline void initialize( void ){
    112 
    113     Interface::copy_matrix(A_ref,A,_size);
    114     Interface::copy_matrix(B_ref,B,_size);
    115     Interface::copy_matrix(X_ref,X,_size);
    116 
    117   }
    118 
    119   inline void calculate( void ) {
    120       Interface::trisolve_lower_matrix(A,B,X,_size);
    121   }
    122 
    123   void check_result( void ){
    124 
    125     // calculation check
    126 
    127 //     Interface::matrix_to_stl(X,resu_stl);
    128 //
    129 //     STL_interface<typename Interface::real_type>::matrix_matrix_product(A_stl,B_stl,X_stl,_size);
    130 //
    131 //     typename Interface::real_type error=
    132 //       STL_interface<typename Interface::real_type>::norm_diff(X_stl,resu_stl);
    133 //
    134 //     if (error>1.e-6){
    135 //       INFOS("WRONG CALCULATION...residual=" << error);
    136 // //       exit(1);
    137 //     }
    138 
    139   }
    140 
    141 private :
    142 
    143   typename Interface::stl_matrix A_stl;
    144   typename Interface::stl_matrix B_stl;
    145   typename Interface::stl_matrix X_stl;
    146   typename Interface::stl_matrix resu_stl;
    147 
    148   typename Interface::gene_matrix A_ref;
    149   typename Interface::gene_matrix B_ref;
    150   typename Interface::gene_matrix X_ref;
    151 
    152   typename Interface::gene_matrix A;
    153   typename Interface::gene_matrix B;
    154   typename Interface::gene_matrix X;
    155 
    156   int _size;
    157   double _cost;
    158 
    159 };
    160 
    161 
    162 #endif
    163 
    164 
    165 
    166