Merge pull request #2626 from mark-mb/jama

Update JAMA to version 1.0.3

libraries/jama/include/jama_cholesky.h

@@ -191,6 +191,9 @@ Array2D<Real> Cholesky<Real>::solve(const Array2D<Real> &B)
 	if (B.dim1() != n)
 		return Array2D<Real>();
 
+	if (!isspd)
+	    return Arran2D<Real>();
+
 
 	Array2D<Real> X = B.copy();
 	int nx = B.dim2();
@@ -224,9 +227,9 @@ Array2D<Real> Cholesky<Real>::solve(const Array2D<Real> &B)
 
 
       // Solve L*y = b;
-  	  for (int j=0; j< nx; j++)
+      for (int k = 0; k < n; k++)
 	  {
-      	for (int k = 0; k < n; k++) 
+  	    for (int j=0; j< nx; j++)
 		{
 			for (int i = 0; i < k; i++) 
                X[k][j] -= X[i][j]*L_[k][i];
@@ -235,9 +238,9 @@ Array2D<Real> Cholesky<Real>::solve(const Array2D<Real> &B)
       }
 
       // Solve L'*X = Y;
-     for (int j=0; j<nx; j++)
+     for (int k = n-1; k >= 0; k--)
 	 {
-      	for (int k = n-1; k >= 0; k--) 
+        for (int j=0; j<nx; j++)
 	  	{
          	for (int i = k+1; i < n; i++) 
                X[k][j] -= X[i][j]*L_[i][k];

libraries/jama/include/jama_eig.h

@@ -663,21 +663,20 @@ class Eigenvalue
             // Double QR step involving rows l:n and columns m:n
    
             for (int k = m; k <= n-1; k++) {
-               int notlast = (k != n-1);
+               bool notlast = (k != n-1);
                if (k != m) {
                   p = H[k][k-1];
                   q = H[k+1][k-1];
                   r = (notlast ? H[k+2][k-1] : 0.0);
                   x = abs(p) + abs(q) + abs(r);
-                  if (x != 0.0) {
-                     p = p / x;
-                     q = q / x;
-                     r = r / x;
+                  if (x == 0.0) {
+                     continue;
                   }
+                  p = p / x;
+                  q = q / x;
+                  r = r / x;
                }
-               if (x == 0.0) {
-                  break;
-               }
+
                s = sqrt(p * p + q * q + r * r);
                if (p < 0) {
                   s = -s;

libraries/jama/include/jama_svd.h

@@ -12,6 +12,8 @@
 // for min(), max() below
 #include <cmath>
 // for abs() below
+#include <limits>
+// for numeric_limits
 
 using namespace TNT;
 using namespace std;
@@ -247,6 +249,7 @@ class SVD
       int pp = p-1;
       int iter = 0;
       Real eps(pow(2.0,-52.0));
+      Real tiny = numeric_limits<Real>::min();
       while (p > 0) {
          int k=0;
 		 int kase=0;
@@ -267,7 +270,7 @@ class SVD
             if (k == -1) {
                break;
             }
-            if (abs(e[k]) <= eps*(abs(s[k]) + abs(s[k+1]))) {
+            if (abs(e[k]) <= tiny + eps*(abs(s[k]) + abs(s[k+1]))) {
                e[k] = 0.0;
                break;
             }
@@ -282,7 +285,7 @@ class SVD
                }
                Real t( (ks != p ? abs(e[ks]) : 0.) + 
                           (ks != k+1 ? abs(e[ks-1]) : 0.));
-               if (abs(s[ks]) <= eps*t)  {
+               if (abs(s[ks]) <= tiny + eps*t)  {
                   s[ks] = 0.0;
                   break;
                }