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:	3494
Committed:	Fri Apr 3 17:44:57 2009 UTC (16 years 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.
#	Content
1	/* -- 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