src/2geom/numeric/linear_system.h

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