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:	3520
Committed:	Mon Sep 7 16:31:51 2009 UTC (14 years, 9 months ago) by cli2
File size:	10806 byte(s)
Log Message:	Added new restraint infrastructure Added MolecularRestraints Added ObjectRestraints Added RestraintStamp Updated thermodynamic integration to use ObjectRestraints Added Quaternion mathematics for twist swing decompositions Significantly updated RestWriter and RestReader to use dump-like files Added selections for x, y, and z coordinates of atoms Removed monolithic Restraints class Fixed a few bugs in gradients of Euler angles in DirectionalAtom and RigidBody Added some rotational capabilities to prinicpalAxisCalculator
#	Content
1	#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	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
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
65	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