IntegerSquareMatrix


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