svd.cpp

有很多的函数库

        Integer p = n_, pp = p-1;
        Integer iter = 0;
        Real eps = std::pow(2.0,-52.0);
        while (p > 0) {
            Integer k;
            Integer kase;

            // Here is where a test for too many iterations would go.

            // This section of the program inspects for
            // negligible elements in the s and e arrays.  On
            // completion the variables kase and k are set as follows.

            // kase = 1     if s(p) and e[k-1] are negligible and k<p
            // kase = 2     if s(k) is negligible and k<p
            // kase = 3     if e[k-1] is negligible, k<p, and
            //              s(k), ..., s(p) are not negligible (qr step).
            // kase = 4     if e(p-1) is negligible (convergence).

            for (k = p-2; k >= -1; --k) {
                if (k == -1) {
                    break;
                }
                if (std::fabs(e[k]) <= eps*(std::fabs(s_[k]) +
                                            std::fabs(s_[k+1]))) {
                    e[k] = 0.0;
                    break;
                }
            }
            if (k == p-2) {
                kase = 4;
            } else {
                Integer ks;
                for (ks = p-1; ks >= k; --ks) {
                    if (ks == k) {
                        break;
                    }
                    Real t = (ks != p ? std::fabs(e[ks]) : 0.) +
                        (ks != k+1 ? std::fabs(e[ks-1]) : 0.);
                    if (std::fabs(s_[ks]) <= eps*t)  {
                        s_[ks] = 0.0;
                        break;
                    }
                }
                if (ks == k) {
                    kase = 3;
                } else if (ks == p-1) {
                    kase = 1;
                } else {
                    kase = 2;
                    k = ks;
                }
            }
            k++;

            // Perform the task indicated by kase.

            switch (kase) {

                // Deflate negligible s(p).

              case 1: {
                  Real f = e[p-2];
                  e[p-2] = 0.0;
                  for (j = p-2; j >= k; --j) {
                      Real t = hypot(s_[j],f);
                      Real cs = s_[j]/t;
                      Real sn = f/t;
                      s_[j] = t;
                      if (j != k) {
                          f = -sn*e[j-1];
                          e[j-1] = cs*e[j-1];
                      }
                      for (i = 0; i < n_; i++) {
                          t = cs*V_[i][j] + sn*V_[i][p-1];
                          V_[i][p-1] = -sn*V_[i][j] + cs*V_[i][p-1];
                          V_[i][j] = t;
                      }
                  }
              }
                break;

                // Split at negligible s(k).

              case 2: {
                  Real f = e[k-1];
                  e[k-1] = 0.0;
                  for (j = k; j < p; j++) {
                      Real t = hypot(s_[j],f);
                      Real cs = s_[j]/t;
                      Real sn = f/t;
                      s_[j] = t;
                      f = -sn*e[j];
                      e[j] = cs*e[j];
                      for (i = 0; i < m_; i++) {
                          t = cs*U_[i][j] + sn*U_[i][k-1];
                          U_[i][k-1] = -sn*U_[i][j] + cs*U_[i][k-1];
                          U_[i][j] = t;
                      }
                  }
              }
                break;

                // Perform one qr step.

              case 3: {

                  // Calculate the shift.
                  Real scale = std::max(
                                     std::max(
                                         std::max(
                                             std::max(std::fabs(s_[p-1]),
                                                    std::fabs(s_[p-2])),
                                             std::fabs(e[p-2])),
                                         std::fabs(s_[k])),
                                     std::fabs(e[k]));
                  Real sp = s_[p-1]/scale;
                  Real spm1 = s_[p-2]/scale;
                  Real epm1 = e[p-2]/scale;
                  Real sk = s_[k]/scale;
                  Real ek = e[k]/scale;
                  Real b = ((spm1 + sp)*(spm1 - sp) + epm1*epm1)/2.0;
                  Real c = (sp*epm1)*(sp*epm1);
                  Real shift = 0.0;
                  if ((b != 0.0) | (c != 0.0)) {
                      shift = std::sqrt(b*b + c);
                      if (b < 0.0) {
                          shift = -shift;
                      }
                      shift = c/(b + shift);
                  }
                  Real f = (sk + sp)*(sk - sp) + shift;
                  Real g = sk*ek;

                  // Chase zeros.

                  for (j = k; j < p-1; j++) {
                      Real t = hypot(f,g);
                      Real cs = f/t;
                      Real sn = g/t;
                      if (j != k) {
                          e[j-1] = t;
                      }
                      f = cs*s_[j] + sn*e[j];
                      e[j] = cs*e[j] - sn*s_[j];
                      g = sn*s_[j+1];
                      s_[j+1] = cs*s_[j+1];
                      for (i = 0; i < n_; i++) {
                          t = cs*V_[i][j] + sn*V_[i][j+1];
                          V_[i][j+1] = -sn*V_[i][j] + cs*V_[i][j+1];
                          V_[i][j] = t;
                      }
                      t = hypot(f,g);
                      cs = f/t;
                      sn = g/t;
                      s_[j] = t;
                      f = cs*e[j] + sn*s_[j+1];
                      s_[j+1] = -sn*e[j] + cs*s_[j+1];
                      g = sn*e[j+1];
                      e[j+1] = cs*e[j+1];
                      if (j < m_-1) {
                          for (i = 0; i < m_; i++) {
                              t = cs*U_[i][j] + sn*U_[i][j+1];
                              U_[i][j+1] = -sn*U_[i][j] + cs*U_[i][j+1];
                              U_[i][j] = t;
                          }
                      }
                  }
                  e[p-2] = f;
                  iter = iter + 1;
              }
                break;

                // Convergence.

              case 4: {

                  // Make the singular values positive.

                  if (s_[k] <= 0.0) {
                      s_[k] = (s_[k] < 0.0 ? -s_[k] : 0.0);
                      for (i = 0; i <= pp; i++) {
                          V_[i][k] = -V_[i][k];
                      }
                  }

                  // Order the singular values.

                  while (k < pp) {
                      if (s_[k] >= s_[k+1]) {
                          break;
                      }
                      swap(s_[k], s_[k+1]);
                      if (k < n_-1) {
                          for (i = 0; i < n_; i++) {
                              swap(V_[i][k], V_[i][k+1]);
                          }
                      }
                      if (k < m_-1) {
                          for (i = 0; i < m_; i++) {
                              swap(U_[i][k], U_[i][k+1]);
                          }
                      }
                      k++;
                  }
                  iter = 0;
                  --p;
              }
                break;
            }
        }
    }

    const Matrix& SVD::U() const {
        return (transpose_ ? V_ : U_);
    }

    const Matrix& SVD::V() const {
        return (transpose_ ? U_ : V_);
    }

    const Array& SVD::singularValues() const {
        return s_;
    }

    Disposable<Matrix> SVD::S() const {
        Matrix S(n_,n_);
        for (Size i = 0; i < Size(n_); i++) {
            for (Size j = 0; j < Size(n_); j++) {
                S[i][j] = 0.0;
            }
            S[i][i] = s_[i];
        }
        return S;
    }

    Real SVD::norm2() {
        return s_[0];
    }

    Real SVD::cond() {
        return s_[0]/s_[n_-1];
    }

    Integer SVD::rank() {
        Real eps = std::pow(2.0,-52.0);
        Real tol = m_*s_[0]*eps;
        Integer r = 0;
        for (Size i = 0; i < s_.size(); i++) {
            if (s_[i] > tol) {
                r++;
            }
        }
        return r;
    }

    Disposable<Array> SVD::solveFor(const Array& b) const{
        Matrix W(n_, n_, 0.0);
        for (Size i=0; i<Size(n_); i++)
            W[i][i] = 1./s_[i];

        Matrix inverse = V()* W * transpose(U());
        Array result = inverse * b;
        return result;
    }

}