src/math/jama_svd.hpp

#ifndef JAMA_SVD_H
#define JAMA_SVD_H


#include "tnt_array1d.hpp"
#include "tnt_array1d_utils.hpp"
#include "tnt_array2d.hpp"
#include "tnt_array2d_utils.hpp"
#include "tnt_math_utils.hpp"

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

using namespace TNT;
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 
{


        Array2D<Real> U, V;
        Array1D<Real> s;
        int m, n;

  public:


   SVD (const Array2D<Real> &Arg) {


      m = Arg.dim1();
      n = Arg.dim2();
      int nu = min(m,n);
      s = Array1D<Real>(min(m+1,n)); 
      U = Array2D<Real>(m, nu, Real(0));
      V = Array2D<Real>(n,n);
      Array1D<Real> e(n);
      Array1D<Real> work(m);
          Array2D<Real> A(Arg.copy());
      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 (Array2D<Real> &A) 
   {
          int minm = min(m+1,n);

          A = Array2D<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 (Array2D<Real> &A) 
   {
          A = V;
   }

   /** Return the one-dimensional array of singular values */

   void getSingularValues (Array1D<Real> &x) 
   {
      x = s;
   }

   /** Return the diagonal matrix of singular values
   @return     S
   */

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