src/math/SVD.hpp

#ifndef JAMA_SVD_H
#define JAMA_SVD_H

#include "math/DynamicRectMatrix.hpp"

#include <algorithm>
// for min(), max() below
#include <cmath>
// for abs() below

using namespace OpenMD;
using namespace std;

namespace JAMA
{
  /** Singular Value Decomposition.
      <P>
      For an m-by-n matrix A with m >= n, the singular value decomposition is
      an m-by-n orthogonal matrix U, an n-by-n diagonal matrix S, and
      an n-by-n orthogonal matrix V so that A = U*S*V'.  
      <P>
      The singular values, sigma(k) = S(k,k), are ordered so that
      sigma(0) >= sigma(1) >= ... >= sigma(n-1).
      <P>
      The singular value decompostion always exists, so the constructor will
      never fail.  The matrix condition number and the effective numerical
      rank can be computed from this decomposition.

      <p>
      (Adapted from JAMA, a Java Matrix Library, developed by jointly 
      by the Mathworks and NIST; see  http://math.nist.gov/javanumerics/jama).
  */
  template <class Real>
  class SVD 
  {
        
    DynamicRectMatrix<Real> U, V;
    DynamicVector<Real> s;
    int m, n;
    
  public:
    
    
    SVD (const DynamicRectMatrix<Real> &Arg) {
      m = Arg.getNRow();
      n = Arg.getNCol();
      int nu = min(m,n);
      s = DynamicVector<Real>(min(m+1,n)); 
      U = DynamicRectMatrix<Real>(m, nu, Real(0));
      V = DynamicRectMatrix<Real>(n,n);
      DynamicVector<Real> e(n);
      DynamicVector<Real> work(m);
      DynamicRectMatrix<Real> A(Arg);

      int wantu = 1;    /* boolean */
      int wantv = 1;    /* boolean */
      int i=0, j=0, k=0;
      
      // Reduce A to bidiagonal form, storing the diagonal elements
      // in s and the super-diagonal elements in e.
      
      int nct = min(m-1,n);
      int nrt = max(0,min(n-2,m));

      for (k = 0; k < 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.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 (wantu & (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);
              }
            }
          }
          if (wantv) {

            // 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 p.

      int p = min(n,m+1);
      if (nct < n) {
        s(nct) = A(nct,nct);
      }
      if (m < p) {
        s(p-1) = 0.0;
      }
      if (nrt+1 < p) {
        e(nrt) = A(nrt,p-1);
      }
      e(p-1) = 0.0;

      // If required, generate U.

      if (wantu) {
        for (j = nct; j < nu; 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 < nu; j++) {
              Real t(0.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;
          }
        }
      }

      // If required, generate V.

      if (wantv) {
        for (k = n-1; k >= 0; k--) {
          if ((k < nrt) & (e(k) != 0.0)) {
            for (j = k+1; j < nu; j++) {
              Real t(0.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.

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

        // 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 (abs(e(k)) <= eps*(abs(s(k)) + abs(s(k+1)))) {
            e(k) = 0.0;
            break;
          }
        }
        if (k == p-2) {
          kase = 4;
        } else {
          int ks;
          for (ks = p-1; ks >= k; ks--) {
            if (ks == k) {
              break;
            }
            Real t( (ks != p ? abs(e(ks)) : 0.) + 
                    (ks != k+1 ? abs(e(ks-1)) : 0.));
            if (abs(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);
            }
            if (wantv) {
              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);
            if (wantu) {
              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 = max(max(max(max(
                                       abs(s(p-1)),abs(s(p-2))),abs(e(p-2))), 
                               abs(s(k))),abs(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 = 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);
            if (wantv) {
              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 (wantu && (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);
            if (wantv) {
              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;
            }
            Real t = s(k);
            s(k) = s(k+1);
            s(k+1) = t;
            if (wantv && (k < n-1)) {
              for (i = 0; i < n; i++) {
                t = V(i,k+1); V(i,k+1) = V(i,k); V(i,k) = t;
              }
            }
            if (wantu && (k < m-1)) {
              for (i = 0; i < m; i++) {
                t = U(i,k+1); U(i,k+1) = U(i,k); U(i,k) = t;
              }
            }
            k++;
          }
          iter = 0;
          p--;
        }
          break;
        }
      }
    }


    void getU (DynamicRectMatrix<Real> &A) {

      int minm = min(m+1,n);
      
      A = DynamicRectMatrix<Real>(m, minm);
      
      for (int i=0; i<m; i++)
        for (int j=0; j<minm; j++)
          A(i,j) = U(i,j);      
    }

    /* Return the right singular vectors */
    void getV (DynamicRectMatrix<Real> &A) {
      A = V;
    }

    /** Return the one-dimensional array of singular values */
    void getSingularValues (DynamicVector<Real> &x) {
      x = s;
    }

    /** Return the diagonal matrix of singular values
        @return     S
    */
    void getS (DynamicRectMatrix<Real> &A) {
      A = DynamicRectMatrix<Real>(n,n);
      for (int i = 0; i < n; i++) {
        for (int j = 0; j < n; j++) {
          A(i,j) = 0.0;
        }
        A(i,i) = s(i);
      }
    }
    
    /** Two norm  (max(S)) */    
    Real norm2 () {
      return s(0);
    }
    
    /** Two norm of condition number (max(S)/min(S)) */
    Real cond () {
      return s(0)/s(min(m,n)-1);
    }
    
    /** Effective numerical matrix rank
        @return     Number of nonnegligible singular values.
    */
    int rank () {
      Real eps = pow(2.0,-52.0);
      Real tol = max(m,n)*s(0)*eps;
      int r = 0;
      for (int i = 0; i < s.dim(); i++) {
        if (s(i) > tol) {
          r++;
        }
      }
      return r;
    }
  };  
}
#endif
// JAMA_SVD_H
Revision:	1390
Committed:	Wed Nov 25 20:02:06 2009 UTC (15 years, 11 months ago) by gezelter
File size:	10807 byte(s)
Log Message:	Almost all of the changes necessary to create OpenMD out of our old project (OOPSE-4)
#	User	Rev	Content
1	cli2	1349	#ifndef JAMA_SVD_H
2			#define JAMA_SVD_H
3
4			#include "math/DynamicRectMatrix.hpp"
5
6			#include <algorithm>
7			// for min(), max() below
8			#include <cmath>
9			// for abs() below
10
11	gezelter	1390	using namespace OpenMD;
12	cli2	1349	using namespace std;
13
14			namespace JAMA
15			{
16			/** Singular Value Decomposition.
17			<P>
18			For an m-by-n matrix A with m >= n, the singular value decomposition is
19			an m-by-n orthogonal matrix U, an n-by-n diagonal matrix S, and
20			an n-by-n orthogonal matrix V so that A = USV'.
21			<P>
22			The singular values, sigma(k) = S(k,k), are ordered so that
23			sigma(0) >= sigma(1) >= ... >= sigma(n-1).
24			<P>
25			The singular value decompostion always exists, so the constructor will
26			never fail. The matrix condition number and the effective numerical
27			rank can be computed from this decomposition.
28
29			<p>
30			(Adapted from JAMA, a Java Matrix Library, developed by jointly
31			by the Mathworks and NIST; see http://math.nist.gov/javanumerics/jama).
32			*/
33			template <class Real>
34			class SVD
35			{
36
37			DynamicRectMatrix<Real> U, V;
38			DynamicVector<Real> s;
39			int m, n;
40
41			public:
42
43
44			SVD (const DynamicRectMatrix<Real> &Arg) {
45			m = Arg.getNRow();
46			n = Arg.getNCol();
47			int nu = min(m,n);
48			s = DynamicVector<Real>(min(m+1,n));
49			U = DynamicRectMatrix<Real>(m, nu, Real(0));
50			V = DynamicRectMatrix<Real>(n,n);
51			DynamicVector<Real> e(n);
52			DynamicVector<Real> work(m);
53			DynamicRectMatrix<Real> A(Arg);
54	cli2	1360
55			int wantu = 1; /* boolean */
56			int wantv = 1; /* boolean */
57	cli2	1349	int i=0, j=0, k=0;
58
59			// Reduce A to bidiagonal form, storing the diagonal elements
60			// in s and the super-diagonal elements in e.
61
62			int nct = min(m-1,n);
63			int nrt = max(0,min(n-2,m));
64	cli2	1360
65	cli2	1349	for (k = 0; k < max(nct,nrt); k++) {
66			if (k < nct) {
67
68			// Compute the transformation for the k-th column and
69			// place the k-th diagonal in s(k).
70			// Compute 2-norm of k-th column without under/overflow.
71			s(k) = 0;
72			for (i = k; i < m; i++) {
73			s(k) = hypot(s(k),A(i,k));
74			}
75			if (s(k) != 0.0) {
76			if (A(k,k) < 0.0) {
77			s(k) = -s(k);
78			}
79			for (i = k; i < m; i++) {
80			A(i,k) /= s(k);
81			}
82			A(k,k) += 1.0;
83			}
84			s(k) = -s(k);
85			}
86			for (j = k+1; j < n; j++) {
87			if ((k < nct) && (s(k) != 0.0)) {
88
89			// Apply the transformation.
90
91			Real t(0.0);
92			for (i = k; i < m; i++) {
93			t += A(i,k)*A(i,j);
94			}
95			t = -t/A(k,k);
96			for (i = k; i < m; i++) {
97			A(i,j) += t*A(i,k);
98			}
99			}
100
101			// Place the k-th row of A into e for the
102			// subsequent calculation of the row transformation.
103
104			e(j) = A(k,j);
105			}
106			if (wantu & (k < nct)) {
107
108			// Place the transformation in U for subsequent back
109			// multiplication.
110
111			for (i = k; i < m; i++) {
112			U(i,k) = A(i,k);
113			}
114			}
115			if (k < nrt) {
116
117			// Compute the k-th row transformation and place the
118			// k-th super-diagonal in e(k).
119			// Compute 2-norm without under/overflow.
120			e(k) = 0;
121			for (i = k+1; i < n; i++) {
122			e(k) = hypot(e(k),e(i));
123			}
124			if (e(k) != 0.0) {
125			if (e(k+1) < 0.0) {
126			e(k) = -e(k);
127			}
128			for (i = k+1; i < n; i++) {
129			e(i) /= e(k);
130			}
131			e(k+1) += 1.0;
132			}
133			e(k) = -e(k);
134			if ((k+1 < m) & (e(k) != 0.0)) {
135
136			// Apply the transformation.
137
138			for (i = k+1; i < m; i++) {
139			work(i) = 0.0;
140			}
141			for (j = k+1; j < n; j++) {
142			for (i = k+1; i < m; i++) {
143			work(i) += e(j)*A(i,j);
144			}
145			}
146			for (j = k+1; j < n; j++) {
147			Real t(-e(j)/e(k+1));
148			for (i = k+1; i < m; i++) {
149			A(i,j) += t*work(i);
150			}
151			}
152			}
153			if (wantv) {
154
155			// Place the transformation in V for subsequent
156			// back multiplication.
157
158			for (i = k+1; i < n; i++) {
159			V(i,k) = e(i);
160			}
161			}
162			}
163			}
164
165			// Set up the final bidiagonal matrix or order p.
166
167			int p = min(n,m+1);
168			if (nct < n) {
169			s(nct) = A(nct,nct);
170			}
171			if (m < p) {
172			s(p-1) = 0.0;
173			}
174			if (nrt+1 < p) {
175			e(nrt) = A(nrt,p-1);
176			}
177			e(p-1) = 0.0;
178
179			// If required, generate U.
180
181			if (wantu) {
182			for (j = nct; j < nu; j++) {
183			for (i = 0; i < m; i++) {
184			U(i,j) = 0.0;
185			}
186			U(j,j) = 1.0;
187			}
188			for (k = nct-1; k >= 0; k--) {
189			if (s(k) != 0.0) {
190			for (j = k+1; j < nu; j++) {
191			Real t(0.0);
192			for (i = k; i < m; i++) {
193			t += U(i,k)*U(i,j);
194			}
195			t = -t/U(k,k);
196			for (i = k; i < m; i++) {
197			U(i,j) += t*U(i,k);
198			}
199			}
200			for (i = k; i < m; i++ ) {
201			U(i,k) = -U(i,k);
202			}
203			U(k,k) = 1.0 + U(k,k);
204			for (i = 0; i < k-1; i++) {
205			U(i,k) = 0.0;
206			}
207			} else {
208			for (i = 0; i < m; i++) {
209			U(i,k) = 0.0;
210			}
211			U(k,k) = 1.0;
212			}
213			}
214			}
215
216			// If required, generate V.
217
218			if (wantv) {
219			for (k = n-1; k >= 0; k--) {
220			if ((k < nrt) & (e(k) != 0.0)) {
221			for (j = k+1; j < nu; j++) {
222			Real t(0.0);
223			for (i = k+1; i < n; i++) {
224			t += V(i,k)*V(i,j);
225			}
226			t = -t/V(k+1,k);
227			for (i = k+1; i < n; i++) {
228			V(i,j) += t*V(i,k);
229			}
230			}
231			}
232			for (i = 0; i < n; i++) {
233			V(i,k) = 0.0;
234			}
235			V(k,k) = 1.0;
236			}
237			}
238
239			// Main iteration loop for the singular values.
240
241			int pp = p-1;
242			int iter = 0;
243			Real eps(pow(2.0,-52.0));
244			while (p > 0) {
245			int k=0;
246			int kase=0;
247
248			// Here is where a test for too many iterations would go.
249
250			// This section of the program inspects for
251			// negligible elements in the s and e arrays. On
252			// completion the variables kase and k are set as follows.
253
254			// kase = 1 if s(p) and e(k-1) are negligible and k<p
255			// kase = 2 if s(k) is negligible and k<p
256			// kase = 3 if e(k-1) is negligible, k<p, and
257			// s(k), ..., s(p) are not negligible (qr step).
258			// kase = 4 if e(p-1) is negligible (convergence).
259
260			for (k = p-2; k >= -1; k--) {
261			if (k == -1) {
262			break;
263			}
264			if (abs(e(k)) <= eps*(abs(s(k)) + abs(s(k+1)))) {
265			e(k) = 0.0;
266			break;
267			}
268			}
269			if (k == p-2) {
270			kase = 4;
271			} else {
272			int ks;
273			for (ks = p-1; ks >= k; ks--) {
274			if (ks == k) {
275			break;
276			}
277			Real t( (ks != p ? abs(e(ks)) : 0.) +
278			(ks != k+1 ? abs(e(ks-1)) : 0.));
279			if (abs(s(ks)) <= eps*t) {
280			s(ks) = 0.0;
281			break;
282			}
283			}
284			if (ks == k) {
285			kase = 3;
286			} else if (ks == p-1) {
287			kase = 1;
288			} else {
289			kase = 2;
290			k = ks;
291			}
292			}
293			k++;
294
295			// Perform the task indicated by kase.
296
297			switch (kase) {
298
299			// Deflate negligible s(p).
300
301			case 1: {
302			Real f(e(p-2));
303			e(p-2) = 0.0;
304			for (j = p-2; j >= k; j--) {
305			Real t( hypot(s(j),f));
306			Real cs(s(j)/t);
307			Real sn(f/t);
308			s(j) = t;
309			if (j != k) {
310			f = -sn*e(j-1);
311			e(j-1) = cs*e(j-1);
312			}
313			if (wantv) {
314			for (i = 0; i < n; i++) {
315			t = csV(i,j) + snV(i,p-1);
316			V(i,p-1) = -snV(i,j) + csV(i,p-1);
317			V(i,j) = t;
318			}
319			}
320			}
321			}
322			break;
323
324			// Split at negligible s(k).
325
326			case 2: {
327			Real f(e(k-1));
328			e(k-1) = 0.0;
329			for (j = k; j < p; j++) {
330			Real t(hypot(s(j),f));
331			Real cs( s(j)/t);
332			Real sn(f/t);
333			s(j) = t;
334			f = -sn*e(j);
335			e(j) = cs*e(j);
336			if (wantu) {
337			for (i = 0; i < m; i++) {
338			t = csU(i,j) + snU(i,k-1);
339			U(i,k-1) = -snU(i,j) + csU(i,k-1);
340			U(i,j) = t;
341			}
342			}
343			}
344			}
345			break;
346
347			// Perform one qr step.
348
349			case 3: {
350
351			// Calculate the shift.
352
353			Real scale = max(max(max(max(
354			abs(s(p-1)),abs(s(p-2))),abs(e(p-2))),
355			abs(s(k))),abs(e(k)));
356			Real sp = s(p-1)/scale;
357			Real spm1 = s(p-2)/scale;
358			Real epm1 = e(p-2)/scale;
359			Real sk = s(k)/scale;
360			Real ek = e(k)/scale;
361			Real b = ((spm1 + sp)(spm1 - sp) + epm1epm1)/2.0;
362			Real c = (spepm1)(sp*epm1);
363			Real shift = 0.0;
364			if ((b != 0.0) \|\| (c != 0.0)) {
365			shift = sqrt(b*b + c);
366			if (b < 0.0) {
367			shift = -shift;
368			}
369			shift = c/(b + shift);
370			}
371			Real f = (sk + sp)*(sk - sp) + shift;
372			Real g = sk*ek;
373
374			// Chase zeros.
375
376			for (j = k; j < p-1; j++) {
377			Real t = hypot(f,g);
378			Real cs = f/t;
379			Real sn = g/t;
380			if (j != k) {
381			e(j-1) = t;
382			}
383			f = css(j) + sne(j);
384			e(j) = cse(j) - sns(j);
385			g = sn*s(j+1);
386			s(j+1) = cs*s(j+1);
387			if (wantv) {
388			for (i = 0; i < n; i++) {
389			t = csV(i,j) + snV(i,j+1);
390			V(i,j+1) = -snV(i,j) + csV(i,j+1);
391			V(i,j) = t;
392			}
393			}
394			t = hypot(f,g);
395			cs = f/t;
396			sn = g/t;
397			s(j) = t;
398			f = cse(j) + sns(j+1);
399			s(j+1) = -sne(j) + css(j+1);
400			g = sn*e(j+1);
401			e(j+1) = cs*e(j+1);
402			if (wantu && (j < m-1)) {
403			for (i = 0; i < m; i++) {
404			t = csU(i,j) + snU(i,j+1);
405			U(i,j+1) = -snU(i,j) + csU(i,j+1);
406			U(i,j) = t;
407			}
408			}
409			}
410			e(p-2) = f;
411			iter = iter + 1;
412			}
413			break;
414
415			// Convergence.
416
417			case 4: {
418
419			// Make the singular values positive.
420
421			if (s(k) <= 0.0) {
422			s(k) = (s(k) < 0.0 ? -s(k) : 0.0);
423			if (wantv) {
424			for (i = 0; i <= pp; i++) {
425			V(i,k) = -V(i,k);
426			}
427			}
428			}
429
430			// Order the singular values.
431
432			while (k < pp) {
433			if (s(k) >= s(k+1)) {
434			break;
435			}
436			Real t = s(k);
437			s(k) = s(k+1);
438			s(k+1) = t;
439			if (wantv && (k < n-1)) {
440			for (i = 0; i < n; i++) {
441			t = V(i,k+1); V(i,k+1) = V(i,k); V(i,k) = t;
442			}
443			}
444			if (wantu && (k < m-1)) {
445			for (i = 0; i < m; i++) {
446			t = U(i,k+1); U(i,k+1) = U(i,k); U(i,k) = t;
447			}
448			}
449			k++;
450			}
451			iter = 0;
452			p--;
453			}
454			break;
455			}
456			}
457			}
458
459
460			void getU (DynamicRectMatrix<Real> &A) {
461
462			int minm = min(m+1,n);
463
464			A = DynamicRectMatrix<Real>(m, minm);
465
466			for (int i=0; i<m; i++)
467			for (int j=0; j<minm; j++)
468			A(i,j) = U(i,j);
469			}
470
471			/* Return the right singular vectors */
472			void getV (DynamicRectMatrix<Real> &A) {
473			A = V;
474			}
475
476			/** Return the one-dimensional array of singular values */
477			void getSingularValues (DynamicVector<Real> &x) {
478			x = s;
479			}
480
481			/** Return the diagonal matrix of singular values
482			@return S
483			*/
484			void getS (DynamicRectMatrix<Real> &A) {
485			A = DynamicRectMatrix<Real>(n,n);
486			for (int i = 0; i < n; i++) {
487			for (int j = 0; j < n; j++) {
488			A(i,j) = 0.0;
489			}
490			A(i,i) = s(i);
491			}
492			}
493
494			/** Two norm (max(S)) */
495			Real norm2 () {
496			return s(0);
497			}
498
499			/** Two norm of condition number (max(S)/min(S)) */
500			Real cond () {
501			return s(0)/s(min(m,n)-1);
502			}
503
504			/** Effective numerical matrix rank
505			@return Number of nonnegligible singular values.
506			*/
507			int rank () {
508			Real eps = pow(2.0,-52.0);
509			Real tol = max(m,n)s(0)eps;
510			int r = 0;
511			for (int i = 0; i < s.dim(); i++) {
512			if (s(i) > tol) {
513			r++;
514			}
515			}
516			return r;
517			}
518			};
519			}
520			#endif
521			// JAMA_SVD_H