1 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 19 public class ComplexSquareMatrix extends AbstractComplexSquareMatrix { 20 23 protected final double matrixRe[][],matrixIm[][]; 24 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 43 public ComplexSquareMatrix(final int size) { 44 this(new double[size][size], new double[size][size]); 45 } 46 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 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 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 99 public String toString() { 100 final StringBuffer buf=new StringBuffer (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 114 public AbstractDoubleMatrix real() { 115 return new DoubleSquareMatrix(matrixRe); 116 } 117 121 public AbstractDoubleMatrix imag() { 122 return new DoubleSquareMatrix(matrixIm); 123 } 124 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 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 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 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 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 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 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 245 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 264 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 307 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 350 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 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 391 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 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 465 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 489 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 511 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 531 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 547 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 567 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 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 633 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 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 double tmp=tmpRe*tmpRe+tmpIm*tmpIm; 678 if(tmp>max) { 679 max=tmp; 680 pivotrow=i; 681 } 682 } 683 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 pivot[N]=-pivot[N]; 696 } 697 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 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 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 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 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 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 809 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 870 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 |