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.

128 lines
3.1 KiB

  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. using namespace std;
  28. template<class Interface>
  29. class Action_cholesky {
  30. public :
  31. // Ctor
  32. Action_cholesky( int size ):_size(size)
  33. {
  34. MESSAGE("Action_cholesky Ctor");
  35. // STL mat/vec initialization
  36. init_matrix_symm<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] = std::abs(X_stl[i][i]) * 1e2 + 100;
  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 = 0;
  46. for (int j=0; j<_size; ++j)
  47. {
  48. double r = std::max(_size - j -1,0);
  49. _cost += 2*(r*j+r+j);
  50. }
  51. }
  52. // invalidate copy ctor
  53. Action_cholesky( const Action_cholesky & )
  54. {
  55. INFOS("illegal call to Action_cholesky Copy Ctor");
  56. exit(1);
  57. }
  58. // Dtor
  59. ~Action_cholesky( void ){
  60. MESSAGE("Action_cholesky Dtor");
  61. // deallocation
  62. Interface::free_matrix(X_ref,_size);
  63. Interface::free_matrix(X,_size);
  64. Interface::free_matrix(C,_size);
  65. }
  66. // action name
  67. static inline std::string name( void )
  68. {
  69. return "cholesky_"+Interface::name();
  70. }
  71. double nb_op_base( void ){
  72. return _cost;
  73. }
  74. inline void initialize( void ){
  75. Interface::copy_matrix(X_ref,X,_size);
  76. }
  77. inline void calculate( void ) {
  78. Interface::cholesky(X,C,_size);
  79. }
  80. void check_result( void ){
  81. // calculation check
  82. // STL_interface<typename Interface::real_type>::cholesky(X_stl,C_stl,_size);
  83. //
  84. // typename Interface::real_type error=
  85. // STL_interface<typename Interface::real_type>::norm_diff(C_stl,resu_stl);
  86. //
  87. // if (error>1.e-6){
  88. // INFOS("WRONG CALCULATION...residual=" << error);
  89. // exit(0);
  90. // }
  91. }
  92. private :
  93. typename Interface::stl_matrix X_stl;
  94. typename Interface::stl_matrix C_stl;
  95. typename Interface::gene_matrix X_ref;
  96. typename Interface::gene_matrix X;
  97. typename Interface::gene_matrix C;
  98. int _size;
  99. double _cost;
  100. };
  101. #endif