001: package JSci.maths;
002:
003: import JSci.maths.matrices.*;
004: import JSci.maths.vectors.*;
005: import JSci.maths.polynomials.RealPolynomial;
006:
007: /**
008: * The linear math library.
009: * This class cannot be subclassed or instantiated because all methods are static.
010: * @version 2.2
011: * @author Mark Hale
012: */
013: public final class LinearMath extends AbstractMath {
014: private LinearMath() {
015: }
016:
017: // LINEAR SYSTEM
018:
019: /**
020: * Solves the linear system Mx=v.
021: * @param M a double square matrix.
022: * @param v a double vector.
023: * @return the double vector x.
024: */
025: public static AbstractDoubleVector solve(
026: final AbstractDoubleSquareMatrix M,
027: final AbstractDoubleVector v) {
028: final int n = v.dimension();
029: final double array[] = new double[n];
030: final int pivot[] = new int[n + 1];
031: final AbstractDoubleSquareMatrix lu[] = M.luDecompose(pivot);
032: int i, j;
033: double sum;
034: // forward substitution
035: for (i = 0; i < n; i++) {
036: sum = v.getComponent(pivot[i]);
037: for (j = 0; j < i; j++)
038: sum -= lu[0].getElement(i, j) * array[j];
039: array[i] = sum / lu[0].getElement(i, i);
040: }
041: // back substitution
042: for (i = n - 1; i >= 0; i--) {
043: sum = array[i];
044: for (j = i + 1; j < n; j++)
045: sum -= lu[1].getElement(i, j) * array[j];
046: array[i] = sum / lu[1].getElement(i, i);
047: }
048: return new DoubleVector(array);
049: }
050:
051: /**
052: * Solves a linear system using cholesky decomposition.
053: */
054: private static double[] solveCholesky(final double m[][],
055: final double v[]) {
056: final int n = v.length;
057: final double array[] = new double[n];
058: final AbstractDoubleSquareMatrix lu[] = new DoubleSquareMatrix(
059: m).choleskyDecompose();
060: int i, j;
061: double sum;
062: // forward substitution
063: for (i = 0; i < n; i++) {
064: sum = v[i];
065: for (j = 0; j < i; j++)
066: sum -= lu[0].getElement(i, j) * array[j];
067: array[i] = sum / lu[0].getElement(i, i);
068: }
069: // back substitution
070: for (i = n - 1; i >= 0; i--) {
071: sum = array[i];
072: for (j = i + 1; j < n; j++)
073: sum -= lu[1].getElement(i, j) * array[j];
074: array[i] = sum / lu[1].getElement(i, i);
075: }
076: return array;
077: }
078:
079: /**
080: * Solves a linear system using QR decomposition.
081: */
082: private static double[] solveQR(final double m[][],
083: final double v[]) {
084: final int n = v.length;
085: final double array[] = new double[n];
086: final AbstractDoubleSquareMatrix lu[] = new DoubleSquareMatrix(
087: m).qrDecompose();
088: int i, j;
089: double sum;
090: // Q.transpose().multiply(v)
091: for (i = 0; i < n; i++) {
092: sum = 0.0;
093: for (j = 0; j < n; j++)
094: sum += lu[0].getElement(j, i) * v[j];
095: array[i] = sum;
096: }
097: // back substitution
098: for (i = n - 1; i >= 0; i--) {
099: sum = array[i];
100: for (j = i + 1; j < n; j++)
101: sum -= lu[1].getElement(i, j) * array[j];
102: array[i] = sum / lu[1].getElement(i, i);
103: }
104: return array;
105: }
106:
107: /**
108: * Solves the unsymmetric linear system Ax=b using the
109: * Generalized Minimum Residual method (doesn't require A
110: * to be nonsingular).
111: * While slower than LU decomposition, it is more
112: * robust and should be used with large matrices.
113: * It is guaranted to converge exactly in N iterations for an
114: * N by N matrix (minus some numerical errors).
115: * @author Alain Beliveau
116: * @author Daniel Lemire
117: * @param max_iter maximum number of iterations.
118: * @param tol tolerance.
119: * @exception IllegalArgumentException If either the tolerance
120: * or the number of iterations is not positive.
121: * Also, if an unexpected error occurs.
122: * @exception MaximumIterationsExceededException If it cannot
123: * converge according to the given parameters.
124: */
125: public static AbstractDoubleVector solveGMRes(
126: final AbstractDoubleMatrix A, final AbstractDoubleVector b,
127: int max_iter, double tol)
128: throws MaximumIterationsExceededException {
129: if (max_iter <= 0)
130: throw new IllegalArgumentException(
131: "Number of allowed iterations must be a positive integer: "
132: + max_iter + " <= 0.");
133: if (tol < 0)
134: throw new IllegalArgumentException(
135: "Tolerance must be positive or zero: " + tol
136: + " < 0.");
137: final int m = A.rows();
138: double resid;
139: int i, j = 1, k;
140: double[] s = new double[m + 1];
141: double[][] cs = new double[m + 1][2];
142: double[] rotmp = new double[2];
143: AbstractDoubleVector w;
144: AbstractDoubleVector x = new DoubleVector(A.rows());
145: double normb = b.norm();
146: AbstractDoubleVector r = b.subtract(A.multiply(x));
147: double beta = r.norm();
148: if (normb == 0.0)
149: normb = 1.0;
150: if ((resid = r.norm() / normb) <= tol) {
151: tol = resid;
152: max_iter = 0;
153: throw new IllegalArgumentException("There is a bug.");
154: }
155: AbstractDoubleVector[] v = new DoubleVector[m + 1];
156: double[][] H = new double[m + 1][m];
157: while (j <= max_iter) {
158: v[0] = r.scalarMultiply(1.0 / beta);
159: for (i = 0; i < m + 1; i++)
160: s[i] = 0.0;
161: s[0] = beta;
162: for (i = 0; i < m && j <= max_iter; i++, j++) {
163: w = A.multiply(v[i]);
164: for (k = 0; k <= i; k++) {
165: H[k][i] = w.scalarProduct(v[k]);
166: w = w.subtract(v[k].scalarMultiply(H[k][i]));
167: }
168: H[i + 1][i] = w.norm();
169: v[i + 1] = w.scalarMultiply(1.0 / H[i + 1][i]);
170: for (k = 0; k < i; k++) {
171: rotmp = applyPlaneRotation(H[k][i], H[k + 1][i],
172: cs[k][0], cs[k][1]);
173: H[k][i] = rotmp[0];
174: H[k + 1][i] = rotmp[1];
175: }
176: cs[i] = generatePlaneRotation(H[i][i], H[i + 1][i]);
177: rotmp = applyPlaneRotation(H[i][i], H[i + 1][i],
178: cs[i][0], cs[i][1]);
179: H[i][i] = rotmp[0];
180: H[i + 1][i] = rotmp[1];
181: rotmp = applyPlaneRotation(s[i], s[i + 1], cs[i][0],
182: cs[i][1]);
183: s[i] = rotmp[0];
184: s[i + 1] = rotmp[1];
185: if ((resid = Math.abs(s[i + 1]) / normb) < tol) {
186: x = update(x, i, H, s, v);
187: tol = resid;
188: max_iter = j;
189: return (x);
190: }
191: }
192: x = update(x, m - 1, H, s, v);
193: r = b.subtract(A.multiply(x));
194: beta = r.norm();
195: if ((resid = beta / normb) < tol) {
196: tol = resid;
197: max_iter = j;
198: return (x);
199: }
200: }
201: tol = resid;
202: throw new MaximumIterationsExceededException(
203: "(tol) "
204: + tol
205: + ". It doesn't converge in "
206: + max_iter
207: + "iterations. Try raising the number of allowed iterations or raising the tolerance.");
208: }
209:
210: private static double[] generatePlaneRotation(double dx, double dy) {
211: double[] cs = new double[2];
212: double temp;
213: if (dy == 0.0) {
214: cs[0] = 1.0;
215: cs[1] = 0.0;
216: } else if (Math.abs(dy) > Math.abs(dx)) {
217: temp = dx / dy;
218: cs[1] = 1.0 / Math.sqrt(1.0 + temp * temp);
219: cs[0] = temp * cs[1];
220: } else {
221: temp = dy / dx;
222: cs[0] = 1.0 / Math.sqrt(1.0 + temp * temp);
223: cs[1] = temp * cs[0];
224: }
225: return (cs);
226: }
227:
228: private static double[] applyPlaneRotation(double dx, double dy,
229: double cs, double sn) {
230: double[] dxy = new double[2];
231: dxy[0] = cs * dx + sn * dy;
232: dxy[1] = -sn * dx + cs * dy;
233: return (dxy);
234: }
235:
236: private static AbstractDoubleVector update(AbstractDoubleVector x,
237: int k, double[][] H, double[] s, AbstractDoubleVector[] v) {
238: /************************************************
239: * Backsolve of a triangular + sub diagonal matrix:
240: ************************************************/
241: for (int i = k; i >= 0; i--) {
242: s[i] = s[i] / H[i][i];
243: for (int j = i - 1; j >= 0; j--)
244: s[j] = s[j] - H[j][i] * s[i];
245: }
246: for (int j = 0; j <= k; j++) {
247: x = x.add(v[j].scalarMultiply(s[j]));
248: }
249: return (x);
250: }
251:
252: // LEAST SQUARES
253:
254: /**
255: * Fits an nth degree polynomial to data using the method of least squares.
256: * @param n the degree of the polynomial (>= 0).
257: * @param data
258: * [0][] contains the x-series,
259: * [1][] contains the y-series.
260: */
261: public static RealPolynomial leastSquaresFit(int n,
262: final double data[][]) {
263: int i, j, nm1 = n++;
264: double xsum, xysum, tmp;
265: final double mArray[][] = new double[n][n];
266: final double vArray[] = new double[n];
267: for (i = 0; i < n; i++) {
268: xsum = xysum = 0.0;
269: for (j = 0; j < data[0].length; j++) {
270: tmp = Math.pow(data[0][j], i);
271: xsum += tmp;
272: xysum += tmp * data[1][j];
273: }
274: mArray[0][i] = xsum;
275: vArray[i] = xysum;
276: }
277: for (i = 1; i < n; i++) {
278: System.arraycopy(mArray[i - 1], 1, mArray[i], 0, nm1);
279: xsum = 0.0;
280: for (j = 0; j < data[0].length; j++)
281: xsum += Math.pow(data[0][j], nm1 + i);
282: mArray[i][nm1] = xsum;
283: }
284: return new RealPolynomial(solveCholesky(mArray, vArray));
285: }
286:
287: // LINEAR REGRESSION
288:
289: /**
290: * Fits a line to multi-dimensional data using the method of least squares.
291: * @param data
292: * [0...n-1][] contains the x-series' (they must be linearly uncorrelated),
293: * [n][] contains the y-series.
294: * @return a vector containing the coefficients (zero component is the intercept, the rest are gradient components).
295: * E.g. y(x1, x2, ...) = coeffs(0) + coeffs(1) * x1 + coeffs(2) * x2 + ...
296: */
297: public static AbstractDoubleVector linearRegression(
298: final double data[][]) {
299: final int y = data.length - 1;
300: int i, j, k;
301: double xsum, xysum;
302: final double mArray[][] = new double[data.length][data.length];
303: final double vArray[] = new double[data.length];
304: mArray[0][0] = data[0].length;
305: for (j = 1; j < data.length; j++) {
306: xsum = 0.0;
307: for (k = 0; k < data[0].length; k++)
308: xsum += data[j - 1][k];
309: mArray[0][j] = mArray[j][0] = xsum;
310: }
311: xysum = 0.0;
312: for (k = 0; k < data[0].length; k++)
313: xysum += data[y][k];
314: vArray[0] = xysum;
315: for (i = 1; i < data.length; i++) {
316: for (j = i; j < data.length; j++) {
317: xsum = 0.0;
318: for (k = 0; k < data[0].length; k++)
319: xsum += data[i - 1][k] * data[j - 1][k];
320: mArray[i][j] = mArray[j][i] = xsum;
321: }
322: xysum = 0.0;
323: for (k = 0; k < data[0].length; k++)
324: xysum += data[i - 1][k] * data[y][k];
325: vArray[i] = xysum;
326: }
327: return new DoubleVector(solveCholesky(mArray, vArray));
328: }
329:
330: // GRAM-SCHMIDT
331:
332: /**
333: * The Gram-Schmidt orthonormalization method.
334: * @param vecs a set of linearly independent vectors.
335: * @return a set of orthonormal vectors.
336: */
337: public static AbstractDoubleVector[] orthonormalize(
338: AbstractDoubleVector vecs[]) {
339: final int N = vecs.length;
340: AbstractDoubleVector orthovecs[] = new DoubleVector[N];
341: for (int i = 0; i < N; i++) {
342: orthovecs[i] = vecs[i];
343: for (int j = 0; j < i; j++)
344: orthovecs[i] = orthovecs[i].subtract(orthovecs[j]
345: .scalarMultiply(orthovecs[j]
346: .scalarProduct(vecs[i])));
347: orthovecs[i].normalize();
348: }
349: return orthovecs;
350: }
351:
352: // EIGENVALUES & EIGENVECTORS
353:
354: /**
355: * This method finds the eigenvalues of a Hermitian matrix.
356: * @param matrix a Hermitian matrix.
357: * @return an array containing the eigenvalues.
358: * @exception MaximumIterationsExceededException If it takes too many iterations to determine an eigenvalue.
359: */
360: public static double[] eigenvalueSolveHermitian(
361: final AbstractComplexSquareMatrix matrix)
362: throws MaximumIterationsExceededException {
363: final int n = matrix.rows();
364: final double matrix2[][] = new double[2 * n][2 * n];
365: double real, imag;
366: for (int j, i = 0; i < n; i++) {
367: for (j = 0; j < n; j++) {
368: real = matrix.getRealElement(i, j);
369: imag = matrix.getImagElement(i, j);
370: matrix2[i][j] = real;
371: matrix2[n + i][n + j] = real;
372: matrix2[n + i][j] = imag;
373: matrix2[i][n + j] = -imag;
374: }
375: }
376: final double eigenvalue2[] = new double[2 * n];
377: final double offdiag[] = new double[2 * n];
378: reduceSymmetric1_SquareToTridiagonal(matrix2, eigenvalue2,
379: offdiag);
380: System.arraycopy(offdiag, 1, offdiag, 0, n - 1);
381: offdiag[n - 1] = 0.0;
382: eigenvalueSolveSymmetricTridiagonalMatrix(eigenvalue2, offdiag);
383: final double eigenvalue[] = new double[n];
384: System.arraycopy(eigenvalue2, 0, eigenvalue, 0, n);
385: return eigenvalue;
386: }
387:
388: /**
389: * This method finds the eigenvalues and eigenvectors of a Hermitian matrix.
390: * @param matrix a Hermitian matrix.
391: * @param eigenvector an empty array of complex vectors to hold the eigenvectors.
392: * All eigenvectors will be orthogonal.
393: * @return an array containing the eigenvalues.
394: * @exception MaximumIterationsExceededException If it takes too many iterations to determine an eigenvalue.
395: */
396: public static double[] eigenSolveHermitian(
397: final AbstractComplexSquareMatrix matrix,
398: final AbstractComplexVector eigenvector[])
399: throws MaximumIterationsExceededException {
400: final int n = matrix.rows();
401: final double matrix2[][] = new double[2 * n][2 * n];
402: int i, j;
403: double real, imag;
404: for (i = 0; i < n; i++) {
405: for (j = 0; j < n; j++) {
406: real = matrix.getRealElement(i, j);
407: imag = matrix.getImagElement(i, j);
408: matrix2[i][j] = real;
409: matrix2[n + i][n + j] = real;
410: matrix2[n + i][j] = imag;
411: matrix2[i][n + j] = -imag;
412: }
413: }
414: final double eigenvalue2[] = new double[2 * n];
415: final double offdiag[] = new double[2 * n];
416: reduceSymmetric2_SquareToTridiagonal(matrix2, eigenvalue2,
417: offdiag);
418: System.arraycopy(offdiag, 1, offdiag, 0, n - 1);
419: offdiag[n - 1] = 0.0;
420: eigenSolveSymmetricTridiagonalMatrix(eigenvalue2, offdiag,
421: matrix2);
422: final double eigenvalue[] = new double[n];
423: double arrayRe[], arrayIm[];
424: for (i = 0; i < n; i++) {
425: eigenvalue[i] = eigenvalue2[i];
426: arrayRe = new double[n];
427: arrayIm = new double[n];
428: for (j = 0; j < n; j++) {
429: arrayRe[j] = matrix2[j][i];
430: arrayIm[j] = matrix2[j + n][i];
431: }
432: eigenvector[i] = new ComplexVector(arrayRe, arrayIm);
433: }
434: return eigenvalue;
435: }
436:
437: /**
438: * This method finds the eigenvalues of a symmetric tridiagonal matrix by the QL method.
439: * It is based on the NETLIB algol/fortran procedure tql1 by Bowdler, Martin, Reinsch and Wilkinson.
440: * @param matrix a double symmetric tridiagonal matrix.
441: * @return an array containing the eigenvalues.
442: * @exception MaximumIterationsExceededException If it takes too many iterations to determine an eigenvalue.
443: */
444: public static double[] eigenvalueSolveSymmetric(
445: final DoubleTridiagonalMatrix matrix)
446: throws MaximumIterationsExceededException {
447: final int n = matrix.rows();
448: final int nm1 = n - 1;
449: final double eigenvalue[] = new double[n];
450: final double offdiag[] = new double[n];
451: for (int i = 0; i < nm1; i++) {
452: eigenvalue[i] = matrix.getElement(i, i);
453: offdiag[i] = matrix.getElement(i, i + 1);
454: }
455: eigenvalue[nm1] = matrix.getElement(nm1, nm1);
456: offdiag[nm1] = 0.0;
457: eigenvalueSolveSymmetricTridiagonalMatrix(eigenvalue, offdiag);
458: return eigenvalue;
459: }
460:
461: /**
462: * This method finds the eigenvalues and eigenvectors of a symmetric tridiagonal matrix by the QL method.
463: * It is based on the NETLIB algol/fortran procedure tql2 by Bowdler, Martin, Reinsch and Wilkinson.
464: * @param matrix a double symmetric tridiagonal matrix.
465: * @param eigenvector an empty array of double vectors to hold the eigenvectors.
466: * All eigenvectors will be orthogonal.
467: * @return an array containing the eigenvalues.
468: * @exception MaximumIterationsExceededException If it takes too many iterations to determine an eigenvalue.
469: */
470: public static double[] eigenSolveSymmetric(
471: final DoubleTridiagonalMatrix matrix,
472: final AbstractDoubleVector eigenvector[])
473: throws MaximumIterationsExceededException {
474: final int n = matrix.rows();
475: final int nm1 = n - 1;
476: final double eigenvalue[] = new double[n];
477: final double offdiag[] = new double[n];
478: final double id[][] = new double[n][n];
479: int i, j;
480: for (i = 0; i < nm1; i++) {
481: id[i][i] = 1.0;
482: eigenvalue[i] = matrix.getElement(i, i);
483: offdiag[i] = matrix.getElement(i, i + 1);
484: }
485: id[nm1][nm1] = 1.0;
486: eigenvalue[nm1] = matrix.getElement(nm1, nm1);
487: offdiag[nm1] = 0.0;
488: eigenSolveSymmetricTridiagonalMatrix(eigenvalue, offdiag, id);
489: for (i = 0; i < n; i++) {
490: DoubleVector evec = new DoubleVector(n);
491: for (j = 0; j < n; j++)
492: evec.setComponent(j, id[j][i]);
493: eigenvector[i] = evec;
494: }
495: return eigenvalue;
496: }
497:
498: /**
499: * This method finds the eigenvalues of a symmetric square matrix.
500: * The matrix is reduced to tridiagonal form and then the QL method is applied.
501: * It is based on the NETLIB algol/fortran procedure tred1/tql1 by Bowdler, Martin, Reinsch and Wilkinson.
502: * @param matrix a double symmetric square matrix.
503: * @return an array containing the eigenvalues.
504: * @exception MaximumIterationsExceededException If it takes too many iterations to determine an eigenvalue.
505: */
506: public static double[] eigenvalueSolveSymmetric(
507: final AbstractDoubleSquareMatrix matrix)
508: throws MaximumIterationsExceededException {
509: final int n = matrix.rows();
510: final double eigenvalue[] = new double[n];
511: final double offdiag[] = new double[n];
512: final double array[][] = new double[n][n];
513: int i, j;
514: for (i = 0; i < n; i++) {
515: for (j = 0; j < n; j++)
516: array[i][j] = matrix.getElement(i, j);
517: }
518: reduceSymmetric1_SquareToTridiagonal(array, eigenvalue, offdiag);
519: System.arraycopy(offdiag, 1, offdiag, 0, n - 1);
520: offdiag[n - 1] = 0.0;
521: eigenvalueSolveSymmetricTridiagonalMatrix(eigenvalue, offdiag);
522: return eigenvalue;
523: }
524:
525: /**
526: * This method finds the eigenvalues and eigenvectors of a symmetric square matrix.
527: * The matrix is reduced to tridiagonal form and then the QL method is applied.
528: * It is based on the NETLIB algol/fortran procedure tred2/tql2 by Bowdler, Martin, Reinsch and Wilkinson.
529: * @param matrix a double symmetric square matrix.
530: * @param eigenvector an empty array of double vectors to hold the eigenvectors.
531: * All eigenvectors will be orthogonal.
532: * @return an array containing the eigenvalues.
533: * @exception MaximumIterationsExceededException If it takes too many iterations to determine an eigenvalue.
534: */
535: public static double[] eigenSolveSymmetric(
536: final AbstractDoubleSquareMatrix matrix,
537: final AbstractDoubleVector eigenvector[])
538: throws MaximumIterationsExceededException {
539: final int n = matrix.rows();
540: final double eigenvalue[] = new double[n];
541: final double offdiag[] = new double[n];
542: final double transf[][] = new double[n][n];
543: int i, j;
544: for (i = 0; i < n; i++) {
545: for (j = 0; j < n; j++)
546: transf[i][j] = matrix.getElement(i, j);
547: }
548: reduceSymmetric2_SquareToTridiagonal(transf, eigenvalue,
549: offdiag);
550: System.arraycopy(offdiag, 1, offdiag, 0, n - 1);
551: offdiag[n - 1] = 0.0;
552: eigenSolveSymmetricTridiagonalMatrix(eigenvalue, offdiag,
553: transf);
554: for (i = 0; i < n; i++) {
555: DoubleVector evec = new DoubleVector(n);
556: for (j = 0; j < n; j++)
557: evec.setComponent(j, transf[j][i]);
558: eigenvector[i] = evec;
559: }
560: return eigenvalue;
561: }
562:
563: private final static int EIGEN_MAX_ITERATIONS = 250;
564:
565: /**
566: * Internal NETLIB tql1 routine.
567: * @param diag output eigenvalues.
568: * @author Richard Cannings
569: */
570: private static void eigenvalueSolveSymmetricTridiagonalMatrix(
571: final double diag[], final double offdiag[])
572: throws MaximumIterationsExceededException {
573: final int n = diag.length;
574: final int nm1 = n - 1;
575: int m, l, iteration, i, k;
576: double s, r, p, g, f, dd, c, b;
577: for (l = 0; l < n; l++) {
578: iteration = 0;
579: do {
580: for (m = l; m < nm1; m++) {
581: dd = Math.abs(diag[m]) + Math.abs(diag[m + 1]);
582: if (Math.abs(offdiag[m]) + dd == dd)
583: break;
584: }
585: if (m != l) {
586: if (iteration++ == EIGEN_MAX_ITERATIONS)
587: throw new MaximumIterationsExceededException(
588: "No convergence after "
589: + EIGEN_MAX_ITERATIONS
590: + " iterations.");
591: g = (diag[l + 1] - diag[l]) / (2.0 * offdiag[l]);
592: r = Math.sqrt(g * g + 1.0);
593: g = diag[m]
594: - diag[l]
595: + offdiag[l]
596: / (g + (g < 0.0 ? -Math.abs(r) : Math
597: .abs(r)));
598: s = c = 1.0;
599: p = 0.0;
600: for (i = m - 1; i >= l; i--) {
601: f = s * offdiag[i];
602: b = c * offdiag[i];
603: if (Math.abs(f) >= Math.abs(g)) {
604: c = g / f;
605: r = Math.sqrt(c * c + 1.0);
606: offdiag[i + 1] = f * r;
607: s = 1 / r;
608: c *= s;
609: } else {
610: s = f / g;
611: r = Math.sqrt(s * s + 1.0);
612: offdiag[i + 1] = g * r;
613: c = 1 / r;
614: s *= c;
615: }
616: g = diag[i + 1] - p;
617: r = (diag[i] - g) * s + 2.0 * c * b;
618: p = s * r;
619: diag[i + 1] = g + p;
620: g = c * r - b;
621: }
622: diag[l] = diag[l] - p;
623: offdiag[l] = g;
624: offdiag[m] = 0.0;
625: }
626: } while (m != l);
627: }
628: }
629:
630: /**
631: * Internal NETLIB tred1 routine.
632: * @author Richard Cannings
633: */
634: private static void reduceSymmetric1_SquareToTridiagonal(
635: final double matrix[][], final double diag[],
636: final double offdiag[]) {
637: final int n = diag.length;
638: int i, j, k, l;
639: double f, g, h, hh, scale;
640: for (i = n - 1; i > 0; i--) {
641: l = i - 1;
642: h = scale = 0.0;
643: if (l > 0) {
644: for (k = 0; k <= l; k++)
645: scale += Math.abs(matrix[i][k]);
646: if (scale == 0.0)
647: offdiag[i] = matrix[i][l];
648: else {
649: for (k = 0; k <= l; k++) {
650: matrix[i][k] /= scale;
651: h += matrix[i][k] * matrix[i][k];
652: }
653: f = matrix[i][l];
654: g = (f >= 0.0 ? -Math.sqrt(h) : Math.sqrt(h));
655: offdiag[i] = scale * g;
656: h -= f * g;
657: matrix[i][l] = f - g;
658: f = 0.0;
659: for (j = 0; j <= l; j++) {
660: g = 0.0;
661: for (k = 0; k <= j; k++)
662: g += matrix[j][k] * matrix[i][k];
663: for (k = j + 1; k <= l; k++)
664: g += matrix[k][j] * matrix[i][k];
665: offdiag[j] = g / h;
666: f += offdiag[j] * matrix[i][j];
667: }
668: hh = f / (h + h);
669: for (j = 0; j <= l; j++) {
670: f = matrix[i][j];
671: offdiag[j] = g = offdiag[j] - hh * f;
672: for (k = 0; k <= j; k++)
673: matrix[j][k] -= f * offdiag[k] + g
674: * matrix[i][k];
675: }
676: }
677: } else
678: offdiag[i] = matrix[i][l];
679: diag[i] = h;
680: }
681: offdiag[0] = 0.0;
682: for (i = 0; i < n; i++)
683: diag[i] = matrix[i][i];
684: }
685:
686: /**
687: * Internal NETLIB tql2 routine.
688: * @param diag output eigenvalues.
689: * @param transf output eigenvectors.
690: * @author Richard Cannings
691: */
692: private static void eigenSolveSymmetricTridiagonalMatrix(
693: final double diag[], final double offdiag[],
694: final double transf[][])
695: throws MaximumIterationsExceededException {
696: final int n = diag.length;
697: final int nm1 = n - 1;
698: int m, l, iteration, i, k;
699: double s, r, p, g, f, dd, c, b;
700: for (l = 0; l < n; l++) {
701: iteration = 0;
702: do {
703: for (m = l; m < nm1; m++) {
704: dd = Math.abs(diag[m]) + Math.abs(diag[m + 1]);
705: if (Math.abs(offdiag[m]) + dd == dd)
706: break;
707: }
708: if (m != l) {
709: if (iteration++ == EIGEN_MAX_ITERATIONS)
710: throw new MaximumIterationsExceededException(
711: "No convergence after "
712: + EIGEN_MAX_ITERATIONS
713: + " iterations.");
714: g = (diag[l + 1] - diag[l]) / (2.0 * offdiag[l]);
715: r = Math.sqrt(g * g + 1.0);
716: g = diag[m]
717: - diag[l]
718: + offdiag[l]
719: / (g + (g < 0.0 ? -Math.abs(r) : Math
720: .abs(r)));
721: s = c = 1.0;
722: p = 0.0;
723: for (i = m - 1; i >= l; i--) {
724: f = s * offdiag[i];
725: b = c * offdiag[i];
726: if (Math.abs(f) >= Math.abs(g)) {
727: c = g / f;
728: r = Math.sqrt(c * c + 1.0);
729: offdiag[i + 1] = f * r;
730: s = 1 / r;
731: c *= s;
732: } else {
733: s = f / g;
734: r = Math.sqrt(s * s + 1.0);
735: offdiag[i + 1] = g * r;
736: c = 1 / r;
737: s *= c;
738: }
739: g = diag[i + 1] - p;
740: r = (diag[i] - g) * s + 2.0 * c * b;
741: p = s * r;
742: diag[i + 1] = g + p;
743: g = c * r - b;
744: for (k = 0; k < n; k++) {
745: f = transf[k][i + 1];
746: transf[k][i + 1] = s * transf[k][i] + c * f;
747: transf[k][i] = c * transf[k][i] - s * f;
748: }
749: }
750: diag[l] = diag[l] - p;
751: offdiag[l] = g;
752: offdiag[m] = 0.0;
753: }
754: } while (m != l);
755: }
756: }
757:
758: /**
759: * Internal NETLIB tred2 routine.
760: * @param matrix output orthogonal transformations.
761: * @author Richard Cannings
762: */
763: private static void reduceSymmetric2_SquareToTridiagonal(
764: final double matrix[][], final double diag[],
765: final double offdiag[]) {
766: final int n = diag.length;
767: int i, j, k, l;
768: double f, g, h, hh, scale;
769: for (i = n - 1; i > 0; i--) {
770: l = i - 1;
771: h = scale = 0.0;
772: if (l > 0) {
773: for (k = 0; k <= l; k++)
774: scale += Math.abs(matrix[i][k]);
775: if (scale == 0.0)
776: offdiag[i] = matrix[i][l];
777: else {
778: for (k = 0; k <= l; k++) {
779: matrix[i][k] /= scale;
780: h += matrix[i][k] * matrix[i][k];
781: }
782: f = matrix[i][l];
783: g = (f >= 0.0 ? -Math.sqrt(h) : Math.sqrt(h));
784: offdiag[i] = scale * g;
785: h -= f * g;
786: matrix[i][l] = f - g;
787: f = 0.0;
788: for (j = 0; j <= l; j++) {
789: matrix[j][i] = matrix[i][j] / h;
790: g = 0.0;
791: for (k = 0; k <= j; k++)
792: g += matrix[j][k] * matrix[i][k];
793: for (k = j + 1; k <= l; k++)
794: g += matrix[k][j] * matrix[i][k];
795: offdiag[j] = g / h;
796: f += offdiag[j] * matrix[i][j];
797: }
798: hh = f / (h + h);
799: for (j = 0; j <= l; j++) {
800: f = matrix[i][j];
801: offdiag[j] = g = offdiag[j] - hh * f;
802: for (k = 0; k <= j; k++)
803: matrix[j][k] -= f * offdiag[k] + g
804: * matrix[i][k];
805: }
806: }
807: } else
808: offdiag[i] = matrix[i][l];
809: diag[i] = h;
810: }
811: diag[0] = offdiag[0] = 0.0;
812: for (i = 0; i < n; i++) {
813: l = i - 1;
814: if (diag[i] != 0.0) {
815: for (j = 0; j <= l; j++) {
816: g = 0.0;
817: for (k = 0; k <= l; k++)
818: g += matrix[i][k] * matrix[k][j];
819: for (k = 0; k <= l; k++)
820: matrix[k][j] -= g * matrix[k][i];
821: }
822: }
823: diag[i] = matrix[i][i];
824: matrix[i][i] = 1.0;
825: for (j = 0; j <= l; j++)
826: matrix[j][i] = matrix[i][j] = 0.0;
827: }
828: }
829: }
|