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, 1 month ago) by gezelter
File size:	14486 byte(s)
Log Message:	adding jama and tnt libraries for various linear algebra routines
#	Content
1	#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