1 /*
2 * LinearSystem class wraps some gsl routines for solving linear systems
3 *
4 * Authors:
5 * Marco Cecchetti <mrcekets at gmail.com>
6 *
7 * Copyright 2008 authors
8 *
9 * This library is free software; you can redistribute it and/or
10 * modify it either under the terms of the GNU Lesser General Public
11 * License version 2.1 as published by the Free Software Foundation
12 * (the "LGPL") or, at your option, under the terms of the Mozilla
13 * Public License Version 1.1 (the "MPL"). If you do not alter this
14 * notice, a recipient may use your version of this file under either
15 * the MPL or the LGPL.
16 *
17 * You should have received a copy of the LGPL along with this library
18 * in the file COPYING-LGPL-2.1; if not, write to the Free Software
19 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
20 * You should have received a copy of the MPL along with this library
21 * in the file COPYING-MPL-1.1
22 *
23 * The contents of this file are subject to the Mozilla Public License
24 * Version 1.1 (the "License"); you may not use this file except in
25 * compliance with the License. You may obtain a copy of the License at
26 * http://www.mozilla.org/MPL/
27 *
28 * This software is distributed on an "AS IS" basis, WITHOUT WARRANTY
29 * OF ANY KIND, either express or implied. See the LGPL or the MPL for
30 * the specific language governing rights and limitations.
31 */
34 #ifndef _NL_LINEAR_SYSTEM_H_
35 #define _NL_LINEAR_SYSTEM_H_
38 #include <cassert>
40 #include <gsl/gsl_linalg.h>
42 #include <2geom/numeric/matrix.h>
43 #include <2geom/numeric/vector.h>
46 namespace Geom { namespace NL {
49 class LinearSystem
50 {
51 public:
52 LinearSystem(MatrixView & _matrix, VectorView & _vector)
53 : m_matrix(_matrix), m_vector(_vector), m_solution(_matrix.columns())
54 {
55 }
57 LinearSystem(Matrix & _matrix, Vector & _vector)
58 : m_matrix(_matrix), m_vector(_vector), m_solution(_matrix.columns())
59 {
60 }
62 const Vector & LU_solve()
63 {
64 assert( matrix().rows() == matrix().columns()
65 && matrix().rows() == vector().size() );
66 int s;
67 gsl_permutation * p = gsl_permutation_alloc(matrix().rows());
68 gsl_linalg_LU_decomp (matrix().get_gsl_matrix(), p, &s);
69 gsl_linalg_LU_solve( matrix().get_gsl_matrix(),
70 p,
71 vector().get_gsl_vector(),
72 m_solution.get_gsl_vector()
73 );
74 gsl_permutation_free(p);
75 return solution();
76 }
78 const Vector & SV_solve()
79 {
80 assert( matrix().rows() >= matrix().columns()
81 && matrix().rows() == vector().size() );
83 gsl_matrix* U = matrix().get_gsl_matrix();
84 gsl_matrix* V = gsl_matrix_alloc(matrix().columns(), matrix().columns());
85 gsl_vector* S = gsl_vector_alloc(matrix().columns());
86 gsl_vector* work = gsl_vector_alloc(matrix().columns());
88 gsl_linalg_SV_decomp( U, V, S, work );
90 gsl_vector* b = vector().get_gsl_vector();
91 gsl_vector* x = m_solution.get_gsl_vector();
93 gsl_linalg_SV_solve( U, V, S, b, x);
95 gsl_matrix_free(V);
96 gsl_vector_free(S);
97 gsl_vector_free(work);
99 return solution();
100 }
102 MatrixView & matrix()
103 {
104 return m_matrix;
105 }
107 VectorView & vector()
108 {
109 return m_vector;
110 }
112 const Vector & solution() const
113 {
114 return m_solution;
115 }
117 private:
118 MatrixView m_matrix;
119 VectorView m_vector;
120 Vector m_solution;
121 };
124 } } // end namespaces
127 #endif /*_NL_LINEAR_SYSTEM_H_*/
129 /*
130 Local Variables:
131 mode:c++
132 c-file-style:"stroustrup"
133 c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +))
134 indent-tabs-mode:nil
135 fill-column:99
136 End:
137 */
138 // vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:encoding=utf-8:textwidth=99 :