Typical usage looks like:
! Array2DA(n,n); ! Array2D L; ... ! Cholesky chol(A); if (chol.is_spd()) L = chol.getL(); --- 21,32 ---- Typical usage looks like:
! Array2DA(n,n); ! Array2D L; ... ! Cholesky chol(A); if (chol.is_spd()) L = chol.getL(); *************** *** 110,119 **** // Main loop. for (int j = 0; j < n; j++) { ! double d = 0.0; for (int k = 0; k < j; k++) { ! Real s = 0.0; for (int i = 0; i < k; i++) { s += L_[k][i]*L_[j][i]; --- 110,119 ---- // Main loop. for (int j = 0; j < n; j++) { ! Real d = 0; for (int k = 0; k < j; k++) { ! Real s = 0; for (int i = 0; i < k; i++) { s += L_[k][i]*L_[j][i]; *************** *** 123,133 **** isspd = isspd && (A[k][j] == A[j][k]); } d = A[j][j] - d; ! isspd = isspd && (d > 0.0); ! L_[j][j] = sqrt(d > 0.0 ? d : 0.0); for (int k = j+1; k < n; k++) { ! L_[j][k] = 0.0; } } } --- 123,133 ---- isspd = isspd && (A[k][j] == A[j][k]); } d = A[j][j] - d; ! isspd = isspd && (d > 0); ! L_[j][j] = sqrt(d > 0 ? d : 0); for (int k = j+1; k < n; k++) { ! L_[j][k] = 0; } } } diff -ucr jama.orig/jama_eig.h jama/jama_eig.h *** jama.orig/jama_eig.h Sat Oct 16 06:07:11 2004 --- jama/jama_eig.h Tue May 10 11:26:54 2005 *************** *** 117,133 **** // Scale to avoid under/overflow. ! Real scale = 0.0; ! Real h = 0.0; for (int k = 0; k < i; k++) { scale = scale + abs(d[k]); } ! if (scale == 0.0) { e[i] = d[i-1]; for (int j = 0; j < i; j++) { d[j] = V[i-1][j]; ! V[i][j] = 0.0; ! V[j][i] = 0.0; } } else { --- 117,133 ---- // Scale to avoid under/overflow. ! Real scale = 0; ! Real h = 0; for (int k = 0; k < i; k++) { scale = scale + abs(d[k]); } ! if (scale == 0) { e[i] = d[i-1]; for (int j = 0; j < i; j++) { d[j] = V[i-1][j]; ! V[i][j] = 0; ! V[j][i] = 0; } } else { *************** *** 146,152 **** h = h - f * g; d[i-1] = f - g; for (int j = 0; j < i; j++) { ! e[j] = 0.0; } // Apply similarity transformation to remaining columns. --- 146,152 ---- h = h - f * g; d[i-1] = f - g; for (int j = 0; j < i; j++) { ! e[j] = 0; } // Apply similarity transformation to remaining columns. *************** *** 161,167 **** } e[j] = g; } ! f = 0.0; for (int j = 0; j < i; j++) { e[j] /= h; f += e[j] * d[j]; --- 161,167 ---- } e[j] = g; } ! f = 0; for (int j = 0; j < i; j++) { e[j] /= h; f += e[j] * d[j]; *************** *** 177,183 **** V[k][j] -= (f * e[k] + g * d[k]); } d[j] = V[i-1][j]; ! V[i][j] = 0.0; } } d[i] = h; --- 177,183 ---- V[k][j] -= (f * e[k] + g * d[k]); } d[j] = V[i-1][j]; ! V[i][j] = 0; } } d[i] = h; *************** *** 187,200 **** for (int i = 0; i < n-1; i++) { V[n-1][i] = V[i][i]; ! V[i][i] = 1.0; Real h = d[i+1]; ! if (h != 0.0) { for (int k = 0; k <= i; k++) { d[k] = V[k][i+1] / h; } for (int j = 0; j <= i; j++) { ! Real g = 0.0; for (int k = 0; k <= i; k++) { g += V[k][i+1] * V[k][j]; } --- 187,200 ---- for (int i = 0; i < n-1; i++) { V[n-1][i] = V[i][i]; ! V[i][i] = 1; Real h = d[i+1]; ! if (h != 0) { for (int k = 0; k <= i; k++) { d[k] = V[k][i+1] / h; } for (int j = 0; j <= i; j++) { ! Real g = 0; for (int k = 0; k <= i; k++) { g += V[k][i+1] * V[k][j]; } *************** *** 204,218 **** } } for (int k = 0; k <= i; k++) { ! V[k][i+1] = 0.0; } } for (int j = 0; j < n; j++) { d[j] = V[n-1][j]; ! V[n-1][j] = 0.0; } ! V[n-1][n-1] = 1.0; ! e[0] = 0.0; } // Symmetric tridiagonal QL algorithm. --- 204,218 ---- } } for (int k = 0; k <= i; k++) { ! V[k][i+1] = 0; } } for (int j = 0; j < n; j++) { d[j] = V[n-1][j]; ! V[n-1][j] = 0; } ! V[n-1][n-1] = 1; ! e[0] = 0; } // Symmetric tridiagonal QL algorithm. *************** *** 227,237 **** for (int i = 1; i < n; i++) { e[i-1] = e[i]; } ! e[n-1] = 0.0; ! Real f = 0.0; ! Real tst1 = 0.0; ! Real eps = pow(2.0,-52.0); for (int l = 0; l < n; l++) { // Find small subdiagonal element --- 227,237 ---- for (int i = 1; i < n; i++) { e[i-1] = e[i]; } ! e[n-1] = 0; ! Real f = 0; ! Real tst1 = 0; ! Real eps = pow(2,-52); for (int l = 0; l < n; l++) { // Find small subdiagonal element *************** *** 259,266 **** // Compute implicit shift Real g = d[l]; ! Real p = (d[l+1] - g) / (2.0 * e[l]); ! Real r = hypot(p,1.0); if (p < 0) { r = -r; } --- 259,266 ---- // Compute implicit shift Real g = d[l]; ! Real p = (d[l+1] - g) / (2 * e[l]); ! Real r = hypot(p,1); if (p < 0) { r = -r; } *************** *** 276,287 **** // Implicit QL transformation. p = d[m]; ! Real c = 1.0; Real c2 = c; Real c3 = c; Real el1 = e[l+1]; ! Real s = 0.0; ! Real s2 = 0.0; for (int i = m-1; i >= l; i--) { c3 = c2; c2 = c; --- 276,287 ---- // Implicit QL transformation. p = d[m]; ! Real c = 1; Real c2 = c; Real c3 = c; Real el1 = e[l+1]; ! Real s = 0; ! Real s2 = 0; for (int i = m-1; i >= l; i--) { c3 = c2; c2 = c; *************** *** 312,318 **** } while (abs(e[l]) > eps*tst1); } d[l] = d[l] + f; ! e[l] = 0.0; } // Sort eigenvalues and corresponding vectors. --- 312,318 ---- } while (abs(e[l]) > eps*tst1); } d[l] = d[l] + f; ! e[l] = 0; } // Sort eigenvalues and corresponding vectors. *************** *** 354,368 **** // Scale column. ! Real scale = 0.0; for (int i = m; i <= high; i++) { scale = scale + abs(H[i][m-1]); } ! if (scale != 0.0) { // Compute Householder transformation. ! Real h = 0.0; for (int i = high; i >= m; i--) { ort[i] = H[i][m-1]/scale; h += ort[i] * ort[i]; --- 354,368 ---- // Scale column. ! Real scale = 0; for (int i = m; i <= high; i++) { scale = scale + abs(H[i][m-1]); } ! if (scale != 0) { // Compute Householder transformation. ! Real h = 0; for (int i = high; i >= m; i--) { ort[i] = H[i][m-1]/scale; h += ort[i] * ort[i]; *************** *** 378,384 **** // H = (I-u*u'/h)*H*(I-u*u')/h) for (int j = m; j < n; j++) { ! Real f = 0.0; for (int i = high; i >= m; i--) { f += ort[i]*H[i][j]; } --- 378,384 ---- // H = (I-u*u'/h)*H*(I-u*u')/h) for (int j = m; j < n; j++) { ! Real f = 0; for (int i = high; i >= m; i--) { f += ort[i]*H[i][j]; } *************** *** 389,395 **** } for (int i = 0; i <= high; i++) { ! Real f = 0.0; for (int j = high; j >= m; j--) { f += ort[j]*H[i][j]; } --- 389,395 ---- } for (int i = 0; i <= high; i++) { ! Real f = 0; for (int j = high; j >= m; j--) { f += ort[j]*H[i][j]; } *************** *** 407,423 **** for (int i = 0; i < n; i++) { for (int j = 0; j < n; j++) { ! V[i][j] = (i == j ? 1.0 : 0.0); } } for (int m = high-1; m >= low+1; m--) { ! if (H[m][m-1] != 0.0) { for (int i = m+1; i <= high; i++) { ort[i] = H[i][m-1]; } for (int j = m; j <= high; j++) { ! Real g = 0.0; for (int i = m; i <= high; i++) { g += ort[i] * V[i][j]; } --- 407,423 ---- for (int i = 0; i < n; i++) { for (int j = 0; j < n; j++) { ! V[i][j] = (i == j ? 1 : 0); } } for (int m = high-1; m >= low+1; m--) { ! if (H[m][m-1] != 0) { for (int i = m+1; i <= high; i++) { ort[i] = H[i][m-1]; } for (int j = m; j <= high; j++) { ! Real g = 0; for (int i = m; i <= high; i++) { g += ort[i] * V[i][j]; } *************** *** 466,482 **** int n = nn-1; int low = 0; int high = nn-1; ! Real eps = pow(2.0,-52.0); ! Real exshift = 0.0; Real p=0,q=0,r=0,s=0,z=0,t,w,x,y; // Store roots isolated by balanc and compute matrix norm ! Real norm = 0.0; for (int i = 0; i < nn; i++) { if ((i < low) || (i > high)) { d[i] = H[i][i]; ! e[i] = 0.0; } for (int j = max(i-1,0); j < nn; j++) { norm = norm + abs(H[i][j]); --- 466,482 ---- int n = nn-1; int low = 0; int high = nn-1; ! Real eps = pow(2,-52); ! Real exshift = 0; Real p=0,q=0,r=0,s=0,z=0,t,w,x,y; // Store roots isolated by balanc and compute matrix norm ! Real norm = 0; for (int i = 0; i < nn; i++) { if ((i < low) || (i > high)) { d[i] = H[i][i]; ! e[i] = 0; } for (int j = max(i-1,0); j < nn; j++) { norm = norm + abs(H[i][j]); *************** *** 493,499 **** int l = n; while (l > low) { s = abs(H[l-1][l-1]) + abs(H[l][l]); ! if (s == 0.0) { s = norm; } if (abs(H[l][l-1]) < eps * s) { --- 493,499 ---- int l = n; while (l > low) { s = abs(H[l-1][l-1]) + abs(H[l][l]); ! if (s == 0) { s = norm; } if (abs(H[l][l-1]) < eps * s) { *************** *** 508,514 **** if (l == n) { H[n][n] = H[n][n] + exshift; d[n] = H[n][n]; ! e[n] = 0.0; n--; iter = 0; --- 508,514 ---- if (l == n) { H[n][n] = H[n][n] + exshift; d[n] = H[n][n]; ! e[n] = 0; n--; iter = 0; *************** *** 516,522 **** } else if (l == n-1) { w = H[n][n-1] * H[n-1][n]; ! p = (H[n-1][n-1] - H[n][n]) / 2.0; q = p * p + w; z = sqrt(abs(q)); H[n][n] = H[n][n] + exshift; --- 516,522 ---- } else if (l == n-1) { w = H[n][n-1] * H[n-1][n]; ! p = (H[n-1][n-1] - H[n][n]) / 2; q = p * p + w; z = sqrt(abs(q)); H[n][n] = H[n][n] + exshift; *************** *** 533,543 **** } d[n-1] = x + z; d[n] = d[n-1]; ! if (z != 0.0) { d[n] = x - w / z; } ! e[n-1] = 0.0; ! e[n] = 0.0; x = H[n][n-1]; s = abs(x) + abs(z); p = x / s; --- 533,543 ---- } d[n-1] = x + z; d[n] = d[n-1]; ! if (z != 0) { d[n] = x - w / z; } ! e[n-1] = 0; ! e[n] = 0; x = H[n][n-1]; s = abs(x) + abs(z); p = x / s; *************** *** 588,595 **** // Form shift x = H[n][n]; ! y = 0.0; ! w = 0.0; if (l < n) { y = H[n-1][n-1]; w = H[n][n-1] * H[n-1][n]; --- 588,595 ---- // Form shift x = H[n][n]; ! y = 0; ! w = 0; if (l < n) { y = H[n-1][n-1]; w = H[n][n-1] * H[n-1][n]; *************** *** 603,616 **** H[i][i] -= x; } s = abs(H[n][n-1]) + abs(H[n-1][n-2]); ! x = y = 0.75 * s; ! w = -0.4375 * s * s; } // MATLAB's new ad hoc shift if (iter == 30) { ! s = (y - x) / 2.0; s = s * s + w; if (s > 0) { s = sqrt(s); --- 603,616 ---- H[i][i] -= x; } s = abs(H[n][n-1]) + abs(H[n-1][n-2]); ! x = y = 3/4 * s; ! w = -4375/10000 * s * s; } // MATLAB's new ad hoc shift if (iter == 30) { ! s = (y - x) / 2; s = s * s + w; if (s > 0) { s = sqrt(s); *************** *** 617,628 **** if (y < x) { s = -s; } ! s = x - w / ((y - x) / 2.0 + s); for (int i = low; i <= n; i++) { H[i][i] -= s; } exshift += s; ! x = y = w = 0.964; } } --- 617,628 ---- if (y < x) { s = -s; } ! s = x - w / ((y - x) / 2 + s); for (int i = low; i <= n; i++) { H[i][i] -= s; } exshift += s; ! x = y = w = 964/1000; } } *************** *** 654,662 **** } for (int i = m+2; i <= n; i++) { ! H[i][i-2] = 0.0; if (i > m+2) { ! H[i][i-3] = 0.0; } } --- 654,662 ---- } for (int i = m+2; i <= n; i++) { ! H[i][i-2] = 0; if (i > m+2) { ! H[i][i-3] = 0; } } *************** *** 667,681 **** 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) { break; } s = sqrt(p * p + q * q + r * r); --- 667,681 ---- if (k != m) { p = H[k][k-1]; q = H[k+1][k-1]; ! r = (notlast ? H[k+2][k-1] : 0); x = abs(p) + abs(q) + abs(r); ! if (x != 0) { p = p / x; q = q / x; r = r / x; } } ! if (x == 0) { break; } s = sqrt(p * p + q * q + r * r); *************** *** 737,743 **** // Backsubstitute to find vectors of upper triangular form ! if (norm == 0.0) { return; } --- 737,743 ---- // Backsubstitute to find vectors of upper triangular form ! if (norm == 0) { return; } *************** *** 749,768 **** if (q == 0) { int l = n; ! H[n][n] = 1.0; for (int i = n-1; i >= 0; i--) { w = H[i][i] - p; ! r = 0.0; for (int j = l; j <= n; j++) { r = r + H[i][j] * H[j][n]; } ! if (e[i] < 0.0) { z = w; s = r; } else { l = i; ! if (e[i] == 0.0) { ! if (w != 0.0) { H[i][n] = -r / w; } else { H[i][n] = -r / (eps * norm); --- 749,768 ---- if (q == 0) { int l = n; ! H[n][n] = 1; for (int i = n-1; i >= 0; i--) { w = H[i][i] - p; ! r = 0; for (int j = l; j <= n; j++) { r = r + H[i][j] * H[j][n]; } ! if (e[i] < 0) { z = w; s = r; } else { l = i; ! if (e[i] == 0) { ! if (w != 0) { H[i][n] = -r / w; } else { H[i][n] = -r / (eps * norm); *************** *** 805,820 **** H[n-1][n-1] = q / H[n][n-1]; H[n-1][n] = -(H[n][n] - p) / H[n][n-1]; } else { ! cdiv(0.0,-H[n-1][n],H[n-1][n-1]-p,q); H[n-1][n-1] = cdivr; H[n-1][n] = cdivi; } ! H[n][n-1] = 0.0; ! H[n][n] = 1.0; for (int i = n-2; i >= 0; i--) { Real ra,sa,vr,vi; ! ra = 0.0; ! sa = 0.0; for (int j = l; j <= n; j++) { ra = ra + H[i][j] * H[j][n-1]; sa = sa + H[i][j] * H[j][n]; --- 805,820 ---- H[n-1][n-1] = q / H[n][n-1]; H[n-1][n] = -(H[n][n] - p) / H[n][n-1]; } else { ! cdiv(0,-H[n-1][n],H[n-1][n-1]-p,q); H[n-1][n-1] = cdivr; H[n-1][n] = cdivi; } ! H[n][n-1] = 0; ! H[n][n] = 1; for (int i = n-2; i >= 0; i--) { Real ra,sa,vr,vi; ! ra = 0; ! sa = 0; for (int j = l; j <= n; j++) { ra = ra + H[i][j] * H[j][n-1]; sa = sa + H[i][j] * H[j][n]; *************** *** 821,827 **** } w = H[i][i] - p; ! if (e[i] < 0.0) { z = w; r = ra; s = sa; --- 821,827 ---- } w = H[i][i] - p; ! if (e[i] < 0) { z = w; r = ra; s = sa; *************** *** 838,845 **** x = H[i][i+1]; y = H[i+1][i]; vr = (d[i] - p) * (d[i] - p) + e[i] * e[i] - q * q; ! vi = (d[i] - p) * 2.0 * q; ! if ((vr == 0.0) && (vi == 0.0)) { vr = eps * norm * (abs(w) + abs(q) + abs(x) + abs(y) + abs(z)); } --- 838,845 ---- x = H[i][i+1]; y = H[i+1][i]; vr = (d[i] - p) * (d[i] - p) + e[i] * e[i] - q * q; ! vi = (d[i] - p) * 2 * q; ! if ((vr == 0) && (vi == 0)) { vr = eps * norm * (abs(w) + abs(q) + abs(x) + abs(y) + abs(z)); } *************** *** 884,890 **** for (int j = nn-1; j >= low; j--) { for (int i = low; i <= high; i++) { ! z = 0.0; for (int k = low; k <= min(j,high); k++) { z = z + V[i][k] * H[k][j]; } --- 884,890 ---- for (int j = nn-1; j >= low; j--) { for (int i = low; i <= high; i++) { ! z = 0; for (int k = low; k <= min(j,high); k++) { z = z + V[i][k] * H[k][j]; } *************** *** 1011,1017 **** D = Array2D (n,n); for (int i = 0; i < n; i++) { for (int j = 0; j < n; j++) { ! D[i][j] = 0.0; } D[i][i] = d[i]; if (e[i] > 0) { --- 1011,1017 ---- D = Array2D (n,n); for (int i = 0; i < n; i++) { for (int j = 0; j < n; j++) { ! D[i][j] = 0; } D[i][i] = d[i]; if (e[i] > 0) { diff -ucr jama.orig/jama_lu.h jama/jama_lu.h *** jama.orig/jama_lu.h Wed Nov 17 04:28:42 2004 --- jama/jama_lu.h Tue May 10 11:26:54 2005 *************** *** 107,113 **** // Most of the time is spent in the following dot product. int kmax = min(i,j); ! double s = 0.0; for (int k = 0; k < kmax; k++) { s += LUrowi[k]*LUcolj[k]; } --- 107,113 ---- // Most of the time is spent in the following dot product. int kmax = min(i,j); ! Real s = 0; for (int k = 0; k < kmax; k++) { s += LUrowi[k]*LUcolj[k]; } *************** *** 126,132 **** if (p != j) { int k=0; for (k = 0; k < n; k++) { ! double t = LU_[p][k]; LU_[p][k] = LU_[j][k]; LU_[j][k] = t; } --- 126,132 ---- if (p != j) { int k=0; for (k = 0; k < n; k++) { ! Real t = LU_[p][k]; LU_[p][k] = LU_[j][k]; LU_[j][k] = t; } *************** *** 138,144 **** // Compute multipliers. ! if ((j < m) && (LU_[j][j] != 0.0)) { for (int i = j+1; i < m; i++) { LU_[i][j] /= LU_[j][j]; } --- 138,144 ---- // Compute multipliers. ! if ((j < m) && (LU_[j][j] != 0)) { for (int i = j+1; i < m; i++) { LU_[i][j] /= LU_[j][j]; } *************** *** 171,179 **** if (i > j) { L_[i][j] = LU_[i][j]; } else if (i == j) { ! L_[i][j] = 1.0; } else { ! L_[i][j] = 0.0; } } } --- 171,179 ---- if (i > j) { L_[i][j] = LU_[i][j]; } else if (i == j) { ! L_[i][j] = 1; } else { ! L_[i][j] = 0; } } } *************** *** 191,197 **** if (i <= j) { U_[i][j] = LU_[i][j]; } else { ! U_[i][j] = 0.0; } } } --- 191,197 ---- if (i <= j) { U_[i][j] = LU_[i][j]; } else { ! U_[i][j] = 0; } } } diff -ucr jama.orig/jama_qr.h jama/jama_qr.h *** jama.orig/jama_qr.h Sat Oct 16 06:07:11 2004 --- jama/jama_qr.h Tue May 10 11:26:54 2005 *************** *** 76,82 **** nrm = hypot(nrm,QR_[i][k]); } ! if (nrm != 0.0) { // Form k-th Householder vector. if (QR_[k][k] < 0) { nrm = -nrm; --- 76,82 ---- nrm = hypot(nrm,QR_[i][k]); } ! if (nrm != 0) { // Form k-th Householder vector. if (QR_[k][k] < 0) { nrm = -nrm; *************** *** 84,94 **** for (i = k; i < m; i++) { QR_[i][k] /= nrm; } ! QR_[k][k] += 1.0; // Apply transformation to remaining columns. for (j = k+1; j < n; j++) { ! Real s = 0.0; for (i = k; i < m; i++) { s += QR_[i][k]*QR_[i][j]; } --- 84,94 ---- for (i = k; i < m; i++) { QR_[i][k] /= nrm; } ! QR_[k][k] += 1; // Apply transformation to remaining columns. for (j = k+1; j < n; j++) { ! Real s = 0; for (i = k; i < m; i++) { s += QR_[i][k]*QR_[i][j]; } *************** *** 141,147 **** if (i >= j) { H[i][j] = QR_[i][j]; } else { ! H[i][j] = 0.0; } } } --- 141,147 ---- if (i >= j) { H[i][j] = QR_[i][j]; } else { ! H[i][j] = 0; } } } *************** *** 164,170 **** } else if (i == j) { R[i][j] = Rdiag[i]; } else { ! R[i][j] = 0.0; } } } --- 164,170 ---- } else if (i == j) { R[i][j] = Rdiag[i]; } else { ! R[i][j] = 0; } } } *************** *** 187,198 **** TNT::Array2D Q(m,n); for (k = n-1; k >= 0; k--) { for (i = 0; i < m; i++) { ! Q[i][k] = 0.0; } ! Q[k][k] = 1.0; for (j = k; j < n; j++) { if (QR_[k][k] != 0) { ! Real s = 0.0; for (i = k; i < m; i++) { s += QR_[i][k]*Q[i][j]; } --- 187,198 ---- TNT::Array2D Q(m,n); for (k = n-1; k >= 0; k--) { for (i = 0; i < m; i++) { ! Q[i][k] = 0; } ! Q[k][k] = 1; for (j = k; j < n; j++) { if (QR_[k][k] != 0) { ! Real s = 0; for (i = k; i < m; i++) { s += QR_[i][k]*Q[i][j]; } *************** *** 229,235 **** // Compute Y = transpose(Q)*b for (int k = 0; k < n; k++) { ! Real s = 0.0; for (int i = k; i < m; i++) { s += QR_[i][k]*x[i]; --- 229,235 ---- // Compute Y = transpose(Q)*b for (int k = 0; k < n; k++) { ! Real s = 0; for (int i = k; i < m; i++) { s += QR_[i][k]*x[i]; *************** *** 282,288 **** // Compute Y = transpose(Q)*B for (k = 0; k < n; k++) { for (j = 0; j < nx; j++) { ! Real s = 0.0; for (i = k; i < m; i++) { s += QR_[i][k]*X[i][j]; } --- 282,288 ---- // Compute Y = transpose(Q)*B for (k = 0; k < n; k++) { for (j = 0; j < nx; j++) { ! Real s = 0; for (i = k; i < m; i++) { s += QR_[i][k]*X[i][j]; } diff -ucr jama.orig/jama_svd.h jama/jama_svd.h *** jama.orig/jama_svd.h Sat Nov 13 03:39:54 2004 --- jama/jama_svd.h Tue May 10 11:26:54 2005 *************** *** 78,100 **** 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. ! double t = 0; for (i = k; i < m; i++) { t += A[i][k]*A[i][j]; } --- 78,100 ---- for (i = k; i < m; i++) { s[k] = hypot(s[k],A[i][k]); } ! if (s[k] != 0) { ! if (A[k][k] < 0) { s[k] = -s[k]; } for (i = k; i < m; i++) { A[i][k] /= s[k]; } ! A[k][k] += 1; } s[k] = -s[k]; } for (j = k+1; j < n; j++) { ! if ((k < nct) && (s[k] != 0)) { // Apply the transformation. ! Real t = 0; for (i = k; i < m; i++) { t += A[i][k]*A[i][j]; } *************** *** 127,148 **** 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++) { --- 127,148 ---- for (i = k+1; i < n; i++) { e[k] = hypot(e[k],e[i]); } ! if (e[k] != 0) { ! if (e[k+1] < 0) { e[k] = -e[k]; } for (i = k+1; i < n; i++) { e[i] /= e[k]; } ! e[k+1] += 1; } e[k] = -e[k]; ! if ((k+1 < m) & (e[k] != 0)) { // Apply the transformation. for (i = k+1; i < m; i++) { ! work[i] = 0; } for (j = k+1; j < n; j++) { for (i = k+1; i < m; i++) { *************** *** 150,156 **** } } for (j = k+1; j < n; j++) { ! double t = -e[j]/e[k+1]; for (i = k+1; i < m; i++) { A[i][j] += t*work[i]; } --- 150,156 ---- } } 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]; } *************** *** 175,186 **** 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. --- 175,186 ---- s[nct] = A[nct][nct]; } if (m < p) { ! s[p-1] = 0; } if (nrt+1 < p) { e[nrt] = A[nrt][p-1]; } ! e[p-1] = 0; // If required, generate U. *************** *** 187,200 **** 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++) { ! double t = 0; for (i = k; i < m; i++) { t += U[i][k]*U[i][j]; } --- 187,200 ---- if (wantu) { for (j = nct; j < nu; j++) { for (i = 0; i < m; i++) { ! U[i][j] = 0; } ! U[j][j] = 1; } for (k = nct-1; k >= 0; k--) { ! if (s[k] != 0) { for (j = k+1; j < nu; j++) { ! Real t = 0; for (i = k; i < m; i++) { t += U[i][k]*U[i][j]; } *************** *** 206,220 **** 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; } } } --- 206,220 ---- for (i = k; i < m; i++ ) { U[i][k] = -U[i][k]; } ! U[k][k] = 1 + U[k][k]; for (i = 0; i < k-1; i++) { ! U[i][k] = 0; } } else { for (i = 0; i < m; i++) { ! U[i][k] = 0; } ! U[k][k] = 1; } } } *************** *** 223,231 **** if (wantv) { for (k = n-1; k >= 0; k--) { ! if ((k < nrt) & (e[k] != 0.0)) { for (j = k+1; j < nu; j++) { ! double t = 0; for (i = k+1; i < n; i++) { t += V[i][k]*V[i][j]; } --- 223,231 ---- if (wantv) { for (k = n-1; k >= 0; k--) { ! if ((k < nrt) & (e[k] != 0)) { for (j = k+1; j < nu; j++) { ! Real t = 0; for (i = k+1; i < n; i++) { t += V[i][k]*V[i][j]; } *************** *** 236,244 **** } } for (i = 0; i < n; i++) { ! V[i][k] = 0.0; } ! V[k][k] = 1.0; } } --- 236,244 ---- } } for (i = 0; i < n; i++) { ! V[i][k] = 0; } ! V[k][k] = 1; } } *************** *** 246,252 **** int pp = p-1; int iter = 0; ! double eps = pow(2.0,-52.0); while (p > 0) { int k=0; int kase=0; --- 246,252 ---- int pp = p-1; int iter = 0; ! Real eps = pow(2,-52); while (p > 0) { int k=0; int kase=0; *************** *** 268,274 **** break; } if (abs(e[k]) <= eps*(abs(s[k]) + abs(s[k+1]))) { ! e[k] = 0.0; break; } } --- 268,274 ---- break; } if (abs(e[k]) <= eps*(abs(s[k]) + abs(s[k+1]))) { ! e[k] = 0; break; } } *************** *** 280,289 **** if (ks == k) { break; } ! double 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; } } --- 280,289 ---- 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; break; } } *************** *** 305,316 **** // Deflate negligible s(p). case 1: { ! double f = e[p-2]; ! e[p-2] = 0.0; for (j = p-2; j >= k; j--) { ! double t = hypot(s[j],f); ! double cs = s[j]/t; ! double sn = f/t; s[j] = t; if (j != k) { f = -sn*e[j-1]; --- 305,316 ---- // Deflate negligible s(p). case 1: { ! Real f = e[p-2]; ! e[p-2] = 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]; *************** *** 330,341 **** // Split at negligible s(k). case 2: { ! double f = e[k-1]; ! e[k-1] = 0.0; for (j = k; j < p; j++) { ! double t = hypot(s[j],f); ! double cs = s[j]/t; ! double sn = f/t; s[j] = t; f = -sn*e[j]; e[j] = cs*e[j]; --- 330,341 ---- // Split at negligible s(k). case 2: { ! Real f = e[k-1]; ! e[k-1] = 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]; *************** *** 356,388 **** // Calculate the shift. ! double scale = max(max(max(max( abs(s[p-1]),abs(s[p-2])),abs(e[p-2])), abs(s[k])),abs(e[k])); ! double sp = s[p-1]/scale; ! double spm1 = s[p-2]/scale; ! double epm1 = e[p-2]/scale; ! double sk = s[k]/scale; ! double ek = e[k]/scale; ! double b = ((spm1 + sp)*(spm1 - sp) + epm1*epm1)/2.0; ! double c = (sp*epm1)*(sp*epm1); ! double 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); } ! double f = (sk + sp)*(sk - sp) + shift; ! double g = sk*ek; // Chase zeros. for (j = k; j < p-1; j++) { ! double t = hypot(f,g); ! double cs = f/t; ! double sn = g/t; if (j != k) { e[j-1] = t; } --- 356,388 ---- // 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; ! Real c = (sp*epm1)*(sp*epm1); ! Real shift = 0; ! if ((b != 0) || (c != 0)) { shift = sqrt(b*b + c); ! if (b < 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; } *************** *** 424,431 **** // 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]; --- 424,431 ---- // Make the singular values positive. ! if (s[k] <= 0) { ! s[k] = (s[k] < 0 ? -s[k] : 0); if (wantv) { for (i = 0; i <= pp; i++) { V[i][k] = -V[i][k]; *************** *** 439,445 **** if (s[k] >= s[k+1]) { break; } ! double t = s[k]; s[k] = s[k+1]; s[k+1] = t; if (wantv && (k < n-1)) { --- 439,445 ---- 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)) { *************** *** 497,503 **** A = Array2D (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]; } --- 497,503 ---- A = Array2D (n,n); for (int i = 0; i < n; i++) { for (int j = 0; j < n; j++) { ! A[i][j] = 0; } A[i][i] = s[i]; } *************** *** 505,517 **** /** Two norm (max(S)) */ ! double norm2 () { return s[0]; } /** Two norm of condition number (max(S)/min(S)) */ ! double cond () { return s[0]/s[min(m,n)-1]; } --- 505,517 ---- /** 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]; } *************** *** 521,528 **** int rank () { ! double eps = pow(2.0,-52.0); ! double tol = max(m,n)*s[0]*eps; int r = 0; for (int i = 0; i < s.dim(); i++) { if (s[i] > tol) { --- 521,528 ---- int rank () { ! Real eps = pow(2,-52); ! Real tol = max(m,n)*s[0]*eps; int r = 0; for (int i = 0; i < s.dim(); i++) { if (s[i] > tol) {