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