Update JAMA to version 1.0.3

Adapted from https://math.nist.gov/javanumerics/jama/

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 j=0; j< nx; j++)
 		{
 			for (int i = 0; i < k; i++) 
                X[k][j] -= X[i][j]*L_[k][i];
@@ -235,10 +238,10 @@ 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 j=0; j<nx; j++)
+	  	{
          	for (int i = k+1; i < n; i++) 
                X[k][j] -= X[i][j]*L_[i][k];
          	X[k][j] /= L_[k][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) {
+                  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;
                }