src/2geom/numeric/matrix.h

   1 /*
   2  * Matrix, MatrixView, ConstMatrixView classes wrap the gsl matrix routines;
   3  * "views" mimic the semantic of C++ references: any operation performed
   4  * on a "view" is actually performed on the "viewed object"
   5  *
   6  * Authors:
   7  *              Marco Cecchetti <mrcekets at gmail.com>
   8  *
   9  * Copyright 2008  authors
  10  *
  11  * This library is free software; you can redistribute it and/or
  12  * modify it either under the terms of the GNU Lesser General Public
  13  * License version 2.1 as published by the Free Software Foundation
  14  * (the "LGPL") or, at your option, under the terms of the Mozilla
  15  * Public License Version 1.1 (the "MPL"). If you do not alter this
  16  * notice, a recipient may use your version of this file under either
  17  * the MPL or the LGPL.
  18  *
  19  * You should have received a copy of the LGPL along with this library
  20  * in the file COPYING-LGPL-2.1; if not, write to the Free Software
  21  * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
  22  * You should have received a copy of the MPL along with this library
  23  * in the file COPYING-MPL-1.1
  24  *
  25  * The contents of this file are subject to the Mozilla Public License
  26  * Version 1.1 (the "License"); you may not use this file except in
  27  * compliance with the License. You may obtain a copy of the License at
  28  * http://www.mozilla.org/MPL/
  29  *
  30  * This software is distributed on an "AS IS" basis, WITHOUT WARRANTY
  31  * OF ANY KIND, either express or implied. See the LGPL or the MPL for
  32  * the specific language governing rights and limitations.
  33  */
  34
  35
  36
  37
  38 #ifndef _NL_MATRIX_H_
  39 #define _NL_MATRIX_H_
  40
  41 #include <2geom/numeric/vector.h>
  42
  43 #include <cassert>
  44 #include <utility>    // for std::pair
  45 #include <algorithm>  // for std::swap
  46 #include <sstream>
  47 #include <string>
  48
  49 #include <gsl/gsl_matrix.h>
  50 #include <gsl/gsl_linalg.h>
  51
  52
  53 namespace Geom { namespace NL {
  54
  55 namespace detail
  56 {
  57
  58 class BaseMatrixImpl
  59 {
  60   public:
  61         virtual ~BaseMatrixImpl()
  62         {
  63         }
  64
  65         ConstVectorView row_const_view(size_t i) const
  66         {
  67                 return ConstVectorView(gsl_matrix_const_row(m_matrix, i));
  68         }
  69
  70         ConstVectorView column_const_view(size_t i) const
  71         {
  72                 return ConstVectorView(gsl_matrix_const_column(m_matrix, i));
  73         }
  74
  75         const double & operator() (size_t i, size_t j) const
  76         {
  77                 return *gsl_matrix_const_ptr(m_matrix, i, j);
  78         }
  79
  80         const gsl_matrix* get_gsl_matrix() const
  81         {
  82                 return m_matrix;
  83         }
  84
  85         bool is_zero() const
  86         {
  87                 return gsl_matrix_isnull(m_matrix);
  88         }
  89
  90         bool is_positive() const
  91         {
  92         for ( unsigned int i = 0; i < rows(); ++i )
  93         {
  94             for ( unsigned int j = 0; j < columns(); ++j )
  95             {
  96                 if ( (*this)(i,j) <= 0 ) return false;
  97             }
  98         }
  99         return true;
 100         }
 101
 102         bool is_negative() const
 103         {
 104         for ( unsigned int i = 0; i < rows(); ++i )
 105         {
 106             for ( unsigned int j = 0; j < columns(); ++j )
 107             {
 108                 if ( (*this)(i,j) >= 0 ) return false;
 109             }
 110         }
 111         return true;
 112         }
 113
 114         bool is_non_negative() const
 115         {
 116                 for ( unsigned int i = 0; i < rows(); ++i )
 117                 {
 118                         for ( unsigned int j = 0; j < columns(); ++j )
 119                         {
 120                                 if ( (*this)(i,j) < 0 ) return false;
 121                         }
 122                 }
 123                 return true;
 124         }
 125
 126         double max() const
 127         {
 128                 return gsl_matrix_max(m_matrix);
 129         }
 130
 131         double min() const
 132         {
 133                 return gsl_matrix_min(m_matrix);
 134         }
 135
 136         std::pair<size_t, size_t>
 137         max_index() const
 138         {
 139                 std::pair<size_t, size_t> indices;
 140                 gsl_matrix_max_index(m_matrix, &(indices.first), &(indices.second));
 141                 return indices;
 142         }
 143
 144         std::pair<size_t, size_t>
 145         min_index() const
 146         {
 147                 std::pair<size_t, size_t> indices;
 148                 gsl_matrix_min_index(m_matrix, &(indices.first), &(indices.second));
 149                 return indices;
 150         }
 151
 152         size_t rows() const
 153         {
 154                 return m_rows;
 155         }
 156
 157         size_t columns() const
 158         {
 159                 return m_columns;
 160         }
 161
 162         std::string str() const;
 163
 164   protected:
 165         size_t m_rows, m_columns;
 166         gsl_matrix* m_matrix;
 167
 168 };  // end class BaseMatrixImpl
 169
 170
 171 inline
 172 bool operator== (BaseMatrixImpl const& m1, BaseMatrixImpl const& m2)
 173 {
 174         if (m1.rows() != m2.rows() || m1.columns() != m2.columns()) return false;
 175
 176         for (size_t i = 0; i < m1.rows(); ++i)
 177                 for (size_t j = 0; j < m1.columns(); ++j)
 178                         if (m1(i,j) != m2(i,j)) return false;
 179
 180         return true;
 181 }
 182
 183 template< class charT >
 184 inline
 185 std::basic_ostream<charT> &
 186 operator<< (std::basic_ostream<charT> & os, const BaseMatrixImpl & _matrix)
 187 {
 188         if (_matrix.rows() == 0 || _matrix.columns() == 0) return os;
 189
 190         os << "[[" << _matrix(0,0);
 191         for (size_t j = 1; j < _matrix.columns(); ++j)
 192         {
 193                 os << ", " << _matrix(0,j);
 194         }
 195         os << "]";
 196
 197         for (size_t i = 1; i < _matrix.rows(); ++i)
 198         {
 199                 os << ", [" << _matrix(i,0);
 200                 for (size_t j = 1; j < _matrix.columns(); ++j)
 201                 {
 202                         os << ", " << _matrix(i,j);
 203                 }
 204                 os << "]";
 205         }
 206         os << "]";
 207         return os;
 208 }
 209
 210 inline
 211 std::string BaseMatrixImpl::str() const
 212 {
 213         std::ostringstream oss;
 214         oss << (*this);
 215         return oss.str();
 216 }
 217
 218
 219 class MatrixImpl : public BaseMatrixImpl
 220 {
 221   public:
 222
 223         typedef BaseMatrixImpl base_type;
 224
 225         void set_all( double x )
 226         {
 227                 gsl_matrix_set_all(m_matrix, x);
 228         }
 229
 230         void set_identity()
 231         {
 232                 gsl_matrix_set_identity(m_matrix);
 233         }
 234
 235         using base_type::operator(); // VSC legacy support
 236         const double & operator() (size_t i, size_t j) const
 237         {
 238             return base_type::operator ()(i, j);
 239         }
 240
 241         double & operator() (size_t i, size_t j)
 242         {
 243             return *gsl_matrix_ptr(m_matrix, i, j);
 244         }
 245
 246         using base_type::get_gsl_matrix;
 247
 248         gsl_matrix* get_gsl_matrix()
 249         {
 250                 return m_matrix;
 251         }
 252
 253         VectorView row_view(size_t i)
 254         {
 255                 return VectorView(gsl_matrix_row(m_matrix, i));
 256         }
 257
 258         VectorView column_view(size_t i)
 259         {
 260                 return VectorView(gsl_matrix_column(m_matrix, i));
 261         }
 262
 263         void swap_rows(size_t i, size_t j)
 264         {
 265                  gsl_matrix_swap_rows(m_matrix, i, j);
 266         }
 267
 268         void swap_columns(size_t i, size_t j)
 269         {
 270                 gsl_matrix_swap_columns(m_matrix, i, j);
 271         }
 272
 273         MatrixImpl & transpose()
 274         {
 275             assert(columns() == rows());
 276                 gsl_matrix_transpose(m_matrix);
 277                 return (*this);
 278         }
 279
 280         MatrixImpl & scale(double x)
 281         {
 282                 gsl_matrix_scale(m_matrix, x);
 283                 return (*this);
 284         }
 285
 286         MatrixImpl & translate(double x)
 287         {
 288                  gsl_matrix_add_constant(m_matrix, x);
 289                  return (*this);
 290         }
 291
 292         MatrixImpl & operator+=(base_type const& _matrix)
 293         {
 294                 gsl_matrix_add(m_matrix, _matrix.get_gsl_matrix());
 295                 return (*this);
 296         }
 297
 298         MatrixImpl & operator-=(base_type const& _matrix)
 299         {
 300                 gsl_matrix_sub(m_matrix, _matrix.get_gsl_matrix());
 301                 return (*this);
 302         }
 303
 304 }; // end class MatrixImpl
 305
 306 }  // end namespace detail
 307
 308
 309 using detail::operator==;
 310 using detail::operator<<;
 311
 312
 313
 314
 315 class Matrix: public detail::MatrixImpl
 316 {
 317   public:
 318         typedef detail::MatrixImpl base_type;
 319
 320   public:
 321         // the matrix is not inizialized
 322         Matrix(size_t n1, size_t n2)
 323         {
 324                 m_rows = n1;
 325                 m_columns = n2;
 326                 m_matrix = gsl_matrix_alloc(n1, n2);
 327         }
 328
 329         Matrix(size_t n1, size_t n2, double x)
 330         {
 331                 m_rows = n1;
 332                 m_columns = n2;
 333                 m_matrix = gsl_matrix_alloc(n1, n2);
 334                 gsl_matrix_set_all(m_matrix, x);
 335         }
 336
 337         Matrix(Matrix const& _matrix)
 338         : base_type()
 339         {
 340                 m_rows = _matrix.rows();
 341                 m_columns = _matrix.columns();
 342                 m_matrix = gsl_matrix_alloc(rows(), columns());
 343                 gsl_matrix_memcpy(m_matrix, _matrix.get_gsl_matrix());
 344         }
 345
 346         explicit
 347         Matrix(base_type::base_type const& _matrix)
 348         {
 349                 m_rows = _matrix.rows();
 350                 m_columns = _matrix.columns();
 351                 m_matrix = gsl_matrix_alloc(rows(), columns());
 352                 gsl_matrix_memcpy(m_matrix, _matrix.get_gsl_matrix());
 353         }
 354
 355         Matrix & operator=(Matrix const& _matrix)
 356         {
 357                 assert( rows() == _matrix.rows() && columns() ==  _matrix.columns() );
 358                 gsl_matrix_memcpy(m_matrix, _matrix.get_gsl_matrix());
 359                 return *this;
 360         }
 361
 362         Matrix & operator=(base_type::base_type const& _matrix)
 363         {
 364                 assert( rows() == _matrix.rows() && columns() ==  _matrix.columns() );
 365                 gsl_matrix_memcpy(m_matrix, _matrix.get_gsl_matrix());
 366                 return *this;
 367         }
 368
 369         virtual ~Matrix()
 370         {
 371                 gsl_matrix_free(m_matrix);
 372         }
 373
 374         Matrix & transpose()
 375         {
 376                 return static_cast<Matrix &>( base_type::transpose() );
 377         }
 378
 379         Matrix & scale(double x)
 380         {
 381                 return static_cast<Matrix &>( base_type::scale(x) );
 382         }
 383
 384         Matrix & translate(double x)
 385         {
 386                 return static_cast<Matrix &>( base_type::translate(x) );
 387         }
 388
 389         Matrix & operator+=(base_type::base_type const& _matrix)
 390         {
 391                 return static_cast<Matrix &>( base_type::operator+=(_matrix) );
 392         }
 393
 394         Matrix & operator-=(base_type::base_type const& _matrix)
 395         {
 396                 return static_cast<Matrix &>( base_type::operator-=(_matrix) );
 397         }
 398
 399         friend
 400         void swap(Matrix & m1, Matrix & m2);
 401         friend
 402         void swap_any(Matrix & m1, Matrix & m2);
 403
 404 };  // end class Matrix
 405
 406
 407 // warning! this operation invalidates any view of the passed matrix objects
 408 inline
 409 void swap(Matrix & m1, Matrix & m2)
 410 {
 411         assert( m1.rows() == m2.rows() && m1.columns() ==  m2.columns() );
 412         std::swap(m1.m_matrix, m2.m_matrix);
 413 }
 414
 415 inline
 416 void swap_any(Matrix & m1, Matrix & m2)
 417 {
 418     std::swap(m1.m_matrix, m2.m_matrix);
 419     std::swap(m1.m_rows, m2.m_rows);
 420     std::swap(m1.m_columns, m2.m_columns);
 421 }
 422
 423
 424
 425 class ConstMatrixView : public detail::BaseMatrixImpl
 426 {
 427   public:
 428         typedef detail::BaseMatrixImpl base_type;
 429
 430   public:
 431         ConstMatrixView(const base_type & _matrix, size_t k1, size_t k2, size_t n1, size_t n2)
 432                 : m_matrix_view( gsl_matrix_const_submatrix(_matrix.get_gsl_matrix(), k1, k2, n1, n2) )
 433         {
 434                 m_rows = n1;
 435                 m_columns = n2;
 436                 m_matrix = const_cast<gsl_matrix*>( &(m_matrix_view.matrix) );
 437         }
 438
 439         ConstMatrixView(const ConstMatrixView & _matrix)
 440                 : base_type(),
 441                   m_matrix_view(_matrix.m_matrix_view)
 442         {
 443                 m_rows = _matrix.rows();
 444                 m_columns = _matrix.columns();
 445                 m_matrix = const_cast<gsl_matrix*>( &(m_matrix_view.matrix) );
 446         }
 447
 448         ConstMatrixView(const base_type & _matrix)
 449                 : m_matrix_view(gsl_matrix_const_submatrix(_matrix.get_gsl_matrix(), 0, 0, _matrix.rows(), _matrix.columns()))
 450         {
 451                 m_rows = _matrix.rows();
 452                 m_columns = _matrix.columns();
 453                 m_matrix = const_cast<gsl_matrix*>( &(m_matrix_view.matrix) );
 454         }
 455
 456   private:
 457         gsl_matrix_const_view m_matrix_view;
 458
 459 };  // end class ConstMatrixView
 460
 461
 462
 463
 464 class MatrixView : public detail::MatrixImpl
 465 {
 466   public:
 467         typedef detail::MatrixImpl base_type;
 468
 469   public:
 470         MatrixView(base_type & _matrix, size_t k1, size_t k2, size_t n1, size_t n2)
 471         {
 472                 m_rows = n1;
 473                 m_columns = n2;
 474                 m_matrix_view
 475                         = gsl_matrix_submatrix(_matrix.get_gsl_matrix(), k1, k2, n1, n2);
 476                 m_matrix = &(m_matrix_view.matrix);
 477         }
 478
 479         MatrixView(const MatrixView & _matrix)
 480         : base_type()
 481         {
 482                 m_rows = _matrix.rows();
 483                 m_columns = _matrix.columns();
 484                 m_matrix_view = _matrix.m_matrix_view;
 485                 m_matrix = &(m_matrix_view.matrix);
 486         }
 487
 488         MatrixView(Matrix & _matrix)
 489         {
 490                 m_rows = _matrix.rows();
 491                 m_columns = _matrix.columns();
 492                 m_matrix_view
 493                         = gsl_matrix_submatrix(_matrix.get_gsl_matrix(), 0, 0, rows(), columns());
 494                 m_matrix = &(m_matrix_view.matrix);
 495         }
 496
 497         MatrixView & operator=(MatrixView const& _matrix)
 498         {
 499                 assert( rows() == _matrix.rows() && columns() ==  _matrix.columns() );
 500                 gsl_matrix_memcpy(m_matrix, _matrix.m_matrix);
 501                 return *this;
 502         }
 503
 504         MatrixView & operator=(base_type::base_type const& _matrix)
 505         {
 506                 assert( rows() == _matrix.rows() && columns() ==  _matrix.columns() );
 507                 gsl_matrix_memcpy(m_matrix, _matrix.get_gsl_matrix());
 508                 return *this;
 509         }
 510
 511         MatrixView & transpose()
 512         {
 513                 return static_cast<MatrixView &>( base_type::transpose() );
 514         }
 515
 516         MatrixView & scale(double x)
 517         {
 518                 return static_cast<MatrixView &>( base_type::scale(x) );
 519         }
 520
 521         MatrixView & translate(double x)
 522         {
 523                 return static_cast<MatrixView &>( base_type::translate(x) );
 524         }
 525
 526         MatrixView & operator+=(base_type::base_type const& _matrix)
 527         {
 528                 return static_cast<MatrixView &>( base_type::operator+=(_matrix) );
 529         }
 530
 531         MatrixView & operator-=(base_type::base_type const& _matrix)
 532         {
 533                 return static_cast<MatrixView &>( base_type::operator-=(_matrix) );
 534         }
 535
 536         friend
 537         void swap_view(MatrixView & m1, MatrixView & m2);
 538
 539   private:
 540         gsl_matrix_view m_matrix_view;
 541
 542 };  // end class MatrixView
 543
 544
 545 inline
 546 void swap_view(MatrixView & m1, MatrixView & m2)
 547 {
 548         assert( m1.rows() == m2.rows() && m1.columns() ==  m2.columns() );
 549         std::swap(m1.m_matrix_view, m2.m_matrix_view);
 550 }
 551
 552 Vector operator*( detail::BaseMatrixImpl const& A,
 553                   detail::BaseVectorImpl const& v );
 554
 555 Matrix operator*( detail::BaseMatrixImpl const& A,
 556                   detail::BaseMatrixImpl const& B );
 557
 558 Matrix pseudo_inverse(detail::BaseMatrixImpl const& A);
 559
 560 } } // end namespaces
 561
 562 #endif /*_NL_MATRIX_H_*/
 563
 564 /*
 565   Local Variables:
 566   mode:c++
 567   c-file-style:"stroustrup"
 568   c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +))
 569   indent-tabs-mode:nil
 570   fill-column:99
 571   End:
 572 */
 573 // vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:fileencoding=utf-8:textwidth=99 :