src/2geom/numeric/linear_system.h

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