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 */
36 #include <numeric/matrix.h>
37 #include <numeric/vector.h>
41 namespace Geom { namespace NL {
44 Vector operator*( detail::BaseMatrixImpl const& A,
45 detail::BaseVectorImpl const& v )
46 {
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];
54 return result;
55 }
57 Matrix operator*( detail::BaseMatrixImpl const& A,
58 detail::BaseMatrixImpl const& B )
59 {
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);
68 return C;
69 }
71 Matrix pseudo_inverse(detail::BaseMatrixImpl const& A)
72 {
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());
79 gsl_linalg_SV_decomp( U.get_gsl_matrix(),
80 V.get_gsl_matrix(),
81 s.get_gsl_vector(),
82 work );
84 Matrix P(A.columns(), A.rows(), 0.0);
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);
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 }
101 return P;
102 }
104 } } // end namespaces
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:encoding=utf-8:textwidth=99 :