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