1 2 package JSci.maths.matrices; 3 4 import JSci.GlobalSettings; 5 import JSci.maths.ArrayMath; 6 import JSci.maths.ExtraMath; 7 import JSci.maths.LinearMath; 8 import JSci.maths.MaximumIterationsExceededException; 9 import JSci.maths.Mapping; 10 import JSci.maths.vectors.AbstractIntegerVector; 11 import JSci.maths.groups.AbelianGroup; 12 import JSci.maths.algebras.*; 13 import JSci.maths.fields.*; 14 15 20 public abstract class AbstractIntegerSquareMatrix extends AbstractIntegerMatrix implements SquareMatrix { 21 protected transient AbstractDoubleSquareMatrix[] LU; 22 protected transient int[] LUpivot; 23 26 protected AbstractIntegerSquareMatrix(final int size) { 27 super(size, size); 28 } 29 33 public AbstractDoubleMatrix toDoubleMatrix() { 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); 38 } 39 return new DoubleSquareMatrix(ans); 40 } 41 45 public AbstractComplexMatrix toComplexMatrix() { 46 ComplexSquareMatrix cm = new ComplexSquareMatrix(numRows); 47 for(int i=0;i<numRows;i++) { 48 for(int j=0;j<numCols;j++) 49 cm.setElement(i, j, getElement(i, j), 0.0); 50 } 51 return cm; 52 } 53 56 public boolean isSymmetric() { 57 return this.equals(this.transpose()); 58 } 59 62 public boolean isUnitary() { 63 return this.multiply(this.transpose()).equals(DoubleDiagonalMatrix.identity(numRows)); 64 } 65 68 public int det() { 69 if(numRows==2) { 70 return getElement(0,0)*getElement(1,1)-getElement(0,1)*getElement(1,0); 71 } else { 72 final AbstractDoubleSquareMatrix lu[]=this.luDecompose(null); 73 double det=lu[1].getElement(0,0); 74 for(int i=1;i<numRows;i++) 75 det*=lu[1].getElement(i,i); 76 return Math.round((float)det)*LUpivot[numRows]; 77 } 78 } 79 82 public int trace() { 83 int result = getElement(0,0); 84 for(int i=1;i<numRows;i++) 85 result += getElement(i,i); 86 return result; 87 } 88 89 93 96 public AbelianGroup.Member negate() { 97 final int array[][]=new int[numRows][numCols]; 98 for(int i=0;i<numRows;i++) { 99 array[i][0] = -getElement(i,0); 100 for(int j=1;j<numCols;j++) 101 array[i][j] = -getElement(i,j); 102 } 103 return new IntegerSquareMatrix(array); 104 } 105 106 108 113 public final AbstractIntegerMatrix add(final AbstractIntegerMatrix m) { 114 if(m instanceof AbstractIntegerSquareMatrix) 115 return add((AbstractIntegerSquareMatrix)m); 116 else if(numRows == m.rows() && numCols == m.columns()) 117 return add(new SquareMatrixAdaptor(m)); 118 else 119 throw new MatrixDimensionException("Matrices are different sizes."); 120 } 121 126 public AbstractIntegerSquareMatrix add(final AbstractIntegerSquareMatrix m) { 127 if(numRows==m.rows() && numCols==m.columns()) { 128 final int array[][]=new int[numRows][numCols]; 129 for(int i=0;i<numRows;i++) { 130 array[i][0] = getElement(i,0)+m.getElement(i,0); 131 for(int j=1;j<numCols;j++) 132 array[i][j] = getElement(i,j)+m.getElement(i,j); 133 } 134 return new IntegerSquareMatrix(array); 135 } else { 136 throw new MatrixDimensionException("Matrices are different sizes."); 137 } 138 } 139 140 142 147 public final AbstractIntegerMatrix subtract(final AbstractIntegerMatrix m) { 148 if(m instanceof AbstractIntegerSquareMatrix) 149 return subtract((AbstractIntegerSquareMatrix)m); 150 else if(numRows == m.rows() && numCols == m.columns()) 151 return subtract(new SquareMatrixAdaptor(m)); 152 else 153 throw new MatrixDimensionException("Matrices are different sizes."); 154 } 155 160 public AbstractIntegerSquareMatrix subtract(final AbstractIntegerSquareMatrix m) { 161 if(numRows==m.rows() && numCols==m.columns()) { 162 final int array[][]=new int[numRows][numCols]; 163 for(int i=0;i<numRows;i++) { 164 array[i][0] = getElement(i,0)-m.getElement(i,0); 165 for(int j=1;j<numCols;j++) 166 array[i][j] = getElement(i,j)-m.getElement(i,j); 167 } 168 return new IntegerSquareMatrix(array); 169 } else { 170 throw new MatrixDimensionException("Matrices are different sizes."); 171 } 172 } 173 174 176 181 public AbstractIntegerMatrix scalarMultiply(final int x) { 182 final int array[][]=new int[numRows][numCols]; 183 for(int i=0;i<numRows;i++) { 184 array[i][0] = x*getElement(i,0); 185 for(int j=1;j<numCols;j++) 186 array[i][j] = x*getElement(i,j); 187 } 188 return new IntegerSquareMatrix(array); 189 } 190 191 193 194 196 201 public final int scalarProduct(final AbstractIntegerMatrix m) { 202 if(m instanceof AbstractIntegerSquareMatrix) 203 return scalarProduct((AbstractIntegerSquareMatrix)m); 204 else if(numRows == m.rows() && numCols == m.columns()) 205 return scalarProduct(new SquareMatrixAdaptor(m)); 206 else 207 throw new MatrixDimensionException("Matrices are different sizes."); 208 } 209 214 public int scalarProduct(final AbstractIntegerSquareMatrix m) { 215 if(numRows==m.rows() && numCols==m.columns()) { 216 int ans = 0; 217 for(int i=0; i<numRows; i++) { 218 ans += getElement(i,0)*m.getElement(i,0); 219 for(int j=1; j<numCols; j++) 220 ans += getElement(i,j)*m.getElement(i,j); 221 } 222 return ans; 223 } else { 224 throw new MatrixDimensionException("Matrices are different sizes."); 225 } 226 } 227 228 230 235 public AbstractIntegerSquareMatrix multiply(final AbstractIntegerSquareMatrix m) { 236 if(numCols==m.rows()) { 237 final int mColumns = m.columns(); 238 final int array[][]=new int[numRows][mColumns]; 239 for(int j=0; j<numRows; j++) { 240 for(int k=0; k<mColumns; k++) { 241 array[j][k] = getElement(j,0)*m.getElement(0,k); 242 for(int n=1; n<numCols; n++) 243 array[j][k] += getElement(j,n)*m.getElement(n,k); 244 } 245 } 246 return new IntegerSquareMatrix(array); 247 } else { 248 throw new MatrixDimensionException("Incompatible matrices."); 249 } 250 } 251 252 254 257 public AbstractIntegerSquareMatrix directSum(final AbstractIntegerSquareMatrix m) { 258 final int array[][]=new int[numRows+m.numRows][numCols+m.numCols]; 259 for(int i=0;i<numRows;i++) { 260 for(int j=0;j<numCols;j++) 261 array[i][j] = getElement(i,j); 262 } 263 for(int i=0;i<m.numRows;i++) { 264 for(int j=0;j<m.numCols;j++) 265 array[i+numRows][j+numCols] = m.getElement(i,j); 266 } 267 return new IntegerSquareMatrix(array); 268 } 269 270 272 275 public AbstractIntegerSquareMatrix tensor(final AbstractIntegerSquareMatrix m) { 276 final int array[][]=new int[numRows*m.numRows][numCols*m.numCols]; 277 for(int i=0;i<numRows;i++) { 278 for(int j=0;j<numCols;j++) { 279 for(int k=0;k<m.numRows;j++) { 280 for(int l=0;l<m.numCols;l++) 281 array[i*m.numRows+k][j*m.numCols+l] = getElement(i,j)*m.getElement(k,l); 282 } 283 } 284 } 285 return new IntegerSquareMatrix(array); 286 } 287 288 290 294 public Matrix transpose() { 295 final int array[][]=new int[numCols][numRows]; 296 for(int i=0;i<numRows;i++) { 297 array[0][i] = getElement(i,0); 298 for(int j=1;j<numCols;j++) 299 array[j][i] = getElement(i,j); 300 } 301 return new IntegerSquareMatrix(array); 302 } 303 304 306 310 public AbstractDoubleSquareMatrix inverse() { 311 final int N=numRows; 312 final double arrayL[][]=new double[N][N]; 313 final double arrayU[][]=new double[N][N]; 314 final AbstractDoubleSquareMatrix lu[]=this.luDecompose(null); 315 arrayL[0][0]=1.0/lu[0].getElement(0,0); 316 arrayU[0][0]=1.0/lu[1].getElement(0,0); 317 for(int i=1;i<N;i++) { 318 arrayL[i][i]=1.0/lu[0].getElement(i,i); 319 arrayU[i][i]=1.0/lu[1].getElement(i,i); 320 } 321 for(int i=0;i<N-1;i++) { 322 for(int j=i+1;j<N;j++) { 323 double tmpL=0.0, tmpU=0.0; 324 for(int k=i;k<j;k++) { 325 tmpL-=lu[0].getElement(j,k)*arrayL[k][i]; 326 tmpU-=arrayU[i][k]*lu[1].getElement(k,j); 327 } 328 arrayL[j][i]=tmpL/lu[0].getElement(j,j); 329 arrayU[i][j]=tmpU/lu[1].getElement(j,j); 330 } 331 } 332 final double inv[][]=new double[N][N]; 334 for(int i=0;i<N;i++) { 335 for(int j=0;j<i;j++) { 336 for(int k=i;k<N;k++) 337 inv[i][LUpivot[j]]+=arrayU[i][k]*arrayL[k][j]; 338 } 339 for(int j=i;j<N;j++) { 340 for(int k=j;k<N;k++) 341 inv[i][LUpivot[j]]+=arrayU[i][k]*arrayL[k][j]; 342 } 343 } 344 return new DoubleSquareMatrix(inv); 345 } 346 347 349 358 public AbstractDoubleSquareMatrix[] luDecompose(int pivot[]) { 359 if(LU!=null) { 360 if(pivot!=null) 361 System.arraycopy(LUpivot,0,pivot,0,pivot.length); 362 return LU; 363 } 364 int pivotrow; 365 final int N=numRows; 366 final double arrayL[][]=new double[N][N]; 367 final double arrayU[][]=new double[N][N]; 368 if(pivot==null) 369 pivot=new int[N+1]; 370 for(int i=0;i<N;i++) 371 pivot[i]=i; 372 pivot[N]=1; 373 for(int j=0;j<N;j++) { 375 for(int i=0;i<j;i++) { 376 double tmp=getElement(pivot[i],j); 377 for(int k=0;k<i;k++) 378 tmp-=arrayU[i][k]*arrayU[k][j]; 379 arrayU[i][j]=tmp; 380 } 381 double max=0.0; 382 pivotrow=j; 383 for(int i=j;i<N;i++) { 384 double tmp=getElement(pivot[i],j); 385 for(int k=0;k<j;k++) 386 tmp-=arrayU[i][k]*arrayU[k][j]; 387 arrayU[i][j]=tmp; 388 tmp=Math.abs(tmp); 390 if(tmp>max) { 391 max=tmp; 392 pivotrow=i; 393 } 394 } 395 if(pivotrow!=j) { 397 double[] tmprow = arrayU[j]; 398 arrayU[j] = arrayU[pivotrow]; 399 arrayU[pivotrow] = tmprow; 400 int k=pivot[j]; 401 pivot[j]=pivot[pivotrow]; 402 pivot[pivotrow]=k; 403 pivot[N]=-pivot[N]; 405 } 406 double tmp=arrayU[j][j]; 408 for(int i=j+1;i<N;i++) 409 arrayU[i][j]/=tmp; 410 } 411 for(int j=0;j<N;j++) { 413 arrayL[j][j]=1.0; 414 for(int i=j+1;i<N;i++) { 415 arrayL[i][j]=arrayU[i][j]; 416 arrayU[i][j]=0.0; 417 } 418 } 419 LU=new AbstractDoubleSquareMatrix[2]; 420 LU[0]=new DoubleSquareMatrix(arrayL); 421 LU[1]=new DoubleSquareMatrix(arrayU); 422 LUpivot=new int[pivot.length]; 423 System.arraycopy(pivot,0,LUpivot,0,pivot.length); 424 return LU; 425 } 426 433 public AbstractDoubleSquareMatrix[] luDecompose() { 434 final int N=numRows; 435 final double arrayL[][]=new double[N][N]; 436 final double arrayU[][]=new double[N][N]; 437 for(int j=0;j<N;j++) { 439 for(int i=0;i<j;i++) { 440 double tmp=getElement(i,j); 441 for(int k=0;k<i;k++) 442 tmp-=arrayU[i][k]*arrayU[k][j]; 443 arrayU[i][j]=tmp; 444 } 445 for(int i=j;i<N;i++) { 446 double tmp=getElement(i,j); 447 for(int k=0;k<j;k++) 448 tmp-=arrayU[i][k]*arrayU[k][j]; 449 arrayU[i][j]=tmp; 450 } 451 double tmp=arrayU[j][j]; 453 for(int i=j+1;i<N;i++) 454 arrayU[i][j]/=tmp; 455 } 456 for(int j=0;j<N;j++) { 458 arrayL[j][j]=1.0; 459 for(int i=j+1;i<N;i++) { 460 arrayL[i][j]=arrayU[i][j]; 461 arrayU[i][j]=0.0; 462 } 463 } 464 AbstractDoubleSquareMatrix[] lu=new AbstractDoubleSquareMatrix[2]; 465 lu[0]=new DoubleSquareMatrix(arrayL); 466 lu[1]=new DoubleSquareMatrix(arrayU); 467 return lu; 468 } 469 470 472 478 public AbstractDoubleSquareMatrix[] choleskyDecompose() { 479 final int N=numRows; 480 final double arrayL[][]=new double[N][N]; 481 final double arrayU[][]=new double[N][N]; 482 double tmp=Math.sqrt(getElement(0,0)); 483 arrayL[0][0]=arrayU[0][0]=tmp; 484 for(int i=1;i<N;i++) 485 arrayL[i][0]=arrayU[0][i]=getElement(i,0)/tmp; 486 for(int j=1;j<N;j++) { 487 tmp=getElement(j,j); 488 for(int i=0;i<j;i++) 489 tmp-=arrayL[j][i]*arrayL[j][i]; 490 arrayL[j][j]=arrayU[j][j]=Math.sqrt(tmp); 491 for(int i=j+1;i<N;i++) { 492 tmp=getElement(i,j); 493 for(int k=0;k<i;k++) 494 tmp-=arrayL[j][k]*arrayU[k][i]; 495 arrayL[i][j]=arrayU[j][i]=tmp/arrayU[j][j]; 496 } 497 } 498 final AbstractDoubleSquareMatrix lu[]=new AbstractDoubleSquareMatrix[2]; 499 lu[0]=new DoubleSquareMatrix(arrayL); 500 lu[1]=new DoubleSquareMatrix(arrayU); 501 return lu; 502 } 503 504 506 512 public AbstractDoubleSquareMatrix[] qrDecompose() { 513 final int N=numRows; 514 final double array[][]=new double[N][N]; 515 final double arrayQ[][]=new double[N][N]; 516 final double arrayR[][]=new double[N][N]; 517 for(int i=0;i<N;i++) { 519 array[i][0]=getElement(i,0); 520 for(int j=1;j<N;j++) 521 array[i][j]=getElement(i,j); 522 } 523 for(int k=0; k<N; k++) { 524 double norm = array[k][k]; 526 for(int i=k+1; i<N; i++) 527 norm = ExtraMath.hypot(norm, array[i][k]); 528 if(norm != 0.0) { 529 if(array[k][k] < 0.0) 531 norm = -norm; 532 for(int i=k; i<N; i++) 533 array[i][k] /= norm; 534 array[k][k] += 1.0; 535 for(int j=k+1; j<N; j++) { 537 double s=array[k][k]*array[k][j]; 538 for(int i=k+1; i<N; i++) 539 s += array[i][k]*array[i][j]; 540 s /= array[k][k]; 541 for(int i=k; i<N; i++) 542 array[i][j] -= s*array[i][k]; 543 } 544 } 545 arrayR[k][k] = -norm; 546 } 547 for(int k=N-1; k>=0; k--) { 548 arrayQ[k][k] = 1.0; 549 for(int j=k; j<N; j++) { 550 if(array[k][k] != 0.0) { 551 double s = array[k][k]*arrayQ[k][j]; 552 for(int i=k+1; i<N; i++) 553 s += array[i][k]*arrayQ[i][j]; 554 s /= array[k][k]; 555 for(int i=k; i<N; i++) 556 arrayQ[i][j] -= s*array[i][k]; 557 } 558 } 559 } 560 for(int i=0; i<N; i++) { 561 for(int j=i+1; j<N; j++) 562 arrayR[i][j] = array[i][j]; 563 } 564 final AbstractDoubleSquareMatrix qr[]=new AbstractDoubleSquareMatrix[2]; 565 qr[0]=new DoubleSquareMatrix(arrayQ); 566 qr[1]=new DoubleSquareMatrix(arrayR); 567 return qr; 568 } 569 570 572 578 public AbstractDoubleSquareMatrix[] singularValueDecompose() { 579 final int N=numRows; 580 final int Nm1=N-1; 581 final double array[][]=new double[N][N]; 582 final double arrayU[][]=new double[N][N]; 583 final double arrayS[]=new double[N]; 584 final double arrayV[][]=new double[N][N]; 585 final double e[]=new double[N]; 586 final double work[]=new double[N]; 587 for(int i=0;i<N;i++) { 589 array[i][0]=getElement(i,0); 590 for(int j=1;j<N;j++) 591 array[i][j]=getElement(i,j); 592 } 593 for(int k=0;k<Nm1;k++) { 595 arrayS[k]=array[k][k]; 598 for(int i=k+1;i<N;i++) 599 arrayS[k]=ExtraMath.hypot(arrayS[k],array[i][k]); 600 if(arrayS[k]!=0.0) { 601 if(array[k][k]<0.0) 602 arrayS[k]=-arrayS[k]; 603 for(int i=k;i<N;i++) 604 array[i][k]/=arrayS[k]; 605 array[k][k]+=1.0; 606 } 607 arrayS[k]=-arrayS[k]; 608 for(int j=k+1;j<N;j++) { 609 if(arrayS[k]!=0.0) { 610 double t=array[k][k]*array[k][j]; 612 for(int i=k+1; i<N; i++) 613 t+=array[i][k]*array[i][j]; 614 t /= array[k][k]; 615 for(int i=k; i<N; i++) 616 array[i][j] -= t*array[i][k]; 617 } 618 e[j]=array[k][j]; 619 } 620 for(int i=k;i<N;i++) 621 arrayU[i][k]=array[i][k]; 622 if(k<N-2) { 623 e[k]=e[k+1]; 626 for(int i=k+2;i<N;i++) 627 e[k]=ExtraMath.hypot(e[k],e[i]); 628 if(e[k]!=0.0) { 629 if(e[k+1]<0.0) 630 e[k]=-e[k]; 631 for(int i=k+1;i<N;i++) 632 e[i]/=e[k]; 633 e[k+1]+=1.0; 634 } 635 e[k]=-e[k]; 636 if(e[k]!=0.0) { 637 for(int i=k+1;i<N;i++) { 639 work[i]=0.0; 640 for(int j=k+1;j<N;j++) 641 work[i]+=e[j]*array[i][j]; 642 } 643 for(int j=k+1;j<N;j++) { 644 double t = e[j]/e[k+1]; 645 for(int i=k+1;i<N;i++) 646 array[i][j] -= t*work[i]; 647 } 648 } 649 for(int i=k+1;i<N;i++) 650 arrayV[i][k]=e[i]; 651 } 652 } 653 int p=N; 655 arrayS[Nm1]=array[Nm1][Nm1]; 656 e[N-2]=array[N-2][Nm1]; 657 e[Nm1]=0.0; 658 arrayU[Nm1][Nm1]=1.0; 659 for(int k=N-2;k>=0;k--) { 660 if(arrayS[k]!=0.0) { 661 for(int j=k+1;j<N;j++) { 662 double t=arrayU[k][k]*arrayU[k][j]; 663 for(int i=k+1;i<N;i++) 664 t+=arrayU[i][k]*arrayU[i][j]; 665 t /= arrayU[k][k]; 666 for(int i=k;i<N;i++) 667 arrayU[i][j] -= t*arrayU[i][k]; 668 } 669 for(int i=k;i<N;i++) 670 arrayU[i][k]=-arrayU[i][k]; 671 arrayU[k][k]+=1.0; 672 for(int i=0;i<k-1;i++) 673 arrayU[i][k]=0.0; 674 } else { 675 for(int i=0;i<N;i++) 676 arrayU[i][k]=0.0; 677 arrayU[k][k]=1.0; 678 } 679 } 680 for(int k=Nm1;k>=0;k--) { 681 if(k<N-2 && e[k]!=0.0) { 682 for(int j=k+1;j<N;j++) { 683 double t=arrayV[k+1][k]*arrayV[k+1][j]; 684 for(int i=k+2;i<N;i++) 685 t+=arrayV[i][k]*arrayV[i][j]; 686 t /= arrayV[k+1][k]; 687 for(int i=k+1;i<N;i++) 688 arrayV[i][j] -= t*arrayV[i][k]; 689 } 690 } 691 for(int i=0;i<N;i++) 692 arrayV[i][k]=0.0; 693 arrayV[k][k]=1.0; 694 } 695 final double eps=Math.pow(2.0,-52.0); 696 int iter=0; 697 while(p>0) { 698 int k,action; 699 for(k=p-2;k>=-1;k--) { 704 if(k==-1) 705 break; 706 if(Math.abs(e[k])<=eps*(Math.abs(arrayS[k])+Math.abs(arrayS[k+1]))) { 707 e[k]=0.0; 708 break; 709 } 710 } 711 if(k==p-2) { 712 action=4; 713 } else { 714 int ks; 715 for(ks=p-1;ks>=k;ks--) { 716 if(ks==k) 717 break; 718 double t=(ks!=p ? Math.abs(e[ks]) : 0.0)+(ks!=k+1 ? Math.abs(e[ks-1]) : 0.0); 719 if(Math.abs(arrayS[ks])<=eps*t) { 720 arrayS[ks]=0.0; 721 break; 722 } 723 } 724 if(ks==k) { 725 action=3; 726 } else if(ks==p-1) { 727 action=1; 728 } else { 729 action=2; 730 k=ks; 731 } 732 } 733 k++; 734 switch(action) { 735 case 1: { 737 double f=e[p-2]; 738 e[p-2]=0.0; 739 for(int j=p-2;j>=k;j--) { 740 double t=ExtraMath.hypot(arrayS[j],f); 741 final double cs=arrayS[j]/t; 742 final double sn=f/t; 743 arrayS[j]=t; 744 if(j!=k) { 745 f=-sn*e[j-1]; 746 e[j-1]*=cs; 747 } 748 for(int i=0;i<N;i++) { 749 t=cs*arrayV[i][j]+sn*arrayV[i][p-1]; 750 arrayV[i][p-1]=-sn*arrayV[i][j]+cs*arrayV[i][p-1]; 751 arrayV[i][j]=t; 752 } 753 } 754 } break; 755 case 2: { 757 double f=e[k-1]; 758 e[k-1]=0.0; 759 for(int j=k;j<p;j++) { 760 double t=ExtraMath.hypot(arrayS[j],f); 761 final double cs=arrayS[j]/t; 762 final double sn=f/t; 763 arrayS[j]=t; 764 f=-sn*e[j]; 765 e[j]*=cs; 766 for(int i=0;i<N;i++) { 767 t=cs*arrayU[i][j]+sn*arrayU[i][k-1]; 768 arrayU[i][k-1]=-sn*arrayU[i][j]+cs*arrayU[i][k-1]; 769 arrayU[i][j]=t; 770 } 771 } 772 } break; 773 case 3: { 775 final double scale=Math.max(Math.max(Math.max(Math.max( 777 Math.abs(arrayS[p-1]),Math.abs(arrayS[p-2])),Math.abs(e[p-2])), 778 Math.abs(arrayS[k])),Math.abs(e[k])); 779 double sp=arrayS[p-1]/scale; 780 double spm1=arrayS[p-2]/scale; 781 double epm1=e[p-2]/scale; 782 double sk=arrayS[k]/scale; 783 double ek=e[k]/scale; 784 double b=((spm1+sp)*(spm1-sp)+epm1*epm1)/2.0; 785 double c=(sp*epm1)*(sp*epm1); 786 double shift=0.0; 787 if(b!=0.0 || c!=0.0) { 788 shift=Math.sqrt(b*b+c); 789 if(b<0.0) 790 shift=-shift; 791 shift=c/(b+shift); 792 } 793 double f=(sk+sp)*(sk-sp)+shift; 794 double g=sk*ek; 795 for(int j=k;j<p-1;j++) { 797 double t=ExtraMath.hypot(f,g); 798 double cs=f/t; 799 double sn=g/t; 800 if(j!=k) 801 e[j-1]=t; 802 f=cs*arrayS[j]+sn*e[j]; 803 e[j]=cs*e[j]-sn*arrayS[j]; 804 g=sn*arrayS[j+1]; 805 arrayS[j+1]*=cs; 806 for(int i=0;i<N;i++) { 807 t=cs*arrayV[i][j]+sn*arrayV[i][j+1]; 808 arrayV[i][j+1]=-sn*arrayV[i][j]+cs*arrayV[i][j+1]; 809 arrayV[i][j]=t; 810 } 811 t=ExtraMath.hypot(f,g); 812 cs=f/t; 813 sn=g/t; 814 arrayS[j]=t; 815 f=cs*e[j]+sn*arrayS[j+1]; 816 arrayS[j+1]=-sn*e[j]+cs*arrayS[j+1]; 817 g=sn*e[j+1]; 818 e[j+1]*=cs; 819 if(j<Nm1) { 820 for(int i=0;i<N;i++) { 821 t=cs*arrayU[i][j]+sn*arrayU[i][j+1]; 822 arrayU[i][j+1]=-sn*arrayU[i][j]+cs*arrayU[i][j+1]; 823 arrayU[i][j]=t; 824 } 825 } 826 } 827 e[p-2]=f; 828 iter++; 829 } break; 830 case 4: { 832 if(arrayS[k]<=0.0) { 834 arrayS[k]=-arrayS[k]; 835 for(int i=0;i<p;i++) 836 arrayV[i][k]=-arrayV[i][k]; 837 } 838 while(k<p-1) { 840 if(arrayS[k]>=arrayS[k+1]) 841 break; 842 double tmp=arrayS[k]; 843 arrayS[k]=arrayS[k+1]; 844 arrayS[k+1]=tmp; 845 if(k<Nm1) { 846 for(int i=0;i<N;i++) { 847 tmp=arrayU[i][k+1]; 848 arrayU[i][k+1]=arrayU[i][k]; 849 arrayU[i][k]=tmp; 850 tmp=arrayV[i][k+1]; 851 arrayV[i][k+1]=arrayV[i][k]; 852 arrayV[i][k]=tmp; 853 } 854 } 855 k++; 856 } 857 iter=0; 858 p--; 859 } break; 860 } 861 } 862 final AbstractDoubleSquareMatrix svd[]=new AbstractDoubleSquareMatrix[3]; 863 svd[0]=new DoubleSquareMatrix(arrayU); 864 svd[1]=new DoubleDiagonalMatrix(arrayS); 865 svd[2]=new DoubleSquareMatrix(arrayV); 866 return svd; 867 } 868 869 870 private static class SquareMatrixAdaptor extends AbstractIntegerSquareMatrix { 871 private final AbstractIntegerMatrix matrix; 872 private SquareMatrixAdaptor(AbstractIntegerMatrix m) { 873 super(m.rows()); 874 matrix = m; 875 } 876 public int getElement(int i, int j) { 877 return matrix.getElement(i,j); 878 } 879 public void setElement(int i, int j, int x) { 880 matrix.setElement(i,j, x); 881 } 882 } 883 } 884 | Popular Tags |