You can not select more than 25 topics Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.

125 lines
3.1 KiB

  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. using namespace std;
  28. template<class Interface>
  29. class Action_partial_lu {
  30. public :
  31. // Ctor
  32. Action_partial_lu( int size ):_size(size)
  33. {
  34. MESSAGE("Action_partial_lu Ctor");
  35. // STL vector initialization
  36. init_matrix<pseudo_random>(X_stl,_size);
  37. init_matrix<null_function>(C_stl,_size);
  38. // make sure X is invertible
  39. for (int i=0; i<_size; ++i)
  40. X_stl[i][i] = X_stl[i][i] * 1e2 + 1;
  41. // generic matrix and vector initialization
  42. Interface::matrix_from_stl(X_ref,X_stl);
  43. Interface::matrix_from_stl(X,X_stl);
  44. Interface::matrix_from_stl(C,C_stl);
  45. _cost = 2.0*size*size*size/3.0 + size*size;
  46. }
  47. // invalidate copy ctor
  48. Action_partial_lu( const Action_partial_lu & )
  49. {
  50. INFOS("illegal call to Action_partial_lu Copy Ctor");
  51. exit(1);
  52. }
  53. // Dtor
  54. ~Action_partial_lu( void ){
  55. MESSAGE("Action_partial_lu Dtor");
  56. // deallocation
  57. Interface::free_matrix(X_ref,_size);
  58. Interface::free_matrix(X,_size);
  59. Interface::free_matrix(C,_size);
  60. }
  61. // action name
  62. static inline std::string name( void )
  63. {
  64. return "partial_lu_decomp_"+Interface::name();
  65. }
  66. double nb_op_base( void ){
  67. return _cost;
  68. }
  69. inline void initialize( void ){
  70. Interface::copy_matrix(X_ref,X,_size);
  71. }
  72. inline void calculate( void ) {
  73. Interface::partial_lu_decomp(X,C,_size);
  74. }
  75. void check_result( void ){
  76. // calculation check
  77. // Interface::matrix_to_stl(C,resu_stl);
  78. // STL_interface<typename Interface::real_type>::lu_decomp(X_stl,C_stl,_size);
  79. //
  80. // typename Interface::real_type error=
  81. // STL_interface<typename Interface::real_type>::norm_diff(C_stl,resu_stl);
  82. //
  83. // if (error>1.e-6){
  84. // INFOS("WRONG CALCULATION...residual=" << error);
  85. // exit(0);
  86. // }
  87. }
  88. private :
  89. typename Interface::stl_matrix X_stl;
  90. typename Interface::stl_matrix C_stl;
  91. typename Interface::gene_matrix X_ref;
  92. typename Interface::gene_matrix X;
  93. typename Interface::gene_matrix C;
  94. int _size;
  95. double _cost;
  96. };
  97. #endif