src/2geom/numeric/linear_system.h

   1 #ifndef _NL_LINEAR_SYSTEM_H_\r
   2 #define _NL_LINEAR_SYSTEM_H_\r
   3 \r
   4 \r
   5 #include <cassert>\r
   6 \r
   7 #include <gsl/gsl_linalg.h>\r
   8 \r
   9 #include "numeric/matrix.h"\r
  10 #include "numeric/vector.h"\r
  11 \r
  12 \r
  13 namespace Geom { namespace NL {\r
  14 \r
  15 \r
  16 class LinearSystem\r
  17 {\r
  18 public:\r
  19         LinearSystem(Matrix & _matrix, Vector & _vector)\r
  20                 : m_matrix(_matrix), m_vector(_vector), m_solution(_matrix.columns())\r
  21         {\r
  22         }\r
  23         \r
  24         const Vector & LU_solve()\r
  25         {\r
  26                 assert( matrix().rows() == matrix().columns() \r
  27                                 && matrix().rows() == vector().size() );\r
  28                 int s;\r
  29                 gsl_permutation * p = gsl_permutation_alloc(matrix().rows());\r
  30                 gsl_linalg_LU_decomp (matrix().get_gsl_matrix(), p, &s);\r
  31                 gsl_linalg_LU_solve( matrix().get_gsl_matrix(), \r
  32                                                          p, \r
  33                                              vector().get_gsl_vector(), \r
  34                                              m_solution.get_gsl_vector()\r
  35                                            );\r
  36                 gsl_permutation_free(p);\r
  37                 return solution();\r
  38         }\r
  39         \r
  40         const Vector & SV_solve()\r
  41         {\r
  42                 assert( matrix().rows() >= matrix().columns()\r
  43                                 && matrix().rows() == vector().size() );\r
  44                 \r
  45                 gsl_matrix* U = matrix().get_gsl_matrix();\r
  46                 gsl_matrix* V = gsl_matrix_alloc(matrix().columns(), matrix().columns());\r
  47                 gsl_vector* S = gsl_vector_alloc(matrix().columns());\r
  48                 gsl_vector* work = gsl_vector_alloc(matrix().columns());\r
  49                 \r
  50                 gsl_linalg_SV_decomp( U, V, S, work );\r
  51                 \r
  52                 gsl_vector* b = vector().get_gsl_vector();\r
  53                 gsl_vector* x = m_solution.get_gsl_vector();\r
  54                 \r
  55                 gsl_linalg_SV_solve( U, V, S, b, x);\r
  56                 \r
  57                 gsl_matrix_free(V);\r
  58                 gsl_vector_free(S);\r
  59                 gsl_vector_free(work);\r
  60                 \r
  61                 return solution();                        \r
  62         }\r
  63         \r
  64         Matrix & matrix()\r
  65         {\r
  66                 return m_matrix;\r
  67         }\r
  68         \r
  69         Vector & vector()\r
  70         {\r
  71                 return m_vector;\r
  72         }\r
  73         \r
  74         const Vector & solution() const\r
  75         {\r
  76                 return m_solution;\r
  77         }\r
  78         \r
  79 private:\r
  80         Matrix & m_matrix;\r
  81         Vector & m_vector;\r
  82         Vector m_solution;\r
  83 };\r
  84 \r
  85 \r
  86 } } // end namespaces\r
  87 \r
  88 \r
  89 #endif /*_NL_LINEAR_SYSTEM_H_*/\r