tensor-0.1.0
 All Data Structures Namespaces Functions Variables Typedefs Enumerations Enumerator Groups Pages
solve_with_svd_z.cc
1 // -*- mode: c++; fill-column: 80; c-basic-offset: 2; indent-tabs-mode: nil -*-
2 /*
3  Copyright (c) 2010 Juan Jose Garcia Ripoll
4 
5  Tensor is free software; you can redistribute it and/or modify it
6  under the terms of the GNU Library General Public License as published
7  by the Free Software Foundation; either version 2 of the License, or
8  (at your option) any later version.
9 
10  This program is distributed in the hope that it will be useful,
11  but WITHOUT ANY WARRANTY; without even the implied warranty of
12  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13  GNU Library General Public License for more details.
14 
15  You should have received a copy of the GNU General Public License along
16  with this program; if not, write to the Free Software Foundation, Inc.,
17  51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
18 */
19 
20 #include <cfloat>
21 #include <tensor/tensor.h>
22 #include <tensor/linalg.h>
23 
24 namespace linalg {
25 
26  using namespace tensor;
27 
36  const CTensor
37  solve_with_svd(const CTensor &A, const CTensor &B, double tol)
38  {
39  CTensor U, VT;
40  RTensor s = svd(A, &U, &VT, SVD_ECONOMIC);
41  if (tol <= 0) {
42  tol = DBL_EPSILON;
43  }
44  for (tensor::index i = 0; i < s.size(); i++) {
45  if (s[i] <= tol) {
46  U = U(range(), range(0,i-1));
47  s = s(range(0,i-1));
48  VT = VT(range(0,i-1), range());
49  }
50  }
51  CTensor X = foldc(U, 0, B, 0);
52  scale_inplace(X, 0, 1.0/s);
53  return foldc(VT, 0, X, 0);
54  }
55 
56 } // namespace linalg
Real Tensor with elements of type "double".
Definition: tensor.h:349
Complex Tensor with elements of type "cdouble".
Definition: tensor.h:435
index size() const
Returns total number of elements in Tensor.
Definition: tensor.h:114
const RTensor solve_with_svd(const RTensor &A, const RTensor &B, double tol=0.0)
Solution of a linear system of equations using Penrose's pseudoinvese.
RTensor svd(RTensor A, RTensor *pU=0, RTensor *pVT=0, bool economic=0)
Singular value decomposition of a real matrix.
Definition: svd_d.cc:50