|
1
2 package JSci.maths.matrices;
3
4 import JSci.GlobalSettings;
5 import JSci.maths.Complex;
6 import JSci.maths.ComplexMapping;
7 import JSci.maths.ArrayMath;
8 import JSci.maths.LinearMath;
9 import JSci.maths.MaximumIterationsExceededException;
10 import JSci.maths.vectors.AbstractComplexVector;
11 import JSci.maths.groups.AbelianGroup;
12 import JSci.maths.algebras.*;
13 import JSci.maths.fields.*;
14
15
20 public abstract class AbstractComplexSquareMatrix extends AbstractComplexMatrix implements CStarAlgebra.Member, SquareMatrix {
21 protected transient AbstractComplexSquareMatrix[] LU;
22 protected transient int[] LUpivot;
23
26 protected AbstractComplexSquareMatrix(final int size) {
27 super(size, size);
28 }
29
33 public AbstractDoubleMatrix real() {
34 final double ans[][]=new double[numRows][numCols];
35 for(int i=0;i<numRows;i++) {
36 for(int j=0;j<numCols;j++)
37 ans[i][j]=getElement(i,j).real();
38 }
39 return new DoubleSquareMatrix(ans);
40 }
41
45 public AbstractDoubleMatrix imag() {
46 final double ans[][]=new double[numRows][numCols];
47 for(int i=0;i<numRows;i++) {
48 for(int j=0;j<numCols;j++)
49 ans[i][j]=getElement(i,j).imag();
50 }
51 return new DoubleSquareMatrix(ans);
52 }
53
56 public boolean isHermitian() {
57 return this.equals(this.hermitianAdjoint());
58 }
59
62 public boolean isUnitary() {
63 return this.multiply(this.hermitianAdjoint()).equals(ComplexDiagonalMatrix.identity(numRows));
64 }
65
68 public Complex det() {
69 if(numRows==2) {
70 return new Complex(
71 getRealElement(0,0)*getRealElement(1,1) - getImagElement(0,0)*getImagElement(1,1)
72 -getRealElement(0,1)*getRealElement(1,0) + getImagElement(0,1)*getImagElement(1,0),
73 getRealElement(0,0)*getImagElement(1,1) + getImagElement(0,0)*getRealElement(1,1)
74 -getRealElement(0,1)*getImagElement(1,0) - getImagElement(0,1)*getRealElement(1,0)
75 );
76 } else {
77 final AbstractComplexSquareMatrix lu[]=this.luDecompose(null);
78 double tmp;
79 double detRe=lu[1].getRealElement(0,0);
80 double detIm=lu[1].getImagElement(0,0);
81 for(int i=1;i<numRows;i++) {
82 tmp=detRe*lu[1].getRealElement(i,i) - detIm*lu[1].getImagElement(i,i);
83 detIm=detRe*lu[1].getImagElement(i,i) + detIm*lu[1].getRealElement(i,i);
84 detRe=tmp;
85 }
86 return new Complex(detRe*LUpivot[numRows],detIm*LUpivot[numRows]);
87 }
88 }
89
92 public Complex trace() {
93 double trRe=getRealElement(0,0);
94 double trIm=getImagElement(0,0);
95 for(int i=1;i<numRows;i++) {
96 trRe+=getRealElement(i,i);
97 trIm+=getImagElement(i,i);
98 }
99 return new Complex(trRe,trIm);
100 }
101
104 public double norm() {
105 try {
106 return operatorNorm();
107 } catch(MaximumIterationsExceededException e) {
108 throw new RuntimeException  (e);
109 }
110 }
111
115 public double operatorNorm() throws MaximumIterationsExceededException {
116 return Math.sqrt(ArrayMath.max(LinearMath.eigenvalueSolveHermitian((AbstractComplexSquareMatrix)(this.hermitianAdjoint().multiply(this)))));
117 }
118
119
123
126 public AbelianGroup.Member negate() {
127 final double arrayRe[][]=new double[numRows][numCols];
128 final double arrayIm[][]=new double[numRows][numCols];
129 for(int i=0;i<numRows;i++) {
130 arrayRe[i][0]=-getRealElement(i,0);
131 arrayIm[i][0]=-getImagElement(i,0);
132 for(int j=1;j<numCols;j++) {
133 arrayRe[i][j]=-getRealElement(i,j);
134 arrayIm[i][j]=-getImagElement(i,j);
135 }
136 }
137 return new ComplexSquareMatrix(arrayRe,arrayIm);
138 }
139
140
142
147 public final AbstractComplexMatrix add(final AbstractComplexMatrix m) {
148 if(m instanceof AbstractComplexSquareMatrix)
149 return add((AbstractComplexSquareMatrix)m);
150 else if(numRows == m.rows() && numCols == m.columns())
151 return add(new SquareMatrixAdaptor(m));
152 else
153 throw new MatrixDimensionException("Matrices are different sizes.");
154 }
155
160 public AbstractComplexSquareMatrix add(final AbstractComplexSquareMatrix m) {
161 if(numRows==m.rows() && numCols==m.columns()) {
162 final double arrayRe[][]=new double[numRows][numCols];
163 final double arrayIm[][]=new double[numRows][numCols];
164 for(int i=0;i<numRows;i++) {
165 arrayRe[i][0] = getRealElement(i,0)+m.getRealElement(i,0);
166 arrayIm[i][0] = getImagElement(i,0)+m.getImagElement(i,0);
167 for(int j=1;j<numCols;j++) {
168 arrayRe[i][j] = getRealElement(i,j)+m.getRealElement(i,j);
169 arrayIm[i][j] = getImagElement(i,j)+m.getImagElement(i,j);
170 }
171 }
172 return new ComplexSquareMatrix(arrayRe,arrayIm);
173 } else
174 throw new MatrixDimensionException("Matrices are different sizes.");
175 }
176
177
179
184 public final AbstractComplexMatrix subtract(final AbstractComplexMatrix m) {
185 if(m instanceof AbstractComplexSquareMatrix)
186 return subtract((AbstractComplexSquareMatrix)m);
187 else if(numRows == m.rows() && numCols == m.columns())
188 return subtract(new SquareMatrixAdaptor(m));
189 else
190 throw new MatrixDimensionException("Matrices are different sizes.");
191 }
192
197 public AbstractComplexSquareMatrix subtract(final AbstractComplexSquareMatrix m) {
198 if(numRows==m.rows() && numCols==m.columns()) {
199 final double arrayRe[][]=new double[numRows][numCols];
200 final double arrayIm[][]=new double[numRows][numCols];
201 for(int i=0;i<numRows;i++) {
202 arrayRe[i][0] = getRealElement(i,0)-m.getRealElement(i,0);
203 arrayIm[i][0] = getImagElement(i,0)-m.getImagElement(i,0);
204 for(int j=1;j<numCols;j++) {
205 arrayRe[i][j] = getRealElement(i,j)-m.getRealElement(i,j);
206 arrayIm[i][j] = getImagElement(i,j)-m.getImagElement(i,j);
207 }
208 }
209 return new ComplexSquareMatrix(arrayRe,arrayIm);
210 } else
211 throw new MatrixDimensionException("Matrices are different sizes.");
212 }
213
214
216
221 public AbstractComplexMatrix scalarMultiply(final Complex z) {
222 final double real=z.real();
223 final double imag=z.imag();
224 final double arrayRe[][]=new double[numRows][numCols];
225 final double arrayIm[][]=new double[numRows][numCols];
226 for(int i=0;i<numRows;i++) {
227 arrayRe[i][0] = real*getRealElement(i,0)-imag*getImagElement(i,0);
228 arrayIm[i][0] = imag*getRealElement(i,0)+real*getImagElement(i,0);
229 for(int j=1;j<numCols;j++) {
230 arrayRe[i][j] = real*getRealElement(i,j)-imag*getImagElement(i,j);
231 arrayIm[i][j] = imag*getRealElement(i,j)+real*getImagElement(i,j);
232 }
233 }
234 return new ComplexSquareMatrix(arrayRe,arrayIm);
235 }
236
241 public AbstractComplexMatrix scalarMultiply(final double x) {
242 final double arrayRe[][]=new double[numRows][numCols];
243 final double arrayIm[][]=new double[numRows][numCols];
244 for(int i=0;i<numRows;i++) {
245 arrayRe[i][0] = x*getRealElement(i,0);
246 arrayIm[i][0] = x*getImagElement(i,0);
247 for(int j=1;j<numCols;j++) {
248 arrayRe[i][j] = x*getRealElement(i,j);
249 arrayIm[i][j] = x*getImagElement(i,j);
250 }
251 }
252 return new ComplexSquareMatrix(arrayRe,arrayIm);
253 }
254
255
257
262 public AbstractComplexMatrix scalarDivide(final Complex z) {
263 final Complex array[][]=new Complex[numRows][numCols];
264 for(int i=0;i<numRows;i++) {
265 array[i][0] = getElement(i,0).divide(z);
266 for(int j=1;j<numCols;j++)
267 array[i][j] = getElement(i,j).divide(z);
268 }
269 return new ComplexSquareMatrix(array);
270 }
271
276 public AbstractComplexMatrix scalarDivide(final double x) {
277 final double arrayRe[][]=new double[numRows][numCols];
278 final double arrayIm[][]=new double[numRows][numCols];
279 for(int i=0;i<numRows;i++) {
280 arrayRe[i][0]=getRealElement(i,0)/x;
281 arrayIm[i][0]=getImagElement(i,0)/x;
282 for(int j=1;j<numCols;j++) {
283 arrayRe[i][j]=getRealElement(i,j)/x;
284 arrayIm[i][j]=getImagElement(i,j)/x;
285 }
286 }
287 return new ComplexSquareMatrix(arrayRe,arrayIm);
288 }
289
290
292
297 public final Complex scalarProduct(final AbstractComplexMatrix m) {
298 if(m instanceof AbstractComplexSquareMatrix)
299 return scalarProduct((AbstractComplexSquareMatrix)m);
300 else if(numRows == m.rows() && numCols == m.columns())
301 return scalarProduct(new SquareMatrixAdaptor(m));
302 else
303 throw new MatrixDimensionException("Matrices are different sizes.");
304 }
305
310 public Complex scalarProduct(final AbstractComplexSquareMatrix m) {
311 if(numRows==m.rows() && numCols==m.columns()) {
312 double real = 0.0, imag = 0.0;
313 for(int i=0; i<numRows; i++) {
314 real += getRealElement(i,0)*m.getRealElement(i,0) + getImagElement(i,0)*m.getImagElement(i,0);
315 imag += getImagElement(i,0)*m.getRealElement(i,0) - getRealElement(i,0)*m.getImagElement(i,0);
316 for(int j=1; j<numCols; j++) {
317 real += getRealElement(i,j)*m.getRealElement(i,j) + getImagElement(i,j)*m.getImagElement(i,j);
318 imag += getImagElement(i,j)*m.getRealElement(i,j) - getRealElement(i,j)*m.getImagElement(i,j);
319 }
320 }
321 return new Complex(real, imag);
322 } else {
323 throw new MatrixDimensionException("Matrices are different sizes.");
324 }
325 }
326
327
329
335 public AbstractComplexSquareMatrix multiply(final AbstractComplexSquareMatrix m) {
336 if(numCols==m.rows()) {
337 final double arrayRe[][]=new double[numRows][m.columns()];
338 final double arrayIm[][]=new double[numRows][m.columns()];
339 Complex tmp;
340 for(int j=0;j<numRows;j++) {
341 for(int k=0;k<m.columns();k++) {
342 tmp=getElement(j,0).multiply(m.getElement(0,k));
343 arrayRe[j][k]=tmp.real();
344 arrayIm[j][k]=tmp.imag();
345 for(int n=1;n<numCols;n++) {
346 tmp=getElement(j,n).multiply(m.getElement(n,k));
347 arrayRe[j][k]+=tmp.real();
348 arrayIm[j][k]+=tmp.imag();
349 }
350 }
351 }
352 return new ComplexSquareMatrix(arrayRe,arrayIm);
353 } else {
354 throw new MatrixDimensionException("Incompatible matrices.");
355 }
356 }
357
358
360
363 public AbstractComplexSquareMatrix directSum(final AbstractComplexSquareMatrix m) {
364 final double arrayRe[][]=new double[numRows+m.numRows][numCols+m.numCols];
365 final double arrayIm[][]=new double[numRows+m.numRows][numCols+m.numCols];
366 for(int j,i=0;i<numRows;i++) {
367 for(j=0;j<numCols;j++) {
368 arrayRe[i][j]=getRealElement(i,j);
369 arrayIm[i][j]=getImagElement(i,j);
370 }
371 }
372 for(int j,i=0;i<m.numRows;i++) {
373 for(j=0;j<m.numCols;j++) {
374 arrayRe[i+numRows][j+numCols]=m.getRealElement(i,j);
375 arrayIm[i+numRows][j+numCols]=m.getImagElement(i,j);
376 }
377 }
378 return new ComplexSquareMatrix(arrayRe,arrayIm);
379 }
380
381
383
386 public AbstractComplexSquareMatrix tensor(final AbstractComplexSquareMatrix m) {
387 final double arrayRe[][]=new double[numRows*m.numRows][numCols*m.numCols];
388 final double arrayIm[][]=new double[numRows*m.numRows][numCols*m.numCols];
389 for(int i=0;i<numRows;i++) {
390 for(int j=0;j<numCols;j++) {
391 for(int k=0;k<m.numRows;j++) {
392 for(int l=0;l<m.numCols;l++) {
393 Complex tmp=getElement(i,j).multiply(m.getElement(k,l));
394 arrayRe[i*m.numRows+k][j*m.numCols+l]=tmp.real();
395 arrayIm[i*m.numRows+k][j*m.numCols+l]=tmp.imag();
396 }
397 }
398 }
399 }
400 return new ComplexSquareMatrix(arrayRe,arrayIm);
401 }
402
403
405
408 public final CStarAlgebra.Member involution() {
409 return (CStarAlgebra.Member) hermitianAdjoint();
410 }
411
415 public AbstractComplexMatrix hermitianAdjoint() {
416 final double arrayRe[][]=new double[numCols][numRows];
417 final double arrayIm[][]=new double[numCols][numRows];
418 for(int i=0;i<numRows;i++) {
419 arrayRe[0][i]=getRealElement(i,0);
420 arrayIm[0][i]=-getImagElement(i,0);
421 for(int j=1;j<numCols;j++) {
422 arrayRe[j][i]=getRealElement(i,j);
423 arrayIm[j][i]=-getImagElement(i,j);
424 }
425 }
426 return new ComplexSquareMatrix(arrayRe,arrayIm);
427 }
428
429
431
435 public AbstractComplexMatrix conjugate() {
436 final double arrayRe[][]=new double[numCols][numRows];
437 final double arrayIm[][]=new double[numCols][numRows];
438 for(int i=0;i<numRows;i++) {
439 arrayRe[i][0]=getRealElement(i,0);
440 arrayIm[i][0]=-getImagElement(i,0);
441 for(int j=1;j<numCols;j++) {
442 arrayRe[i][j]=getRealElement(i,j);
443 arrayIm[i][j]=-getImagElement(i,j);
444 }
445 }
446 return new ComplexSquareMatrix(arrayRe,arrayIm);
447 }
448
449
451
455 public Matrix transpose() {
456 final double arrayRe[][]=new double[numCols][numRows];
457 final double arrayIm[][]=new double[numCols][numRows];
458 for(int i=0;i<numRows;i++) {
459 arrayRe[0][i]=getRealElement(i,0);
460 arrayIm[0][i]=getImagElement(i,0);
461 for(int j=1;j<numCols;j++) {
462 arrayRe[j][i]=getRealElement(i,j);
463 arrayIm[j][i]=getImagElement(i,j);
464 }
465 }
466 return new ComplexSquareMatrix(arrayRe,arrayIm);
467 }
468
469
471
475 public AbstractComplexSquareMatrix inverse() {
476 final int N=numRows;
477 final double arrayLRe[][]=new double[N][N];
478 final double arrayLIm[][]=new double[N][N];
479 final double arrayURe[][]=new double[N][N];
480 final double arrayUIm[][]=new double[N][N];
481 final AbstractComplexSquareMatrix lu[]=this.luDecompose(null);
482 double denom;
483 denom = (lu[0].getRealElement(0,0)*lu[0].getRealElement(0,0) + lu[0].getImagElement(0,0)*lu[0].getImagElement(0,0));
484 arrayLRe[0][0] = lu[0].getRealElement(0,0)/denom;
485 arrayLIm[0][0] = -lu[0].getImagElement(0,0)/denom;
486 denom = (lu[1].getRealElement(0,0)*lu[1].getRealElement(0,0) + lu[1].getImagElement(0,0)*lu[1].getImagElement(0,0));
487 arrayURe[0][0] = lu[1].getRealElement(0,0)/denom;
488 arrayUIm[0][0] = -lu[1].getImagElement(0,0)/denom;
489 for(int i=1;i<N;i++) {
490 denom = (lu[0].getRealElement(i,i)*lu[0].getRealElement(i,i) + lu[0].getImagElement(i,i)*lu[0].getImagElement(i,i));
491 arrayLRe[i][i] = lu[0].getRealElement(i,i)/denom;
492 arrayLIm[i][i] = -lu[0].getImagElement(i,i)/denom;
493 denom = (lu[1].getRealElement(i,i)*lu[1].getRealElement(i,i) + lu[1].getImagElement(i,i)*lu[1].getImagElement(i,i));
494 arrayURe[i][i] = lu[1].getRealElement(i,i)/denom;
495 arrayUIm[i][i] = -lu[1].getImagElement(i,i)/denom;
496 }
497 for(int i=0;i<N-1;i++) {
498 for(int j=i+1;j<N;j++) {
499 double tmpLRe=0.0, tmpLIm=0.0;
500 double tmpURe=0.0, tmpUIm=0.0;
501 for(int k=i;k<j;k++) {
502 tmpLRe-=(lu[0].getRealElement(j,k)*arrayLRe[k][i] - lu[0].getImagElement(j,k)*arrayLIm[k][i]);
503 tmpLIm-=(lu[0].getImagElement(j,k)*arrayLRe[k][i] + lu[0].getRealElement(j,k)*arrayLIm[k][i]);
504 tmpURe-=(arrayURe[i][k]*lu[1].getRealElement(k,j) - arrayUIm[i][k]*lu[1].getImagElement(k,j));
505 tmpUIm-=(arrayUIm[i][k]*lu[1].getRealElement(k,j) + arrayURe[i][k]*lu[1].getImagElement(k,j));
506 }
507 denom = (lu[0].getRealElement(j,j)*lu[0].getRealElement(j,j) + lu[0].getImagElement(j,j)*lu[0].getImagElement(j,j));
508 arrayLRe[j][i]=(tmpLRe*lu[0].getRealElement(j,j)+tmpLIm*lu[0].getImagElement(j,j))/denom;
509 arrayLIm[j][i]=(tmpLIm*lu[0].getRealElement(j,j)-tmpLRe*lu[0].getImagElement(j,j))/denom;
510 denom = (lu[1].getRealElement(j,j)*lu[1].getRealElement(j,j) + lu[1].getImagElement(j,j)*lu[1].getImagElement(j,j));
511 arrayURe[i][j]=(tmpURe*lu[1].getRealElement(j,j)+tmpUIm*lu[1].getImagElement(j,j))/denom;
512 arrayUIm[i][j]=(tmpUIm*lu[1].getRealElement(j,j)-tmpURe*lu[1].getImagElement(j,j))/denom;
513 }
514 }
515 final double invRe[][]=new double[N][N];
517 final double invIm[][]=new double[N][N];
518 for(int i=0;i<N;i++) {
519 for(int j=0;j<i;j++) {
520 for(int k=i;k<N;k++) {
521 invRe[i][LUpivot[j]]+=(arrayURe[i][k]*arrayLRe[k][j] - arrayUIm[i][k]*arrayLIm[k][j]);
522 invIm[i][LUpivot[j]]+=(arrayUIm[i][k]*arrayLRe[k][j] + arrayURe[i][k]*arrayLIm[k][j]);
523 }
524 }
525 for(int j=i;j<N;j++) {
526 for(int k=j;k<N;k++) {
527 invRe[i][LUpivot[j]]+=(arrayURe[i][k]*arrayLRe[k][j] - arrayUIm[i][k]*arrayLIm[k][j]);
528 invIm[i][LUpivot[j]]+=(arrayUIm[i][k]*arrayLRe[k][j] + arrayURe[i][k]*arrayLIm[k][j]);
529 }
530 }
531 }
532 return new ComplexSquareMatrix(invRe,invIm);
533 }
534
535
537
541 public AbstractComplexSquareMatrix[] luDecompose(int pivot[]) {
542 if(LU!=null) {
543 if(pivot!=null)
544 System.arraycopy(LUpivot,0,pivot,0,pivot.length);
545 return LU;
546 }
547 final int N=numRows;
548 final double arrayLRe[][]=new double[N][N];
549 final double arrayLIm[][]=new double[N][N];
550 final double arrayURe[][]=new double[N][N];
551 final double arrayUIm[][]=new double[N][N];
552 if(pivot==null)
553 pivot=new int[N+1];
554 for(int i=0;i<N;i++)
555 pivot[i]=i;
556 pivot[N]=1;
557 for(int j=0;j<N;j++) {
559 for(int i=0;i<j;i++) {
560 double tmpRe = getRealElement(pivot[i],j);
561 double tmpIm = getImagElement(pivot[i],j);
562 for(int k=0;k<i;k++) {
563 tmpRe-=(arrayURe[i][k]*arrayURe[k][j] - arrayUIm[i][k]*arrayUIm[k][j]);
564 tmpIm-=(arrayUIm[i][k]*arrayURe[k][j] + arrayURe[i][k]*arrayUIm[k][j]);
565 }
566 arrayURe[i][j]=tmpRe;
567 arrayUIm[i][j]=tmpIm;
568 }
569 double max=0.0;
570 int pivotrow=j;
571 for(int i=j;i<N;i++) {
572 double tmpRe = getRealElement(pivot[i],j);
573 double tmpIm = getImagElement(pivot[i],j);
574 for(int k=0;k<j;k++) {
575 tmpRe-=(arrayURe[i][k]*arrayURe[k][j] - arrayUIm[i][k]*arrayUIm[k][j]);
576 tmpIm-=(arrayUIm[i][k]*arrayURe[k][j] + arrayURe[i][k]*arrayUIm[k][j]);
577 }
578 arrayURe[i][j]=tmpRe;
579 arrayUIm[i][j]=tmpIm;
580 double tmp=tmpRe*tmpRe+tmpIm*tmpIm;
582 if(tmp>max) {
583 max=tmp;
584 pivotrow=i;
585 }
586 }
587 if(pivotrow!=j) {
589 double[] tmprow = arrayURe[j];
590 arrayURe[j] = arrayURe[pivotrow];
591 arrayURe[pivotrow] = tmprow;
592 tmprow = arrayUIm[j];
593 arrayUIm[j] = arrayUIm[pivotrow];
594 arrayUIm[pivotrow] = tmprow;
595 int k=pivot[j];
596 pivot[j]=pivot[pivotrow];
597 pivot[pivotrow]=k;
598 pivot[N]=-pivot[N];
600 }
601 double tmpRe=arrayURe[j][j];
603 double tmpIm=arrayUIm[j][j];
604 if(Math.abs(tmpRe)<Math.abs(tmpIm)) {
605 double a=tmpRe/tmpIm;
606 double denom=tmpRe*a+tmpIm;
607 for(int i=j+1;i<N;i++) {
608 double tmp=(arrayURe[i][j]*a+arrayUIm[i][j])/denom;
609 arrayUIm[i][j]=(arrayUIm[i][j]*a-arrayURe[i][j])/denom;
610 arrayURe[i][j]=tmp;
611 }
612 } else {
613 double a=tmpIm/tmpRe;
614 double denom=tmpRe+tmpIm*a;
615 for(int i=j+1;i<N;i++) {
616 double tmp=(arrayURe[i][j]+arrayUIm[i][j]*a)/denom;
617 arrayUIm[i][j]=(arrayUIm[i][j]-arrayURe[i][j]*a)/denom;
618 arrayURe[i][j]=tmp;
619 }
620 }
621 }
622 for(int j=0;j<N;j++) {
624 arrayLRe[j][j]=1.0;
625 for(int i=j+1;i<N;i++) {
626 arrayLRe[i][j]=arrayURe[i][j];
627 arrayLIm[i][j]=arrayUIm[i][j];
628 arrayURe[i][j]=0.0;
629 arrayUIm[i][j]=0.0;
630 }
631 }
632 LU=new AbstractComplexSquareMatrix[2];
633 LU[0]=new ComplexSquareMatrix(arrayLRe,arrayLIm);
634 LU[1]=new ComplexSquareMatrix(arrayURe,arrayUIm);
635 LUpivot=new int[pivot.length];
636 System.arraycopy(pivot,0,LUpivot,0,pivot.length);
637 return LU;
638 }
639
646 public AbstractComplexSquareMatrix[] luDecompose() {
647 final int N=numRows;
648 final double arrayLRe[][]=new double[N][N];
649 final double arrayLIm[][]=new double[N][N];
650 final double arrayURe[][]=new double[N][N];
651 final double arrayUIm[][]=new double[N][N];
652 for(int j=0;j<N;j++) {
654 for(int i=0;i<j;i++) {
655 double tmpRe = getRealElement(i,j);
656 double tmpIm = getImagElement(i,j);
657 for(int k=0;k<i;k++) {
658 tmpRe-=(arrayURe[i][k]*arrayURe[k][j] - arrayUIm[i][k]*arrayUIm[k][j]);
659 tmpIm-=(arrayUIm[i][k]*arrayURe[k][j] + arrayURe[i][k]*arrayUIm[k][j]);
660 }
661 arrayURe[i][j]=tmpRe;
662 arrayUIm[i][j]=tmpIm;
663 }
664 for(int i=j;i<N;i++) {
665 double tmpRe = getRealElement(i,j);
666 double tmpIm = getImagElement(i,j);
667 for(int k=0;k<j;k++) {
668 tmpRe-=(arrayURe[i][k]*arrayURe[k][j] - arrayUIm[i][k]*arrayUIm[k][j]);
669 tmpIm-=(arrayUIm[i][k]*arrayURe[k][j] + arrayURe[i][k]*arrayUIm[k][j]);
670 }
671 arrayURe[i][j]=tmpRe;
672 arrayUIm[i][j]=tmpIm;
673 }
674 double tmpRe=arrayURe[j][j];
676 double tmpIm=arrayUIm[j][j];
677 if(Math.abs(tmpRe)<Math.abs(tmpIm)) {
678 double a=tmpRe/tmpIm;
679 double denom=tmpRe*a+tmpIm;
680 for(int i=j+1;i<N;i++) {
681 double tmp=(arrayURe[i][j]*a+arrayUIm[i][j])/denom;
682 arrayUIm[i][j]=(arrayUIm[i][j]*a-arrayURe[i][j])/denom;
683 arrayURe[i][j]=tmp;
684 }
685 } else {
686 double a=tmpIm/tmpRe;
687 double denom=tmpRe+tmpIm*a;
688 for(int i=j+1;i<N;i++) {
689 double tmp=(arrayURe[i][j]+arrayUIm[i][j]*a)/denom;
690 arrayUIm[i][j]=(arrayUIm[i][j]-arrayURe[i][j]*a)/denom;
691 arrayURe[i][j]=tmp;
692 }
693 }
694 }
695 for(int j=0;j<N;j++) {
697 arrayLRe[j][j]=1.0;
698 for(int i=j+1;i<N;i++) {
699 arrayLRe[i][j]=arrayURe[i][j];
700 arrayLIm[i][j]=arrayUIm[i][j];
701 arrayURe[i][j]=0.0;
702 arrayUIm[i][j]=0.0;
703 }
704 }
705 AbstractComplexSquareMatrix[] lu=new AbstractComplexSquareMatrix[2];
706 lu[0]=new ComplexSquareMatrix(arrayLRe,arrayLIm);
707 lu[1]=new ComplexSquareMatrix(arrayURe,arrayUIm);
708 return lu;
709 }
710
711
713
716 public AbstractComplexSquareMatrix[] polarDecompose() {
717 final int N=numRows;
718 final AbstractComplexVector evec[]=new AbstractComplexVector[N];
719 double eval[];
720 try {
721 eval=LinearMath.eigenSolveHermitian(this,evec);
722 } catch(MaximumIterationsExceededException e) {
723 return null;
724 }
725 final double tmpaRe[][]=new double[N][N];
726 final double tmpaIm[][]=new double[N][N];
727 final double tmpmRe[][]=new double[N][N];
728 final double tmpmIm[][]=new double[N][N];
729 double abs;
730 Complex comp;
731 for(int i=0;i<N;i++) {
732 abs=Math.abs(eval[i]);
733 comp=evec[i].getComponent(0).conjugate();
734 tmpaRe[i][0]=eval[i]*comp.real()/abs;
735 tmpaIm[i][0]=eval[i]*comp.imag()/abs;
736 tmpmRe[i][0]=abs*comp.real();
737 tmpmIm[i][0]=abs*comp.imag();
738 for(int j=1;j<N;j++) {
739 comp=evec[i].getComponent(j).conjugate();
740 tmpaRe[i][j]=eval[i]*comp.real()/abs;
741 tmpaIm[i][j]=eval[i]*comp.imag()/abs;
742 tmpmRe[i][j]=abs*comp.real();
743 tmpmIm[i][j]=abs*comp.imag();
744 }
745 }
746 final double argRe[][]=new double[N][N];
747 final double argIm[][]=new double[N][N];
748 final double modRe[][]=new double[N][N];
749 final double modIm[][]=new double[N][N];
750 for(int i=0;i<N;i++) {
751 for(int j=0;j<N;j++) {
752 comp=evec[0].getComponent(i);
753 argRe[i][j]=(tmpaRe[0][j]*comp.real() - tmpaIm[0][j]*comp.imag());
754 argIm[i][j]=(tmpaIm[0][j]*comp.real() + tmpaRe[0][j]*comp.imag());
755 modRe[i][j]=(tmpmRe[0][j]*comp.real() - tmpmIm[0][j]*comp.imag());
756 modIm[i][j]=(tmpmIm[0][j]*comp.real() + tmpmRe[0][j]*comp.imag());
757 for(int k=1;k<N;k++) {
758 comp=evec[k].getComponent(i);
759 argRe[i][j]+=(tmpaRe[k][j]*comp.real() - tmpaIm[k][j]*comp.imag());
760 argIm[i][j]+=(tmpaIm[k][j]*comp.real() + tmpaRe[k][j]*comp.imag());
761 modRe[i][j]+=(tmpmRe[k][j]*comp.real() - tmpmIm[k][j]*comp.imag());
762 modIm[i][j]+=(tmpmIm[k][j]*comp.real() + tmpmRe[k][j]*comp.imag());
763 }
764 }
765 }
766 final AbstractComplexSquareMatrix us[]=new AbstractComplexSquareMatrix[2];
767 us[0]=new ComplexSquareMatrix(argRe,argIm);
768 us[1]=new ComplexSquareMatrix(modRe,modIm);
769 return us;
770 }
771
772
774
779 public AbstractComplexMatrix mapElements(final ComplexMapping f) {
780 final Complex array[][]=new Complex[numRows][numCols];
781 for(int i=0;i<numRows;i++) {
782 array[i][0]=f.map(getElement(i,0));
783 for(int j=1;j<numCols;j++)
784 array[i][j]=f.map(getElement(i,j));
785 }
786 return new ComplexSquareMatrix(array);
787 }
788
789 private static class SquareMatrixAdaptor extends AbstractComplexSquareMatrix {
790 private final AbstractComplexMatrix matrix;
791 private SquareMatrixAdaptor(AbstractComplexMatrix m) {
792 super(m.rows());
793 matrix = m;
794 }
795 public Complex getElement(int i, int j) {
796 return matrix.getElement(i,j);
797 }
798 public double getRealElement(int i, int j) {
799 return matrix.getRealElement(i,j);
800 }
801 public double getImagElement(int i, int j) {
802 return matrix.getImagElement(i,j);
803 }
804 public void setElement(int i, int j, Complex z) {
805 matrix.setElement(i,j, z);
806 }
807 public void setElement(int i, int j, double x, double y) {
808 matrix.setElement(i,j, x,y);
809 }
810 }
811 }
812 |
Popular Tags |
Java API By Example, From Geeks To Geeks.