KickJava   Java API By Example, From Geeks To Geeks.

Java > Open Source Codes > JSci > maths > matrices > AbstractDoubleSquareMatrix


1 /* AUTO-GENERATED */
2 package JSci.maths.matrices;
3
4 import JSci.GlobalSettings;
5 import JSci.maths.ArrayMath;
6 import JSci.maths.ExtraMath;
7 import JSci.maths.LinearMath;
8 import JSci.maths.MaximumIterationsExceededException;
9 import JSci.maths.Mapping;
10 import JSci.maths.vectors.AbstractDoubleVector;
11 import JSci.maths.groups.AbelianGroup;
12 import JSci.maths.algebras.*;
13 import JSci.maths.fields.*;
14
15 /**
16 * The AbstractDoubleSquareMatrix class provides an object for encapsulating double square matrices.
17 * @version 2.2
18 * @author Mark Hale
19 */
20 public abstract class AbstractDoubleSquareMatrix extends AbstractDoubleMatrix implements SquareMatrix {
21         protected transient AbstractDoubleSquareMatrix[] LU;
22         protected transient int[] LUpivot;
23         /**
24         * Constructs a matrix.
25         */
26         protected AbstractDoubleSquareMatrix(final int size) {
27                 super(size, size);
28         }
29         /**
30         * Converts this matrix to an integer matrix.
31         * @return an integer square matrix
32         */
33         public AbstractIntegerMatrix toIntegerMatrix() {
34                 final int ans[][]=new int[numRows][numCols];
35                 for(int i=0;i<numRows;i++) {
36                         for(int j=0;j<numCols;j++)
37                                 ans[i][j]=Math.round((float)getElement(i,j));
38                 }
39                 return new IntegerSquareMatrix(ans);
40         }
41         /**
42         * Converts this matrix to a complex matrix.
43         * @return a complex square matrix
44         */
45         public AbstractComplexMatrix toComplexMatrix() {
46                 ComplexSquareMatrix cm = new ComplexSquareMatrix(numRows);
47                 for(int i=0;i<numRows;i++) {
48                         for(int j=0;j<numCols;j++)
49                                 cm.setElement(i, j, getElement(i, j), 0.0);
50                 }
51                 return cm;
52         }
53         /**
54         * Returns true if this matrix is symmetric.
55         */
56         public boolean isSymmetric() {
57                 return this.equals(this.transpose());
58         }
59         /**
60         * Returns true if this matrix is unitary.
61         */
62         public boolean isUnitary() {
63                 return this.multiply(this.transpose()).equals(DoubleDiagonalMatrix.identity(numRows));
64         }
65         /**
66         * Returns the determinant.
67         */
68         public double det() {
69                 if(numRows==2) {
70                         return getElement(0,0)*getElement(1,1)-getElement(0,1)*getElement(1,0);
71                 } else {
72                         final AbstractDoubleSquareMatrix lu[]=this.luDecompose(null);
73                         double det=lu[1].getElement(0,0);
74                         for(int i=1;i<numRows;i++)
75                                 det*=lu[1].getElement(i,i);
76                         return det*LUpivot[numRows];
77                 }
78         }
79         /**
80         * Returns the trace.
81         */
82         public double trace() {
83                 double result = getElement(0,0);
84                 for(int i=1;i<numRows;i++)
85                         result += getElement(i,i);
86                 return result;
87         }
88         /**
89         * Returns the operator norm.
90         * @exception MaximumIterationsExceededException If it takes more than 50 iterations to determine an eigenvalue.
91         */
92         public double operatorNorm() throws MaximumIterationsExceededException {
93                 return Math.sqrt(ArrayMath.max(LinearMath.eigenvalueSolveSymmetric((AbstractDoubleSquareMatrix)(this.transpose().multiply(this)))));
94         }
95
96 //============
97 // OPERATIONS
98 //============
99
100         /**
101         * Returns the negative of this matrix.
102         */
103         public AbelianGroup.Member negate() {
104                 final double array[][]=new double[numRows][numCols];
105                 for(int i=0;i<numRows;i++) {
106                         array[i][0] = -getElement(i,0);
107                         for(int j=1;j<numCols;j++)
108                                 array[i][j] = -getElement(i,j);
109                 }
110                 return new DoubleSquareMatrix(array);
111         }
112
113 // ADDITION
114
115         /**
116         * Returns the addition of this matrix and another.
117         * @param m a double square matrix
118         * @exception MatrixDimensionException If the matrices are not square or different sizes.
119         */
120         public final AbstractDoubleMatrix add(final AbstractDoubleMatrix m) {
121                 if(m instanceof AbstractDoubleSquareMatrix)
122                         return add((AbstractDoubleSquareMatrix)m);
123                 else if(numRows == m.rows() && numCols == m.columns())
124                         return add(new SquareMatrixAdaptor(m));
125                 else
126                         throw new MatrixDimensionException("Matrices are different sizes.");
127         }
128         /**
129         * Returns the addition of this matrix and another.
130         * @param m a double square matrix
131         * @exception MatrixDimensionException If the matrices are different sizes.
132         */
133         public AbstractDoubleSquareMatrix add(final AbstractDoubleSquareMatrix m) {
134                 if(numRows==m.rows() && numCols==m.columns()) {
135                         final double array[][]=new double[numRows][numCols];
136                         for(int i=0;i<numRows;i++) {
137                                 array[i][0] = getElement(i,0)+m.getElement(i,0);
138                                 for(int j=1;j<numCols;j++)
139                                         array[i][j] = getElement(i,j)+m.getElement(i,j);
140                         }
141                         return new DoubleSquareMatrix(array);
142                 } else {
143                         throw new MatrixDimensionException("Matrices are different sizes.");
144                 }
145         }
146
147 // SUBTRACTION
148
149         /**
150         * Returns the subtraction of this matrix and another.
151         * @param m a double square matrix
152         * @exception MatrixDimensionException If the matrices are not square or different sizes.
153         */
154         public final AbstractDoubleMatrix subtract(final AbstractDoubleMatrix m) {
155                 if(m instanceof AbstractDoubleSquareMatrix)
156                         return subtract((AbstractDoubleSquareMatrix)m);
157                 else if(numRows == m.rows() && numCols == m.columns())
158                         return subtract(new SquareMatrixAdaptor(m));
159                 else
160                         throw new MatrixDimensionException("Matrices are different sizes.");
161         }
162         /**
163         * Returns the subtraction of this matrix by another.
164         * @param m a double square matrix
165         * @exception MatrixDimensionException If the matrices are different sizes.
166         */
167         public AbstractDoubleSquareMatrix subtract(final AbstractDoubleSquareMatrix m) {
168                 if(numRows==m.rows() && numCols==m.columns()) {
169                         final double array[][]=new double[numRows][numCols];
170                         for(int i=0;i<numRows;i++) {
171                                 array[i][0] = getElement(i,0)-m.getElement(i,0);
172                                 for(int j=1;j<numCols;j++)
173                                         array[i][j] = getElement(i,j)-m.getElement(i,j);
174                         }
175                         return new DoubleSquareMatrix(array);
176                 } else {
177                         throw new MatrixDimensionException("Matrices are different sizes.");
178                 }
179         }
180
181 // SCALAR MULTIPLICATION
182
183         /**
184         * Returns the multiplication of this matrix by a scalar.
185         * @param x a double.
186         * @return a double square matrix.
187         */
188         public AbstractDoubleMatrix scalarMultiply(final double x) {
189                 final double array[][]=new double[numRows][numCols];
190                 for(int i=0;i<numRows;i++) {
191                         array[i][0] = x*getElement(i,0);
192                         for(int j=1;j<numCols;j++)
193                                 array[i][j] = x*getElement(i,j);
194                 }
195                 return new DoubleSquareMatrix(array);
196         }
197
198 // SCALAR DIVISON
199
200         /**
201         * Returns the division of this matrix by a scalar.
202         * @param x a double.
203         * @return a double square matrix.
204         */
205         public AbstractDoubleMatrix scalarDivide(final double x) {
206                 final double array[][]=new double[numRows][numCols];
207                 for(int i=0;i<numRows;i++) {
208                         array[i][0] = getElement(i,0)/x;
209                         for(int j=1;j<numCols;j++)
210                                 array[i][j] = getElement(i,j)/x;
211                 }
212                 return new DoubleSquareMatrix(array);
213         }
214
215 // SCALAR PRODUCT
216
217         /**
218         * Returns the scalar product of this matrix and another.
219         * @param m a double square matrix.
220         * @exception MatrixDimensionException If the matrices are not square or different sizes.
221         */
222         public final double scalarProduct(final AbstractDoubleMatrix m) {
223                 if(m instanceof AbstractDoubleSquareMatrix)
224                         return scalarProduct((AbstractDoubleSquareMatrix)m);
225                 else if(numRows == m.rows() && numCols == m.columns())
226                         return scalarProduct(new SquareMatrixAdaptor(m));
227                 else
228                         throw new MatrixDimensionException("Matrices are different sizes.");
229         }
230         /**
231         * Returns the scalar product of this matrix and another.
232         * @param m a double square matrix.
233         * @exception MatrixDimensionException If the matrices are different sizes.
234         */
235         public double scalarProduct(final AbstractDoubleSquareMatrix m) {
236                 if(numRows==m.rows() && numCols==m.columns()) {
237                         double ans = 0;
238                         for(int i=0; i<numRows; i++) {
239                                 ans += getElement(i,0)*m.getElement(i,0);
240                                 for(int j=1; j<numCols; j++)
241                                         ans += getElement(i,j)*m.getElement(i,j);
242                         }
243                         return ans;
244                 } else {
245                        throw new MatrixDimensionException("Matrices are different sizes.");
246                 }
247         }
248
249 // MATRIX MULTIPLICATION
250
251         /**
252         * Returns the multiplication of this matrix and another.
253         * @param m a double square matrix
254         * @exception MatrixDimensionException If the matrices are different sizes.
255         */
256         public AbstractDoubleSquareMatrix multiply(final AbstractDoubleSquareMatrix m) {
257                 if(numCols==m.rows()) {
258                         final int mColumns = m.columns();
259                         final double array[][]=new double[numRows][mColumns];
260                         for(int j=0; j<numRows; j++) {
261                                 for(int k=0; k<mColumns; k++) {
262                                         array[j][k] = getElement(j,0)*m.getElement(0,k);
263                                         for(int n=1; n<numCols; n++)
264                                                 array[j][k] += getElement(j,n)*m.getElement(n,k);
265                                 }
266                         }
267                         return new DoubleSquareMatrix(array);
268                 } else {
269                         throw new MatrixDimensionException("Incompatible matrices.");
270                 }
271         }
272
273 // DIRECT SUM
274
275         /**
276         * Returns the direct sum of this matrix and another.
277         */
278         public AbstractDoubleSquareMatrix directSum(final AbstractDoubleSquareMatrix m) {
279                 final double array[][]=new double[numRows+m.numRows][numCols+m.numCols];
280                 for(int i=0;i<numRows;i++) {
281                         for(int j=0;j<numCols;j++)
282                                 array[i][j] = getElement(i,j);
283                 }
284                 for(int i=0;i<m.numRows;i++) {
285                         for(int j=0;j<m.numCols;j++)
286                                 array[i+numRows][j+numCols] = m.getElement(i,j);
287                 }
288                 return new DoubleSquareMatrix(array);
289         }
290
291 // TENSOR PRODUCT
292
293         /**
294         * Returns the tensor product of this matrix and another.
295         */
296         public AbstractDoubleSquareMatrix tensor(final AbstractDoubleSquareMatrix m) {
297                 final double array[][]=new double[numRows*m.numRows][numCols*m.numCols];
298                 for(int i=0;i<numRows;i++) {
299                         for(int j=0;j<numCols;j++) {
300                                 for(int k=0;k<m.numRows;j++) {
301                                         for(int l=0;l<m.numCols;l++)
302                                                 array[i*m.numRows+k][j*m.numCols+l] = getElement(i,j)*m.getElement(k,l);
303                                 }
304                         }
305                 }
306                 return new DoubleSquareMatrix(array);
307         }
308
309 // TRANSPOSE
310
311         /**
312         * Returns the transpose of this matrix.
313         * @return a double square matrix
314         */
315         public Matrix transpose() {
316                 final double array[][]=new double[numCols][numRows];
317                 for(int i=0;i<numRows;i++) {
318                         array[0][i] = getElement(i,0);
319                         for(int j=1;j<numCols;j++)
320                                 array[j][i] = getElement(i,j);
321                 }
322                 return new DoubleSquareMatrix(array);
323         }
324
325 // INVERSE
326
327         /**
328         * Returns the inverse of this matrix.
329         * @return a double square matrix.
330         */
331         public AbstractDoubleSquareMatrix inverse() {
332                 final int N=numRows;
333                 final double arrayL[][]=new double[N][N];
334                 final double arrayU[][]=new double[N][N];
335                 final AbstractDoubleSquareMatrix lu[]=this.luDecompose(null);
336                 arrayL[0][0]=1.0/lu[0].getElement(0,0);
337                 arrayU[0][0]=1.0/lu[1].getElement(0,0);
338                 for(int i=1;i<N;i++) {
339                         arrayL[i][i]=1.0/lu[0].getElement(i,i);
340                         arrayU[i][i]=1.0/lu[1].getElement(i,i);
341                 }
342                 for(int i=0;i<N-1;i++) {
343                         for(int j=i+1;j<N;j++) {
344                                 double tmpL=0.0, tmpU=0.0;
345                                 for(int k=i;k<j;k++) {
346                                         tmpL-=lu[0].getElement(j,k)*arrayL[k][i];
347                                         tmpU-=arrayU[i][k]*lu[1].getElement(k,j);
348                                 }
349                                 arrayL[j][i]=tmpL/lu[0].getElement(j,j);
350                                 arrayU[i][j]=tmpU/lu[1].getElement(j,j);
351                         }
352                 }
353                 // matrix multiply arrayU x arrayL
354 final double inv[][]=new double[N][N];
355                 for(int i=0;i<N;i++) {
356                         for(int j=0;j<i;j++) {
357                                 for(int k=i;k<N;k++)
358                                         inv[i][LUpivot[j]]+=arrayU[i][k]*arrayL[k][j];
359                         }
360                         for(int j=i;j<N;j++) {
361                                 for(int k=j;k<N;k++)
362                                         inv[i][LUpivot[j]]+=arrayU[i][k]*arrayL[k][j];
363                         }
364                 }
365                 return new DoubleSquareMatrix(inv);
366         }
367
368 // LU DECOMPOSITION
369
370         /**
371         * Returns the LU decomposition of this matrix.
372         * @param pivot an empty array of length <code>rows()+1</code>
373         * to hold the pivot information (null if not interested).
374         * The last array element will contain the parity.
375         * @return an array with [0] containing the L-matrix
376         * and [1] containing the U-matrix.
377         * @jsci.planetmath LUDecomposition
378         */
379         public AbstractDoubleSquareMatrix[] luDecompose(int pivot[]) {
380                 if(LU!=null) {
381                         if(pivot!=null)
382                                 System.arraycopy(LUpivot,0,pivot,0,pivot.length);
383                         return LU;
384                 }
385                 int pivotrow;
386                 final int N=numRows;
387                 final double arrayL[][]=new double[N][N];
388                 final double arrayU[][]=new double[N][N];
389                 if(pivot==null)
390                         pivot=new int[N+1];
391                 for(int i=0;i<N;i++)
392                         pivot[i]=i;
393                 pivot[N]=1;
394         // LU decomposition to arrayU
395 for(int j=0;j<N;j++) {
396                         for(int i=0;i<j;i++) {
397                                 double tmp=getElement(pivot[i],j);
398                                 for(int k=0;k<i;k++)
399                                         tmp-=arrayU[i][k]*arrayU[k][j];
400                                 arrayU[i][j]=tmp;
401                         }
402                         double max=0.0;
403                         pivotrow=j;
404                         for(int i=j;i<N;i++) {
405                                 double tmp=getElement(pivot[i],j);
406                                 for(int k=0;k<j;k++)
407                                         tmp-=arrayU[i][k]*arrayU[k][j];
408                                 arrayU[i][j]=tmp;
409                         // while we're here search for a pivot for arrayU[j][j]
410 tmp=Math.abs(tmp);
411                                 if(tmp>max) {
412                                         max=tmp;
413                                         pivotrow=i;
414                                 }
415                         }
416                 // swap row j with pivotrow
417 if(pivotrow!=j) {
418                                 double[] tmprow = arrayU[j];
419                                 arrayU[j] = arrayU[pivotrow];
420                                 arrayU[pivotrow] = tmprow;
421                                 int k=pivot[j];
422                                 pivot[j]=pivot[pivotrow];
423                                 pivot[pivotrow]=k;
424                                 // update parity
425 pivot[N]=-pivot[N];
426                         }
427                 // divide by pivot
428 double tmp=arrayU[j][j];
429                         for(int i=j+1;i<N;i++)
430                                 arrayU[i][j]/=tmp;
431                 }
432                 // move lower triangular part to arrayL
433 for(int j=0;j<N;j++) {
434                         arrayL[j][j]=1.0;
435                         for(int i=j+1;i<N;i++) {
436                                 arrayL[i][j]=arrayU[i][j];
437                                 arrayU[i][j]=0.0;
438                         }
439                 }
440                 LU=new AbstractDoubleSquareMatrix[2];
441                 LU[0]=new DoubleSquareMatrix(arrayL);
442                 LU[1]=new DoubleSquareMatrix(arrayU);
443                 LUpivot=new int[pivot.length];
444                 System.arraycopy(pivot,0,LUpivot,0,pivot.length);
445                 return LU;
446         }
447         /**
448         * Returns the LU decomposition of this matrix.
449         * Warning: no pivoting.
450         * @return an array with [0] containing the L-matrix
451         * and [1] containing the U-matrix.
452         * @jsci.planetmath LUDecomposition
453         */
454         public AbstractDoubleSquareMatrix[] luDecompose() {
455                 final int N=numRows;
456                 final double arrayL[][]=new double[N][N];
457                 final double arrayU[][]=new double[N][N];
458         // LU decomposition to arrayU
459 for(int j=0;j<N;j++) {
460                         for(int i=0;i<j;i++) {
461                                 double tmp=getElement(i,j);
462                                 for(int k=0;k<i;k++)
463                                         tmp-=arrayU[i][k]*arrayU[k][j];
464                                 arrayU[i][j]=tmp;
465                         }
466                         for(int i=j;i<N;i++) {
467                                 double tmp=getElement(i,j);
468                                 for(int k=0;k<j;k++)
469                                         tmp-=arrayU[i][k]*arrayU[k][j];
470                                 arrayU[i][j]=tmp;
471                         }
472                 // divide
473 double tmp=arrayU[j][j];
474                         for(int i=j+1;i<N;i++)
475                                 arrayU[i][j]/=tmp;
476                 }
477                 // move lower triangular part to arrayL
478 for(int j=0;j<N;j++) {
479                         arrayL[j][j]=1.0;
480                         for(int i=j+1;i<N;i++) {
481                                 arrayL[i][j]=arrayU[i][j];
482                                 arrayU[i][j]=0.0;
483                         }
484                 }
485                 AbstractDoubleSquareMatrix[] lu=new AbstractDoubleSquareMatrix[2];
486                 lu[0]=new DoubleSquareMatrix(arrayL);
487                 lu[1]=new DoubleSquareMatrix(arrayU);
488                 return lu;
489         }
490
491 // CHOLESKY DECOMPOSITION
492
493         /**
494         * Returns the Cholesky decomposition of this matrix.
495         * Matrix must be symmetric and positive definite.
496         * @return an array with [0] containing the L-matrix and [1] containing the U-matrix.
497         * @jsci.planetmath CholeskyDecomposition
498         */
499         public AbstractDoubleSquareMatrix[] choleskyDecompose() {
500                 final int N=numRows;
501                 final double arrayL[][]=new double[N][N];
502                 final double arrayU[][]=new double[N][N];
503                 double tmp=Math.sqrt(getElement(0,0));
504                 arrayL[0][0]=arrayU[0][0]=tmp;
505                 for(int i=1;i<N;i++)
506                         arrayL[i][0]=arrayU[0][i]=getElement(i,0)/tmp;
507                 for(int j=1;j<N;j++) {
508                         tmp=getElement(j,j);
509                         for(int i=0;i<j;i++)
510                                 tmp-=arrayL[j][i]*arrayL[j][i];
511                         arrayL[j][j]=arrayU[j][j]=Math.sqrt(tmp);
512                         for(int i=j+1;i<N;i++) {
513                                 tmp=getElement(i,j);
514                                 for(int k=0;k<i;k++)
515                                         tmp-=arrayL[j][k]*arrayU[k][i];
516                                 arrayL[i][j]=arrayU[j][i]=tmp/arrayU[j][j];
517                         }
518                 }
519                 final AbstractDoubleSquareMatrix lu[]=new AbstractDoubleSquareMatrix[2];
520                 lu[0]=new DoubleSquareMatrix(arrayL);
521                 lu[1]=new DoubleSquareMatrix(arrayU);
522                 return lu;
523         }
524
525 // QR DECOMPOSITION
526
527         /**
528         * Returns the QR decomposition of this matrix.
529         * Based on the code from <a HREF="http://math.nist.gov/javanumerics/jama/">JAMA</a> (public domain).
530         * @return an array with [0] containing the Q-matrix and [1] containing the R-matrix.
531         * @jsci.planetmath QRDecomposition
532         */
533         public AbstractDoubleSquareMatrix[] qrDecompose() {
534                 final int N=numRows;
535                 final double array[][]=new double[N][N];
536                 final double arrayQ[][]=new double[N][N];
537                 final double arrayR[][]=new double[N][N];
538                 // copy matrix
539 for(int i=0;i<N;i++) {
540                         array[i][0]=getElement(i,0);
541                         for(int j=1;j<N;j++)
542                                 array[i][j]=getElement(i,j);
543                 }
544                 for(int k=0; k<N; k++) {
545                         // compute l2-norm of kth column
546 double norm = array[k][k];
547                         for(int i=k+1; i<N; i++)
548                                 norm = ExtraMath.hypot(norm, array[i][k]);
549                         if(norm != 0.0) {
550                                 // form kth Householder vector
551 if(array[k][k] < 0.0)
552                                         norm = -norm;
553                                 for(int i=k; i<N; i++)
554                                         array[i][k] /= norm;
555                                 array[k][k] += 1.0;
556                                 // apply transformation to remaining columns
557 for(int j=k+1; j<N; j++) {
558                                         double s=array[k][k]*array[k][j];
559                                         for(int i=k+1; i<N; i++)
560                                                 s += array[i][k]*array[i][j];
561                                         s /= array[k][k];
562                                         for(int i=k; i<N; i++)
563                                                 array[i][j] -= s*array[i][k];
564                                 }
565                         }
566                         arrayR[k][k] = -norm;
567                 }
568                 for(int k=N-1; k>=0; k--) {
569                         arrayQ[k][k] = 1.0;
570                         for(int j=k; j<N; j++) {
571                                 if(array[k][k] != 0.0) {
572                                         double s = array[k][k]*arrayQ[k][j];
573                                         for(int i=k+1; i<N; i++)
574                                                 s += array[i][k]*arrayQ[i][j];
575                                         s /= array[k][k];
576                                         for(int i=k; i<N; i++)
577                                                 arrayQ[i][j] -= s*array[i][k];
578                                 }
579                         }
580                 }
581                 for(int i=0; i<N; i++) {
582                         for(int j=i+1; j<N; j++)
583                                 arrayR[i][j] = array[i][j];
584                 }
585                 final AbstractDoubleSquareMatrix qr[]=new AbstractDoubleSquareMatrix[2];
586                 qr[0]=new DoubleSquareMatrix(arrayQ);
587                 qr[1]=new DoubleSquareMatrix(arrayR);
588                 return qr;
589         }
590
591 // SINGULAR VALUE DECOMPOSITION
592
593         /**
594         * Returns the singular value decomposition of this matrix.
595         * Based on the code from <a HREF="http://math.nist.gov/javanumerics/jama/">JAMA</a> (public domain).
596         * @return an array with [0] containing the U-matrix, [1] containing the S-matrix and [2] containing the V-matrix.
597         * @jsci.planetmath SingularValueDecomposition
598         */
599         public AbstractDoubleSquareMatrix[] singularValueDecompose() {
600                 final int N=numRows;
601                 final int Nm1=N-1;
602                 final double array[][]=new double[N][N];
603                 final double arrayU[][]=new double[N][N];
604                 final double arrayS[]=new double[N];
605                 final double arrayV[][]=new double[N][N];
606                 final double e[]=new double[N];
607                 final double work[]=new double[N];
608                 // copy matrix
609 for(int i=0;i<N;i++) {
610                         array[i][0]=getElement(i,0);
611                         for(int j=1;j<N;j++)
612                                 array[i][j]=getElement(i,j);
613                 }
614                 // reduce matrix to bidiagonal form
615 for(int k=0;k<Nm1;k++) {
616                         // compute the transformation for the kth column
617 // compute l2-norm of kth column
618 arrayS[k]=array[k][k];
619                         for(int i=k+1;i<N;i++)
620                                 arrayS[k]=ExtraMath.hypot(arrayS[k],array[i][k]);
621                         if(arrayS[k]!=0.0) {
622                                 if(array[k][k]<0.0)
623                                         arrayS[k]=-arrayS[k];
624                                 for(int i=k;i<N;i++)
625                                         array[i][k]/=arrayS[k];
626                                 array[k][k]+=1.0;
627                         }
628                         arrayS[k]=-arrayS[k];
629                         for(int j=k+1;j<N;j++) {
630                                 if(arrayS[k]!=0.0) {
631                                         // apply the transformation
632 double t=array[k][k]*array[k][j];
633                                         for(int i=k+1; i<N; i++)
634                                                 t+=array[i][k]*array[i][j];
635                                         t /= array[k][k];
636                                         for(int i=k; i<N; i++)
637                                                 array[i][j] -= t*array[i][k];
638                                 }
639                                 e[j]=array[k][j];
640                         }
641                         for(int i=k;i<N;i++)
642                                 arrayU[i][k]=array[i][k];
643                         if(k<N-2) {
644                                 // compute the kth row transformation
645 // compute l2-norm of kth column
646 e[k]=e[k+1];
647                                 for(int i=k+2;i<N;i++)
648                                         e[k]=ExtraMath.hypot(e[k],e[i]);
649                                 if(e[k]!=0.0) {
650                                         if(e[k+1]<0.0)
651                                                 e[k]=-e[k];
652                                         for(int i=k+1;i<N;i++)
653                                                 e[i]/=e[k];
654                                         e[k+1]+=1.0;
655                                 }
656                                 e[k]=-e[k];
657                                 if(e[k]!=0.0) {
658                                         // apply the transformation
659 for(int i=k+1;i<N;i++) {
660                                                 work[i]=0.0;
661                                                 for(int j=k+1;j<N;j++)
662                                                         work[i]+=e[j]*array[i][j];
663                                         }
664                                         for(int j=k+1;j<N;j++) {
665                                                 double t = e[j]/e[k+1];
666                                                 for(int i=k+1;i<N;i++)
667                                                         array[i][j] -= t*work[i];
668                                         }
669                                 }
670                                 for(int i=k+1;i<N;i++)
671                                         arrayV[i][k]=e[i];
672                         }
673                 }
674                 // setup the final bidiagonal matrix of order p
675 int p=N;
676                 arrayS[Nm1]=array[Nm1][Nm1];
677                 e[N-2]=array[N-2][Nm1];
678                 e[Nm1]=0.0;
679                 arrayU[Nm1][Nm1]=1.0;
680                 for(int k=N-2;k>=0;k--) {
681                         if(arrayS[k]!=0.0) {
682                                 for(int j=k+1;j<N;j++) {
683                                         double t=arrayU[k][k]*arrayU[k][j];
684                                         for(int i=k+1;i<N;i++)
685                                                 t+=arrayU[i][k]*arrayU[i][j];
686                                         t /= arrayU[k][k];
687                                         for(int i=k;i<N;i++)
688                                                 arrayU[i][j] -= t*arrayU[i][k];
689                                 }
690                                 for(int i=k;i<N;i++)
691                                         arrayU[i][k]=-arrayU[i][k];
692                                 arrayU[k][k]+=1.0;
693                                 for(int i=0;i<k-1;i++)
694                                         arrayU[i][k]=0.0;
695                         } else {
696                                 for(int i=0;i<N;i++)
697                                         arrayU[i][k]=0.0;
698                                 arrayU[k][k]=1.0;
699                         }
700                 }
701                 for(int k=Nm1;k>=0;k--) {
702                         if(k<N-2 && e[k]!=0.0) {
703                                 for(int j=k+1;j<N;j++) {
704                                         double t=arrayV[k+1][k]*arrayV[k+1][j];
705                                         for(int i=k+2;i<N;i++)
706                                                 t+=arrayV[i][k]*arrayV[i][j];
707                                         t /= arrayV[k+1][k];
708                                         for(int i=k+1;i<N;i++)
709                                                 arrayV[i][j] -= t*arrayV[i][k];
710                                 }
711                         }
712                         for(int i=0;i<N;i++)
713                                 arrayV[i][k]=0.0;
714                         arrayV[k][k]=1.0;
715                 }
716                 final double eps=Math.pow(2.0,-52.0);
717                 int iter=0;
718                 while(p>0) {
719                         int k,action;
720                         // action = 1 if arrayS[p] and e[k-1] are negligible and k<p
721 // action = 2 if arrayS[k] is negligible and k<p
722 // action = 3 if e[k-1] is negligible, k<p, and arrayS[k], ..., arrayS[p] are not negligible (QR step)
723 // action = 4 if e[p-1] is negligible (convergence)
724 for(k=p-2;k>=-1;k--) {
725                                 if(k==-1)
726                                         break;
727                                 if(Math.abs(e[k])<=eps*(Math.abs(arrayS[k])+Math.abs(arrayS[k+1]))) {
728                                         e[k]=0.0;
729                                         break;
730                                 }
731                         }
732                         if(k==p-2) {
733                                 action=4;
734                         } else {
735                                 int ks;
736                                 for(ks=p-1;ks>=k;ks--) {
737                                         if(ks==k)
738                                                 break;
739                                         double t=(ks!=p ? Math.abs(e[ks]) : 0.0)+(ks!=k+1 ? Math.abs(e[ks-1]) : 0.0);
740                                         if(Math.abs(arrayS[ks])<=eps*t) {
741                                                 arrayS[ks]=0.0;
742                                                 break;
743                                         }
744                                 }
745                                 if(ks==k) {
746                                         action=3;
747                                 } else if(ks==p-1) {
748                                         action=1;
749                                 } else {
750                                         action=2;
751                                         k=ks;
752                                 }
753                         }
754                         k++;
755                         switch(action) {
756                                 // deflate negligible arrayS[p]
757 case 1: {
758                                         double f=e[p-2];
759                                         e[p-2]=0.0;
760                                         for(int j=p-2;j>=k;j--) {
761                                                 double t=ExtraMath.hypot(arrayS[j],f);
762                                                 final double cs=arrayS[j]/t;
763                                                 final double sn=f/t;
764                                                 arrayS[j]=t;
765                                                 if(j!=k) {
766                                                         f=-sn*e[j-1];
767                                                         e[j-1]*=cs;
768                                                 }
769                                                 for(int i=0;i<N;i++) {
770                                                         t=cs*arrayV[i][j]+sn*arrayV[i][p-1];
771                                                         arrayV[i][p-1]=-sn*arrayV[i][j]+cs*arrayV[i][p-1];
772                                                         arrayV[i][j]=t;
773                                                 }
774                                         }
775                                         } break;
776                                 // split at negligible arrayS[k]
777 case 2: {
778                                         double f=e[k-1];
779                                         e[k-1]=0.0;
780                                         for(int j=k;j<p;j++) {
781                                                 double t=ExtraMath.hypot(arrayS[j],f);
782                                                 final double cs=arrayS[j]/t;
783                                                 final double sn=f/t;
784                                                 arrayS[j]=t;
785                                                 f=-sn*e[j];
786                                                 e[j]*=cs;
787                                                 for(int i=0;i<N;i++) {
788                                                         t=cs*arrayU[i][j]+sn*arrayU[i][k-1];
789                                                         arrayU[i][k-1]=-sn*arrayU[i][j]+cs*arrayU[i][k-1];
790                                                         arrayU[i][j]=t;
791                                                 }
792                                         }
793                                         } break;
794                                 // perform one QR step
795 case 3: {
796                                         // calculate the shift
797 final double scale=Math.max(Math.max(Math.max(Math.max(
798                                                 Math.abs(arrayS[p-1]),Math.abs(arrayS[p-2])),Math.abs(e[p-2])),
799                                                 Math.abs(arrayS[k])),Math.abs(e[k]));
800                                         double sp=arrayS[p-1]/scale;
801                                         double spm1=arrayS[p-2]/scale;
802                                         double epm1=e[p-2]/scale;
803                                         double sk=arrayS[k]/scale;
804                                         double ek=e[k]/scale;
805                                         double b=((spm1+sp)*(spm1-sp)+epm1*epm1)/2.0;
806                                         double c=(sp*epm1)*(sp*epm1);
807                                         double shift=0.0;
808                                         if(b!=0.0 || c!=0.0) {
809                                                 shift=Math.sqrt(b*b+c);
810                                                 if(b<0.0)
811                                                         shift=-shift;
812                                                 shift=c/(b+shift);
813                                         }
814                                         double f=(sk+sp)*(sk-sp)+shift;
815                                         double g=sk*ek;
816                                         // chase zeros
817 for(int j=k;j<p-1;j++) {
818                                                 double t=ExtraMath.hypot(f,g);
819                                                 double cs=f/t;
820                                                 double sn=g/t;
821                                                 if(j!=k)
822                                                         e[j-1]=t;
823                                                 f=cs*arrayS[j]+sn*e[j];
824                                                 e[j]=cs*e[j]-sn*arrayS[j];
825                                                 g=sn*arrayS[j+1];
826                                                 arrayS[j+1]*=cs;
827                                                 for(int i=0;i<N;i++) {
828                                                         t=cs*arrayV[i][j]+sn*arrayV[i][j+1];
829                                                         arrayV[i][j+1]=-sn*arrayV[i][j]+cs*arrayV[i][j+1];
830                                                         arrayV[i][j]=t;
831                                                 }
832                                                 t=ExtraMath.hypot(f,g);
833                                                 cs=f/t;
834                                                 sn=g/t;
835                                                 arrayS[j]=t;
836                                                 f=cs*e[j]+sn*arrayS[j+1];
837                                                 arrayS[j+1]=-sn*e[j]+cs*arrayS[j+1];
838                                                 g=sn*e[j+1];
839                                                 e[j+1]*=cs;
840                                                 if(j<Nm1) {
841                                                         for(int i=0;i<N;i++) {
842                                                                 t=cs*arrayU[i][j]+sn*arrayU[i][j+1];
843                                                                 arrayU[i][j+1]=-sn*arrayU[i][j]+cs*arrayU[i][j+1];
844                                                                 arrayU[i][j]=t;
845                                                         }
846                                                 }
847                                         }
848                                         e[p-2]=f;
849                                         iter++;
850                                         } break;
851                                 // convergence
852 case 4: {
853                                         // make the singular values positive
854 if(arrayS[k]<=0.0) {
855                                                 arrayS[k]=-arrayS[k];
856                                                 for(int i=0;i<p;i++)
857                                                         arrayV[i][k]=-arrayV[i][k];
858                                         }
859                                         // order the singular values
860 while(k<p-1) {
861                                                 if(arrayS[k]>=arrayS[k+1])
862                                                         break;
863                                                 double tmp=arrayS[k];
864                                                 arrayS[k]=arrayS[k+1];
865                                                 arrayS[k+1]=tmp;
866                                                 if(k<Nm1) {
867                                                         for(int i=0;i<N;i++) {
868                                                                 tmp=arrayU[i][k+1];
869                                                                 arrayU[i][k+1]=arrayU[i][k];
870                                                                 arrayU[i][k]=tmp;
871                                                                 tmp=arrayV[i][k+1];
872                                                                 arrayV[i][k+1]=arrayV[i][k];
873                                                                 arrayV[i][k]=tmp;
874                                                         }
875                                                 }
876                                                 k++;
877                                         }
878                                         iter=0;
879                                         p--;
880                                         } break;
881                         }
882                 }
883                 final AbstractDoubleSquareMatrix svd[]=new AbstractDoubleSquareMatrix[3];
884                 svd[0]=new DoubleSquareMatrix(arrayU);
885                 svd[1]=new DoubleDiagonalMatrix(arrayS);
886                 svd[2]=new DoubleSquareMatrix(arrayV);
887                 return svd;
888         }
889
890 // POLAR DECOMPOSITION
891
892         /**
893         * Returns the polar decomposition of this matrix.
894         */
895         public AbstractDoubleSquareMatrix[] polarDecompose() {
896                 final int N=numRows;
897                 final AbstractDoubleVector evec[]=new AbstractDoubleVector[N];
898                 double eval[];
899                 try {
900                         eval=LinearMath.eigenSolveSymmetric(this,evec);
901                 } catch(MaximumIterationsExceededException e) {
902                         return null;
903                 }
904                 final double tmpa[][]=new double[N][N];
905                 final double tmpm[][]=new double[N][N];
906                 double abs;
907                 for(int i=0;i<N;i++) {
908                         abs=Math.abs(eval[i]);
909                         tmpa[i][0]=eval[i]*evec[i].getComponent(0)/abs;
910                         tmpm[i][0]=abs*evec[i].getComponent(0);
911                         for(int j=1;j<N;j++) {
912                                 tmpa[i][j]=eval[i]*evec[i].getComponent(j)/abs;
913                                 tmpm[i][j]=abs*evec[i].getComponent(j);
914                         }
915                 }
916                 final double arg[][]=new double[N][N];
917                 final double mod[][]=new double[N][N];
918                 for(int i=0;i<N;i++) {
919                         for(int j=0;j<N;j++) {
920                                 arg[i][j]=evec[0].getComponent(i)*tmpa[0][j];
921                                 mod[i][j]=evec[0].getComponent(i)*tmpm[0][j];
922                                 for(int k=1;k<N;k++) {
923                                         arg[i][j]+=evec[k].getComponent(i)*tmpa[k][j];
924                                         mod[i][j]+=evec[k].getComponent(i)*tmpm[k][j];
925                                 }
926                         }
927                 }
928                 final AbstractDoubleSquareMatrix us[]=new AbstractDoubleSquareMatrix[2];
929                 us[0]=new DoubleSquareMatrix(arg);
930                 us[1]=new DoubleSquareMatrix(mod);
931                 return us;
932         }
933
934 // MAP ELEMENTS
935
936         /**
937         * Applies a function on all the matrix elements.
938         * @param f a user-defined function
939         * @return a double square matrix
940         */
941         public AbstractDoubleMatrix mapElements(final Mapping f) {
942                 final double array[][]=new double[numRows][numCols];
943                 for(int i=0;i<numRows;i++) {
944                         array[i][0] = f.map(getElement(i,0));
945                         for(int j=1;j<numCols;j++)
946                                 array[i][j] = f.map(getElement(i,j));
947                 }
948                 return new DoubleSquareMatrix(array);
949         }
950
951         private static class SquareMatrixAdaptor extends AbstractDoubleSquareMatrix {
952                 private final AbstractDoubleMatrix matrix;
953                 private SquareMatrixAdaptor(AbstractDoubleMatrix m) {
954                         super(m.rows());
955                         matrix = m;
956                 }
957                 public double getElement(int i, int j) {
958                         return matrix.getElement(i,j);
959                 }
960                 public void setElement(int i, int j, double x) {
961                         matrix.setElement(i,j, x);
962                 }
963         }
964 }
965
Popular Tags