src/2geom/numeric/matrix.cpp

   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 #include <2geom/numeric/matrix.h>
  37 #include <2geom/numeric/vector.h>
  38
  39
  40
  41 namespace Geom { namespace NL {
  42
  43
  44 Vector operator*( detail::BaseMatrixImpl const& A,
  45                   detail::BaseVectorImpl const& v )
  46 {
  47     assert(A.columns() == v.size());
  48
  49     Vector result(A.rows(), 0.0);
  50     for (size_t i = 0; i < A.rows(); ++i)
  51         for (size_t j = 0; j < A.columns(); ++j)
  52             result[i] += A(i,j) * v[j];
  53
  54     return result;
  55 }
  56
  57 Matrix operator*( detail::BaseMatrixImpl const& A,
  58                   detail::BaseMatrixImpl const& B )
  59 {
  60     assert(A.columns() == B.rows());
  61
  62     Matrix C(A.rows(), B.columns(), 0.0);
  63     for (size_t i = 0; i < C.rows(); ++i)
  64         for (size_t j = 0; j < C.columns(); ++j)
  65             for (size_t k = 0; k < A.columns(); ++k)
  66                 C(i,j) += A(i,k) * B(k, j);
  67
  68     return C;
  69 }
  70
  71 Matrix pseudo_inverse(detail::BaseMatrixImpl const& A)
  72 {
  73
  74     Matrix U(A);
  75     Matrix V(A.columns(), A.columns());
  76     Vector s(A.columns());
  77     gsl_vector* work = gsl_vector_alloc(A.columns());
  78
  79     gsl_linalg_SV_decomp( U.get_gsl_matrix(),
  80                           V.get_gsl_matrix(),
  81                           s.get_gsl_vector(),
  82                           work );
  83
  84     Matrix P(A.columns(), A.rows(), 0.0);
  85
  86     int sz = s.size();
  87     while ( sz-- > 0 && s[sz] == 0 ) {}
  88     ++sz;
  89     if (sz == 0) return P;
  90     VectorView sv(s, sz);
  91
  92     for (size_t i = 0; i < sv.size(); ++i)
  93     {
  94         VectorView v = V.column_view(i);
  95         v.scale(1/sv[i]);
  96         for (size_t h = 0; h < P.rows(); ++h)
  97             for (size_t k = 0; k < P.columns(); ++k)
  98                 P(h,k) += V(h,i) * U(k,i);
  99     }
 100
 101     return P;
 102 }
 103
 104 } }  // end namespaces
 105
 106 /*
 107   Local Variables:
 108   mode:c++
 109   c-file-style:"stroustrup"
 110   c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +))
 111   indent-tabs-mode:nil
 112   fill-column:99
 113   End:
 114 */
 115 // vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:fileencoding=utf-8:textwidth=99 :