KickJava   Java API By Example, From Geeks To Geeks.

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


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.AbstractIntegerVector;
11 import JSci.maths.groups.AbelianGroup;
12 import JSci.maths.algebras.*;
13 import JSci.maths.fields.*;
14
15 /**
16 * The AbstractIntegerSquareMatrix class provides an object for encapsulating integer square matrices.
17 * @version 2.2
18 * @author Mark Hale
19 */

20 public abstract class AbstractIntegerSquareMatrix extends AbstractIntegerMatrix implements SquareMatrix {
21         protected transient AbstractDoubleSquareMatrix[] LU;
22         protected transient int[] LUpivot;
23         /**
24         * Constructs a matrix.
25         */

26         protected AbstractIntegerSquareMatrix(final int size) {
27                 super(size, size);
28         }
29         /**
30         * Converts this matrix to a double matrix.
31         * @return a double square matrix
32         */

33         public AbstractDoubleMatrix toDoubleMatrix() {
34                 final double ans[][]=new double[numRows][numCols];
35                 for(int i=0;i<numRows;i++) {
36                         for(int j=0;j<numCols;j++)
37                                 ans[i][j]=getElement(i,j);
38                 }
39                 return new DoubleSquareMatrix(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 int 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 Math.round((float)det)*LUpivot[numRows];
77                 }
78         }
79         /**
80         * Returns the trace.
81         */

82         public int trace() {
83                 int result = getElement(0,0);
84                 for(int i=1;i<numRows;i++)
85                         result += getElement(i,i);
86                 return result;
87         }
88
89 //============
90
// OPERATIONS
91
//============
92

93         /**
94         * Returns the negative of this matrix.
95         */

96         public AbelianGroup.Member negate() {
97                 final int array[][]=new int[numRows][numCols];
98                 for(int i=0;i<numRows;i++) {
99                         array[i][0] = -getElement(i,0);
100                         for(int j=1;j<numCols;j++)
101                                 array[i][j] = -getElement(i,j);
102                 }
103                 return new IntegerSquareMatrix(array);
104         }
105
106 // ADDITION
107

108         /**
109         * Returns the addition of this matrix and another.
110         * @param m a int square matrix
111         * @exception MatrixDimensionException If the matrices are not square or different sizes.
112         */

113         public final AbstractIntegerMatrix add(final AbstractIntegerMatrix m) {
114                 if(m instanceof AbstractIntegerSquareMatrix)
115                         return add((AbstractIntegerSquareMatrix)m);
116                 else if(numRows == m.rows() && numCols == m.columns())
117                         return add(new SquareMatrixAdaptor(m));
118                 else
119                         throw new MatrixDimensionException("Matrices are different sizes.");
120         }
121         /**
122         * Returns the addition of this matrix and another.
123         * @param m a int square matrix
124         * @exception MatrixDimensionException If the matrices are different sizes.
125         */

126         public AbstractIntegerSquareMatrix add(final AbstractIntegerSquareMatrix m) {
127                 if(numRows==m.rows() && numCols==m.columns()) {
128                         final int array[][]=new int[numRows][numCols];
129                         for(int i=0;i<numRows;i++) {
130                                 array[i][0] = getElement(i,0)+m.getElement(i,0);
131                                 for(int j=1;j<numCols;j++)
132                                         array[i][j] = getElement(i,j)+m.getElement(i,j);
133                         }
134                         return new IntegerSquareMatrix(array);
135                 } else {
136                         throw new MatrixDimensionException("Matrices are different sizes.");
137                 }
138         }
139
140 // SUBTRACTION
141

142         /**
143         * Returns the subtraction of this matrix and another.
144         * @param m a int square matrix
145         * @exception MatrixDimensionException If the matrices are not square or different sizes.
146         */

147         public final AbstractIntegerMatrix subtract(final AbstractIntegerMatrix m) {
148                 if(m instanceof AbstractIntegerSquareMatrix)
149                         return subtract((AbstractIntegerSquareMatrix)m);
150                 else if(numRows == m.rows() && numCols == m.columns())
151                         return subtract(new SquareMatrixAdaptor(m));
152                 else
153                         throw new MatrixDimensionException("Matrices are different sizes.");
154         }
155         /**
156         * Returns the subtraction of this matrix by another.
157         * @param m a int square matrix
158         * @exception MatrixDimensionException If the matrices are different sizes.
159         */

160         public AbstractIntegerSquareMatrix subtract(final AbstractIntegerSquareMatrix m) {
161                 if(numRows==m.rows() && numCols==m.columns()) {
162                         final int array[][]=new int[numRows][numCols];
163                         for(int i=0;i<numRows;i++) {
164                                 array[i][0] = getElement(i,0)-m.getElement(i,0);
165                                 for(int j=1;j<numCols;j++)
166                                         array[i][j] = getElement(i,j)-m.getElement(i,j);
167                         }
168                         return new IntegerSquareMatrix(array);
169                 } else {
170                         throw new MatrixDimensionException("Matrices are different sizes.");
171                 }
172         }
173
174 // SCALAR MULTIPLICATION
175

176         /**
177         * Returns the multiplication of this matrix by a scalar.
178         * @param x a int.
179         * @return a int square matrix.
180         */

181         public AbstractIntegerMatrix scalarMultiply(final int x) {
182                 final int array[][]=new int[numRows][numCols];
183                 for(int i=0;i<numRows;i++) {
184                         array[i][0] = x*getElement(i,0);
185                         for(int j=1;j<numCols;j++)
186                                 array[i][j] = x*getElement(i,j);
187                 }
188                 return new IntegerSquareMatrix(array);
189         }
190
191 // SCALAR DIVISON
192

193
194 // SCALAR PRODUCT
195

196         /**
197         * Returns the scalar product of this matrix and another.
198         * @param m a int square matrix.
199         * @exception MatrixDimensionException If the matrices are not square or different sizes.
200         */

201         public final int scalarProduct(final AbstractIntegerMatrix m) {
202                 if(m instanceof AbstractIntegerSquareMatrix)
203                         return scalarProduct((AbstractIntegerSquareMatrix)m);
204                 else if(numRows == m.rows() && numCols == m.columns())
205                         return scalarProduct(new SquareMatrixAdaptor(m));
206                 else
207                         throw new MatrixDimensionException("Matrices are different sizes.");
208         }
209         /**
210         * Returns the scalar product of this matrix and another.
211         * @param m a int square matrix.
212         * @exception MatrixDimensionException If the matrices are different sizes.
213         */

214         public int scalarProduct(final AbstractIntegerSquareMatrix m) {
215                 if(numRows==m.rows() && numCols==m.columns()) {
216                         int ans = 0;
217                         for(int i=0; i<numRows; i++) {
218                                 ans += getElement(i,0)*m.getElement(i,0);
219                                 for(int j=1; j<numCols; j++)
220                                         ans += getElement(i,j)*m.getElement(i,j);
221                         }
222                         return ans;
223                 } else {
224                        throw new MatrixDimensionException("Matrices are different sizes.");
225                 }
226         }
227
228 // MATRIX MULTIPLICATION
229

230         /**
231         * Returns the multiplication of this matrix and another.
232         * @param m a int square matrix
233         * @exception MatrixDimensionException If the matrices are different sizes.
234         */

235         public AbstractIntegerSquareMatrix multiply(final AbstractIntegerSquareMatrix m) {
236                 if(numCols==m.rows()) {
237                         final int mColumns = m.columns();
238                         final int array[][]=new int[numRows][mColumns];
239                         for(int j=0; j<numRows; j++) {
240                                 for(int k=0; k<mColumns; k++) {
241                                         array[j][k] = getElement(j,0)*m.getElement(0,k);
242                                         for(int n=1; n<numCols; n++)
243                                                 array[j][k] += getElement(j,n)*m.getElement(n,k);
244                                 }
245                         }
246                         return new IntegerSquareMatrix(array);
247                 } else {
248                         throw new MatrixDimensionException("Incompatible matrices.");
249                 }
250         }
251
252 // DIRECT SUM
253

254         /**
255         * Returns the direct sum of this matrix and another.
256         */

257         public AbstractIntegerSquareMatrix directSum(final AbstractIntegerSquareMatrix m) {
258                 final int array[][]=new int[numRows+m.numRows][numCols+m.numCols];
259                 for(int i=0;i<numRows;i++) {
260                         for(int j=0;j<numCols;j++)
261                                 array[i][j] = getElement(i,j);
262                 }
263                 for(int i=0;i<m.numRows;i++) {
264                         for(int j=0;j<m.numCols;j++)
265                                 array[i+numRows][j+numCols] = m.getElement(i,j);
266                 }
267                 return new IntegerSquareMatrix(array);
268         }
269
270 // TENSOR PRODUCT
271

272         /**
273         * Returns the tensor product of this matrix and another.
274         */

275         public AbstractIntegerSquareMatrix tensor(final AbstractIntegerSquareMatrix m) {
276                 final int array[][]=new int[numRows*m.numRows][numCols*m.numCols];
277                 for(int i=0;i<numRows;i++) {
278                         for(int j=0;j<numCols;j++) {
279                                 for(int k=0;k<m.numRows;j++) {
280                                         for(int l=0;l<m.numCols;l++)
281                                                 array[i*m.numRows+k][j*m.numCols+l] = getElement(i,j)*m.getElement(k,l);
282                                 }
283                         }
284                 }
285                 return new IntegerSquareMatrix(array);
286         }
287
288 // TRANSPOSE
289

290         /**
291         * Returns the transpose of this matrix.
292         * @return a int square matrix
293         */

294         public Matrix transpose() {
295                 final int array[][]=new int[numCols][numRows];
296                 for(int i=0;i<numRows;i++) {
297                         array[0][i] = getElement(i,0);
298                         for(int j=1;j<numCols;j++)
299                                 array[j][i] = getElement(i,j);
300                 }
301                 return new IntegerSquareMatrix(array);
302         }
303
304 // INVERSE
305

306         /**
307         * Returns the inverse of this matrix.
308         * @return a double square matrix.
309         */

310         public AbstractDoubleSquareMatrix inverse() {
311                 final int N=numRows;
312                 final double arrayL[][]=new double[N][N];
313                 final double arrayU[][]=new double[N][N];
314                 final AbstractDoubleSquareMatrix lu[]=this.luDecompose(null);
315                 arrayL[0][0]=1.0/lu[0].getElement(0,0);
316                 arrayU[0][0]=1.0/lu[1].getElement(0,0);
317                 for(int i=1;i<N;i++) {
318                         arrayL[i][i]=1.0/lu[0].getElement(i,i);
319                         arrayU[i][i]=1.0/lu[1].getElement(i,i);
320                 }
321                 for(int i=0;i<N-1;i++) {
322                         for(int j=i+1;j<N;j++) {
323                                 double tmpL=0.0, tmpU=0.0;
324                                 for(int k=i;k<j;k++) {
325                                         tmpL-=lu[0].getElement(j,k)*arrayL[k][i];
326                                         tmpU-=arrayU[i][k]*lu[1].getElement(k,j);
327                                 }
328                                 arrayL[j][i]=tmpL/lu[0].getElement(j,j);
329                                 arrayU[i][j]=tmpU/lu[1].getElement(j,j);
330                         }
331                 }
332                 // matrix multiply arrayU x arrayL
333
final double inv[][]=new double[N][N];
334                 for(int i=0;i<N;i++) {
335                         for(int j=0;j<i;j++) {
336                                 for(int k=i;k<N;k++)
337                                         inv[i][LUpivot[j]]+=arrayU[i][k]*arrayL[k][j];
338                         }
339                         for(int j=i;j<N;j++) {
340                                 for(int k=j;k<N;k++)
341                                         inv[i][LUpivot[j]]+=arrayU[i][k]*arrayL[k][j];
342                         }
343                 }
344                 return new DoubleSquareMatrix(inv);
345         }
346
347 // LU DECOMPOSITION
348

349         /**
350         * Returns the LU decomposition of this matrix.
351         * @param pivot an empty array of length <code>rows()+1</code>
352         * to hold the pivot information (null if not interested).
353         * The last array element will contain the parity.
354         * @return an array with [0] containing the L-matrix
355         * and [1] containing the U-matrix.
356         * @jsci.planetmath LUDecomposition
357         */

358         public AbstractDoubleSquareMatrix[] luDecompose(int pivot[]) {
359                 if(LU!=null) {
360                         if(pivot!=null)
361                                 System.arraycopy(LUpivot,0,pivot,0,pivot.length);
362                         return LU;
363                 }
364                 int pivotrow;
365                 final int N=numRows;
366                 final double arrayL[][]=new double[N][N];
367                 final double arrayU[][]=new double[N][N];
368                 if(pivot==null)
369                         pivot=new int[N+1];
370                 for(int i=0;i<N;i++)
371                         pivot[i]=i;
372                 pivot[N]=1;
373         // LU decomposition to arrayU
374
for(int j=0;j<N;j++) {
375                         for(int i=0;i<j;i++) {
376                                 double tmp=getElement(pivot[i],j);
377                                 for(int k=0;k<i;k++)
378                                         tmp-=arrayU[i][k]*arrayU[k][j];
379                                 arrayU[i][j]=tmp;
380                         }
381                         double max=0.0;
382                         pivotrow=j;
383                         for(int i=j;i<N;i++) {
384                                 double tmp=getElement(pivot[i],j);
385                                 for(int k=0;k<j;k++)
386                                         tmp-=arrayU[i][k]*arrayU[k][j];
387                                 arrayU[i][j]=tmp;
388                         // while we're here search for a pivot for arrayU[j][j]
389
tmp=Math.abs(tmp);
390                                 if(tmp>max) {
391                                         max=tmp;
392                                         pivotrow=i;
393                                 }
394                         }
395                 // swap row j with pivotrow
396
if(pivotrow!=j) {
397                                 double[] tmprow = arrayU[j];
398                                 arrayU[j] = arrayU[pivotrow];
399                                 arrayU[pivotrow] = tmprow;
400                                 int k=pivot[j];
401                                 pivot[j]=pivot[pivotrow];
402                                 pivot[pivotrow]=k;
403                                 // update parity
404
pivot[N]=-pivot[N];
405                         }
406                 // divide by pivot
407
double tmp=arrayU[j][j];
408                         for(int i=j+1;i<N;i++)
409                                 arrayU[i][j]/=tmp;
410                 }
411                 // move lower triangular part to arrayL
412
for(int j=0;j<N;j++) {
413                         arrayL[j][j]=1.0;
414                         for(int i=j+1;i<N;i++) {
415                                 arrayL[i][j]=arrayU[i][j];
416                                 arrayU[i][j]=0.0;
417                         }
418                 }
419                 LU=new AbstractDoubleSquareMatrix[2];
420                 LU[0]=new DoubleSquareMatrix(arrayL);
421                 LU[1]=new DoubleSquareMatrix(arrayU);
422                 LUpivot=new int[pivot.length];
423                 System.arraycopy(pivot,0,LUpivot,0,pivot.length);
424                 return LU;
425         }
426         /**
427         * Returns the LU decomposition of this matrix.
428         * Warning: no pivoting.
429         * @return an array with [0] containing the L-matrix
430         * and [1] containing the U-matrix.
431         * @jsci.planetmath LUDecomposition
432         */

433         public AbstractDoubleSquareMatrix[] luDecompose() {
434                 final int N=numRows;
435                 final double arrayL[][]=new double[N][N];
436                 final double arrayU[][]=new double[N][N];
437         // LU decomposition to arrayU
438
for(int j=0;j<N;j++) {
439                         for(int i=0;i<j;i++) {
440                                 double tmp=getElement(i,j);
441                                 for(int k=0;k<i;k++)
442                                         tmp-=arrayU[i][k]*arrayU[k][j];
443                                 arrayU[i][j]=tmp;
444                         }
445                         for(int i=j;i<N;i++) {
446                                 double tmp=getElement(i,j);
447                                 for(int k=0;k<j;k++)
448                                         tmp-=arrayU[i][k]*arrayU[k][j];
449                                 arrayU[i][j]=tmp;
450                         }
451                 // divide
452
double tmp=arrayU[j][j];
453                         for(int i=j+1;i<N;i++)
454                                 arrayU[i][j]/=tmp;
455                 }
456                 // move lower triangular part to arrayL
457
for(int j=0;j<N;j++) {
458                         arrayL[j][j]=1.0;
459                         for(int i=j+1;i<N;i++) {
460                                 arrayL[i][j]=arrayU[i][j];
461                                 arrayU[i][j]=0.0;
462                         }
463                 }
464                 AbstractDoubleSquareMatrix[] lu=new AbstractDoubleSquareMatrix[2];
465                 lu[0]=new DoubleSquareMatrix(arrayL);
466                 lu[1]=new DoubleSquareMatrix(arrayU);
467                 return lu;
468         }
469
470 // CHOLESKY DECOMPOSITION
471

472         /**
473         * Returns the Cholesky decomposition of this matrix.
474         * Matrix must be symmetric and positive definite.
475         * @return an array with [0] containing the L-matrix and [1] containing the U-matrix.
476         * @jsci.planetmath CholeskyDecomposition
477         */

478         public AbstractDoubleSquareMatrix[] choleskyDecompose() {
479                 final int N=numRows;
480                 final double arrayL[][]=new double[N][N];
481                 final double arrayU[][]=new double[N][N];
482                 double tmp=Math.sqrt(getElement(0,0));
483                 arrayL[0][0]=arrayU[0][0]=tmp;
484                 for(int i=1;i<N;i++)
485                         arrayL[i][0]=arrayU[0][i]=getElement(i,0)/tmp;
486                 for(int j=1;j<N;j++) {
487                         tmp=getElement(j,j);
488                         for(int i=0;i<j;i++)
489                                 tmp-=arrayL[j][i]*arrayL[j][i];
490                         arrayL[j][j]=arrayU[j][j]=Math.sqrt(tmp);
491                         for(int i=j+1;i<N;i++) {
492                                 tmp=getElement(i,j);
493                                 for(int k=0;k<i;k++)
494                                         tmp-=arrayL[j][k]*arrayU[k][i];
495                                 arrayL[i][j]=arrayU[j][i]=tmp/arrayU[j][j];
496                         }
497                 }
498                 final AbstractDoubleSquareMatrix lu[]=new AbstractDoubleSquareMatrix[2];
499                 lu[0]=new DoubleSquareMatrix(arrayL);
500                 lu[1]=new DoubleSquareMatrix(arrayU);
501                 return lu;
502         }
503
504 // QR DECOMPOSITION
505

506         /**
507         * Returns the QR decomposition of this matrix.
508         * Based on the code from <a HREF="http://math.nist.gov/javanumerics/jama/">JAMA</a> (public domain).
509         * @return an array with [0] containing the Q-matrix and [1] containing the R-matrix.
510         * @jsci.planetmath QRDecomposition
511         */

512         public AbstractDoubleSquareMatrix[] qrDecompose() {
513                 final int N=numRows;
514                 final double array[][]=new double[N][N];
515                 final double arrayQ[][]=new double[N][N];
516                 final double arrayR[][]=new double[N][N];
517                 // copy matrix
518
for(int i=0;i<N;i++) {
519                         array[i][0]=getElement(i,0);
520                         for(int j=1;j<N;j++)
521                                 array[i][j]=getElement(i,j);
522                 }
523                 for(int k=0; k<N; k++) {
524                         // compute l2-norm of kth column
525
double norm = array[k][k];
526                         for(int i=k+1; i<N; i++)
527                                 norm = ExtraMath.hypot(norm, array[i][k]);
528                         if(norm != 0.0) {
529                                 // form kth Householder vector
530
if(array[k][k] < 0.0)
531                                         norm = -norm;
532                                 for(int i=k; i<N; i++)
533                                         array[i][k] /= norm;
534                                 array[k][k] += 1.0;
535                                 // apply transformation to remaining columns
536
for(int j=k+1; j<N; j++) {
537                                         double s=array[k][k]*array[k][j];
538                                         for(int i=k+1; i<N; i++)
539                                                 s += array[i][k]*array[i][j];
540                                         s /= array[k][k];
541                                         for(int i=k; i<N; i++)
542                                                 array[i][j] -= s*array[i][k];
543                                 }
544                         }
545                         arrayR[k][k] = -norm;
546                 }
547                 for(int k=N-1; k>=0; k--) {
548                         arrayQ[k][k] = 1.0;
549                         for(int j=k; j<N; j++) {
550                                 if(array[k][k] != 0.0) {
551                                         double s = array[k][k]*arrayQ[k][j];
552                                         for(int i=k+1; i<N; i++)
553                                                 s += array[i][k]*arrayQ[i][j];
554                                         s /= array[k][k];
555                                         for(int i=k; i<N; i++)
556                                                 arrayQ[i][j] -= s*array[i][k];
557                                 }
558                         }
559                 }
560                 for(int i=0; i<N; i++) {
561                         for(int j=i+1; j<N; j++)
562                                 arrayR[i][j] = array[i][j];
563                 }
564                 final AbstractDoubleSquareMatrix qr[]=new AbstractDoubleSquareMatrix[2];
565                 qr[0]=new DoubleSquareMatrix(arrayQ);
566                 qr[1]=new DoubleSquareMatrix(arrayR);
567                 return qr;
568         }
569
570 // SINGULAR VALUE DECOMPOSITION
571

572         /**
573         * Returns the singular value decomposition of this matrix.
574         * Based on the code from <a HREF="http://math.nist.gov/javanumerics/jama/">JAMA</a> (public domain).
575         * @return an array with [0] containing the U-matrix, [1] containing the S-matrix and [2] containing the V-matrix.
576         * @jsci.planetmath SingularValueDecomposition
577         */

578         public AbstractDoubleSquareMatrix[] singularValueDecompose() {
579                 final int N=numRows;
580                 final int Nm1=N-1;
581                 final double array[][]=new double[N][N];
582                 final double arrayU[][]=new double[N][N];
583                 final double arrayS[]=new double[N];
584                 final double arrayV[][]=new double[N][N];
585                 final double e[]=new double[N];
586                 final double work[]=new double[N];
587                 // copy matrix
588
for(int i=0;i<N;i++) {
589                         array[i][0]=getElement(i,0);
590                         for(int j=1;j<N;j++)
591                                 array[i][j]=getElement(i,j);
592                 }
593                 // reduce matrix to bidiagonal form
594
for(int k=0;k<Nm1;k++) {
595                         // compute the transformation for the kth column
596
// compute l2-norm of kth column
597
arrayS[k]=array[k][k];
598                         for(int i=k+1;i<N;i++)
599                                 arrayS[k]=ExtraMath.hypot(arrayS[k],array[i][k]);
600                         if(arrayS[k]!=0.0) {
601                                 if(array[k][k]<0.0)
602                                         arrayS[k]=-arrayS[k];
603                                 for(int i=k;i<N;i++)
604                                         array[i][k]/=arrayS[k];
605                                 array[k][k]+=1.0;
606                         }
607                         arrayS[k]=-arrayS[k];
608                         for(int j=k+1;j<N;j++) {
609                                 if(arrayS[k]!=0.0) {
610                                         // apply the transformation
611
double t=array[k][k]*array[k][j];
612                                         for(int i=k+1; i<N; i++)
613                                                 t+=array[i][k]*array[i][j];
614                                         t /= array[k][k];
615                                         for(int i=k; i<N; i++)
616                                                 array[i][j] -= t*array[i][k];
617                                 }
618                                 e[j]=array[k][j];
619                         }
620                         for(int i=k;i<N;i++)
621                                 arrayU[i][k]=array[i][k];
622                         if(k<N-2) {
623                                 // compute the kth row transformation
624
// compute l2-norm of kth column
625
e[k]=e[k+1];
626                                 for(int i=k+2;i<N;i++)
627                                         e[k]=ExtraMath.hypot(e[k],e[i]);
628                                 if(e[k]!=0.0) {
629                                         if(e[k+1]<0.0)
630                                                 e[k]=-e[k];
631                                         for(int i=k+1;i<N;i++)
632                                                 e[i]/=e[k];
633                                         e[k+1]+=1.0;
634                                 }
635                                 e[k]=-e[k];
636                                 if(e[k]!=0.0) {
637                                         // apply the transformation
638
for(int i=k+1;i<N;i++) {
639                                                 work[i]=0.0;
640                                                 for(int j=k+1;j<N;j++)
641                                                         work[i]+=e[j]*array[i][j];
642                                         }
643                                         for(int j=k+1;j<N;j++) {
644                                                 double t = e[j]/e[k+1];
645                                                 for(int i=k+1;i<N;i++)
646                                                         array[i][j] -= t*work[i];
647                                         }
648                                 }
649                                 for(int i=k+1;i<N;i++)
650                                         arrayV[i][k]=e[i];
651                         }
652                 }
653                 // setup the final bidiagonal matrix of order p
654
int p=N;
655                 arrayS[Nm1]=array[Nm1][Nm1];
656                 e[N-2]=array[N-2][Nm1];
657                 e[Nm1]=0.0;
658                 arrayU[Nm1][Nm1]=1.0;
659                 for(int k=N-2;k>=0;k--) {
660                         if(arrayS[k]!=0.0) {
661                                 for(int j=k+1;j<N;j++) {
662                                         double t=arrayU[k][k]*arrayU[k][j];
663                                         for(int i=k+1;i<N;i++)
664                                                 t+=arrayU[i][k]*arrayU[i][j];
665                                         t /= arrayU[k][k];
666                                         for(int i=k;i<N;i++)
667                                                 arrayU[i][j] -= t*arrayU[i][k];
668                                 }
669                                 for(int i=k;i<N;i++)
670                                         arrayU[i][k]=-arrayU[i][k];
671                                 arrayU[k][k]+=1.0;
672                                 for(int i=0;i<k-1;i++)
673                                         arrayU[i][k]=0.0;
674                         } else {
675                                 for(int i=0;i<N;i++)
676                                         arrayU[i][k]=0.0;
677                                 arrayU[k][k]=1.0;
678                         }
679                 }
680                 for(int k=Nm1;k>=0;k--) {
681                         if(k<N-2 && e[k]!=0.0) {
682                                 for(int j=k+1;j<N;j++) {
683                                         double t=arrayV[k+1][k]*arrayV[k+1][j];
684                                         for(int i=k+2;i<N;i++)
685                                                 t+=arrayV[i][k]*arrayV[i][j];
686                                         t /= arrayV[k+1][k];
687                                         for(int i=k+1;i<N;i++)
688                                                 arrayV[i][j] -= t*arrayV[i][k];
689                                 }
690                         }
691                         for(int i=0;i<N;i++)
692                                 arrayV[i][k]=0.0;
693                         arrayV[k][k]=1.0;
694                 }
695                 final double eps=Math.pow(2.0,-52.0);
696                 int iter=0;
697                 while(p>0) {
698                         int k,action;
699                         // action = 1 if arrayS[p] and e[k-1] are negligible and k<p
700
// action = 2 if arrayS[k] is negligible and k<p
701
// action = 3 if e[k-1] is negligible, k<p, and arrayS[k], ..., arrayS[p] are not negligible (QR step)
702
// action = 4 if e[p-1] is negligible (convergence)
703
for(k=p-2;k>=-1;k--) {
704                                 if(k==-1)
705                                         break;
706                                 if(Math.abs(e[k])<=eps*(Math.abs(arrayS[k])+Math.abs(arrayS[k+1]))) {
707                                         e[k]=0.0;
708                                         break;
709                                 }
710                         }
711                         if(k==p-2) {
712                                 action=4;
713                         } else {
714                                 int ks;
715                                 for(ks=p-1;ks>=k;ks--) {
716                                         if(ks==k)
717                                                 break;
718                                         double t=(ks!=p ? Math.abs(e[ks]) : 0.0)+(ks!=k+1 ? Math.abs(e[ks-1]) : 0.0);
719                                         if(Math.abs(arrayS[ks])<=eps*t) {
720                                                 arrayS[ks]=0.0;
721                                                 break;
722                                         }
723                                 }
724                                 if(ks==k) {
725                                         action=3;
726                                 } else if(ks==p-1) {
727                                         action=1;
728                                 } else {
729                                         action=2;
730                                         k=ks;
731                                 }
732                         }
733                         k++;
734                         switch(action) {
735                                 // deflate negligible arrayS[p]
736
case 1: {
737                                         double f=e[p-2];
738                                         e[p-2]=0.0;
739                                         for(int j=p-2;j>=k;j--) {
740                                                 double t=ExtraMath.hypot(arrayS[j],f);
741                                                 final double cs=arrayS[j]/t;
742                                                 final double sn=f/t;
743                                                 arrayS[j]=t;
744                                                 if(j!=k) {
745                                                         f=-sn*e[j-1];
746                                                         e[j-1]*=cs;
747                                                 }
748                                                 for(int i=0;i<N;i++) {
749                                                         t=cs*arrayV[i][j]+sn*arrayV[i][p-1];
750                                                         arrayV[i][p-1]=-sn*arrayV[i][j]+cs*arrayV[i][p-1];
751                                                         arrayV[i][j]=t;
752                                                 }
753                                         }
754                                         } break;
755                                 // split at negligible arrayS[k]
756
case 2: {
757                                         double f=e[k-1];
758                                         e[k-1]=0.0;
759                                         for(int j=k;j<p;j++) {
760                                                 double t=ExtraMath.hypot(arrayS[j],f);
761                                                 final double cs=arrayS[j]/t;
762                                                 final double sn=f/t;
763                                                 arrayS[j]=t;
764                                                 f=-sn*e[j];
765                                                 e[j]*=cs;
766                                                 for(int i=0;i<N;i++) {
767                                                         t=cs*arrayU[i][j]+sn*arrayU[i][k-1];
768                                                         arrayU[i][k-1]=-sn*arrayU[i][j]+cs*arrayU[i][k-1];
769                                                         arrayU[i][j]=t;
770                                                 }
771                                         }
772                                         } break;
773                                 // perform one QR step
774
case 3: {
775                                         // calculate the shift
776
final double scale=Math.max(Math.max(Math.max(Math.max(
777                                                 Math.abs(arrayS[p-1]),Math.abs(arrayS[p-2])),Math.abs(e[p-2])),
778                                                 Math.abs(arrayS[k])),Math.abs(e[k]));
779                                         double sp=arrayS[p-1]/scale;
780                                         double spm1=arrayS[p-2]/scale;
781                                         double epm1=e[p-2]/scale;
782                                         double sk=arrayS[k]/scale;
783                                         double ek=e[k]/scale;
784                                         double b=((spm1+sp)*(spm1-sp)+epm1*epm1)/2.0;
785                                         double c=(sp*epm1)*(sp*epm1);
786                                         double shift=0.0;
787                                         if(b!=0.0 || c!=0.0) {
788                                                 shift=Math.sqrt(b*b+c);
789                                                 if(b<0.0)
790                                                         shift=-shift;
791                                                 shift=c/(b+shift);
792                                         }
793                                         double f=(sk+sp)*(sk-sp)+shift;
794                                         double g=sk*ek;
795                                         // chase zeros
796
for(int j=k;j<p-1;j++) {
797                                                 double t=ExtraMath.hypot(f,g);
798                                                 double cs=f/t;
799                                                 double sn=g/t;
800                                                 if(j!=k)
801                                                         e[j-1]=t;
802                                                 f=cs*arrayS[j]+sn*e[j];
803                                                 e[j]=cs*e[j]-sn*arrayS[j];
804                                                 g=sn*arrayS[j+1];
805                                                 arrayS[j+1]*=cs;
806                                                 for(int i=0;i<N;i++) {
807                                                         t=cs*arrayV[i][j]+sn*arrayV[i][j+1];
808                                                         arrayV[i][j+1]=-sn*arrayV[i][j]+cs*arrayV[i][j+1];
809                                                         arrayV[i][j]=t;
810                                                 }
811                                                 t=ExtraMath.hypot(f,g);
812                                                 cs=f/t;
813                                                 sn=g/t;
814                                                 arrayS[j]=t;
815                                                 f=cs*e[j]+sn*arrayS[j+1];
816                                                 arrayS[j+1]=-sn*e[j]+cs*arrayS[j+1];
817                                                 g=sn*e[j+1];
818                                                 e[j+1]*=cs;
819                                                 if(j<Nm1) {
820                                                         for(int i=0;i<N;i++) {
821                                                                 t=cs*arrayU[i][j]+sn*arrayU[i][j+1];
822                                                                 arrayU[i][j+1]=-sn*arrayU[i][j]+cs*arrayU[i][j+1];
823                                                                 arrayU[i][j]=t;
824                                                         }
825                                                 }
826                                         }
827                                         e[p-2]=f;
828                                         iter++;
829                                         } break;
830                                 // convergence
831
case 4: {
832                                         // make the singular values positive
833
if(arrayS[k]<=0.0) {
834                                                 arrayS[k]=-arrayS[k];
835                                                 for(int i=0;i<p;i++)
836                                                         arrayV[i][k]=-arrayV[i][k];
837                                         }
838                                         // order the singular values
839
while(k<p-1) {
840                                                 if(arrayS[k]>=arrayS[k+1])
841                                                         break;
842                                                 double tmp=arrayS[k];
843                                                 arrayS[k]=arrayS[k+1];
844                                                 arrayS[k+1]=tmp;
845                                                 if(k<Nm1) {
846                                                         for(int i=0;i<N;i++) {
847                                                                 tmp=arrayU[i][k+1];
848                                                                 arrayU[i][k+1]=arrayU[i][k];
849                                                                 arrayU[i][k]=tmp;
850                                                                 tmp=arrayV[i][k+1];
851                                                                 arrayV[i][k+1]=arrayV[i][k];
852                                                                 arrayV[i][k]=tmp;
853                                                         }
854                                                 }
855                                                 k++;
856                                         }
857                                         iter=0;
858                                         p--;
859                                         } break;
860                         }
861                 }
862                 final AbstractDoubleSquareMatrix svd[]=new AbstractDoubleSquareMatrix[3];
863                 svd[0]=new DoubleSquareMatrix(arrayU);
864                 svd[1]=new DoubleDiagonalMatrix(arrayS);
865                 svd[2]=new DoubleSquareMatrix(arrayV);
866                 return svd;
867         }
868
869
870         private static class SquareMatrixAdaptor extends AbstractIntegerSquareMatrix {
871                 private final AbstractIntegerMatrix matrix;
872                 private SquareMatrixAdaptor(AbstractIntegerMatrix m) {
873                         super(m.rows());
874                         matrix = m;
875                 }
876                 public int getElement(int i, int j) {
877                         return matrix.getElement(i,j);
878                 }
879                 public void setElement(int i, int j, int x) {
880                         matrix.setElement(i,j, x);
881                 }
882         }
883 }
884
Popular Tags