KickJava   Java API By Example, From Geeks To Geeks.

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


1 /* AUTO-GENERATED */
2 package JSci.maths.matrices;
3
4 import JSci.maths.ArrayMath;
5 import JSci.maths.Complex;
6 import JSci.maths.ComplexMapping;
7 import JSci.maths.LinearMath;
8 import JSci.maths.DimensionException;
9 import JSci.maths.MaximumIterationsExceededException;
10 import JSci.maths.vectors.AbstractComplexVector;
11 import JSci.maths.vectors.ComplexVector;
12 import JSci.maths.groups.AbelianGroup;
13
14 /**
15 * The ComplexSquareMatrix class provides an object for encapsulating square matrices containing complex numbers.
16 * @version 2.2
17 * @author Mark Hale
18 */
19 public class ComplexSquareMatrix extends AbstractComplexSquareMatrix {
20         /**
21         * Arrays containing the elements of the matrix.
22         */
23         protected final double matrixRe[][],matrixIm[][];
24         /**
25         * Constructs a matrix by wrapping two arrays.
26         * @param arrayRe an array of real values
27         * @param arrayIm an array of imaginary values
28         * @exception MatrixDimensionException If the array is not square.
29         */
30         public ComplexSquareMatrix(final double arrayRe[][],final double arrayIm[][]) {
31                 super(arrayRe.length);
32                 if(!ArrayMath.isSquare(arrayRe))
33                         throw new MatrixDimensionException("Array is not square.");
34                 if(!ArrayMath.isSquare(arrayIm))
35                         throw new MatrixDimensionException("Array is not square.");
36                 matrixRe=arrayRe;
37                 matrixIm=arrayIm;
38         }
39         /**
40         * Constructs an empty matrix.
41         * @param size the number of rows/columns
42         */
43         public ComplexSquareMatrix(final int size) {
44                 this(new double[size][size], new double[size][size]);
45         }
46         /**
47         * Constructs a matrix from an array.
48         * @param array an assigned value
49         * @exception MatrixDimensionException If the array is not square.
50         */
51         public ComplexSquareMatrix(final Complex array[][]) {
52                 this(array.length);
53                 for(int i=0;i<numRows;i++) {
54                         if(array[i].length != array.length)
55                                 throw new MatrixDimensionException("Array is not square.");
56                         for(int j=0;j<numCols;j++) {
57                                 matrixRe[i][j]=array[i][j].real();
58                                 matrixIm[i][j]=array[i][j].imag();
59                         }
60                 }
61         }
62         /**
63         * Constructs a matrix from an array of vectors.
64         * The vectors form columns in the matrix.
65         * @param array an assigned value.
66         * @exception MatrixDimensionException If the array is not square.
67         */
68         public ComplexSquareMatrix(final AbstractComplexVector array[]) {
69                 this(array.length);
70                 for(int i=0;i<numRows;i++) {
71                         for(int j=0;j<numCols;j++) {
72                                 matrixRe[i][j]=array[j].getComponent(i).real();
73                                 matrixIm[i][j]=array[j].getComponent(i).imag();
74                         }
75                 }
76         }
77         /**
78         * Compares two complex matrices for equality.
79         * @param m a complex matrix
80         */
81         public boolean equals(AbstractComplexMatrix m, double tol) {
82                 if(m != null && numRows == m.rows() && numCols == m.columns()) {
83             double sumSqr = 0.0;
84                         for(int i=0;i<numRows;i++) {
85                                 for(int j=0;j<numCols;j++) {
86                     double deltaRe = matrixRe[i][j]-m.getRealElement(i,j);
87                     double deltaIm = matrixIm[i][j]-m.getImagElement(i,j);
88                     sumSqr += deltaRe*deltaRe + deltaIm*deltaIm;
89                                 }
90                         }
91                         return (sumSqr <= tol*tol);
92                 } else {
93                         return false;
94                 }
95         }
96         /**
97         * Returns a string representing this matrix.
98         */
99         public String JavaDoc toString() {
100                 final StringBuffer JavaDoc buf=new StringBuffer JavaDoc(5*numRows*numCols);
101                 for(int j,i=0;i<numRows;i++) {
102                         for(j=0;j<numCols;j++) {
103                                 buf.append(Complex.toString(matrixRe[i][j],matrixIm[i][j]));
104                                 buf.append(' ');
105                         }
106                         buf.append('\n');
107                 }
108                 return buf.toString();
109         }
110         /**
111         * Returns the real part of this complex matrix.
112         * @return a double square matrix
113         */
114         public AbstractDoubleMatrix real() {
115                 return new DoubleSquareMatrix(matrixRe);
116         }
117         /**
118         * Returns the imaginary part of this complex matrix.
119         * @return a double square matrix
120         */
121         public AbstractDoubleMatrix imag() {
122                 return new DoubleSquareMatrix(matrixIm);
123         }
124         /**
125         * Returns an element of the matrix.
126         * @param i row index of the element
127         * @param j column index of the element
128         * @exception MatrixDimensionException If attempting to access an invalid element.
129         */
130         public Complex getElement(final int i, final int j) {
131                 if(i>=0 && i<numRows && j>=0 && j<numCols)
132                         return new Complex(matrixRe[i][j],matrixIm[i][j]);
133                 else
134                         throw new MatrixDimensionException(getInvalidElementMsg(i,j));
135         }
136         public double getRealElement(final int i, final int j) {
137                 if(i>=0 && i<numRows && j>=0 && j<numCols)
138                         return matrixRe[i][j];
139                 else
140                         throw new MatrixDimensionException(getInvalidElementMsg(i,j));
141         }
142         public double getImagElement(final int i, final int j) {
143                 if(i>=0 && i<numRows && j>=0 && j<numCols)
144                         return matrixIm[i][j];
145                 else
146                         throw new MatrixDimensionException(getInvalidElementMsg(i,j));
147         }
148         /**
149         * Sets the value of an element of the matrix.
150         * Should only be used to initialise this matrix.
151         * @param i row index of the element
152         * @param j column index of the element
153         * @param z a complex number
154         * @exception MatrixDimensionException If attempting to access an invalid element.
155         */
156         public void setElement(final int i, final int j, final Complex z) {
157                 if(i>=0 && i<numRows && j>=0 && j<numCols) {
158                         matrixRe[i][j]=z.real();
159                         matrixIm[i][j]=z.imag();
160                 } else
161                         throw new MatrixDimensionException(getInvalidElementMsg(i,j));
162         }
163         /**
164         * Sets the value of an element of the matrix.
165         * Should only be used to initialise this matrix.
166         * @param i row index of the element
167         * @param j column index of the element
168         * @param x the real part of a complex number
169         * @param y the imaginary part of a complex number
170         * @exception MatrixDimensionException If attempting to access an invalid element.
171         */
172         public void setElement(final int i, final int j, final double x, final double y) {
173                 if(i>=0 && i<numRows && j>=0 && j<numCols) {
174                         matrixRe[i][j]=x;
175                         matrixIm[i][j]=y;
176                 } else
177                         throw new MatrixDimensionException(getInvalidElementMsg(i,j));
178         }
179         /**
180         * Returns the l<sup><img border=0 alt="infinity" SRC="doc-files/infinity.gif"></sup>-norm.
181         * @author Taber Smith
182         */
183         public double infNorm() {
184                 double result=0.0,tmpResult;
185                 for(int i=0;i<numRows;i++) {
186                         tmpResult=0.0;
187                         for(int j=0;j<numCols;j++)
188                                 tmpResult += Math.sqrt((matrixRe[i][j]*matrixRe[i][j] + matrixIm[i][j]*matrixIm[i][j]));
189                         if(tmpResult>result)
190                                 result=tmpResult;
191                 }
192                 return result;
193         }
194         /**
195         * Returns the Frobenius or Hilbert-Schmidt (l<sup>2</sup>) norm.
196         * @jsci.planetmath FrobeniusMatrixNorm
197         * @author Taber Smith
198         */
199         public double frobeniusNorm() {
200                 double result=0.0;
201                 for(int j,i=0;i<numRows;i++)
202                         for(j=0;j<numCols;j++)
203                                 result+=matrixRe[i][j]*matrixRe[i][j]+matrixIm[i][j]*matrixIm[i][j];
204                 return Math.sqrt(result);
205         }
206         /**
207         * Returns the determinant.
208         */
209         public Complex det() {
210                 if(numRows==2) {
211                         return new Complex(
212                                 matrixRe[0][0]*matrixRe[1][1]-matrixIm[0][0]*matrixIm[1][1]-matrixRe[0][1]*matrixRe[1][0]+matrixIm[0][1]*matrixIm[1][0],
213                                 matrixRe[0][0]*matrixIm[1][1]+matrixIm[0][0]*matrixRe[1][1]-matrixRe[0][1]*matrixIm[1][0]-matrixIm[0][1]*matrixRe[1][0]
214                         );
215                 } else {
216                         final ComplexSquareMatrix lu[] = (ComplexSquareMatrix[]) this.luDecompose(null);
217                         double tmp;
218                         double detRe=lu[1].matrixRe[0][0];
219                         double detIm=lu[1].matrixIm[0][0];
220                         for(int i=1;i<numRows;i++) {
221                                 tmp=detRe*lu[1].matrixRe[i][i]-detIm*lu[1].matrixIm[i][i];
222                                 detIm=detRe*lu[1].matrixIm[i][i]+detIm*lu[1].matrixRe[i][i];
223                                 detRe=tmp;
224                         }
225                         return new Complex(detRe*LUpivot[numRows],detIm*LUpivot[numRows]);
226                 }
227         }
228         /**
229         * Returns the trace.
230         */
231         public Complex trace() {
232                 double trRe=matrixRe[0][0];
233                 double trIm=matrixIm[0][0];
234                 for(int i=1;i<numRows;i++) {
235                         trRe+=matrixRe[i][i];
236                         trIm+=matrixIm[i][i];
237                 }
238                 return new Complex(trRe,trIm);
239         }
240
241 //============
242 // OPERATIONS
243 //============
244
245         /**
246         * Returns the negative of this matrix.
247         */
248         public AbelianGroup.Member negate() {
249                 final double arrayRe[][]=new double[numRows][numCols];
250                 final double arrayIm[][]=new double[numRows][numCols];
251                 for(int j,i=0;i<numRows;i++) {
252                         arrayRe[i][0]=-matrixRe[i][0];
253                         arrayIm[i][0]=-matrixIm[i][0];
254                         for(j=1;j<numCols;j++) {
255                                 arrayRe[i][j]=-matrixRe[i][j];
256                                 arrayIm[i][j]=-matrixIm[i][j];
257                         }
258                 }
259                 return new ComplexSquareMatrix(arrayRe,arrayIm);
260         }
261
262 // ADDITION
263
264         /**
265         * Returns the addition of this matrix and another.
266         * @param m a complex matrix
267         * @exception MatrixDimensionException If the matrices are different sizes.
268         */
269         public AbstractComplexSquareMatrix add(final AbstractComplexSquareMatrix m) {
270                 if(m instanceof ComplexSquareMatrix)
271                         return add((ComplexSquareMatrix)m);
272
273                 if(numRows==m.rows() && numCols==m.columns()) {
274                         final double arrayRe[][]=new double[numRows][numCols];
275                         final double arrayIm[][]=new double[numRows][numCols];
276                         for(int j,i=0;i<numRows;i++) {
277                                 arrayRe[i][0]=matrixRe[i][0]+m.getElement(i,0).real();
278                                 arrayIm[i][0]=matrixIm[i][0]+m.getElement(i,0).imag();
279                                 for(j=1;j<numCols;j++) {
280                                         arrayRe[i][j]=matrixRe[i][j]+m.getElement(i,j).real();
281                                         arrayIm[i][j]=matrixIm[i][j]+m.getElement(i,j).imag();
282                                 }
283                         }
284                         return new ComplexSquareMatrix(arrayRe,arrayIm);
285                 } else
286                         throw new MatrixDimensionException("Matrices are different sizes.");
287         }
288         public ComplexSquareMatrix add(final ComplexSquareMatrix m) {
289                 if(numRows==m.numRows && numCols==m.numCols) {
290                         final double arrayRe[][]=new double[numRows][numCols];
291                         final double arrayIm[][]=new double[numRows][numCols];
292                         for(int j,i=0;i<numRows;i++) {
293                                 arrayRe[i][0]=matrixRe[i][0]+m.matrixRe[i][0];
294                                 arrayIm[i][0]=matrixIm[i][0]+m.matrixIm[i][0];
295                                 for(j=1;j<numCols;j++) {
296                                         arrayRe[i][j]=matrixRe[i][j]+m.matrixRe[i][j];
297                                         arrayIm[i][j]=matrixIm[i][j]+m.matrixIm[i][j];
298                                 }
299                         }
300                         return new ComplexSquareMatrix(arrayRe,arrayIm);
301                 } else
302                         throw new MatrixDimensionException("Matrices are different sizes.");
303         }
304
305 // SUBTRACTION
306
307         /**
308         * Returns the subtraction of this matrix by another.
309         * @param m a complex matrix
310         * @exception MatrixDimensionException If the matrices are different sizes.
311         */
312         public AbstractComplexSquareMatrix subtract(final AbstractComplexSquareMatrix m) {
313                 if(m instanceof ComplexSquareMatrix)
314                         return subtract((ComplexSquareMatrix)m);
315
316                 if(numRows==m.rows() && numCols==m.columns()) {
317                         final double arrayRe[][]=new double[numRows][numCols];
318                         final double arrayIm[][]=new double[numRows][numCols];
319                         for(int i=0;i<numRows;i++) {
320                                 arrayRe[i][0]=matrixRe[i][0]-m.getElement(i,0).real();
321                                 arrayIm[i][0]=matrixIm[i][0]-m.getElement(i,0).imag();
322                                 for(int j=1;j<numCols;j++) {
323                                         arrayRe[i][j]=matrixRe[i][j]-m.getElement(i,j).real();
324                                         arrayIm[i][j]=matrixIm[i][j]-m.getElement(i,j).imag();
325                                 }
326                         }
327                         return new ComplexSquareMatrix(arrayRe,arrayIm);
328                 } else
329                         throw new MatrixDimensionException("Matrices are different sizes.");
330         }
331         public ComplexSquareMatrix subtract(final ComplexSquareMatrix m) {
332                 if(numRows==m.numRows && numCols==m.numCols) {
333                         final double arrayRe[][]=new double[numRows][numCols];
334                         final double arrayIm[][]=new double[numRows][numCols];
335                         for(int j,i=0;i<numRows;i++) {
336                                 arrayRe[i][0]=matrixRe[i][0]-m.matrixRe[i][0];
337                                 arrayIm[i][0]=matrixIm[i][0]-m.matrixIm[i][0];
338                                 for(j=1;j<numCols;j++) {
339                                         arrayRe[i][j]=matrixRe[i][j]-m.matrixRe[i][j];
340                                         arrayIm[i][j]=matrixIm[i][j]-m.matrixIm[i][j];
341                                 }
342                         }
343                         return new ComplexSquareMatrix(arrayRe,arrayIm);
344                 } else
345                         throw new MatrixDimensionException("Matrices are different sizes.");
346         }
347
348 // SCALAR MULTIPLICATION
349
350         /**
351         * Returns the multiplication of this matrix by a scalar.
352         * @param z a complex number
353         * @return a complex square matrix
354         */
355         public AbstractComplexMatrix scalarMultiply(final Complex z) {
356                 final double real=z.real();
357                 final double imag=z.imag();
358                 final double arrayRe[][]=new double[numRows][numCols];
359                 final double arrayIm[][]=new double[numRows][numCols];
360                 for(int i=0;i<numRows;i++) {
361                         arrayRe[i][0]=matrixRe[i][0]*real-matrixIm[i][0]*imag;
362                         arrayIm[i][0]=matrixRe[i][0]*imag+matrixIm[i][0]*real;
363                         for(int j=1;j<numCols;j++) {
364                                 arrayRe[i][j]=matrixRe[i][j]*real-matrixIm[i][j]*imag;
365                                 arrayIm[i][j]=matrixRe[i][j]*imag+matrixIm[i][j]*real;
366                         }
367                 }
368                 return new ComplexSquareMatrix(arrayRe,arrayIm);
369         }
370         /**
371         * Returns the multiplication of this matrix by a scalar.
372         * @param x a double
373         * @return a complex square matrix
374         */
375         public AbstractComplexMatrix scalarMultiply(final double x) {
376                 final double arrayRe[][]=new double[numRows][numCols];
377                 final double arrayIm[][]=new double[numRows][numCols];
378                 for(int i=0;i<numRows;i++) {
379                         arrayRe[i][0]=x*matrixRe[i][0];
380                         arrayIm[i][0]=x*matrixIm[i][0];
381                         for(int j=1;j<numCols;j++) {
382                                 arrayRe[i][j]=x*matrixRe[i][j];
383                                 arrayIm[i][j]=x*matrixIm[i][j];
384                         }
385                 }
386                 return new ComplexSquareMatrix(arrayRe,arrayIm);
387         }
388
389 // MATRIX MULTIPLICATION
390
391         /**
392         * Returns the multiplication of a vector by this matrix.
393         * @param v a complex vector
394         * @exception DimensionException If the matrix and vector are incompatible.
395         */
396         public AbstractComplexVector multiply(final AbstractComplexVector v) {
397                 if(numCols==v.dimension()) {
398                         final double arrayRe[]=new double[numRows];
399                         final double arrayIm[]=new double[numRows];
400                         Complex comp;
401                         for(int j,i=0;i<numRows;i++) {
402                                 comp=v.getComponent(0);
403                                 arrayRe[i]=(matrixRe[i][0]*comp.real() - matrixIm[i][0]*comp.imag());
404                                 arrayIm[i]=(matrixIm[i][0]*comp.real() + matrixRe[i][0]*comp.imag());
405                                 for(j=1;j<numCols;j++) {
406                                         comp=v.getComponent(j);
407                                         arrayRe[i]+=(matrixRe[i][j]*comp.real() - matrixIm[i][j]*comp.imag());
408                                         arrayIm[i]+=(matrixIm[i][j]*comp.real() + matrixRe[i][j]*comp.imag());
409                                 }
410                         }
411                         return new ComplexVector(arrayRe,arrayIm);
412                 } else
413                         throw new DimensionException("Matrix and vector are incompatible.");
414         }
415         /**
416         * Returns the multiplication of this matrix and another.
417         * @param m a complex square matrix
418         * @exception MatrixDimensionException If the matrices are incompatible.
419         */
420         public AbstractComplexSquareMatrix multiply(final AbstractComplexSquareMatrix m) {
421                 if(m instanceof ComplexSquareMatrix)
422                         return multiply((ComplexSquareMatrix)m);
423
424                 if(numCols==m.rows()) {
425                         final double arrayRe[][]=new double[numRows][numCols];
426                         final double arrayIm[][]=new double[numRows][numCols];
427                         Complex elem;
428                         for(int j=0;j<numRows;j++) {
429                                 for(int k=0;k<numCols;k++) {
430                                         elem=m.getElement(0,k);
431                                         arrayRe[j][k]=(matrixRe[j][0]*elem.real() - matrixIm[j][0]*elem.imag());
432                                         arrayIm[j][k]=(matrixIm[j][0]*elem.real() + matrixRe[j][0]*elem.imag());
433                                         for(int n=1;n<numCols;n++) {
434                                                 elem=m.getElement(n,k);
435                                                 arrayRe[j][k]+=(matrixRe[j][n]*elem.real() - matrixIm[j][n]*elem.imag());
436                                                 arrayIm[j][k]+=(matrixIm[j][n]*elem.real() + matrixRe[j][n]*elem.imag());
437                                         }
438                                 }
439                         }
440                         return new ComplexSquareMatrix(arrayRe,arrayIm);
441                 } else
442                         throw new MatrixDimensionException("Incompatible matrices.");
443         }
444         public ComplexSquareMatrix multiply(final ComplexSquareMatrix m) {
445                 if(numCols==m.numRows) {
446                         final double arrayRe[][]=new double[numRows][numCols];
447                         final double arrayIm[][]=new double[numRows][numCols];
448                         for(int j=0;j<numRows;j++) {
449                                 for(int k=0;k<numCols;k++) {
450                                         arrayRe[j][k]=(matrixRe[j][0]*m.matrixRe[0][k] - matrixIm[j][0]*m.matrixIm[0][k]);
451                                         arrayIm[j][k]=(matrixIm[j][0]*m.matrixRe[0][k] + matrixRe[j][0]*m.matrixIm[0][k]);
452                                         for(int n=1;n<numCols;n++) {
453                                                 arrayRe[j][k]+=(matrixRe[j][n]*m.matrixRe[n][k] - matrixIm[j][n]*m.matrixIm[n][k]);
454                                                 arrayIm[j][k]+=(matrixIm[j][n]*m.matrixRe[n][k] + matrixRe[j][n]*m.matrixIm[n][k]);
455                                         }
456                                 }
457                         }
458                         return new ComplexSquareMatrix(arrayRe,arrayIm);
459                 } else
460                         throw new MatrixDimensionException("Incompatible matrices.");
461         }
462
463 // DIRECT SUM
464
465         /**
466         * Returns the direct sum of this matrix and another.
467         */
468         public AbstractComplexSquareMatrix directSum(final AbstractComplexSquareMatrix m) {
469                 final double arrayRe[][]=new double[numRows+m.numRows][numCols+m.numCols];
470                 final double arrayIm[][]=new double[numRows+m.numRows][numCols+m.numCols];
471                 for(int i=0;i<numRows;i++) {
472                         for(int j=0;j<numCols;j++) {
473                                 arrayRe[i][j]=matrixRe[i][j];
474                                 arrayIm[i][j]=matrixIm[i][j];
475                         }
476                 }
477                 for(int i=0;i<m.numRows;i++) {
478                         for(int j=0;j<m.numCols;j++) {
479                                 Complex elem=m.getElement(i,j);
480                                 arrayRe[i+numRows][j+numCols]=elem.real();
481                                 arrayIm[i+numRows][j+numCols]=elem.imag();
482                         }
483                 }
484                 return new ComplexSquareMatrix(arrayRe,arrayIm);
485         }
486
487 // TENSOR PRODUCT
488
489         /**
490         * Returns the tensor product of this matrix and another.
491         */
492         public AbstractComplexSquareMatrix tensor(final AbstractComplexSquareMatrix m) {
493                 final double arrayRe[][]=new double[numRows*m.numRows][numCols*m.numCols];
494                 final double arrayIm[][]=new double[numRows*m.numRows][numCols*m.numCols];
495                 for(int i=0;i<numRows;i++) {
496                         for(int j=0;j<numCols;j++) {
497                                 for(int k=0;k<m.numRows;j++) {
498                                         for(int l=0;l<m.numCols;l++) {
499                                                 Complex elem=m.getElement(k,l);
500                                                 arrayRe[i*m.numRows+k][j*m.numCols+l]=(matrixRe[i][j]*elem.real() - matrixIm[i][j]*elem.imag());
501                                                 arrayIm[i*m.numRows+k][j*m.numCols+l]=(matrixIm[i][j]*elem.real() + matrixRe[i][j]*elem.imag());
502                                         }
503                                 }
504                         }
505                 }
506                 return new ComplexSquareMatrix(arrayRe,arrayIm);
507         }
508
509 // HERMITIAN ADJOINT
510
511         /**
512         * Returns the hermitian adjoint of this matrix.
513         * @return a complex square matrix
514         */
515         public AbstractComplexMatrix hermitianAdjoint() {
516                 final double arrayRe[][]=new double[numCols][numRows];
517                 final double arrayIm[][]=new double[numCols][numRows];
518                 for(int j,i=0;i<numRows;i++) {
519                         arrayRe[0][i]=matrixRe[i][0];
520                         arrayIm[0][i]=-matrixIm[i][0];
521                         for(j=1;j<numCols;j++) {
522                                 arrayRe[j][i]=matrixRe[i][j];
523                                 arrayIm[j][i]=-matrixIm[i][j];
524                         }
525                 }
526                 return new ComplexSquareMatrix(arrayRe,arrayIm);
527         }
528
529 // CONJUGATE
530
531         /**
532         * Returns the complex conjugate of this matrix.
533         * @return a complex square matrix
534         */
535         public AbstractComplexMatrix conjugate() {
536                 final double arrayIm[][]=new double[numRows][numCols];
537                 for(int j,i=0;i<numRows;i++) {
538                         arrayIm[i][0]=-matrixIm[i][0];
539                         for(j=1;j<numCols;j++)
540                                 arrayIm[i][j]=-matrixIm[i][j];
541                 }
542                 return new ComplexSquareMatrix(matrixRe,arrayIm);
543         }
544
545 // TRANSPOSE
546
547         /**
548         * Returns the transpose of this matrix.
549         * @return a complex square matrix
550         */
551         public Matrix transpose() {
552                 final double arrayRe[][]=new double[numCols][numRows];
553                 final double arrayIm[][]=new double[numCols][numRows];
554                 for(int j,i=0;i<numRows;i++) {
555                         arrayRe[0][i]=matrixRe[i][0];
556                         arrayIm[0][i]=matrixIm[i][0];
557                         for(j=1;j<numCols;j++) {
558                                 arrayRe[j][i]=matrixRe[i][j];
559                                 arrayIm[j][i]=matrixIm[i][j];
560                         }
561                 }
562                 return new ComplexSquareMatrix(arrayRe,arrayIm);
563         }
564
565 // INVERSE
566
567         /**
568         * Returns the inverse of this matrix.
569         * @return a complex square matrix
570         */
571         public AbstractComplexSquareMatrix inverse() {
572                 final int N=numRows;
573                 final double arrayLRe[][]=new double[N][N];
574                 final double arrayLIm[][]=new double[N][N];
575                 final double arrayURe[][]=new double[N][N];
576                 final double arrayUIm[][]=new double[N][N];
577                 final ComplexSquareMatrix lu[] = (ComplexSquareMatrix[]) this.luDecompose(null);
578                 double denom;
579                 denom=lu[0].matrixRe[0][0]*lu[0].matrixRe[0][0]+lu[0].matrixIm[0][0]*lu[0].matrixIm[0][0];
580                 arrayLRe[0][0]=lu[0].matrixRe[0][0]/denom;
581                 arrayLIm[0][0]=-lu[0].matrixIm[0][0]/denom;
582                 denom=lu[1].matrixRe[0][0]*lu[1].matrixRe[0][0]+lu[1].matrixIm[0][0]*lu[1].matrixIm[0][0];
583                 arrayURe[0][0]=lu[1].matrixRe[0][0]/denom;
584                 arrayUIm[0][0]=-lu[1].matrixIm[0][0]/denom;
585                 for(int i=1;i<N;i++) {
586                         denom=lu[0].matrixRe[i][i]*lu[0].matrixRe[i][i]+lu[0].matrixIm[i][i]*lu[0].matrixIm[i][i];
587                         arrayLRe[i][i]=lu[0].matrixRe[i][i]/denom;
588                         arrayLIm[i][i]=-lu[0].matrixIm[i][i]/denom;
589                         denom=lu[1].matrixRe[i][i]*lu[1].matrixRe[i][i]+lu[1].matrixIm[i][i]*lu[1].matrixIm[i][i];
590                         arrayURe[i][i]=lu[1].matrixRe[i][i]/denom;
591                         arrayUIm[i][i]=-lu[1].matrixIm[i][i]/denom;
592                 }
593                 for(int i=0;i<N-1;i++) {
594                         for(int j=i+1;j<N;j++) {
595                                 double tmpLRe=0.0, tmpLIm=0.0;
596                                 double tmpURe=0.0, tmpUIm=0.0;
597                                 for(int k=i;k<j;k++) {
598                                         tmpLRe-=(lu[0].matrixRe[j][k]*arrayLRe[k][i] - lu[0].matrixIm[j][k]*arrayLIm[k][i]);
599                                         tmpLIm-=(lu[0].matrixIm[j][k]*arrayLRe[k][i] + lu[0].matrixRe[j][k]*arrayLIm[k][i]);
600                                         tmpURe-=(arrayURe[i][k]*lu[1].matrixRe[k][j] - arrayUIm[i][k]*lu[1].matrixIm[k][j]);
601                                         tmpUIm-=(arrayUIm[i][k]*lu[1].matrixRe[k][j] + arrayURe[i][k]*lu[1].matrixIm[k][j]);
602                                 }
603                                 denom=lu[0].matrixRe[j][j]*lu[0].matrixRe[j][j]+lu[0].matrixIm[j][j]*lu[0].matrixIm[j][j];
604                                 arrayLRe[j][i]=(tmpLRe*lu[0].matrixRe[j][j]+tmpLIm*lu[0].matrixIm[j][j])/denom;
605                                 arrayLIm[j][i]=(tmpLIm*lu[0].matrixRe[j][j]-tmpLRe*lu[0].matrixIm[j][j])/denom;
606                                 denom=lu[1].matrixRe[j][j]*lu[1].matrixRe[j][j]+lu[1].matrixIm[j][j]*lu[1].matrixIm[j][j];
607                                 arrayURe[i][j]=(tmpURe*lu[1].matrixRe[j][j]+tmpUIm*lu[1].matrixIm[j][j])/denom;
608                                 arrayUIm[i][j]=(tmpUIm*lu[1].matrixRe[j][j]-tmpURe*lu[1].matrixIm[j][j])/denom;
609                         }
610                 }
611                 // matrix multiply arrayU x arrayL
612 final double invRe[][]=new double[N][N];
613                 final double invIm[][]=new double[N][N];
614                 for(int i=0;i<N;i++) {
615                         for(int j=0;j<i;j++) {
616                                 for(int k=i;k<N;k++) {
617                                         invRe[i][LUpivot[j]]+=(arrayURe[i][k]*arrayLRe[k][j] - arrayUIm[i][k]*arrayLIm[k][j]);
618                                         invIm[i][LUpivot[j]]+=(arrayUIm[i][k]*arrayLRe[k][j] + arrayURe[i][k]*arrayLIm[k][j]);
619                                 }
620                         }
621                         for(int j=i;j<N;j++) {
622                                 for(int k=j;k<N;k++) {
623                                         invRe[i][LUpivot[j]]+=(arrayURe[i][k]*arrayLRe[k][j] - arrayUIm[i][k]*arrayLIm[k][j]);
624                                         invIm[i][LUpivot[j]]+=(arrayUIm[i][k]*arrayLRe[k][j] + arrayURe[i][k]*arrayLIm[k][j]);
625                                 }
626                         }
627                 }
628                 return new ComplexSquareMatrix(invRe,invIm);
629         }
630
631 // LU DECOMPOSITION
632
633         /**
634         * Returns the LU decomposition of this matrix.
635         * @return an array with [0] containing the L-matrix and [1] containing the U-matrix.
636         */
637         public final AbstractComplexSquareMatrix[] luDecompose(int pivot[]) {
638                 if(LU!=null) {
639                         if(pivot!=null)
640                                 System.arraycopy(LUpivot,0,pivot,0,pivot.length);
641                         return LU;
642                 }
643                 final int N=numRows;
644                 final double arrayLRe[][]=new double[N][N];
645                 final double arrayLIm[][]=new double[N][N];
646                 final double arrayURe[][]=new double[N][N];
647                 final double arrayUIm[][]=new double[N][N];
648                 if(pivot==null)
649                         pivot=new int[N+1];
650                 for(int i=0;i<N;i++)
651                         pivot[i]=i;
652                 pivot[N]=1;
653         // LU decomposition to arrayU
654 for(int j=0;j<N;j++) {
655                         for(int i=0;i<j;i++) {
656                                 double tmpRe=matrixRe[pivot[i]][j];
657                                 double tmpIm=matrixIm[pivot[i]][j];
658                                 for(int k=0;k<i;k++) {
659                                         tmpRe-=(arrayURe[i][k]*arrayURe[k][j] - arrayUIm[i][k]*arrayUIm[k][j]);
660                                         tmpIm-=(arrayUIm[i][k]*arrayURe[k][j] + arrayURe[i][k]*arrayUIm[k][j]);
661                                 }
662                                 arrayURe[i][j]=tmpRe;
663                                 arrayUIm[i][j]=tmpIm;
664                         }
665                         double max=0.0;
666                         int pivotrow=j;
667                         for(int i=j;i<N;i++) {
668                                 double tmpRe=matrixRe[pivot[i]][j];
669                                 double tmpIm=matrixIm[pivot[i]][j];
670                                 for(int k=0;k<j;k++) {
671                                         tmpRe-=(arrayURe[i][k]*arrayURe[k][j] - arrayUIm[i][k]*arrayUIm[k][j]);
672                                         tmpIm-=(arrayUIm[i][k]*arrayURe[k][j] + arrayURe[i][k]*arrayUIm[k][j]);
673                                 }
674                                 arrayURe[i][j]=tmpRe;
675                                 arrayUIm[i][j]=tmpIm;
676                         // while we're here search for a pivot for arrayU[j][j]
677 double tmp=tmpRe*tmpRe+tmpIm*tmpIm;
678                                 if(tmp>max) {
679                                         max=tmp;
680                                         pivotrow=i;
681                                 }
682                         }
683                 // swap row j with pivotrow
684 if(pivotrow!=j) {
685                                 double[] tmprow = arrayURe[j];
686                                 arrayURe[j] = arrayURe[pivotrow];
687                                 arrayURe[pivotrow] = tmprow;
688                                 tmprow = arrayUIm[j];
689                                 arrayUIm[j] = arrayUIm[pivotrow];
690                                 arrayUIm[pivotrow] = tmprow;
691                                 int k=pivot[j];
692                                 pivot[j]=pivot[pivotrow];
693                                 pivot[pivotrow]=k;
694                                 // update parity
695 pivot[N]=-pivot[N];
696                         }
697                 // divide by pivot
698 double tmpRe=arrayURe[j][j];
699                         double tmpIm=arrayUIm[j][j];
700                         if(Math.abs(tmpRe)<Math.abs(tmpIm)) {
701                                 double a=tmpRe/tmpIm;
702                                 double denom=tmpRe*a+tmpIm;
703                                 for(int i=j+1;i<N;i++) {
704                                         double tmp=(arrayURe[i][j]*a+arrayUIm[i][j])/denom;
705                                         arrayUIm[i][j]=(arrayUIm[i][j]*a-arrayURe[i][j])/denom;
706                                         arrayURe[i][j]=tmp;
707                                 }
708                         } else {
709                                 double a=tmpIm/tmpRe;
710                                 double denom=tmpRe+tmpIm*a;
711                                 for(int i=j+1;i<N;i++) {
712                                         double tmp=(arrayURe[i][j]+arrayUIm[i][j]*a)/denom;
713                                         arrayUIm[i][j]=(arrayUIm[i][j]-arrayURe[i][j]*a)/denom;
714                                         arrayURe[i][j]=tmp;
715                                 }
716                         }
717                 }
718                 // move lower triangular part to arrayL
719 for(int j=0;j<N;j++) {
720                         arrayLRe[j][j]=1.0;
721                         for(int i=j+1;i<N;i++) {
722                                 arrayLRe[i][j]=arrayURe[i][j];
723                                 arrayLIm[i][j]=arrayUIm[i][j];
724                                 arrayURe[i][j]=0.0;
725                                 arrayUIm[i][j]=0.0;
726                         }
727                 }
728                 LU=new ComplexSquareMatrix[2];
729                 LU[0]=new ComplexSquareMatrix(arrayLRe,arrayLIm);
730                 LU[1]=new ComplexSquareMatrix(arrayURe,arrayUIm);
731                 LUpivot=new int[pivot.length];
732                 System.arraycopy(pivot,0,LUpivot,0,pivot.length);
733                 return LU;
734         }
735         /**
736         * Returns the LU decomposition of this matrix.
737         * Warning: no pivoting.
738         * @return an array with [0] containing the L-matrix
739         * and [1] containing the U-matrix.
740         * @jsci.planetmath LUDecomposition
741         */
742         public AbstractComplexSquareMatrix[] luDecompose() {
743                 final int N=numRows;
744                 final double arrayLRe[][]=new double[N][N];
745                 final double arrayLIm[][]=new double[N][N];
746                 final double arrayURe[][]=new double[N][N];
747                 final double arrayUIm[][]=new double[N][N];
748         // LU decomposition to arrayU
749 for(int j=0;j<N;j++) {
750                         for(int i=0;i<j;i++) {
751                                 double tmpRe=matrixRe[i][j];
752                                 double tmpIm=matrixIm[i][j];
753                                 for(int k=0;k<i;k++) {
754                                         tmpRe-=(arrayURe[i][k]*arrayURe[k][j] - arrayUIm[i][k]*arrayUIm[k][j]);
755                                         tmpIm-=(arrayUIm[i][k]*arrayURe[k][j] + arrayURe[i][k]*arrayUIm[k][j]);
756                                 }
757                                 arrayURe[i][j]=tmpRe;
758                                 arrayUIm[i][j]=tmpIm;
759                         }
760                         for(int i=j;i<N;i++) {
761                                 double tmpRe=matrixRe[i][j];
762                                 double tmpIm=matrixIm[i][j];
763                                 for(int k=0;k<j;k++) {
764                                         tmpRe-=(arrayURe[i][k]*arrayURe[k][j] - arrayUIm[i][k]*arrayUIm[k][j]);
765                                         tmpIm-=(arrayUIm[i][k]*arrayURe[k][j] + arrayURe[i][k]*arrayUIm[k][j]);
766                                 }
767                                 arrayURe[i][j]=tmpRe;
768                                 arrayUIm[i][j]=tmpIm;
769                         }
770                 // divide
771 double tmpRe=arrayURe[j][j];
772                         double tmpIm=arrayUIm[j][j];
773                         if(Math.abs(tmpRe)<Math.abs(tmpIm)) {
774                                 double a=tmpRe/tmpIm;
775                                 double denom=tmpRe*a+tmpIm;
776                                 for(int i=j+1;i<N;i++) {
777                                         double tmp=(arrayURe[i][j]*a+arrayUIm[i][j])/denom;
778                                         arrayUIm[i][j]=(arrayUIm[i][j]*a-arrayURe[i][j])/denom;
779                                         arrayURe[i][j]=tmp;
780                                 }
781                         } else {
782                                 double a=tmpIm/tmpRe;
783                                 double denom=tmpRe+tmpIm*a;
784                                 for(int i=j+1;i<N;i++) {
785                                         double tmp=(arrayURe[i][j]+arrayUIm[i][j]*a)/denom;
786                                         arrayUIm[i][j]=(arrayUIm[i][j]-arrayURe[i][j]*a)/denom;
787                                         arrayURe[i][j]=tmp;
788                                 }
789                         }
790                 }
791                 // move lower triangular part to arrayL
792 for(int j=0;j<N;j++) {
793                         arrayLRe[j][j]=1.0;
794                         for(int i=j+1;i<N;i++) {
795                                 arrayLRe[i][j]=arrayURe[i][j];
796                                 arrayLIm[i][j]=arrayUIm[i][j];
797                                 arrayURe[i][j]=0.0;
798                                 arrayUIm[i][j]=0.0;
799                         }
800                 }
801                 ComplexSquareMatrix[] lu=new ComplexSquareMatrix[2];
802                 lu[0]=new ComplexSquareMatrix(arrayLRe,arrayLIm);
803                 lu[1]=new ComplexSquareMatrix(arrayURe,arrayUIm);
804                 return lu;
805         }
806
807 // POLAR DECOMPOSITION
808
809         /**
810         * Returns the polar decomposition of this matrix.
811         */
812         public AbstractComplexSquareMatrix[] polarDecompose() {
813                 final int N=numRows;
814                 final AbstractComplexVector evec[]=new AbstractComplexVector[N];
815                 double eval[];
816                 try {
817                         eval=LinearMath.eigenSolveHermitian(this,evec);
818                 } catch(MaximumIterationsExceededException e) {
819                         return null;
820                 }
821                 final double tmpaRe[][]=new double[N][N];
822                 final double tmpaIm[][]=new double[N][N];
823                 final double tmpmRe[][]=new double[N][N];
824                 final double tmpmIm[][]=new double[N][N];
825                 double abs;
826                 Complex comp;
827                 for(int i=0;i<N;i++) {
828                         abs=Math.abs(eval[i]);
829                         comp=evec[i].getComponent(0).conjugate();
830                         tmpaRe[i][0]=eval[i]*comp.real()/abs;
831                         tmpaIm[i][0]=eval[i]*comp.imag()/abs;
832                         tmpmRe[i][0]=abs*comp.real();
833                         tmpmIm[i][0]=abs*comp.imag();
834                         for(int j=1;j<N;j++) {
835                                 comp=evec[i].getComponent(j).conjugate();
836                                 tmpaRe[i][j]=eval[i]*comp.real()/abs;
837                                 tmpaIm[i][j]=eval[i]*comp.imag()/abs;
838                                 tmpmRe[i][j]=abs*comp.real();
839                                 tmpmIm[i][j]=abs*comp.imag();
840                         }
841                 }
842                 final double argRe[][]=new double[N][N];
843                 final double argIm[][]=new double[N][N];
844                 final double modRe[][]=new double[N][N];
845                 final double modIm[][]=new double[N][N];
846                 for(int i=0;i<N;i++) {
847                         for(int j=0;j<N;j++) {
848                                 comp=evec[0].getComponent(i);
849                                 argRe[i][j]=(tmpaRe[0][j]*comp.real() - tmpaIm[0][j]*comp.imag());
850                                 argIm[i][j]=(tmpaIm[0][j]*comp.real() + tmpaRe[0][j]*comp.imag());
851                                 modRe[i][j]=(tmpmRe[0][j]*comp.real() - tmpmIm[0][j]*comp.imag());
852                                 modIm[i][j]=(tmpmIm[0][j]*comp.real() + tmpmRe[0][j]*comp.imag());
853                                 for(int k=1;k<N;k++) {
854                                         comp=evec[k].getComponent(i);
855                                         argRe[i][j]+=(tmpaRe[k][j]*comp.real() - tmpaIm[k][j]*comp.imag());
856                                         argIm[i][j]+=(tmpaIm[k][j]*comp.real() + tmpaRe[k][j]*comp.imag());
857                                         modRe[i][j]+=(tmpmRe[k][j]*comp.real() - tmpmIm[k][j]*comp.imag());
858                                         modIm[i][j]+=(tmpmIm[k][j]*comp.real() + tmpmRe[k][j]*comp.imag());
859                                 }
860                         }
861                 }
862                 final ComplexSquareMatrix us[]=new ComplexSquareMatrix[2];
863                 us[0]=new ComplexSquareMatrix(argRe,argIm);
864                 us[1]=new ComplexSquareMatrix(modRe,modIm);
865                 return us;
866         }
867
868 // MAP ELEMENTS
869
870         /**
871         * Applies a function on all the matrix elements.
872         * @param f a user-defined function
873         * @return a complex square matrix
874         */
875         public AbstractComplexMatrix mapElements(final ComplexMapping f) {
876                 final Complex array[][]=new Complex[numRows][numCols];
877                 for(int j,i=0;i<numRows;i++) {
878                         array[i][0]=f.map(matrixRe[i][0],matrixIm[i][0]);
879                         for(j=1;j<numCols;j++)
880                                 array[i][j]=f.map(matrixRe[i][j],matrixIm[i][j]);
881                 }
882                 return new ComplexSquareMatrix(array);
883         }
884 }
885
Popular Tags