src/2geom/numeric/linear_system.h

   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  */
  32
  33
  34 #ifndef _NL_LINEAR_SYSTEM_H_
  35 #define _NL_LINEAR_SYSTEM_H_
  36
  37
  38 #include <cassert>
  39
  40 #include <gsl/gsl_linalg.h>
  41
  42 #include <2geom/numeric/matrix.h>
  43 #include <2geom/numeric/vector.h>
  44
  45
  46 namespace Geom { namespace NL {
  47
  48
  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         }
  56
  57         LinearSystem(Matrix & _matrix, Vector & _vector)
  58                 : m_matrix(_matrix), m_vector(_vector), m_solution(_matrix.columns())
  59         {
  60         }
  61
  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         }
  77
  78         const Vector & SV_solve()
  79         {
  80                 assert( matrix().rows() >= matrix().columns()
  81                                 && matrix().rows() == vector().size() );
  82
  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());
  87
  88                 gsl_linalg_SV_decomp( U, V, S, work );
  89
  90                 gsl_vector* b = vector().get_gsl_vector();
  91                 gsl_vector* x = m_solution.get_gsl_vector();
  92
  93                 gsl_linalg_SV_solve( U, V, S, b, x);
  94
  95                 gsl_matrix_free(V);
  96                 gsl_vector_free(S);
  97                 gsl_vector_free(work);
  98
  99                 return solution();
 100         }
 101
 102         MatrixView & matrix()
 103         {
 104                 return m_matrix;
 105         }
 106
 107         VectorView & vector()
 108         {
 109                 return m_vector;
 110         }
 111
 112         const Vector & solution() const
 113         {
 114                 return m_solution;
 115         }
 116
 117 private:
 118         MatrixView m_matrix;
 119         VectorView m_vector;
 120         Vector m_solution;
 121 };
 122
 123
 124 } } // end namespaces
 125
 126
 127 #endif /*_NL_LINEAR_SYSTEM_H_*/
 128
 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:fileencoding=utf-8:textwidth=99 :