DoubleSquareMatrix


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.AbstractDoubleVector;
11  import JSci.maths.vectors.DoubleVector;
12  import JSci.maths.groups.AbelianGroup;
13  import JSci.maths.algebras.*;
14  import JSci.maths.fields.*;
15  
16  /**
17  * The DoubleSquareMatrix class provides an object for encapsulating double square matrices.
18  * @version 2.2
19  * @author Mark Hale
20  */
21  public class DoubleSquareMatrix extends AbstractDoubleSquareMatrix {
22          /**
23          * Array containing the elements of the matrix.
24          */
25          protected final double matrix[][];
26          /**
27          * Constructs a matrix by wrapping an array.
28          * @param array an assigned value
29          */
30          public DoubleSquareMatrix(final double 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 DoubleSquareMatrix(final int size) {
40                  this(new double[size][size]);
41          }
42          /**
43          * Constructs a matrix from an array of vectors (columns).
44          * @param array an assigned value
45          */
46          public DoubleSquareMatrix(final AbstractDoubleVector 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 double matrix
56          */
57          public boolean equals(AbstractDoubleMatrix m, double tol) {
58                  if(m != null && numRows == m.rows() && numCols == m.columns()) {
59              double sumSqr = 0;
60                          for(int i=0;i<numRows;i++) {
61                                  for(int j=0;j<numCols;j++) {
62                      double 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 an integer matrix.
87          * @return an integer matrix
88          */
89          public AbstractIntegerMatrix toIntegerMatrix() {
90                  final int ans[][]=new int[numRows][numCols];
91                  for(int i=0;i<numRows;i++) {
92                          for(int j=0;j<numCols;j++)
93                                  ans[i][j]=Math.round((float)matrix[i][j]);
94                  }
95                  return new IntegerSquareMatrix(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 double 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, double 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 double infNorm() {
140                 double result=0;
141                 for(int i=0;i<numRows;i++) {
142                         double 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 double 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 det*LUpivot[numRows];
174                 }
175         }
176         /**
177         * Returns the trace.
178         */
179         public double trace() {
180                 double 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 double array[][]=new double[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 DoubleSquareMatrix(array);
201         }
202 
203 // ADDITION
204 
205         /**
206         * Returns the addition of this matrix and another.
207         * @param m a double matrix
208         * @exception MatrixDimensionException If the matrices are different sizes.
209         */
210         public AbstractDoubleSquareMatrix add(final AbstractDoubleSquareMatrix m) {
211                 if(m instanceof DoubleSquareMatrix)
212                         return add((DoubleSquareMatrix)m);
213 
214                 if(numRows==m.rows() && numCols==m.columns()) {
215                         final double array[][]=new double[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 DoubleSquareMatrix(array);
222                 } else {
223                         throw new MatrixDimensionException("Matrices are different sizes.");
224                 }
225         }
226         public DoubleSquareMatrix add(final DoubleSquareMatrix m) {
227                 if(numRows==m.numRows && numCols==m.numCols) {
228                         final double array[][]=new double[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 DoubleSquareMatrix(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 double matrix
244         * @exception MatrixDimensionException If the matrices are different sizes.
245         */
246         public AbstractDoubleSquareMatrix subtract(final AbstractDoubleSquareMatrix m) {
247                 if(m instanceof DoubleSquareMatrix)
248                         return subtract((DoubleSquareMatrix)m);
249 
250                 if(numRows==m.rows() && numCols==m.columns()) {
251                         final double array[][]=new double[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 DoubleSquareMatrix(array);
258                 } else {
259                         throw new MatrixDimensionException("Matrices are different sizes.");
260                 }
261         }
262         public DoubleSquareMatrix subtract(final DoubleSquareMatrix m) {
263                 if(numRows==m.numRows && numCols==m.numCols) {
264                         final double array[][]=new double[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 DoubleSquareMatrix(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 double.
280         * @return a double square matrix.
281         */
282         public AbstractDoubleMatrix scalarMultiply(final double x) {
283                 final double array[][]=new double[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 DoubleSquareMatrix(array);
290         }
291 
292 // SCALAR DIVISON
293 
294         /**
295         * Returns the division of this matrix by a scalar.
296         * @param x a double.
297         * @return a double square matrix.
298         */
299         public AbstractDoubleMatrix scalarDivide(final double x) {
300                 final double array[][]=new double[numRows][numCols];
301                 for(int i=0;i<numRows;i++) {
302                         array[i][0] = matrix[i][0]/x;
303                         for(int j=1;j<numCols;j++)
304                                 array[i][j] = matrix[i][j]/x;
305                 }
306                 return new DoubleSquareMatrix(array);
307         }
308 
309 // SCALAR PRODUCT
310 
311         /**
312         * Returns the scalar product of this matrix and another.
313         * @param m a double matrix.
314         * @exception MatrixDimensionException If the matrices are different sizes.
315         */
316         public double scalarProduct(final AbstractDoubleSquareMatrix m) {
317                 if(m instanceof DoubleSquareMatrix)
318                         return scalarProduct((DoubleSquareMatrix)m);
319 
320                 if(numRows==m.rows() && numCols==m.columns()) {
321                         double ans=0;
322                         for(int i=0;i<numRows;i++) {
323                                 ans += matrix[i][0]*m.getElement(i,0);
324                                 for(int j=1;j<numCols;j++)
325                                         ans += matrix[i][j]*m.getElement(i,j);
326                         }
327                         return ans;
328                 } else {
329                        throw new MatrixDimensionException("Matrices are different sizes.");
330                 }
331         }
332         public double scalarProduct(final DoubleSquareMatrix m) {
333                 if(numRows==m.numRows && numCols==m.numCols) {
334                         double ans=0;
335                         for(int i=0;i<numRows;i++) {
336                                 ans+=matrix[i][0]*m.matrix[i][0];
337                                 for(int j=1;j<numCols;j++)
338                                         ans+=matrix[i][j]*m.matrix[i][j];
339                         }
340                         return ans;
341                 } else
342                         throw new MatrixDimensionException("Matrices are different sizes.");
343         }
344 
345 // MATRIX MULTIPLICATION
346 
347         /**
348         * Returns the multiplication of a vector by this matrix.
349         * @param v a double vector.
350         * @exception DimensionException If the matrix and vector are incompatible.
351         */
352         public AbstractDoubleVector multiply(final AbstractDoubleVector v) {
353                 if(numCols==v.dimension()) {
354                         final double array[]=new double[numRows];
355                         for(int i=0;i<numRows;i++) {
356                                 array[i]=matrix[i][0]*v.getComponent(0);
357                                 for(int j=1;j<numCols;j++)
358                                         array[i]+=matrix[i][j]*v.getComponent(j);
359                         }
360                         return new DoubleVector(array);
361                 } else {
362                         throw new DimensionException("Matrix and vector are incompatible.");
363                 }
364         }
365         /**
366         * Returns the multiplication of this matrix and another.
367         * @param m a double matrix
368         * @return a AbstractDoubleMatrix or a AbstractDoubleSquareMatrix as appropriate
369         * @exception MatrixDimensionException If the matrices are incompatible.
370         */
371         public AbstractDoubleSquareMatrix multiply(final AbstractDoubleSquareMatrix m) {
372                 if(m instanceof DoubleSquareMatrix)
373                         return multiply((DoubleSquareMatrix)m);
374 
375                 if(numCols==m.rows()) {
376                         final int mColumns = m.columns();
377                         final double array[][]=new double[numRows][mColumns];
378                         for(int j=0; j<numRows; j++) {
379                                 for(int k=0; k<mColumns; k++) {
380                                         array[j][k] = matrix[j][0]*m.getElement(0,k);
381                                         for(int n=1; n<numCols; n++)
382                                                 array[j][k] += matrix[j][n]*m.getElement(n,k);
383                                 }
384                         }
385                         return new DoubleSquareMatrix(array);
386                 } else {
387                         throw new MatrixDimensionException("Incompatible matrices.");
388                 }
389         }
390         public DoubleSquareMatrix multiply(final DoubleSquareMatrix m) {
391                 if(numCols==m.numRows) {
392                         final double array[][]=new double[numRows][m.numCols];
393                         for(int j=0;j<numRows;j++) {
394                                 for(int k=0;k<m.numCols;k++) {
395                                         array[j][k]=matrix[j][0]*m.matrix[0][k];
396                                         for(int n=1;n<numCols;n++)
397                                                 array[j][k]+=matrix[j][n]*m.matrix[n][k];
398                                 }
399                         }
400                         return new DoubleSquareMatrix(array);
401                 } else
402                         throw new MatrixDimensionException("Incompatible matrices.");
403         }
404 
405 // DIRECT SUM
406 
407         /**
408         * Returns the direct sum of this matrix and another.
409         */
410         public AbstractDoubleSquareMatrix directSum(final AbstractDoubleSquareMatrix m) {
411                 final double array[][]=new double[numRows+m.numRows][numCols+m.numCols];
412                 for(int i=0;i<numRows;i++) {
413                         for(int j=0;j<numCols;j++)
414                                 array[i][j] = matrix[i][j];
415                 }
416                 for(int i=0;i<m.numRows;i++) {
417                         for(int j=0;j<m.numCols;j++)
418                                 array[i+numRows][j+numCols] = m.getElement(i,j);
419                 }
420                 return new DoubleSquareMatrix(array);
421         }
422 
423 // TENSOR PRODUCT
424 
425         /**
426         * Returns the tensor product of this matrix and another.
427         */
428         public AbstractDoubleSquareMatrix tensor(final AbstractDoubleSquareMatrix m) {
429                 final double array[][]=new double[numRows*m.numRows][numCols*m.numCols];
430                 for(int i=0;i<numRows;i++) {
431                         for(int j=0;j<numCols;j++) {
432                                 for(int k=0;k<m.numRows;j++) {
433                                         for(int l=0;l<m.numCols;l++)
434                                                 array[i*m.numRows+k][j*m.numCols+l] = matrix[i][j]*m.getElement(k,l);
435                                 }
436                         }
437                 }
438                 return new DoubleSquareMatrix(array);
439         }
440 
441 // TRANSPOSE
442 
443         /**
444         * Returns the transpose of this matrix.
445         * @return a double matrix
446         */
447         public Matrix transpose() {
448                 final double array[][]=new double[numCols][numRows];
449                 for(int i=0;i<numRows;i++) {
450                         array[0][i] = matrix[i][0];
451                         for(int j=1;j<numCols;j++)
452                                 array[j][i] = matrix[i][j];
453                 }
454                 return new DoubleSquareMatrix(array);
455         }
456 
457 // INVERSE
458 
459         /**
460         * Returns the inverse of this matrix.
461         * @return a double square matrix.
462         */
463         public AbstractDoubleSquareMatrix inverse() {
464                 final int N=numRows;
465                 final double arrayL[][]=new double[N][N];
466                 final double arrayU[][]=new double[N][N];
467                 final DoubleSquareMatrix lu[] = (DoubleSquareMatrix[]) this.luDecompose(null);
468                 arrayL[0][0]=1.0/lu[0].matrix[0][0];
469                 arrayU[0][0]=1.0/lu[1].matrix[0][0];
470                 for(int i=1;i<N;i++) {
471                         arrayL[i][i]=1.0/lu[0].matrix[i][i];
472                         arrayU[i][i]=1.0/lu[1].matrix[i][i];
473                 }
474                 for(int i=0;i<N-1;i++) {
475                         for(int j=i+1;j<N;j++) {
476                                 double tmpL=0.0, tmpU=0.0;
477                                 for(int k=i;k<j;k++) {
478                                         tmpL-=lu[0].matrix[j][k]*arrayL[k][i];
479                                         tmpU-=arrayU[i][k]*lu[1].matrix[k][j];
480                                 }
481                                 arrayL[j][i]=tmpL/lu[0].matrix[j][j];
482                                 arrayU[i][j]=tmpU/lu[1].matrix[j][j];
483                         }
484                 }
485                 // matrix multiply arrayU x arrayL
486                 final double inv[][]=new double[N][N];
487                 for(int i=0;i<N;i++) {
488                         for(int j=0;j<i;j++) {
489                                 for(int k=i;k<N;k++)
490                                         inv[i][LUpivot[j]]+=arrayU[i][k]*arrayL[k][j];
491                         }
492                         for(int j=i;j<N;j++) {
493                                 for(int k=j;k<N;k++)
494                                         inv[i][LUpivot[j]]+=arrayU[i][k]*arrayL[k][j];
495                         }
496                 }
497                 return new DoubleSquareMatrix(inv);
498         }
499 
500 // LU DECOMPOSITION
501 
502         /**
503         * Returns the LU decomposition of this matrix.
504         * @param pivot an empty array of length <code>rows()+1</code>
505         * to hold the pivot information (null if not interested).
506         * The last array element will contain the parity.
507         * @return an array with [0] containing the L-matrix
508         * and [1] containing the U-matrix.
509         * @jsci.planetmath LUDecomposition
510         */
511         public final AbstractDoubleSquareMatrix[] luDecompose(int pivot[]) {
512                 if(LU!=null) {
513                         if(pivot!=null)
514                                 System.arraycopy(LUpivot,0,pivot,0,pivot.length);
515                         return LU;
516                 }
517                 int pivotrow;
518                 final int N=numRows;
519                 final double arrayL[][]=new double[N][N];
520                 final double arrayU[][]=new double[N][N];
521                 double tmp;
522                 if(pivot==null)
523                         pivot=new int[N+1];
524                 for(int i=0;i<N;i++)
525                         pivot[i]=i;
526                 pivot[N]=1;
527         // LU decomposition to arrayU
528                 for(int j=0;j<N;j++) {
529                         for(int i=0;i<j;i++) {
530                                 tmp=matrix[pivot[i]][j];
531                                 for(int k=0;k<i;k++)
532                                         tmp-=arrayU[i][k]*arrayU[k][j];
533                                 arrayU[i][j]=tmp;
534                         }
535                         double max=0.0;
536                         pivotrow=j;
537                         for(int i=j;i<N;i++) {
538                                 tmp=matrix[pivot[i]][j];
539                                 for(int k=0;k<j;k++)
540                                         tmp-=arrayU[i][k]*arrayU[k][j];
541                                 arrayU[i][j]=tmp;
542                         // while we're here search for a pivot for arrayU[j][j]
543                                 tmp=Math.abs(tmp);
544                                 if(tmp>max) {
545                                         max=tmp;
546                                         pivotrow=i;
547                                 }
548                         }
549                 // swap row j with pivotrow
550                         if(pivotrow!=j) {
551                                 double[] tmprow = arrayU[j];
552                                 arrayU[j] = arrayU[pivotrow];
553                                 arrayU[pivotrow] = tmprow;
554                                 int k=pivot[j];
555                                 pivot[j]=pivot[pivotrow];
556                                 pivot[pivotrow]=k;
557                                 // update parity
558                                 pivot[N]=-pivot[N];
559                         }
560                 // divide by pivot
561                         tmp=arrayU[j][j];
562                         for(int i=j+1;i<N;i++)
563                                 arrayU[i][j]/=tmp;
564                 }
565                 // move lower triangular part to arrayL
566                 for(int j=0;j<N;j++) {
567                         arrayL[j][j]=1.0;
568                         for(int i=j+1;i<N;i++) {
569                                 arrayL[i][j]=arrayU[i][j];
570                                 arrayU[i][j]=0.0;
571                         }
572                 }
573                 LU=new DoubleSquareMatrix[2];
574                 LU[0]=new DoubleSquareMatrix(arrayL);
575                 LU[1]=new DoubleSquareMatrix(arrayU);
576                 LUpivot=new int[pivot.length];
577                 System.arraycopy(pivot,0,LUpivot,0,pivot.length);
578                 return LU;
579         }
580         /**
581         * Returns the LU decomposition of this matrix.
582         * Warning: no pivoting.
583         * @return an array with [0] containing the L-matrix
584         * and [1] containing the U-matrix.
585         * @jsci.planetmath LUDecomposition
586         */
587         public final AbstractDoubleSquareMatrix[] luDecompose() {
588                 final int N=numRows;
589                 final double arrayL[][]=new double[N][N];
590                 final double arrayU[][]=new double[N][N];
591                 double tmp;
592         // LU decomposition to arrayU
593                 for(int j=0;j<N;j++) {
594                         for(int i=0;i<j;i++) {
595                                 tmp=matrix[i][j];
596                                 for(int k=0;k<i;k++)
597                                         tmp-=arrayU[i][k]*arrayU[k][j];
598                                 arrayU[i][j]=tmp;
599                         }
600                         for(int i=j;i<N;i++) {
601                                 tmp=matrix[i][j];
602                                 for(int k=0;k<j;k++)
603                                         tmp-=arrayU[i][k]*arrayU[k][j];
604                                 arrayU[i][j]=tmp;
605                         }
606                 // divide
607                         tmp=arrayU[j][j];
608                         for(int i=j+1;i<N;i++)
609                                 arrayU[i][j]/=tmp;
610                 }
611                 // move lower triangular part to arrayL
612                 for(int j=0;j<N;j++) {
613                         arrayL[j][j]=1.0;
614                         for(int i=j+1;i<N;i++) {
615                                 arrayL[i][j]=arrayU[i][j];
616                                 arrayU[i][j]=0.0;
617                         }
618                 }
619                 DoubleSquareMatrix[] lu=new DoubleSquareMatrix[2];
620                 lu[0]=new DoubleSquareMatrix(arrayL);
621                 lu[1]=new DoubleSquareMatrix(arrayU);
622                 return lu;
623         }
624 
625 // CHOLESKY DECOMPOSITION
626 
627         /**
628         * Returns the Cholesky decomposition of this matrix.
629         * Matrix must be symmetric and positive definite.
630         * @return an array with [0] containing the L-matrix and [1] containing the U-matrix.
631         * @jsci.planetmath CholeskyDecomposition
632         */
633         public AbstractDoubleSquareMatrix[] choleskyDecompose() {
634                 final int N=numRows;
635                 final double arrayL[][]=new double[N][N];
636                 final double arrayU[][]=new double[N][N];
637                 double tmp=Math.sqrt(matrix[0][0]);
638                 arrayL[0][0]=arrayU[0][0]=tmp;
639                 for(int i=1;i<N;i++)
640                         arrayL[i][0]=arrayU[0][i]=matrix[i][0]/tmp;
641                 for(int j=1;j<N;j++) {
642                         tmp=matrix[j][j];
643                         for(int i=0;i<j;i++)
644                                 tmp-=arrayL[j][i]*arrayL[j][i];
645                         arrayL[j][j]=arrayU[j][j]=Math.sqrt(tmp);
646                         for(int i=j+1;i<N;i++) {
647                                 tmp=matrix[i][j];
648                                 for(int k=0;k<i;k++)
649                                         tmp-=arrayL[j][k]*arrayU[k][i];
650                                 arrayL[i][j]=arrayU[j][i]=tmp/arrayU[j][j];
651                         }
652                 }
653                 final AbstractDoubleSquareMatrix lu[]=new AbstractDoubleSquareMatrix[2];
654                 lu[0]=new DoubleSquareMatrix(arrayL);
655                 lu[1]=new DoubleSquareMatrix(arrayU);
656                 return lu;
657         }
658 
659 // QR DECOMPOSITION
660 
661         /**
662         * Returns the QR decomposition of this matrix.
663         * Based on the code from <a HREF="http://math.nist.gov/javanumerics/jama/">JAMA</a> (public domain).
664         * @return an array with [0] containing the Q-matrix and [1] containing the R-matrix.
665         * @jsci.planetmath QRDecomposition
666         */
667         public AbstractDoubleSquareMatrix[] qrDecompose() {
668                 final int N=numRows;
669                 final double array[][]=new double[N][N];
670                 final double arrayQ[][]=new double[N][N];
671                 final double arrayR[][]=new double[N][N];
672                 // copy matrix
673                 for(int i=0;i<N;i++) {
674                         array[i][0]=matrix[i][0];
675                         for(int j=1;j<N;j++)
676                                 array[i][j]=matrix[i][j];
677                 }
678                 for(int k=0; k<N; k++) {
679                         // compute l2-norm of kth column
680                         double norm = array[k][k];
681                         for(int i=k+1; i<N; i++)
682                                 norm = ExtraMath.hypot(norm, array[i][k]);
683                         if(norm != 0.0) {
684                                 // form kth Householder vector
685                                 if(array[k][k] < 0.0)
686                                         norm = -norm;
687                                 for(int i=k; i<N; i++)
688                                         array[i][k] /= norm;
689                                 array[k][k] += 1.0;
690                                 // apply transformation to remaining columns
691                                 for(int j=k+1; j<N; j++) {
692                                         double s=array[k][k]*array[k][j];
693                                         for(int i=k+1; i<N; i++)
694                                                 s += array[i][k]*array[i][j];
695                                         s /= array[k][k];
696                                         for(int i=k; i<N; i++)
697                                                 array[i][j] -= s*array[i][k];
698                                 }
699                         }
700                         arrayR[k][k] = -norm;
701                 }
702                 for(int k=N-1; k>=0; k--) {
703                         arrayQ[k][k] = 1.0;
704                         for(int j=k; j<N; j++) {
705                                 if(array[k][k] != 0.0) {
706                                         double s = array[k][k]*arrayQ[k][j];
707                                         for(int i=k+1; i<N; i++)
708                                                 s += array[i][k]*arrayQ[i][j];
709                                         s /= array[k][k];
710                                         for(int i=k; i<N; i++)
711                                                 arrayQ[i][j] -= s*array[i][k];
712                                 }
713                         }
714                 }
715                 for(int i=0; i<N; i++) {
716                         for(int j=i+1; j<N; j++)
717                                 arrayR[i][j] = array[i][j];
718                 }
719                 final AbstractDoubleSquareMatrix qr[]=new AbstractDoubleSquareMatrix[2];
720                 qr[0]=new DoubleSquareMatrix(arrayQ);
721                 qr[1]=new DoubleSquareMatrix(arrayR);
722                 return qr;
723         }
724 
725 // SINGULAR VALUE DECOMPOSITION
726 
727         /**
728         * Returns the singular value decomposition of this matrix.
729         * Based on the code from <a HREF="http://math.nist.gov/javanumerics/jama/">JAMA</a> (public domain).
730         * @return an array with [0] containing the U-matrix, [1] containing the S-matrix and [2] containing the V-matrix.
731         * @jsci.planetmath SingularValueDecomposition
732         */
733         public AbstractDoubleSquareMatrix[] singularValueDecompose() {
734                 final int N=numRows;
735                 final int Nm1=N-1;
736                 final double array[][]=new double[N][N];
737                 final double arrayU[][]=new double[N][N];
738                 final double arrayS[]=new double[N];
739                 final double arrayV[][]=new double[N][N];
740                 final double e[]=new double[N];
741                 final double work[]=new double[N];
742                 // copy matrix
743                 for(int i=0;i<N;i++) {
744                         array[i][0]=matrix[i][0];
745                         for(int j=1;j<N;j++)
746                                 array[i][j]=matrix[i][j];
747                 }
748                 // reduce matrix to bidiagonal form
749                 for(int k=0;k<Nm1;k++) {
750                         // compute the transformation for the kth column
751                         // compute l2-norm of kth column
752                         arrayS[k]=array[k][k];
753                         for(int i=k+1;i<N;i++)
754                                 arrayS[k]=ExtraMath.hypot(arrayS[k],array[i][k]);
755                         if(arrayS[k]!=0.0) {
756                                 if(array[k][k]<0.0)
757                                         arrayS[k]=-arrayS[k];
758                                 for(int i=k;i<N;i++)
759                                         array[i][k]/=arrayS[k];
760                                 array[k][k]+=1.0;
761                         }
762                         arrayS[k]=-arrayS[k];
763                         for(int j=k+1;j<N;j++) {
764                                 if(arrayS[k]!=0.0) {
765                                         // apply the transformation
766                                         double t=array[k][k]*array[k][j];
767                                         for(int i=k+1; i<N; i++)
768                                                 t+=array[i][k]*array[i][j];
769                                         t /= array[k][k];
770                                         for(int i=k; i<N; i++)
771                                                 array[i][j] -= t*array[i][k];
772                                 }
773                                 e[j]=array[k][j];
774                         }
775                         for(int i=k;i<N;i++)
776                                 arrayU[i][k]=array[i][k];
777                         if(k<N-2) {
778                                 // compute the kth row transformation
779                                 // compute l2-norm of kth column
780                                 e[k]=e[k+1];
781                                 for(int i=k+2;i<N;i++)
782                                         e[k]=ExtraMath.hypot(e[k],e[i]);
783                                 if(e[k]!=0.0) {
784                                         if(e[k+1]<0.0)
785                                                 e[k]=-e[k];
786                                         for(int i=k+1;i<N;i++)
787                                                 e[i]/=e[k];
788                                         e[k+1]+=1.0;
789                                 }
790                                 e[k]=-e[k];
791                                 if(e[k]!=0.0) {
792                                         // apply the transformation
793                                         for(int i=k+1;i<N;i++) {
794                                                 work[i]=0.0;
795                                                 for(int j=k+1;j<N;j++)
796                                                         work[i]+=e[j]*array[i][j];
797                                         }
798                                         for(int j=k+1;j<N;j++) {
799                                                 double t = e[j]/e[k+1];
800                                                 for(int i=k+1;i<N;i++)
801                                                         array[i][j] -= t*work[i];
802                                         }
803                                 }
804                                 for(int i=k+1;i<N;i++)
805                                         arrayV[i][k]=e[i];
806                         }
807                 }
808                 // setup the final bidiagonal matrix of order p
809                 int p=N;
810                 arrayS[Nm1]=array[Nm1][Nm1];
811                 e[N-2]=array[N-2][Nm1];
812                 e[Nm1]=0.0;
813                 arrayU[Nm1][Nm1]=1.0;
814                 for(int k=N-2;k>=0;k--) {
815                         if(arrayS[k]!=0.0) {
816                                 for(int j=k+1;j<N;j++) {
817                                         double t=arrayU[k][k]*arrayU[k][j];
818                                         for(int i=k+1;i<N;i++)
819                                                 t+=arrayU[i][k]*arrayU[i][j];
820                                         t /= arrayU[k][k];
821                                         for(int i=k;i<N;i++)
822                                                 arrayU[i][j] -= t*arrayU[i][k];
823                                 }
824                                 for(int i=k;i<N;i++)
825                                         arrayU[i][k]=-arrayU[i][k];
826                                 arrayU[k][k]+=1.0;
827                                 for(int i=0;i<k-1;i++)
828                                         arrayU[i][k]=0.0;
829                         } else {
830                                 for(int i=0;i<N;i++)
831                                         arrayU[i][k]=0.0;
832                                 arrayU[k][k]=1.0;
833                         }
834                 }
835                 for(int k=Nm1;k>=0;k--) {
836                         if(k<N-2 && e[k]!=0.0) {
837                                 for(int j=k+1;j<N;j++) {
838                                         double t=arrayV[k+1][k]*arrayV[k+1][j];
839                                         for(int i=k+2;i<N;i++)
840                                                 t+=arrayV[i][k]*arrayV[i][j];
841                                         t /= arrayV[k+1][k];
842                                         for(int i=k+1;i<N;i++)
843                                                 arrayV[i][j] -= t*arrayV[i][k];
844                                 }
845                         }
846                         for(int i=0;i<N;i++)
847                                 arrayV[i][k]=0.0;
848                         arrayV[k][k]=1.0;
849                 }
850                 final double eps=Math.pow(2.0,-52.0);
851                 int iter=0;
852                 while(p>0) {
853                         int k, action;
854                         // action = 1 if arrayS[p] and e[k-1] are negligible and k<p
855                         // action = 2 if arrayS[k] is negligible and k<p
856                         // action = 3 if e[k-1] is negligible, k<p, and arrayS[k], ..., arrayS[p] are not negligible (QR step)
857                         // action = 4 if e[p-1] is negligible (convergence)
858                         for(k=p-2;k>=-1;k--) {
859                                 if(k==-1)
860                                         break;
861                                 if(Math.abs(e[k])<=eps*(Math.abs(arrayS[k])+Math.abs(arrayS[k+1]))) {
862                                         e[k]=0.0;
863                                         break;
864                                 }
865                         }
866                         if(k==p-2) {
867                                 action=4;
868                         } else {
869                                 int ks;
870                                 for(ks=p-1;ks>=k;ks--) {
871                                         if(ks==k)
872                                                 break;
873                                         double t=(ks!=p ? Math.abs(e[ks]) : 0.0)+(ks!=k+1 ? Math.abs(e[ks-1]) : 0.0);
874                                         if(Math.abs(arrayS[ks])<=eps*t) {
875                                                 arrayS[ks]=0.0;
876                                                 break;
877                                         }
878                                 }
879                                 if(ks==k) {
880                                         action=3;
881                                 } else if(ks==p-1) {
882                                         action=1;
883                                 } else {
884                                         action=2;
885                                         k=ks;
886                                 }
887                         }
888                         k++;
889                         switch(action) {
890                                 // deflate negligible arrayS[p]
891                                 case 1: {
892                                         double f=e[p-2];
893                                         e[p-2]=0.0;
894                                         for(int j=p-2;j>=k;j--) {
895                                                 double t=ExtraMath.hypot(arrayS[j],f);
896                                                 final double cs=arrayS[j]/t;
897                                                 final double sn=f/t;
898                                                 arrayS[j]=t;
899                                                 if(j!=k) {
900                                                         f=-sn*e[j-1];
901                                                         e[j-1]*=cs;
902                                                 }
903                                                 for(int i=0;i<N;i++) {
904                                                         t=cs*arrayV[i][j]+sn*arrayV[i][p-1];
905                                                         arrayV[i][p-1]=-sn*arrayV[i][j]+cs*arrayV[i][p-1];
906                                                         arrayV[i][j]=t;
907                                                 }
908                                         }
909                                         } break;
910                                 // split at negligible arrayS[k]
911                                 case 2: {
912                                         double f=e[k-1];
913                                         e[k-1]=0.0;
914                                         for(int j=k;j<p;j++) {
915                                                 double t=ExtraMath.hypot(arrayS[j],f);
916                                                 final double cs=arrayS[j]/t;
917                                                 final double sn=f/t;
918                                                 arrayS[j]=t;
919                                                 f=-sn*e[j];
920                                                 e[j]*=cs;
921                                                 for(int i=0;i<N;i++) {
922                                                         t=cs*arrayU[i][j]+sn*arrayU[i][k-1];
923                                                         arrayU[i][k-1]=-sn*arrayU[i][j]+cs*arrayU[i][k-1];
924                                                         arrayU[i][j]=t;
925                                                 }
926                                         }
927                                         } break;
928                                 // perform one QR step
929                                 case 3: {
930                                         // calculate the shift
931                                         final double scale=Math.max(Math.max(Math.max(Math.max(
932                                                 Math.abs(arrayS[p-1]),Math.abs(arrayS[p-2])),Math.abs(e[p-2])),
933                                                 Math.abs(arrayS[k])),Math.abs(e[k]));
934                                         double sp=arrayS[p-1]/scale;
935                                         double spm1=arrayS[p-2]/scale;
936                                         double epm1=e[p-2]/scale;
937                                         double sk=arrayS[k]/scale;
938                                         double ek=e[k]/scale;
939                                         double b=((spm1+sp)*(spm1-sp)+epm1*epm1)/2.0;
940                                         double c=(sp*epm1)*(sp*epm1);
941                                         double shift=0.0;
942                                         if(b!=0.0 || c!=0.0) {
943                                                 shift=Math.sqrt(b*b+c);
944                                                 if(b<0.0)
945                                                         shift=-shift;
946                                                 shift=c/(b+shift);
947                                         }
948                                         double f=(sk+sp)*(sk-sp)+shift;
949                                         double g=sk*ek;
950                                         // chase zeros
951                                         for(int j=k;j<p-1;j++) {
952                                                 double t=ExtraMath.hypot(f,g);
953                                                 double cs=f/t;
954                                                 double sn=g/t;
955                                                 if(j!=k)
956                                                         e[j-1]=t;
957                                                 f=cs*arrayS[j]+sn*e[j];
958                                                 e[j]=cs*e[j]-sn*arrayS[j];
959                                                 g=sn*arrayS[j+1];
960                                                 arrayS[j+1]*=cs;
961                                                 for(int i=0;i<N;i++) {
962                                                         t=cs*arrayV[i][j]+sn*arrayV[i][j+1];
963                                                         arrayV[i][j+1]=-sn*arrayV[i][j]+cs*arrayV[i][j+1];
964                                                         arrayV[i][j]=t;
965                                                 }
966                                                 t=ExtraMath.hypot(f,g);
967                                                 cs=f/t;
968                                                 sn=g/t;
969                                                 arrayS[j]=t;
970                                                 f=cs*e[j]+sn*arrayS[j+1];
971                                                 arrayS[j+1]=-sn*e[j]+cs*arrayS[j+1];
972                                                 g=sn*e[j+1];
973                                                 e[j+1]*=cs;
974                                                 if(j<Nm1) {
975                                                         for(int i=0;i<N;i++) {
976                                                                 t=cs*arrayU[i][j]+sn*arrayU[i][j+1];
977                                                                 arrayU[i][j+1]=-sn*arrayU[i][j]+cs*arrayU[i][j+1];
978                                                                 arrayU[i][j]=t;
979                                                         }
980                                                 }
981                                         }
982                                         e[p-2]=f;
983                                         iter++;
984                                         } break;
985                                 // convergence
986                                 case 4: {
987                                         // make the singular values positive
988                                         if(arrayS[k]<=0.0) {
989                                                 arrayS[k]=-arrayS[k];
990                                                 for(int i=0;i<p;i++)
991                                                         arrayV[i][k]=-arrayV[i][k];
992                                         }
993                                         // order the singular values
994                                         while(k<p-1) {
995                                                 if(arrayS[k]>=arrayS[k+1])
996                                                         break;
997                                                 double tmp=arrayS[k];
998                                                 arrayS[k]=arrayS[k+1];
999                                                 arrayS[k+1]=tmp;
1000                                                if(k<Nm1) {
1001                                                        for(int i=0;i<N;i++) {
1002                                                                tmp=arrayU[i][k+1];
1003                                                                arrayU[i][k+1]=arrayU[i][k];
1004                                                                arrayU[i][k]=tmp;
1005                                                                tmp=arrayV[i][k+1];
1006                                                                arrayV[i][k+1]=arrayV[i][k];
1007                                                                arrayV[i][k]=tmp;
1008                                                        }
1009                                                }
1010                                                k++;
1011                                        }
1012                                        iter=0;
1013                                        p--;
1014                                        } break;
1015                        }
1016                }
1017                final AbstractDoubleSquareMatrix svd[]=new AbstractDoubleSquareMatrix[3];
1018                svd[0]=new DoubleSquareMatrix(arrayU);
1019                svd[1]=new DoubleDiagonalMatrix(arrayS);
1020                svd[2]=new DoubleSquareMatrix(arrayV);
1021                return svd;
1022        }
1023
1024// POLAR DECOMPOSITION
1025
1026        /**
1027        * Returns the polar decomposition of this matrix.
1028        */
1029        public AbstractDoubleSquareMatrix[] polarDecompose() {
1030                final int N=numRows;
1031                final AbstractDoubleVector evec[]=new AbstractDoubleVector[N];
1032                double eval[];
1033                try {
1034                        eval=LinearMath.eigenSolveSymmetric(this,evec);
1035                } catch(MaximumIterationsExceededException e) {
1036                        return null;
1037                }
1038                final double tmpa[][]=new double[N][N];
1039                final double tmpm[][]=new double[N][N];
1040                double abs;
1041                for(int i=0;i<N;i++) {
1042                        abs=Math.abs(eval[i]);
1043                        tmpa[i][0]=eval[i]*evec[i].getComponent(0)/abs;
1044                        tmpm[i][0]=abs*evec[i].getComponent(0);
1045                        for(int j=1;j<N;j++) {
1046                                tmpa[i][j]=eval[i]*evec[i].getComponent(j)/abs;
1047                                tmpm[i][j]=abs*evec[i].getComponent(j);
1048                        }
1049                }
1050                final double arg[][]=new double[N][N];
1051                final double mod[][]=new double[N][N];
1052                for(int i=0;i<N;i++) {
1053                        for(int j=0;j<N;j++) {
1054                                arg[i][j]=evec[0].getComponent(i)*tmpa[0][j];
1055                                mod[i][j]=evec[0].getComponent(i)*tmpm[0][j];
1056                                for(int k=1;k<N;k++) {
1057                                        arg[i][j]+=evec[k].getComponent(i)*tmpa[k][j];
1058                                        mod[i][j]+=evec[k].getComponent(i)*tmpm[k][j];
1059                                }
1060                        }
1061                }
1062                final AbstractDoubleSquareMatrix us[]=new AbstractDoubleSquareMatrix[2];
1063                us[0]=new DoubleSquareMatrix(arg);
1064                us[1]=new DoubleSquareMatrix(mod);
1065                return us;
1066        }
1067
1068// MAP ELEMENTS
1069
1070        /**
1071        * Applies a function on all the matrix elements.
1072        * @param f a user-defined function
1073        * @return a double matrix
1074        */
1075        public AbstractDoubleMatrix mapElements(final Mapping f) {
1076                final double array[][]=new double[numRows][numCols];
1077                for(int i=0;i<numRows;i++) {
1078                        array[i][0] = f.map(matrix[i][0]);
1079                        for(int j=1;j<numCols;j++)
1080                                array[i][j] = f.map(matrix[i][j]);
1081                }
1082                return new DoubleSquareMatrix(array);
1083        }
1084}
1085
A to Z: JavaDoc & Examples Daily Java News & Articles Open Source Projects Open Source Codes Free Computer Books Remove Frame
Popular Tags