9cb521eb2a2f919a8d7982d46117109362c5643f
1 #ifndef _NL_LINEAR_SYSTEM_H_
2 #define _NL_LINEAR_SYSTEM_H_
5 #include <cassert>
7 #include <gsl/gsl_linalg.h>
9 #include "numeric/matrix.h"
10 #include "numeric/vector.h"
13 namespace Geom { namespace NL {
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 }
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 }
40 const Vector & SV_solve()
41 {
42 assert( matrix().rows() >= matrix().columns()
43 && matrix().rows() == vector().size() );
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());
50 gsl_linalg_SV_decomp( U, V, S, work );
52 gsl_vector* b = vector().get_gsl_vector();
53 gsl_vector* x = m_solution.get_gsl_vector();
55 gsl_linalg_SV_solve( U, V, S, b, x);
57 gsl_matrix_free(V);
58 gsl_vector_free(S);
59 gsl_vector_free(work);
61 return solution();
62 }
64 Matrix & matrix()
65 {
66 return m_matrix;
67 }
69 Vector & vector()
70 {
71 return m_vector;
72 }
74 const Vector & solution() const
75 {
76 return m_solution;
77 }
79 private:
80 Matrix & m_matrix;
81 Vector & m_vector;
82 Vector m_solution;
83 };
86 } } // end namespaces
89 #endif /*_NL_LINEAR_SYSTEM_H_*/