src/math/svd.cpp

/* -*- mode: c++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- */

/*
 Copyright (C) 2003 Neil Firth

 This file is part of QuantLib, a free-software/open-source library
 for financial quantitative analysts and developers - http://quantlib.org/

 QuantLib is free software: you can redistribute it and/or modify it
 under the terms of the QuantLib license.  You should have received a
 copy of the license along with this program; if not, please email
 <quantlib-dev@lists.sf.net>. The license is also available online at
 <http://quantlib.org/license.shtml>.

 This program is distributed in the hope that it will be useful, but WITHOUT
 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
 FOR A PARTICULAR PURPOSE.  See the license for more details.

        Adapted from the TNT project
        http://math.nist.gov/tnt/download.html

        This software was developed at the National Institute of Standards
        and Technology (NIST) by employees of the Federal Government in the
        course of their official duties. Pursuant to title 17 Section 105
        of the United States Code this software is not subject to copyright
        protection and is in the public domain. NIST assumes no responsibility
        whatsoever for its use by other parties, and makes no guarantees,
        expressed or implied, about its quality, reliability, or any other
        characteristic.

        We would appreciate acknowledgement if the software is incorporated in
        redistributable libraries or applications.
*/


#include "math/svd.hpp"

namespace oopse {

    namespace {

        /*  returns hypotenuse of real (non-complex) scalars a and b by
            avoiding underflow/overflow
            using (a * sqrt( 1 + (b/a) * (b/a))), rather than
            sqrt(a*a + b*b).
        */
        Real hypot(const Real &a, const Real &b) {
            if (a == 0) {
                return std::fabs(b);
            } else {
                Real c = b/a;
                return std::fabs(a) * std::sqrt(1 + c*c);
            }
        }

    }


    SVD::SVD(const RectMatrix& M) {

        using std::swap;

        RectMatrix A;

        /* The implementation requires that rows > columns.
           If this is not the case, we decompose M^T instead.
           Swapping the resulting U and V gives the desired
           result for M as

           M^T = U S V^T           (decomposition of M^T)

           M = (U S V^T)^T         (transpose)

           M = (V^T^T S^T U^T)     ((AB)^T = B^T A^T)

           M = V S U^T             (idempotence of transposition,
                                    symmetry of diagonal matrix S)

        */

        if (M.rows() >= M.columns()) {
            A = M;
            transpose_ = false;
        } else {
            A = transpose(M);
            transpose_ = true;
        }

        m_ = A.rows();
        n_ = A.columns();

        // we're sure that m_ >= n_

        s_ = Vector(n_);
        U_ = RectMatrix(m_,n_, 0.0);
        V_ = RectMatrix(n_,n_);
        Vector e(n_);
        Vector work(m_);
        int i, j, k;

        // Reduce A to bidiagonal form, storing the diagonal elements
        // in s and the super-diagonal elements in e.

        int nct = std::min(m_-1,n_);
        int nrt = std::max(0,n_-2);
        for (k = 0; k < std::max(nct,nrt); k++) {
            if (k < nct) {

                // Compute the transformation for the k-th column and
                // place the k-th diagonal in s[k].
                // Compute 2-norm of k-th column without under/overflow.
                s_[k] = 0;
                for (i = k; i < m_; i++) {
                    s_[k] = hypot(s_[k],A[i][k]);
                }
                if (s_[k] != 0.0) {
                    if (A[k][k] < 0.0) {
                        s_[k] = -s_[k];
                    }
                    for (i = k; i < m_; i++) {
                        A[i][k] /= s_[k];
                    }
                    A[k][k] += 1.0;
                }
                s_[k] = -s_[k];
            }
            for (j = k+1; j < n_; j++) {
                if ((k < nct) && (s_[k] != 0.0))  {

                    // Apply the transformation.

                    Real t = 0;
                    for (i = k; i < m_; i++) {
                        t += A[i][k]*A[i][j];
                    }
                    t = -t/A[k][k];
                    for (i = k; i < m_; i++) {
                        A[i][j] += t*A[i][k];
                    }
                }

                // Place the k-th row of A into e for the
                // subsequent calculation of the row transformation.

                e[j] = A[k][j];
            }
            if (k < nct) {

                // Place the transformation in U for subsequent back
                // multiplication.

                for (i = k; i < m_; i++) {
                    U_[i][k] = A[i][k];
                }
            }
            if (k < nrt) {

                // Compute the k-th row transformation and place the
                // k-th super-diagonal in e[k].
                // Compute 2-norm without under/overflow.
                e[k] = 0;
                for (i = k+1; i < n_; i++) {
                    e[k] = hypot(e[k],e[i]);
                }
                if (e[k] != 0.0) {
                    if (e[k+1] < 0.0) {
                        e[k] = -e[k];
                    }
                    for (i = k+1; i < n_; i++) {
                        e[i] /= e[k];
                    }
                    e[k+1] += 1.0;
                }
                e[k] = -e[k];
                if ((k+1 < m_) & (e[k] != 0.0)) {

                    // Apply the transformation.

                    for (i = k+1; i < m_; i++) {
                        work[i] = 0.0;
                    }
                    for (j = k+1; j < n_; j++) {
                        for (i = k+1; i < m_; i++) {
                            work[i] += e[j]*A[i][j];
                        }
                    }
                    for (j = k+1; j < n_; j++) {
                        Real t = -e[j]/e[k+1];
                        for (i = k+1; i < m_; i++) {
                            A[i][j] += t*work[i];
                        }
                    }
                }

                // Place the transformation in V for subsequent
                // back multiplication.

                for (i = k+1; i < n_; i++) {
                    V_[i][k] = e[i];
                }
            }
        }

        // Set up the final bidiagonal matrix or order n.

        if (nct < n_) {
            s_[nct] = A[nct][nct];
        }
        if (nrt+1 < n_) {
            e[nrt] = A[nrt][n_-1];
        }
        e[n_-1] = 0.0;

        // generate U

        for (j = nct; j < n_; j++) {
            for (i = 0; i < m_; i++) {
                U_[i][j] = 0.0;
            }
            U_[j][j] = 1.0;
        }
        for (k = nct-1; k >= 0; --k) {
            if (s_[k] != 0.0) {
                for (j = k+1; j < n_; ++j) {
                    Real t = 0;
                    for (i = k; i < m_; i++) {
                        t += U_[i][k]*U_[i][j];
                    }
                    t = -t/U_[k][k];
                    for (i = k; i < m_; i++) {
                        U_[i][j] += t*U_[i][k];
                    }
                }
                for (i = k; i < m_; i++ ) {
                    U_[i][k] = -U_[i][k];
                }
                U_[k][k] = 1.0 + U_[k][k];
                for (i = 0; i < k-1; i++) {
                    U_[i][k] = 0.0;
                }
            } else {
                for (i = 0; i < m_; i++) {
                    U_[i][k] = 0.0;
                }
                U_[k][k] = 1.0;
            }
        }

        // generate V

        for (k = n_-1; k >= 0; --k) {
            if ((k < nrt) & (e[k] != 0.0)) {
                for (j = k+1; j < n_; ++j) {
                    Real t = 0;
                    for (i = k+1; i < n_; i++) {
                        t += V_[i][k]*V_[i][j];
                    }
                    t = -t/V_[k+1][k];
                    for (i = k+1; i < n_; i++) {
                        V_[i][j] += t*V_[i][k];
                    }
                }
            }
            for (i = 0; i < n_; i++) {
                V_[i][k] = 0.0;
            }
            V_[k][k] = 1.0;
        }

        // Main iteration loop for the singular values.

        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 RectMatrix& SVD::U() const {
        return (transpose_ ? V_ : U_);
    }

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

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

    RectMatrix 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;
    }

    Vector SVD::solveFor(const Array& b) const{
        RectMatrix W(n_, n_, 0.0);
        for (Size i=0; i<Size(n_); i++)
            W[i][i] = 1./s_[i];
        
        RectMatrix inverse = V()* W * transpose(U());
        Vector result = inverse * b;
        return result;
    }

}

Revision:	1335
Committed:	Fri Apr 3 17:44:57 2009 UTC (16 years, 6 months ago) by gezelter
File size:	17083 byte(s)
Log Message:	Beginning import of SVD and RMSD code to handle orientational restraints of flexible molecules. Work in progress, doesn't compile yet.
#	User	Rev	Content
1	gezelter	1335	/* -- mode: c++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -- */
2
3			/*
4			Copyright (C) 2003 Neil Firth
5
6			This file is part of QuantLib, a free-software/open-source library
7			for financial quantitative analysts and developers - http://quantlib.org/
8
9			QuantLib is free software: you can redistribute it and/or modify it
10			under the terms of the QuantLib license. You should have received a
11			copy of the license along with this program; if not, please email
12			<quantlib-dev@lists.sf.net>. The license is also available online at
13			<http://quantlib.org/license.shtml>.
14
15			This program is distributed in the hope that it will be useful, but WITHOUT
16			ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
17			FOR A PARTICULAR PURPOSE. See the license for more details.
18
19			Adapted from the TNT project
20			http://math.nist.gov/tnt/download.html
21
22			This software was developed at the National Institute of Standards
23			and Technology (NIST) by employees of the Federal Government in the
24			course of their official duties. Pursuant to title 17 Section 105
25			of the United States Code this software is not subject to copyright
26			protection and is in the public domain. NIST assumes no responsibility
27			whatsoever for its use by other parties, and makes no guarantees,
28			expressed or implied, about its quality, reliability, or any other
29			characteristic.
30
31			We would appreciate acknowledgement if the software is incorporated in
32			redistributable libraries or applications.
33			*/
34
35
36			#include "math/svd.hpp"
37
38			namespace oopse {
39
40			namespace {
41
42			/* returns hypotenuse of real (non-complex) scalars a and b by
43			avoiding underflow/overflow
44			using (a * sqrt( 1 + (b/a) * (b/a))), rather than
45			sqrt(aa + bb).
46			*/
47			Real hypot(const Real &a, const Real &b) {
48			if (a == 0) {
49			return std::fabs(b);
50			} else {
51			Real c = b/a;
52			return std::fabs(a) * std::sqrt(1 + c*c);
53			}
54			}
55
56			}
57
58
59			SVD::SVD(const RectMatrix& M) {
60
61			using std::swap;
62
63			RectMatrix A;
64
65			/* The implementation requires that rows > columns.
66			If this is not the case, we decompose M^T instead.
67			Swapping the resulting U and V gives the desired
68			result for M as
69
70			M^T = U S V^T (decomposition of M^T)
71
72			M = (U S V^T)^T (transpose)
73
74			M = (V^T^T S^T U^T) ((AB)^T = B^T A^T)
75
76			M = V S U^T (idempotence of transposition,
77			symmetry of diagonal matrix S)
78
79			*/
80
81			if (M.rows() >= M.columns()) {
82			A = M;
83			transpose_ = false;
84			} else {
85			A = transpose(M);
86			transpose_ = true;
87			}
88
89			m_ = A.rows();
90			n_ = A.columns();
91
92			// we're sure that m_ >= n_
93
94			s_ = Vector(n_);
95			U_ = RectMatrix(m_,n_, 0.0);
96			V_ = RectMatrix(n_,n_);
97			Vector e(n_);
98			Vector work(m_);
99			int i, j, k;
100
101			// Reduce A to bidiagonal form, storing the diagonal elements
102			// in s and the super-diagonal elements in e.
103
104			int nct = std::min(m_-1,n_);
105			int nrt = std::max(0,n_-2);
106			for (k = 0; k < std::max(nct,nrt); k++) {
107			if (k < nct) {
108
109			// Compute the transformation for the k-th column and
110			// place the k-th diagonal in s[k].
111			// Compute 2-norm of k-th column without under/overflow.
112			s_[k] = 0;
113			for (i = k; i < m_; i++) {
114			s_[k] = hypot(s_[k],A[i][k]);
115			}
116			if (s_[k] != 0.0) {
117			if (A[k][k] < 0.0) {
118			s_[k] = -s_[k];
119			}
120			for (i = k; i < m_; i++) {
121			A[i][k] /= s_[k];
122			}
123			A[k][k] += 1.0;
124			}
125			s_[k] = -s_[k];
126			}
127			for (j = k+1; j < n_; j++) {
128			if ((k < nct) && (s_[k] != 0.0)) {
129
130			// Apply the transformation.
131
132			Real t = 0;
133			for (i = k; i < m_; i++) {
134			t += A[i][k]*A[i][j];
135			}
136			t = -t/A[k][k];
137			for (i = k; i < m_; i++) {
138			A[i][j] += t*A[i][k];
139			}
140			}
141
142			// Place the k-th row of A into e for the
143			// subsequent calculation of the row transformation.
144
145			e[j] = A[k][j];
146			}
147			if (k < nct) {
148
149			// Place the transformation in U for subsequent back
150			// multiplication.
151
152			for (i = k; i < m_; i++) {
153			U_[i][k] = A[i][k];
154			}
155			}
156			if (k < nrt) {
157
158			// Compute the k-th row transformation and place the
159			// k-th super-diagonal in e[k].
160			// Compute 2-norm without under/overflow.
161			e[k] = 0;
162			for (i = k+1; i < n_; i++) {
163			e[k] = hypot(e[k],e[i]);
164			}
165			if (e[k] != 0.0) {
166			if (e[k+1] < 0.0) {
167			e[k] = -e[k];
168			}
169			for (i = k+1; i < n_; i++) {
170			e[i] /= e[k];
171			}
172			e[k+1] += 1.0;
173			}
174			e[k] = -e[k];
175			if ((k+1 < m_) & (e[k] != 0.0)) {
176
177			// Apply the transformation.
178
179			for (i = k+1; i < m_; i++) {
180			work[i] = 0.0;
181			}
182			for (j = k+1; j < n_; j++) {
183			for (i = k+1; i < m_; i++) {
184			work[i] += e[j]*A[i][j];
185			}
186			}
187			for (j = k+1; j < n_; j++) {
188			Real t = -e[j]/e[k+1];
189			for (i = k+1; i < m_; i++) {
190			A[i][j] += t*work[i];
191			}
192			}
193			}
194
195			// Place the transformation in V for subsequent
196			// back multiplication.
197
198			for (i = k+1; i < n_; i++) {
199			V_[i][k] = e[i];
200			}
201			}
202			}
203
204			// Set up the final bidiagonal matrix or order n.
205
206			if (nct < n_) {
207			s_[nct] = A[nct][nct];
208			}
209			if (nrt+1 < n_) {
210			e[nrt] = A[nrt][n_-1];
211			}
212			e[n_-1] = 0.0;
213
214			// generate U
215
216			for (j = nct; j < n_; j++) {
217			for (i = 0; i < m_; i++) {
218			U_[i][j] = 0.0;
219			}
220			U_[j][j] = 1.0;
221			}
222			for (k = nct-1; k >= 0; --k) {
223			if (s_[k] != 0.0) {
224			for (j = k+1; j < n_; ++j) {
225			Real t = 0;
226			for (i = k; i < m_; i++) {
227			t += U_[i][k]*U_[i][j];
228			}
229			t = -t/U_[k][k];
230			for (i = k; i < m_; i++) {
231			U_[i][j] += t*U_[i][k];
232			}
233			}
234			for (i = k; i < m_; i++ ) {
235			U_[i][k] = -U_[i][k];
236			}
237			U_[k][k] = 1.0 + U_[k][k];
238			for (i = 0; i < k-1; i++) {
239			U_[i][k] = 0.0;
240			}
241			} else {
242			for (i = 0; i < m_; i++) {
243			U_[i][k] = 0.0;
244			}
245			U_[k][k] = 1.0;
246			}
247			}
248
249			// generate V
250
251			for (k = n_-1; k >= 0; --k) {
252			if ((k < nrt) & (e[k] != 0.0)) {
253			for (j = k+1; j < n_; ++j) {
254			Real t = 0;
255			for (i = k+1; i < n_; i++) {
256			t += V_[i][k]*V_[i][j];
257			}
258			t = -t/V_[k+1][k];
259			for (i = k+1; i < n_; i++) {
260			V_[i][j] += t*V_[i][k];
261			}
262			}
263			}
264			for (i = 0; i < n_; i++) {
265			V_[i][k] = 0.0;
266			}
267			V_[k][k] = 1.0;
268			}
269
270			// Main iteration loop for the singular values.
271
272			Integer p = n_, pp = p-1;
273			Integer iter = 0;
274			Real eps = std::pow(2.0,-52.0);
275			while (p > 0) {
276			Integer k;
277			Integer kase;
278
279			// Here is where a test for too many iterations would go.
280
281			// This section of the program inspects for
282			// negligible elements in the s and e arrays. On
283			// completion the variables kase and k are set as follows.
284
285			// kase = 1 if s(p) and e[k-1] are negligible and k<p
286			// kase = 2 if s(k) is negligible and k<p
287			// kase = 3 if e[k-1] is negligible, k<p, and
288			// s(k), ..., s(p) are not negligible (qr step).
289			// kase = 4 if e(p-1) is negligible (convergence).
290
291			for (k = p-2; k >= -1; --k) {
292			if (k == -1) {
293			break;
294			}
295			if (std::fabs(e[k]) <= eps*(std::fabs(s_[k]) +
296			std::fabs(s_[k+1]))) {
297			e[k] = 0.0;
298			break;
299			}
300			}
301			if (k == p-2) {
302			kase = 4;
303			} else {
304			Integer ks;
305			for (ks = p-1; ks >= k; --ks) {
306			if (ks == k) {
307			break;
308			}
309			Real t = (ks != p ? std::fabs(e[ks]) : 0.) +
310			(ks != k+1 ? std::fabs(e[ks-1]) : 0.);
311			if (std::fabs(s_[ks]) <= eps*t) {
312			s_[ks] = 0.0;
313			break;
314			}
315			}
316			if (ks == k) {
317			kase = 3;
318			} else if (ks == p-1) {
319			kase = 1;
320			} else {
321			kase = 2;
322			k = ks;
323			}
324			}
325			k++;
326
327			// Perform the task indicated by kase.
328
329			switch (kase) {
330
331			// Deflate negligible s(p).
332
333			case 1: {
334			Real f = e[p-2];
335			e[p-2] = 0.0;
336			for (j = p-2; j >= k; --j) {
337			Real t = hypot(s_[j],f);
338			Real cs = s_[j]/t;
339			Real sn = f/t;
340			s_[j] = t;
341			if (j != k) {
342			f = -sn*e[j-1];
343			e[j-1] = cs*e[j-1];
344			}
345			for (i = 0; i < n_; i++) {
346			t = csV_[i][j] + snV_[i][p-1];
347			V_[i][p-1] = -snV_[i][j] + csV_[i][p-1];
348			V_[i][j] = t;
349			}
350			}
351			}
352			break;
353
354			// Split at negligible s(k).
355
356			case 2: {
357			Real f = e[k-1];
358			e[k-1] = 0.0;
359			for (j = k; j < p; j++) {
360			Real t = hypot(s_[j],f);
361			Real cs = s_[j]/t;
362			Real sn = f/t;
363			s_[j] = t;
364			f = -sn*e[j];
365			e[j] = cs*e[j];
366			for (i = 0; i < m_; i++) {
367			t = csU_[i][j] + snU_[i][k-1];
368			U_[i][k-1] = -snU_[i][j] + csU_[i][k-1];
369			U_[i][j] = t;
370			}
371			}
372			}
373			break;
374
375			// Perform one qr step.
376
377			case 3: {
378
379			// Calculate the shift.
380			Real scale = std::max(
381			std::max(
382			std::max(
383			std::max(std::fabs(s_[p-1]),
384			std::fabs(s_[p-2])),
385			std::fabs(e[p-2])),
386			std::fabs(s_[k])),
387			std::fabs(e[k]));
388			Real sp = s_[p-1]/scale;
389			Real spm1 = s_[p-2]/scale;
390			Real epm1 = e[p-2]/scale;
391			Real sk = s_[k]/scale;
392			Real ek = e[k]/scale;
393			Real b = ((spm1 + sp)(spm1 - sp) + epm1epm1)/2.0;
394			Real c = (spepm1)(sp*epm1);
395			Real shift = 0.0;
396			if ((b != 0.0) \| (c != 0.0)) {
397			shift = std::sqrt(b*b + c);
398			if (b < 0.0) {
399			shift = -shift;
400			}
401			shift = c/(b + shift);
402			}
403			Real f = (sk + sp)*(sk - sp) + shift;
404			Real g = sk*ek;
405
406			// Chase zeros.
407
408			for (j = k; j < p-1; j++) {
409			Real t = hypot(f,g);
410			Real cs = f/t;
411			Real sn = g/t;
412			if (j != k) {
413			e[j-1] = t;
414			}
415			f = css_[j] + sne[j];
416			e[j] = cse[j] - sns_[j];
417			g = sn*s_[j+1];
418			s_[j+1] = cs*s_[j+1];
419			for (i = 0; i < n_; i++) {
420			t = csV_[i][j] + snV_[i][j+1];
421			V_[i][j+1] = -snV_[i][j] + csV_[i][j+1];
422			V_[i][j] = t;
423			}
424			t = hypot(f,g);
425			cs = f/t;
426			sn = g/t;
427			s_[j] = t;
428			f = cse[j] + sns_[j+1];
429			s_[j+1] = -sne[j] + css_[j+1];
430			g = sn*e[j+1];
431			e[j+1] = cs*e[j+1];
432			if (j < m_-1) {
433			for (i = 0; i < m_; i++) {
434			t = csU_[i][j] + snU_[i][j+1];
435			U_[i][j+1] = -snU_[i][j] + csU_[i][j+1];
436			U_[i][j] = t;
437			}
438			}
439			}
440			e[p-2] = f;
441			iter = iter + 1;
442			}
443			break;
444
445			// Convergence.
446
447			case 4: {
448
449			// Make the singular values positive.
450
451			if (s_[k] <= 0.0) {
452			s_[k] = (s_[k] < 0.0 ? -s_[k] : 0.0);
453			for (i = 0; i <= pp; i++) {
454			V_[i][k] = -V_[i][k];
455			}
456			}
457
458			// Order the singular values.
459
460			while (k < pp) {
461			if (s_[k] >= s_[k+1]) {
462			break;
463			}
464			swap(s_[k], s_[k+1]);
465			if (k < n_-1) {
466			for (i = 0; i < n_; i++) {
467			swap(V_[i][k], V_[i][k+1]);
468			}
469			}
470			if (k < m_-1) {
471			for (i = 0; i < m_; i++) {
472			swap(U_[i][k], U_[i][k+1]);
473			}
474			}
475			k++;
476			}
477			iter = 0;
478			--p;
479			}
480			break;
481			}
482			}
483			}
484
485			const RectMatrix& SVD::U() const {
486			return (transpose_ ? V_ : U_);
487			}
488
489			const RectMatrix& SVD::V() const {
490			return (transpose_ ? U_ : V_);
491			}
492
493			const Vector& SVD::singularValues() const {
494			return s_;
495			}
496
497			RectMatrix SVD::S() const {
498			Matrix S(n_,n_);
499			for (Size i = 0; i < Size(n_); i++) {
500			for (Size j = 0; j < Size(n_); j++) {
501			S[i][j] = 0.0;
502			}
503			S[i][i] = s_[i];
504			}
505			return S;
506			}
507
508			Real SVD::norm2() {
509			return s_[0];
510			}
511
512			Real SVD::cond() {
513			return s_[0]/s_[n_-1];
514			}
515
516			integer SVD::rank() {
517			Real eps = std::pow(2.0,-52.0);
518			Real tol = m_s_[0]eps;
519			Integer r = 0;
520			for (Size i = 0; i < s_.size(); i++) {
521			if (s_[i] > tol) {
522			r++;
523			}
524			}
525			return r;
526			}
527
528			Vector SVD::solveFor(const Array& b) const{
529			RectMatrix W(n_, n_, 0.0);
530			for (Size i=0; i<Size(n_); i++)
531			W[i][i] = 1./s_[i];
532
533			RectMatrix inverse = V()* W * transpose(U());
534			Vector result = inverse * b;
535			return result;
536			}
537
538			}
539