KickJava   Java API By Example, From Geeks To Geeks.

Java > Open Source Codes > JSci > maths > LinearMath


1 package JSci.maths;
2
3 import JSci.maths.matrices.*;
4 import JSci.maths.vectors.*;
5 import JSci.maths.polynomials.RealPolynomial;
6
7 /**
8 * The linear math library.
9 * This class cannot be subclassed or instantiated because all methods are static.
10 * @version 2.2
11 * @author Mark Hale
12 */

13 public final class LinearMath extends AbstractMath {
14         private LinearMath() {}
15
16 // LINEAR SYSTEM
17

18         /**
19         * Solves the linear system Mx=v.
20         * @param M a double square matrix.
21         * @param v a double vector.
22         * @return the double vector x.
23         */

24         public static AbstractDoubleVector solve(final AbstractDoubleSquareMatrix M,final AbstractDoubleVector v) {
25                 final int n=v.dimension();
26                 final double array[]=new double[n];
27                 final int pivot[]=new int[n+1];
28                 final AbstractDoubleSquareMatrix lu[]=M.luDecompose(pivot);
29                 int i,j;
30                 double sum;
31         // forward substitution
32
for(i=0;i<n;i++) {
33                         sum=v.getComponent(pivot[i]);
34                         for(j=0;j<i;j++)
35                                 sum-=lu[0].getElement(i,j)*array[j];
36                         array[i]=sum/lu[0].getElement(i,i);
37                 }
38         // back substitution
39
for(i=n-1;i>=0;i--) {
40                         sum=array[i];
41                         for(j=i+1;j<n;j++)
42                                 sum-=lu[1].getElement(i,j)*array[j];
43                         array[i]=sum/lu[1].getElement(i,i);
44                 }
45                 return new DoubleVector(array);
46         }
47         /**
48         * Solves a linear system using cholesky decomposition.
49         */

50         private static double[] solveCholesky(final double m[][], final double v[]) {
51                 final int n=v.length;
52                 final double array[]=new double[n];
53                 final AbstractDoubleSquareMatrix lu[]=new DoubleSquareMatrix(m).choleskyDecompose();
54                 int i,j;
55                 double sum;
56         // forward substitution
57
for(i=0;i<n;i++) {
58                         sum=v[i];
59                         for(j=0;j<i;j++)
60                                 sum-=lu[0].getElement(i,j)*array[j];
61                         array[i]=sum/lu[0].getElement(i,i);
62                 }
63         // back substitution
64
for(i=n-1;i>=0;i--) {
65                         sum=array[i];
66                         for(j=i+1;j<n;j++)
67                                 sum-=lu[1].getElement(i,j)*array[j];
68                         array[i]=sum/lu[1].getElement(i,i);
69                 }
70                 return array;
71         }
72         /**
73         * Solves a linear system using QR decomposition.
74         */

75         private static double[] solveQR(final double m[][], final double v[]) {
76                 final int n=v.length;
77                 final double array[]=new double[n];
78                 final AbstractDoubleSquareMatrix lu[]=new DoubleSquareMatrix(m).qrDecompose();
79                 int i,j;
80                 double sum;
81         // Q.transpose().multiply(v)
82
for(i=0;i<n;i++) {
83                         sum=0.0;
84                         for(j=0;j<n;j++)
85                                 sum+=lu[0].getElement(j,i)*v[j];
86                         array[i]=sum;
87                 }
88         // back substitution
89
for(i=n-1;i>=0;i--) {
90                         sum=array[i];
91                         for(j=i+1;j<n;j++)
92                                 sum-=lu[1].getElement(i,j)*array[j];
93                         array[i]=sum/lu[1].getElement(i,i);
94                 }
95                 return array;
96         }
97         /**
98         * Solves the unsymmetric linear system Ax=b using the
99         * Generalized Minimum Residual method (doesn't require A
100         * to be nonsingular).
101         * While slower than LU decomposition, it is more
102         * robust and should be used with large matrices.
103         * It is guaranted to converge exactly in N iterations for an
104         * N by N matrix (minus some numerical errors).
105         * @author Alain Beliveau
106         * @author Daniel Lemire
107         * @param max_iter maximum number of iterations.
108     * @param tol tolerance.
109         * @exception IllegalArgumentException If either the tolerance
110         * or the number of iterations is not positive.
111         * Also, if an unexpected error occurs.
112         * @exception MaximumIterationsExceededException If it cannot
113         * converge according to the given parameters.
114         */

115         public static AbstractDoubleVector solveGMRes(final AbstractDoubleMatrix A,final AbstractDoubleVector b,int max_iter,double tol) throws MaximumIterationsExceededException {
116                 if(max_iter<=0)
117                         throw new IllegalArgumentException JavaDoc("Number of allowed iterations must be a positive integer: "+max_iter+" <= 0.");
118                 if(tol<0)
119                         throw new IllegalArgumentException JavaDoc("Tolerance must be positive or zero: "+tol+" < 0.");
120                 final int m=A.rows();
121         double resid;
122         int i,j=1,k;
123         double[] s=new double[m+1];
124         double[][] cs=new double[m+1][2];
125         double[] rotmp=new double[2];
126         AbstractDoubleVector w;
127                 AbstractDoubleVector x=new DoubleVector(A.rows());
128         double normb=b.norm();
129                 AbstractDoubleVector r=b.subtract(A.multiply(x));
130         double beta=r.norm();
131         if(normb==0.0)
132                         normb=1.0;
133         if((resid=r.norm()/normb)<= tol) {
134                         tol = resid;
135                 max_iter = 0;
136                 throw new IllegalArgumentException JavaDoc("There is a bug.");
137                 }
138         AbstractDoubleVector[] v = new DoubleVector[m+1];
139                 double[][] H = new double[m+1][m];
140         while (j <= max_iter) {
141                 v[0] = r.scalarMultiply(1.0 / beta);
142                 for (i = 0; i < m+1; i++)
143                     s[i] = 0.0;
144                 s[0] = beta;
145                 for (i = 0; i < m && j <= max_iter; i++, j++) {
146                                 w = A.multiply(v[i]);
147                         for (k = 0; k <= i; k++) {
148                                         H[k][i]=w.scalarProduct(v[k]);
149                         w = w.subtract(v[k].scalarMultiply(H[k][i]));
150                             }
151                         H[i+1][i]=w.norm();
152                         v[i+1] = w.scalarMultiply(1.0 / H[i+1][i]);
153                 for (k = 0; k < i; k++) {
154                             rotmp=applyPlaneRotation(H[k][i], H[k+1][i], cs[k][0], cs[k][1]);
155                             H[k][i]=rotmp[0];
156                             H[k+1][i]=rotmp[1];
157                             }
158                         cs[i]=generatePlaneRotation(H[i][i], H[i+1][i]);
159                         rotmp=applyPlaneRotation(H[i][i], H[i+1][i], cs[i][0], cs[i][1]);
160                         H[i][i]=rotmp[0];
161                         H[i+1][i]=rotmp[1];
162                         rotmp=applyPlaneRotation(s[i], s[i+1], cs[i][0], cs[i][1]);
163                         s[i]=rotmp[0];
164                         s[i+1]=rotmp[1];
165                         if ((resid = Math.abs(s[i+1]) / normb) < tol) {
166                         x = update(x, i, H, s, v);
167                         tol = resid;
168                         max_iter = j;
169                         return (x);
170                         }
171                         }
172                         x = update(x, m - 1, H, s, v);
173                         r = b.subtract(A.multiply(x));
174                         beta = r.norm();
175                         if ((resid = beta / normb) < tol) {
176                         tol = resid;
177                         max_iter = j;
178                         return (x);
179             }
180         }
181         tol = resid;
182         throw new MaximumIterationsExceededException("(tol) "+tol+". It doesn't converge in "+max_iter+"iterations. Try raising the number of allowed iterations or raising the tolerance.");
183     }
184     private static double[] generatePlaneRotation(double dx, double dy) {
185         double [] cs = new double[2];
186         double temp;
187         if (dy == 0.0) {
188                 cs[0] = 1.0;
189                 cs[1] = 0.0;
190             } else if (Math.abs(dy) > Math.abs(dx)) {
191                 temp = dx / dy;
192                 cs[1] = 1.0 / Math.sqrt( 1.0 + temp*temp );
193                 cs[0] = temp * cs[1];
194                 } else {
195                 temp = dy / dx;
196                 cs[0] = 1.0 / Math.sqrt( 1.0 + temp*temp );
197                 cs[1] = temp * cs[0];
198                 }
199             return (cs);
200     }
201     private static double[] applyPlaneRotation(double dx, double dy, double cs, double sn) {
202         double [] dxy = new double[2];
203         dxy[0] = cs * dx + sn * dy;
204         dxy[1] = -sn * dx + cs * dy;
205             return (dxy);
206     }
207     private static AbstractDoubleVector update(AbstractDoubleVector x, int k, double[][] H, double [] s, AbstractDoubleVector [] v) {
208         /************************************************
209             * Backsolve of a triangular + sub diagonal matrix:
210             ************************************************/

211         for (int i = k; i >= 0; i--) {
212                 s[i]=s[i]/H[i][i];
213                 for (int j = i - 1; j >= 0; j--)
214                         s[j]=s[j]-H[j][i]*s[i];
215                 }
216         for (int j = 0; j <= k; j++){
217                 x = x.add(v[j].scalarMultiply(s[j]));
218                 }
219             return(x);
220     }
221
222 // LEAST SQUARES
223

224         /**
225         * Fits an nth degree polynomial to data using the method of least squares.
226         * @param n the degree of the polynomial (>= 0).
227         * @param data
228         * [0][] contains the x-series,
229         * [1][] contains the y-series.
230         */

231         public static RealPolynomial leastSquaresFit(int n, final double data[][]) {
232                 int i,j,nm1=n++;
233                 double xsum,xysum,tmp;
234                 final double mArray[][]=new double[n][n];
235                 final double vArray[]=new double[n];
236                 for(i=0;i<n;i++) {
237                         xsum=xysum=0.0;
238                         for(j=0;j<data[0].length;j++) {
239                                 tmp=Math.pow(data[0][j],i);
240                                 xsum+=tmp;
241                                 xysum+=tmp*data[1][j];
242                         }
243                         mArray[0][i]=xsum;
244                         vArray[i]=xysum;
245                 }
246                 for(i=1;i<n;i++) {
247                         System.arraycopy(mArray[i-1],1,mArray[i],0,nm1);
248                         xsum=0.0;
249                         for(j=0;j<data[0].length;j++)
250                                 xsum+=Math.pow(data[0][j],nm1+i);
251                         mArray[i][nm1]=xsum;
252                 }
253                 return new RealPolynomial(solveCholesky(mArray,vArray));
254         }
255
256 // LINEAR REGRESSION
257

258         /**
259         * Fits a line to multi-dimensional data using the method of least squares.
260         * @param data
261         * [0...n-1][] contains the x-series' (they must be linearly uncorrelated),
262         * [n][] contains the y-series.
263         * @return a vector containing the coefficients (zero component is the intercept, the rest are gradient components).
264         * E.g. y(x1, x2, ...) = coeffs(0) + coeffs(1) * x1 + coeffs(2) * x2 + ...
265         */

266         public static AbstractDoubleVector linearRegression(final double data[][]) {
267                 final int y=data.length-1;
268                 int i,j,k;
269                 double xsum,xysum;
270                 final double mArray[][]=new double[data.length][data.length];
271                 final double vArray[]=new double[data.length];
272                 mArray[0][0]=data[0].length;
273                 for(j=1;j<data.length;j++) {
274                         xsum=0.0;
275                         for(k=0;k<data[0].length;k++)
276                                 xsum+=data[j-1][k];
277                         mArray[0][j]=mArray[j][0]=xsum;
278                 }
279                 xysum=0.0;
280                 for(k=0;k<data[0].length;k++)
281                         xysum+=data[y][k];
282                 vArray[0]=xysum;
283                 for(i=1;i<data.length;i++) {
284                         for(j=i;j<data.length;j++) {
285                                 xsum=0.0;
286                                 for(k=0;k<data[0].length;k++)
287                                         xsum+=data[i-1][k]*data[j-1][k];
288                                 mArray[i][j]=mArray[j][i]=xsum;
289                         }
290                         xysum=0.0;
291                         for(k=0;k<data[0].length;k++)
292                                 xysum+=data[i-1][k]*data[y][k];
293                         vArray[i]=xysum;
294                 }
295                 return new DoubleVector(solveCholesky(mArray,vArray));
296         }
297
298 // GRAM-SCHMIDT
299

300         /**
301         * The Gram-Schmidt orthonormalization method.
302         * @param vecs a set of linearly independent vectors.
303         * @return a set of orthonormal vectors.
304         */

305         public static AbstractDoubleVector[] orthonormalize(AbstractDoubleVector vecs[]) {
306                 final int N=vecs.length;
307                 AbstractDoubleVector orthovecs[]=new DoubleVector[N];
308                 for(int i=0;i<N;i++) {
309                         orthovecs[i]=vecs[i];
310                         for(int j=0;j<i;j++)
311                                 orthovecs[i]=orthovecs[i].subtract(orthovecs[j].scalarMultiply(orthovecs[j].scalarProduct(vecs[i])));
312                         orthovecs[i].normalize();
313                 }
314                 return orthovecs;
315         }
316
317 // EIGENVALUES & EIGENVECTORS
318

319         /**
320         * This method finds the eigenvalues of a Hermitian matrix.
321         * @param matrix a Hermitian matrix.
322         * @return an array containing the eigenvalues.
323         * @exception MaximumIterationsExceededException If it takes too many iterations to determine an eigenvalue.
324         */

325         public static double[] eigenvalueSolveHermitian(final AbstractComplexSquareMatrix matrix) throws MaximumIterationsExceededException {
326                 final int n=matrix.rows();
327                 final double matrix2[][]=new double[2*n][2*n];
328                 double real,imag;
329                 for(int j,i=0;i<n;i++) {
330                         for(j=0;j<n;j++) {
331                                 real=matrix.getRealElement(i,j);
332                                 imag=matrix.getImagElement(i,j);
333                                 matrix2[i][j]=real;
334                                 matrix2[n+i][n+j]=real;
335                                 matrix2[n+i][j]=imag;
336                                 matrix2[i][n+j]=-imag;
337                         }
338                 }
339                 final double eigenvalue2[]=new double[2*n];
340                 final double offdiag[]=new double[2*n];
341                 reduceSymmetric1_SquareToTridiagonal(matrix2,eigenvalue2,offdiag);
342                 System.arraycopy(offdiag,1,offdiag,0,n-1);
343                 offdiag[n-1]=0.0;
344                 eigenvalueSolveSymmetricTridiagonalMatrix(eigenvalue2,offdiag);
345                 final double eigenvalue[]=new double[n];
346                 System.arraycopy(eigenvalue2,0,eigenvalue,0,n);
347                 return eigenvalue;
348         }
349         /**
350         * This method finds the eigenvalues and eigenvectors of a Hermitian matrix.
351         * @param matrix a Hermitian matrix.
352         * @param eigenvector an empty array of complex vectors to hold the eigenvectors.
353         * All eigenvectors will be orthogonal.
354         * @return an array containing the eigenvalues.
355         * @exception MaximumIterationsExceededException If it takes too many iterations to determine an eigenvalue.
356         */

357         public static double[] eigenSolveHermitian(final AbstractComplexSquareMatrix matrix,final AbstractComplexVector eigenvector[]) throws MaximumIterationsExceededException {
358                 final int n=matrix.rows();
359                 final double matrix2[][]=new double[2*n][2*n];
360                 int i,j;
361                 double real,imag;
362                 for(i=0;i<n;i++) {
363                         for(j=0;j<n;j++) {
364                                 real=matrix.getRealElement(i,j);
365                                 imag=matrix.getImagElement(i,j);
366                                 matrix2[i][j]=real;
367                                 matrix2[n+i][n+j]=real;
368                                 matrix2[n+i][j]=imag;
369                                 matrix2[i][n+j]=-imag;
370                         }
371                 }
372                 final double eigenvalue2[]=new double[2*n];
373                 final double offdiag[]=new double[2*n];
374                 reduceSymmetric2_SquareToTridiagonal(matrix2,eigenvalue2,offdiag);
375                 System.arraycopy(offdiag,1,offdiag,0,n-1);
376                 offdiag[n-1]=0.0;
377                 eigenSolveSymmetricTridiagonalMatrix(eigenvalue2,offdiag,matrix2);
378                 final double eigenvalue[]=new double[n];
379                 double arrayRe[],arrayIm[];
380                 for(i=0;i<n;i++) {
381                         eigenvalue[i]=eigenvalue2[i];
382                         arrayRe=new double[n];
383                         arrayIm=new double[n];
384                         for(j=0;j<n;j++) {
385                                 arrayRe[j]=matrix2[j][i];
386                                 arrayIm[j]=matrix2[j+n][i];
387                         }
388                         eigenvector[i]=new ComplexVector(arrayRe,arrayIm);
389                 }
390                 return eigenvalue;
391         }
392         /**
393         * This method finds the eigenvalues of a symmetric tridiagonal matrix by the QL method.
394         * It is based on the NETLIB algol/fortran procedure tql1 by Bowdler, Martin, Reinsch and Wilkinson.
395         * @param matrix a double symmetric tridiagonal matrix.
396         * @return an array containing the eigenvalues.
397         * @exception MaximumIterationsExceededException If it takes too many iterations to determine an eigenvalue.
398         */

399         public static double[] eigenvalueSolveSymmetric(final DoubleTridiagonalMatrix matrix) throws MaximumIterationsExceededException {
400                 final int n=matrix.rows();
401                 final int nm1=n-1;
402                 final double eigenvalue[]=new double[n];
403                 final double offdiag[]=new double[n];
404                 for(int i=0;i<nm1;i++) {
405                         eigenvalue[i]=matrix.getElement(i,i);
406                         offdiag[i]=matrix.getElement(i,i+1);
407                 }
408                 eigenvalue[nm1]=matrix.getElement(nm1,nm1);
409                 offdiag[nm1]=0.0;
410                 eigenvalueSolveSymmetricTridiagonalMatrix(eigenvalue,offdiag);
411                 return eigenvalue;
412         }
413         /**
414         * This method finds the eigenvalues and eigenvectors of a symmetric tridiagonal matrix by the QL method.
415         * It is based on the NETLIB algol/fortran procedure tql2 by Bowdler, Martin, Reinsch and Wilkinson.
416         * @param matrix a double symmetric tridiagonal matrix.
417         * @param eigenvector an empty array of double vectors to hold the eigenvectors.
418         * All eigenvectors will be orthogonal.
419         * @return an array containing the eigenvalues.
420         * @exception MaximumIterationsExceededException If it takes too many iterations to determine an eigenvalue.
421         */

422         public static double[] eigenSolveSymmetric(final DoubleTridiagonalMatrix matrix,final AbstractDoubleVector eigenvector[]) throws MaximumIterationsExceededException {
423                 final int n=matrix.rows();
424                 final int nm1=n-1;
425                 final double eigenvalue[]=new double[n];
426                 final double offdiag[]=new double[n];
427                 final double id[][]=new double[n][n];
428                 int i,j;
429                 for(i=0;i<nm1;i++) {
430                         id[i][i]=1.0;
431                         eigenvalue[i]=matrix.getElement(i,i);
432                         offdiag[i]=matrix.getElement(i,i+1);
433                 }
434                 id[nm1][nm1]=1.0;
435                 eigenvalue[nm1]=matrix.getElement(nm1,nm1);
436                 offdiag[nm1]=0.0;
437                 eigenSolveSymmetricTridiagonalMatrix(eigenvalue,offdiag,id);
438                 for(i=0;i<n;i++) {
439                         DoubleVector evec = new DoubleVector(n);
440                         for(j=0; j<n; j++)
441                                 evec.setComponent(j, id[j][i]);
442                         eigenvector[i] = evec;
443                 }
444                 return eigenvalue;
445         }
446         /**
447         * This method finds the eigenvalues of a symmetric square matrix.
448         * The matrix is reduced to tridiagonal form and then the QL method is applied.
449         * It is based on the NETLIB algol/fortran procedure tred1/tql1 by Bowdler, Martin, Reinsch and Wilkinson.
450         * @param matrix a double symmetric square matrix.
451         * @return an array containing the eigenvalues.
452         * @exception MaximumIterationsExceededException If it takes too many iterations to determine an eigenvalue.
453         */

454         public static double[] eigenvalueSolveSymmetric(final AbstractDoubleSquareMatrix matrix) throws MaximumIterationsExceededException {
455                 final int n=matrix.rows();
456                 final double eigenvalue[]=new double[n];
457                 final double offdiag[]=new double[n];
458                 final double array[][]=new double[n][n];
459                 int i,j;
460                 for(i=0;i<n;i++) {
461                         for(j=0;j<n;j++)
462                                 array[i][j]=matrix.getElement(i,j);
463                 }
464                 reduceSymmetric1_SquareToTridiagonal(array,eigenvalue,offdiag);
465                 System.arraycopy(offdiag,1,offdiag,0,n-1);
466                 offdiag[n-1]=0.0;
467                 eigenvalueSolveSymmetricTridiagonalMatrix(eigenvalue,offdiag);
468                 return eigenvalue;
469         }
470         /**
471         * This method finds the eigenvalues and eigenvectors of a symmetric square matrix.
472         * The matrix is reduced to tridiagonal form and then the QL method is applied.
473         * It is based on the NETLIB algol/fortran procedure tred2/tql2 by Bowdler, Martin, Reinsch and Wilkinson.
474         * @param matrix a double symmetric square matrix.
475         * @param eigenvector an empty array of double vectors to hold the eigenvectors.
476         * All eigenvectors will be orthogonal.
477         * @return an array containing the eigenvalues.
478         * @exception MaximumIterationsExceededException If it takes too many iterations to determine an eigenvalue.
479         */

480         public static double[] eigenSolveSymmetric(final AbstractDoubleSquareMatrix matrix,final AbstractDoubleVector eigenvector[]) throws MaximumIterationsExceededException {
481                 final int n=matrix.rows();
482                 final double eigenvalue[]=new double[n];
483                 final double offdiag[]=new double[n];
484                 final double transf[][]=new double[n][n];
485                 int i,j;
486                 for(i=0;i<n;i++) {
487                         for(j=0;j<n;j++)
488                                 transf[i][j]=matrix.getElement(i,j);
489                 }
490                 reduceSymmetric2_SquareToTridiagonal(transf,eigenvalue,offdiag);
491                 System.arraycopy(offdiag,1,offdiag,0,n-1);
492                 offdiag[n-1]=0.0;
493                 eigenSolveSymmetricTridiagonalMatrix(eigenvalue,offdiag,transf);
494                 for(i=0;i<n;i++) {
495                         DoubleVector evec = new DoubleVector(n);
496                         for(j=0; j<n; j++)
497                                 evec.setComponent(j, transf[j][i]);
498                         eigenvector[i] = evec;
499                 }
500                 return eigenvalue;
501         }
502
503         private final static int EIGEN_MAX_ITERATIONS = 250;
504
505         /**
506         * Internal NETLIB tql1 routine.
507         * @param diag output eigenvalues.
508         * @author Richard Cannings
509         */

510         private static void eigenvalueSolveSymmetricTridiagonalMatrix(final double diag[],final double offdiag[]) throws MaximumIterationsExceededException {
511                 final int n=diag.length;
512                 final int nm1=n-1;
513                 int m,l,iteration,i,k;
514                 double s,r,p,g,f,dd,c,b;
515                 for(l=0;l<n;l++) {
516                         iteration=0;
517                         do {
518                                 for(m=l;m<nm1;m++) {
519                                         dd=Math.abs(diag[m])+Math.abs(diag[m+1]);
520                                         if(Math.abs(offdiag[m])+dd==dd)
521                                                 break;
522                                 }
523                                 if(m!=l) {
524                                         if(iteration++ == EIGEN_MAX_ITERATIONS)
525                                                 throw new MaximumIterationsExceededException("No convergence after "+EIGEN_MAX_ITERATIONS+" iterations.");
526                                         g=(diag[l+1]-diag[l])/(2.0*offdiag[l]);
527                                         r=Math.sqrt(g*g+1.0);
528                                         g=diag[m]-diag[l]+offdiag[l]/(g+(g<0.0?-Math.abs(r):Math.abs(r)));
529                                         s=c=1.0;
530                                         p=0.0;
531                                         for(i=m-1;i>=l;i--) {
532                                                 f=s*offdiag[i];
533                                                 b=c*offdiag[i];
534                                                 if(Math.abs(f)>=Math.abs(g)) {
535                                                         c=g/f;
536                                                         r=Math.sqrt(c*c+1.0);
537                                                         offdiag[i+1]=f*r;
538                                                         s=1/r;
539                                                         c*=s;
540                                                 } else {
541                                                         s=f/g;
542                                                         r=Math.sqrt(s*s+1.0);
543                                                         offdiag[i+1]=g*r;
544                                                         c=1/r;
545                                                         s*=c;
546                                                 }
547                                                 g=diag[i+1]-p;
548                                                 r=(diag[i]-g)*s+2.0*c*b;
549                                                 p=s*r;
550                                                 diag[i+1]=g+p;
551                                                 g=c*r-b;
552                                         }
553                                         diag[l]=diag[l]-p;
554                                         offdiag[l]=g;
555                                         offdiag[m]=0.0;
556                                 }
557                         } while(m!=l);
558                 }
559         }
560         /**
561         * Internal NETLIB tred1 routine.
562         * @author Richard Cannings
563         */

564         private static void reduceSymmetric1_SquareToTridiagonal(final double matrix[][],final double diag[],final double offdiag[]) {
565                 final int n=diag.length;
566                 int i,j,k,l;
567                 double f,g,h,hh,scale;
568                 for(i=n-1;i>0;i--) {
569                         l=i-1;
570                         h=scale=0.0;
571                         if(l>0) {
572                                 for(k=0;k<=l;k++)
573                                         scale+=Math.abs(matrix[i][k]);
574                                 if(scale==0.0)
575                                         offdiag[i]=matrix[i][l];
576                                 else {
577                                         for(k=0;k<=l;k++) {
578                                                 matrix[i][k]/=scale;
579                                                 h+=matrix[i][k]*matrix[i][k];
580                                         }
581                                         f=matrix[i][l];
582                                         g=(f>=0.0?-Math.sqrt(h):Math.sqrt(h));
583                                         offdiag[i]=scale*g;
584                                         h-=f*g;
585                                         matrix[i][l]=f-g;
586                                         f=0.0;
587                                         for(j=0;j<=l;j++) {
588                                                 g=0.0;
589                                                 for(k=0;k<=j;k++)
590                                                         g+=matrix[j][k]*matrix[i][k];
591                                                 for(k=j+1;k<=l;k++)
592                                                         g+=matrix[k][j]*matrix[i][k];
593                                                 offdiag[j]=g/h;
594                                                 f+=offdiag[j]*matrix[i][j];
595                                         }
596                                         hh=f/(h+h);
597                                         for(j=0;j<=l;j++) {
598                                                 f=matrix[i][j];
599                                                 offdiag[j]=g=offdiag[j]-hh*f;
600                                                 for(k=0;k<=j;k++)
601                                                         matrix[j][k]-=f*offdiag[k]+g*matrix[i][k];
602                                         }
603                                 }
604                         } else
605                                 offdiag[i]=matrix[i][l];
606                         diag[i]=h;
607                 }
608                 offdiag[0]=0.0;
609                 for(i=0;i<n;i++)
610                         diag[i]=matrix[i][i];
611         }
612         /**
613         * Internal NETLIB tql2 routine.
614         * @param diag output eigenvalues.
615         * @param transf output eigenvectors.
616         * @author Richard Cannings
617         */

618         private static void eigenSolveSymmetricTridiagonalMatrix(final double diag[],final double offdiag[],final double transf[][]) throws MaximumIterationsExceededException {
619                 final int n=diag.length;
620                 final int nm1=n-1;
621                 int m,l,iteration,i,k;
622                 double s,r,p,g,f,dd,c,b;
623                 for(l=0;l<n;l++) {
624                         iteration=0;
625                         do {
626                                 for(m=l;m<nm1;m++) {
627                                         dd=Math.abs(diag[m])+Math.abs(diag[m+1]);
628                                         if(Math.abs(offdiag[m])+dd==dd)
629                                                 break;
630                                 }
631                                 if(m!=l) {
632                                         if(iteration++ == EIGEN_MAX_ITERATIONS)
633                                                 throw new MaximumIterationsExceededException("No convergence after "+EIGEN_MAX_ITERATIONS+" iterations.");
634                                         g=(diag[l+1]-diag[l])/(2.0*offdiag[l]);
635                                         r=Math.sqrt(g*g+1.0);
636                                         g=diag[m]-diag[l]+offdiag[l]/(g+(g<0.0?-Math.abs(r):Math.abs(r)));
637                                         s=c=1.0;
638                                         p=0.0;
639                                         for(i=m-1;i>=l;i--) {
640                                                 f=s*offdiag[i];
641                                                 b=c*offdiag[i];
642                                                 if(Math.abs(f)>=Math.abs(g)) {
643                                                         c=g/f;
644                                                         r=Math.sqrt(c*c+1.0);
645                                                         offdiag[i+1]=f*r;
646                                                         s=1/r;
647                                                         c*=s;
648                                                 } else {
649                                                         s=f/g;
650                                                         r=Math.sqrt(s*s+1.0);
651                                                         offdiag[i+1]=g*r;
652                                                         c=1/r;
653                                                         s*=c;
654                                                 }
655                                                 g=diag[i+1]-p;
656                                                 r=(diag[i]-g)*s+2.0*c*b;
657                                                 p=s*r;
658                                                 diag[i+1]=g+p;
659                                                 g=c*r-b;
660                                                 for(k=0;k<n;k++) {
661                                                         f=transf[k][i+1];
662                                                         transf[k][i+1]=s*transf[k][i]+c*f;
663                                                         transf[k][i]=c*transf[k][i]-s*f;
664                                                 }
665                                         }
666                                         diag[l]=diag[l]-p;
667                                         offdiag[l]=g;
668                                         offdiag[m]=0.0;
669                                 }
670                         } while(m!=l);
671                 }
672         }
673         /**
674         * Internal NETLIB tred2 routine.
675         * @param matrix output orthogonal transformations.
676         * @author Richard Cannings
677         */

678         private static void reduceSymmetric2_SquareToTridiagonal(final double matrix[][],final double diag[],final double offdiag[]) {
679                 final int n=diag.length;
680                 int i,j,k,l;
681                 double f,g,h,hh,scale;
682                 for(i=n-1;i>0;i--) {
683                         l=i-1;
684                         h=scale=0.0;
685                         if(l>0) {
686                                 for(k=0;k<=l;k++)
687                                         scale+=Math.abs(matrix[i][k]);
688                                 if(scale==0.0)
689                                         offdiag[i]=matrix[i][l];
690                                 else {
691                                         for(k=0;k<=l;k++) {
692                                                 matrix[i][k]/=scale;
693                                                 h+=matrix[i][k]*matrix[i][k];
694                                         }
695                                         f=matrix[i][l];
696                                         g=(f>=0.0?-Math.sqrt(h):Math.sqrt(h));
697                                         offdiag[i]=scale*g;
698                                         h-=f*g;
699                                         matrix[i][l]=f-g;
700                                         f=0.0;
701                                         for(j=0;j<=l;j++) {
702                                                 matrix[j][i]=matrix[i][j]/h;
703                                                 g=0.0;
704                                                 for(k=0;k<=j;k++)
705                                                         g+=matrix[j][k]*matrix[i][k];
706                                                 for(k=j+1;k<=l;k++)
707                                                         g+=matrix[k][j]*matrix[i][k];
708                                                 offdiag[j]=g/h;
709                                                 f+=offdiag[j]*matrix[i][j];
710                                         }
711                                         hh=f/(h+h);
712                                         for(j=0;j<=l;j++) {
713                                                 f=matrix[i][j];
714                                                 offdiag[j]=g=offdiag[j]-hh*f;
715                                                 for(k=0;k<=j;k++)
716                                                         matrix[j][k]-=f*offdiag[k]+g*matrix[i][k];
717                                         }
718                                 }
719                         } else
720                                 offdiag[i]=matrix[i][l];
721                         diag[i]=h;
722                 }
723                 diag[0]=offdiag[0]=0.0;
724                 for(i=0;i<n;i++) {
725                         l=i-1;
726                         if(diag[i]!=0.0) {
727                                 for(j=0;j<=l;j++) {
728                                         g=0.0;
729                                         for(k=0;k<=l;k++)
730                                                 g+=matrix[i][k]*matrix[k][j];
731                                         for(k=0;k<=l;k++)
732                                                 matrix[k][j]-=g*matrix[k][i];
733                                 }
734                         }
735                         diag[i]=matrix[i][i];
736                         matrix[i][i]=1.0;
737                         for(j=0;j<=l;j++)
738                                 matrix[j][i]=matrix[i][j]=0.0;
739                 }
740         }
741 }
742
Popular Tags