function_svd.hpp 1.48 KB
Newer Older
Paul McCarthy's avatar
Paul McCarthy committed
1
2
3
#ifndef __FUNCTION_SVD_HPP__
#define __FUNCTION_SVD_HPP__

Paul McCarthy's avatar
Paul McCarthy committed
4
5
6
/*
 * Singular value decomposition.
 */
Paul McCarthy's avatar
Paul McCarthy committed
7

8
namespace armawrap {
9

10
  template<typename T1, typename T2, typename T3, typename T4>
11
  inline  void SVD(const T1  &A,
12
                   T2        &D,
13
                   T3        &U,
14
15
16
                   T4        &V,
                   bool withU=true,
                   bool withV=true) {
Paul McCarthy's avatar
Paul McCarthy committed
17

18

19
20
21
    arma::Col<typename T1::elem_type> s;
    arma::Mat<typename T1::elem_type> Ut;
    arma::Mat<typename T1::elem_type> At(A);
Paul McCarthy's avatar
Paul McCarthy committed
22

23
24
25
26
27
28
29
30
31
32
33
    /*
     * NEWMAT::SVD does not support decomposition on
     * matrices with more columns than rows, and
     * raises a ProgramException. In armawrap/newmat.h,
     * ProgramException is aliased to runtime_error,
     * so that's what we raise here.
     */
    if (At.n_rows < At.n_cols) {
      throw std::runtime_error("SVD requires that m >= n for a m*n matrix");
    }

Matthew Webster's avatar
Matthew Webster committed
34
    arma::svd_econ(Ut, s, V, At);
Paul McCarthy's avatar
Paul McCarthy committed
35

36
    if (!withV) V = 0;
37
    if ( withU) U = Ut;
38
    D = s;
39
  }
Paul McCarthy's avatar
Paul McCarthy committed
40

41
  template<typename T1, typename T2>
42
  inline void SVD(const T1 &A, T2 &D) {
Paul McCarthy's avatar
Paul McCarthy committed
43

44
45
    AWMatrix<typename T1::elem_type> U;
    AWMatrix<typename T1::elem_type> V;
Paul McCarthy's avatar
Paul McCarthy committed
46

47
48
    SVD(A, D, U, V, false, false);
  }
Paul McCarthy's avatar
Paul McCarthy committed
49

50
  template<typename T1, typename T2, typename T3>
51
  inline void SVD(const T1 &A,
52
53
54
                  T2       &D,
                  T3       &U,
                  bool withU=true) {
Paul McCarthy's avatar
Paul McCarthy committed
55

56
57
58
    AWMatrix<typename T1::elem_type> V(D.Nrows(), D.Ncols());
    SVD(A, D, U, V, withU, false);
  }
59
}
Paul McCarthy's avatar
Paul McCarthy committed
60
61

#endif /* __FUNCTION_SVD_HPP__ */