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"
7
* Marco Cecchetti <mrcekets at gmail.com>
9
* Copyright 2008 authors
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.
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
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/
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.
36
#include <2geom/numeric/matrix.h>
37
#include <2geom/numeric/vector.h>
41
namespace Geom { namespace NL {
44
Vector operator*( detail::BaseMatrixImpl const& A,
45
detail::BaseVectorImpl const& v )
47
assert(A.columns() == v.size());
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];
57
Matrix operator*( detail::BaseMatrixImpl const& A,
58
detail::BaseMatrixImpl const& B )
60
assert(A.columns() == B.rows());
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);
71
Matrix pseudo_inverse(detail::BaseMatrixImpl const& A)
75
Matrix V(A.columns(), A.columns());
76
Vector s(A.columns());
77
gsl_vector* work = gsl_vector_alloc(A.columns());
79
gsl_linalg_SV_decomp( U.get_gsl_matrix(),
84
Matrix P(A.columns(), A.rows(), 0.0);
87
while ( sz-- > 0 && s[sz] == 0 ) {}
89
if (sz == 0) return P;
92
for (size_t i = 0; i < sv.size(); ++i)
94
VectorView v = V.column_view(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);
104
} } // end namespaces
109
c-file-style:"stroustrup"
110
c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +))
115
// vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:encoding=utf-8:textwidth=99 :