1 package JSci.maths; 2 3 import JSci.maths.matrices.*; 4 import JSci.maths.vectors.*; 5 import JSci.maths.polynomials.RealPolynomial; 6 7 13 public final class LinearMath extends AbstractMath { 14 private LinearMath() {} 15 16 18 24 public static AbstractDoubleVector solve(final AbstractDoubleSquareMatrix M,final AbstractDoubleVector v) { 25 final int n=v.dimension(); 26 final double array[]=new double[n]; 27 final int pivot[]=new int[n+1]; 28 final AbstractDoubleSquareMatrix lu[]=M.luDecompose(pivot); 29 int i,j; 30 double sum; 31 for(i=0;i<n;i++) { 33 sum=v.getComponent(pivot[i]); 34 for(j=0;j<i;j++) 35 sum-=lu[0].getElement(i,j)*array[j]; 36 array[i]=sum/lu[0].getElement(i,i); 37 } 38 for(i=n-1;i>=0;i--) { 40 sum=array[i]; 41 for(j=i+1;j<n;j++) 42 sum-=lu[1].getElement(i,j)*array[j]; 43 array[i]=sum/lu[1].getElement(i,i); 44 } 45 return new DoubleVector(array); 46 } 47 50 private static double[] solveCholesky(final double m[][], final double v[]) { 51 final int n=v.length; 52 final double array[]=new double[n]; 53 final AbstractDoubleSquareMatrix lu[]=new DoubleSquareMatrix(m).choleskyDecompose(); 54 int i,j; 55 double sum; 56 for(i=0;i<n;i++) { 58 sum=v[i]; 59 for(j=0;j<i;j++) 60 sum-=lu[0].getElement(i,j)*array[j]; 61 array[i]=sum/lu[0].getElement(i,i); 62 } 63 for(i=n-1;i>=0;i--) { 65 sum=array[i]; 66 for(j=i+1;j<n;j++) 67 sum-=lu[1].getElement(i,j)*array[j]; 68 array[i]=sum/lu[1].getElement(i,i); 69 } 70 return array; 71 } 72 75 private static double[] solveQR(final double m[][], final double v[]) { 76 final int n=v.length; 77 final double array[]=new double[n]; 78 final AbstractDoubleSquareMatrix lu[]=new DoubleSquareMatrix(m).qrDecompose(); 79 int i,j; 80 double sum; 81 for(i=0;i<n;i++) { 83 sum=0.0; 84 for(j=0;j<n;j++) 85 sum+=lu[0].getElement(j,i)*v[j]; 86 array[i]=sum; 87 } 88 for(i=n-1;i>=0;i--) { 90 sum=array[i]; 91 for(j=i+1;j<n;j++) 92 sum-=lu[1].getElement(i,j)*array[j]; 93 array[i]=sum/lu[1].getElement(i,i); 94 } 95 return array; 96 } 97 115 public static AbstractDoubleVector solveGMRes(final AbstractDoubleMatrix A,final AbstractDoubleVector b,int max_iter,double tol) throws MaximumIterationsExceededException { 116 if(max_iter<=0) 117 throw new IllegalArgumentException ("Number of allowed iterations must be a positive integer: "+max_iter+" <= 0."); 118 if(tol<0) 119 throw new IllegalArgumentException ("Tolerance must be positive or zero: "+tol+" < 0."); 120 final int m=A.rows(); 121 double resid; 122 int i,j=1,k; 123 double[] s=new double[m+1]; 124 double[][] cs=new double[m+1][2]; 125 double[] rotmp=new double[2]; 126 AbstractDoubleVector w; 127 AbstractDoubleVector x=new DoubleVector(A.rows()); 128 double normb=b.norm(); 129 AbstractDoubleVector r=b.subtract(A.multiply(x)); 130 double beta=r.norm(); 131 if(normb==0.0) 132 normb=1.0; 133 if((resid=r.norm()/normb)<= tol) { 134 tol = resid; 135 max_iter = 0; 136 throw new IllegalArgumentException ("There is a bug."); 137 } 138 AbstractDoubleVector[] v = new DoubleVector[m+1]; 139 double[][] H = new double[m+1][m]; 140 while (j <= max_iter) { 141 v[0] = r.scalarMultiply(1.0 / beta); 142 for (i = 0; i < m+1; i++) 143 s[i] = 0.0; 144 s[0] = beta; 145 for (i = 0; i < m && j <= max_iter; i++, j++) { 146 w = A.multiply(v[i]); 147 for (k = 0; k <= i; k++) { 148 H[k][i]=w.scalarProduct(v[k]); 149 w = w.subtract(v[k].scalarMultiply(H[k][i])); 150 } 151 H[i+1][i]=w.norm(); 152 v[i+1] = w.scalarMultiply(1.0 / H[i+1][i]); 153 for (k = 0; k < i; k++) { 154 rotmp=applyPlaneRotation(H[k][i], H[k+1][i], cs[k][0], cs[k][1]); 155 H[k][i]=rotmp[0]; 156 H[k+1][i]=rotmp[1]; 157 } 158 cs[i]=generatePlaneRotation(H[i][i], H[i+1][i]); 159 rotmp=applyPlaneRotation(H[i][i], H[i+1][i], cs[i][0], cs[i][1]); 160 H[i][i]=rotmp[0]; 161 H[i+1][i]=rotmp[1]; 162 rotmp=applyPlaneRotation(s[i], s[i+1], cs[i][0], cs[i][1]); 163 s[i]=rotmp[0]; 164 s[i+1]=rotmp[1]; 165 if ((resid = Math.abs(s[i+1]) / normb) < tol) { 166 x = update(x, i, H, s, v); 167 tol = resid; 168 max_iter = j; 169 return (x); 170 } 171 } 172 x = update(x, m - 1, H, s, v); 173 r = b.subtract(A.multiply(x)); 174 beta = r.norm(); 175 if ((resid = beta / normb) < tol) { 176 tol = resid; 177 max_iter = j; 178 return (x); 179 } 180 } 181 tol = resid; 182 throw new MaximumIterationsExceededException("(tol) "+tol+". It doesn't converge in "+max_iter+"iterations. Try raising the number of allowed iterations or raising the tolerance."); 183 } 184 private static double[] generatePlaneRotation(double dx, double dy) { 185 double [] cs = new double[2]; 186 double temp; 187 if (dy == 0.0) { 188 cs[0] = 1.0; 189 cs[1] = 0.0; 190 } else if (Math.abs(dy) > Math.abs(dx)) { 191 temp = dx / dy; 192 cs[1] = 1.0 / Math.sqrt( 1.0 + temp*temp ); 193 cs[0] = temp * cs[1]; 194 } else { 195 temp = dy / dx; 196 cs[0] = 1.0 / Math.sqrt( 1.0 + temp*temp ); 197 cs[1] = temp * cs[0]; 198 } 199 return (cs); 200 } 201 private static double[] applyPlaneRotation(double dx, double dy, double cs, double sn) { 202 double [] dxy = new double[2]; 203 dxy[0] = cs * dx + sn * dy; 204 dxy[1] = -sn * dx + cs * dy; 205 return (dxy); 206 } 207 private static AbstractDoubleVector update(AbstractDoubleVector x, int k, double[][] H, double [] s, AbstractDoubleVector [] v) { 208 211 for (int i = k; i >= 0; i--) { 212 s[i]=s[i]/H[i][i]; 213 for (int j = i - 1; j >= 0; j--) 214 s[j]=s[j]-H[j][i]*s[i]; 215 } 216 for (int j = 0; j <= k; j++){ 217 x = x.add(v[j].scalarMultiply(s[j])); 218 } 219 return(x); 220 } 221 222 224 231 public static RealPolynomial leastSquaresFit(int n, final double data[][]) { 232 int i,j,nm1=n++; 233 double xsum,xysum,tmp; 234 final double mArray[][]=new double[n][n]; 235 final double vArray[]=new double[n]; 236 for(i=0;i<n;i++) { 237 xsum=xysum=0.0; 238 for(j=0;j<data[0].length;j++) { 239 tmp=Math.pow(data[0][j],i); 240 xsum+=tmp; 241 xysum+=tmp*data[1][j]; 242 } 243 mArray[0][i]=xsum; 244 vArray[i]=xysum; 245 } 246 for(i=1;i<n;i++) { 247 System.arraycopy(mArray[i-1],1,mArray[i],0,nm1); 248 xsum=0.0; 249 for(j=0;j<data[0].length;j++) 250 xsum+=Math.pow(data[0][j],nm1+i); 251 mArray[i][nm1]=xsum; 252 } 253 return new RealPolynomial(solveCholesky(mArray,vArray)); 254 } 255 256 258 266 public static AbstractDoubleVector linearRegression(final double data[][]) { 267 final int y=data.length-1; 268 int i,j,k; 269 double xsum,xysum; 270 final double mArray[][]=new double[data.length][data.length]; 271 final double vArray[]=new double[data.length]; 272 mArray[0][0]=data[0].length; 273 for(j=1;j<data.length;j++) { 274 xsum=0.0; 275 for(k=0;k<data[0].length;k++) 276 xsum+=data[j-1][k]; 277 mArray[0][j]=mArray[j][0]=xsum; 278 } 279 xysum=0.0; 280 for(k=0;k<data[0].length;k++) 281 xysum+=data[y][k]; 282 vArray[0]=xysum; 283 for(i=1;i<data.length;i++) { 284 for(j=i;j<data.length;j++) { 285 xsum=0.0; 286 for(k=0;k<data[0].length;k++) 287 xsum+=data[i-1][k]*data[j-1][k]; 288 mArray[i][j]=mArray[j][i]=xsum; 289 } 290 xysum=0.0; 291 for(k=0;k<data[0].length;k++) 292 xysum+=data[i-1][k]*data[y][k]; 293 vArray[i]=xysum; 294 } 295 return new DoubleVector(solveCholesky(mArray,vArray)); 296 } 297 298 300 305 public static AbstractDoubleVector[] orthonormalize(AbstractDoubleVector vecs[]) { 306 final int N=vecs.length; 307 AbstractDoubleVector orthovecs[]=new DoubleVector[N]; 308 for(int i=0;i<N;i++) { 309 orthovecs[i]=vecs[i]; 310 for(int j=0;j<i;j++) 311 orthovecs[i]=orthovecs[i].subtract(orthovecs[j].scalarMultiply(orthovecs[j].scalarProduct(vecs[i]))); 312 orthovecs[i].normalize(); 313 } 314 return orthovecs; 315 } 316 317 319 325 public static double[] eigenvalueSolveHermitian(final AbstractComplexSquareMatrix matrix) throws MaximumIterationsExceededException { 326 final int n=matrix.rows(); 327 final double matrix2[][]=new double[2*n][2*n]; 328 double real,imag; 329 for(int j,i=0;i<n;i++) { 330 for(j=0;j<n;j++) { 331 real=matrix.getRealElement(i,j); 332 imag=matrix.getImagElement(i,j); 333 matrix2[i][j]=real; 334 matrix2[n+i][n+j]=real; 335 matrix2[n+i][j]=imag; 336 matrix2[i][n+j]=-imag; 337 } 338 } 339 final double eigenvalue2[]=new double[2*n]; 340 final double offdiag[]=new double[2*n]; 341 reduceSymmetric1_SquareToTridiagonal(matrix2,eigenvalue2,offdiag); 342 System.arraycopy(offdiag,1,offdiag,0,n-1); 343 offdiag[n-1]=0.0; 344 eigenvalueSolveSymmetricTridiagonalMatrix(eigenvalue2,offdiag); 345 final double eigenvalue[]=new double[n]; 346 System.arraycopy(eigenvalue2,0,eigenvalue,0,n); 347 return eigenvalue; 348 } 349 357 public static double[] eigenSolveHermitian(final AbstractComplexSquareMatrix matrix,final AbstractComplexVector eigenvector[]) throws MaximumIterationsExceededException { 358 final int n=matrix.rows(); 359 final double matrix2[][]=new double[2*n][2*n]; 360 int i,j; 361 double real,imag; 362 for(i=0;i<n;i++) { 363 for(j=0;j<n;j++) { 364 real=matrix.getRealElement(i,j); 365 imag=matrix.getImagElement(i,j); 366 matrix2[i][j]=real; 367 matrix2[n+i][n+j]=real; 368 matrix2[n+i][j]=imag; 369 matrix2[i][n+j]=-imag; 370 } 371 } 372 final double eigenvalue2[]=new double[2*n]; 373 final double offdiag[]=new double[2*n]; 374 reduceSymmetric2_SquareToTridiagonal(matrix2,eigenvalue2,offdiag); 375 System.arraycopy(offdiag,1,offdiag,0,n-1); 376 offdiag[n-1]=0.0; 377 eigenSolveSymmetricTridiagonalMatrix(eigenvalue2,offdiag,matrix2); 378 final double eigenvalue[]=new double[n]; 379 double arrayRe[],arrayIm[]; 380 for(i=0;i<n;i++) { 381 eigenvalue[i]=eigenvalue2[i]; 382 arrayRe=new double[n]; 383 arrayIm=new double[n]; 384 for(j=0;j<n;j++) { 385 arrayRe[j]=matrix2[j][i]; 386 arrayIm[j]=matrix2[j+n][i]; 387 } 388 eigenvector[i]=new ComplexVector(arrayRe,arrayIm); 389 } 390 return eigenvalue; 391 } 392 399 public static double[] eigenvalueSolveSymmetric(final DoubleTridiagonalMatrix matrix) throws MaximumIterationsExceededException { 400 final int n=matrix.rows(); 401 final int nm1=n-1; 402 final double eigenvalue[]=new double[n]; 403 final double offdiag[]=new double[n]; 404 for(int i=0;i<nm1;i++) { 405 eigenvalue[i]=matrix.getElement(i,i); 406 offdiag[i]=matrix.getElement(i,i+1); 407 } 408 eigenvalue[nm1]=matrix.getElement(nm1,nm1); 409 offdiag[nm1]=0.0; 410 eigenvalueSolveSymmetricTridiagonalMatrix(eigenvalue,offdiag); 411 return eigenvalue; 412 } 413 422 public static double[] eigenSolveSymmetric(final DoubleTridiagonalMatrix matrix,final AbstractDoubleVector eigenvector[]) throws MaximumIterationsExceededException { 423 final int n=matrix.rows(); 424 final int nm1=n-1; 425 final double eigenvalue[]=new double[n]; 426 final double offdiag[]=new double[n]; 427 final double id[][]=new double[n][n]; 428 int i,j; 429 for(i=0;i<nm1;i++) { 430 id[i][i]=1.0; 431 eigenvalue[i]=matrix.getElement(i,i); 432 offdiag[i]=matrix.getElement(i,i+1); 433 } 434 id[nm1][nm1]=1.0; 435 eigenvalue[nm1]=matrix.getElement(nm1,nm1); 436 offdiag[nm1]=0.0; 437 eigenSolveSymmetricTridiagonalMatrix(eigenvalue,offdiag,id); 438 for(i=0;i<n;i++) { 439 DoubleVector evec = new DoubleVector(n); 440 for(j=0; j<n; j++) 441 evec.setComponent(j, id[j][i]); 442 eigenvector[i] = evec; 443 } 444 return eigenvalue; 445 } 446 454 public static double[] eigenvalueSolveSymmetric(final AbstractDoubleSquareMatrix matrix) throws MaximumIterationsExceededException { 455 final int n=matrix.rows(); 456 final double eigenvalue[]=new double[n]; 457 final double offdiag[]=new double[n]; 458 final double array[][]=new double[n][n]; 459 int i,j; 460 for(i=0;i<n;i++) { 461 for(j=0;j<n;j++) 462 array[i][j]=matrix.getElement(i,j); 463 } 464 reduceSymmetric1_SquareToTridiagonal(array,eigenvalue,offdiag); 465 System.arraycopy(offdiag,1,offdiag,0,n-1); 466 offdiag[n-1]=0.0; 467 eigenvalueSolveSymmetricTridiagonalMatrix(eigenvalue,offdiag); 468 return eigenvalue; 469 } 470 480 public static double[] eigenSolveSymmetric(final AbstractDoubleSquareMatrix matrix,final AbstractDoubleVector eigenvector[]) throws MaximumIterationsExceededException { 481 final int n=matrix.rows(); 482 final double eigenvalue[]=new double[n]; 483 final double offdiag[]=new double[n]; 484 final double transf[][]=new double[n][n]; 485 int i,j; 486 for(i=0;i<n;i++) { 487 for(j=0;j<n;j++) 488 transf[i][j]=matrix.getElement(i,j); 489 } 490 reduceSymmetric2_SquareToTridiagonal(transf,eigenvalue,offdiag); 491 System.arraycopy(offdiag,1,offdiag,0,n-1); 492 offdiag[n-1]=0.0; 493 eigenSolveSymmetricTridiagonalMatrix(eigenvalue,offdiag,transf); 494 for(i=0;i<n;i++) { 495 DoubleVector evec = new DoubleVector(n); 496 for(j=0; j<n; j++) 497 evec.setComponent(j, transf[j][i]); 498 eigenvector[i] = evec; 499 } 500 return eigenvalue; 501 } 502 503 private final static int EIGEN_MAX_ITERATIONS = 250; 504 505 510 private static void eigenvalueSolveSymmetricTridiagonalMatrix(final double diag[],final double offdiag[]) throws MaximumIterationsExceededException { 511 final int n=diag.length; 512 final int nm1=n-1; 513 int m,l,iteration,i,k; 514 double s,r,p,g,f,dd,c,b; 515 for(l=0;l<n;l++) { 516 iteration=0; 517 do { 518 for(m=l;m<nm1;m++) { 519 dd=Math.abs(diag[m])+Math.abs(diag[m+1]); 520 if(Math.abs(offdiag[m])+dd==dd) 521 break; 522 } 523 if(m!=l) { 524 if(iteration++ == EIGEN_MAX_ITERATIONS) 525 throw new MaximumIterationsExceededException("No convergence after "+EIGEN_MAX_ITERATIONS+" iterations."); 526 g=(diag[l+1]-diag[l])/(2.0*offdiag[l]); 527 r=Math.sqrt(g*g+1.0); 528 g=diag[m]-diag[l]+offdiag[l]/(g+(g<0.0?-Math.abs(r):Math.abs(r))); 529 s=c=1.0; 530 p=0.0; 531 for(i=m-1;i>=l;i--) { 532 f=s*offdiag[i]; 533 b=c*offdiag[i]; 534 if(Math.abs(f)>=Math.abs(g)) { 535 c=g/f; 536 r=Math.sqrt(c*c+1.0); 537 offdiag[i+1]=f*r; 538 s=1/r; 539 c*=s; 540 } else { 541 s=f/g; 542 r=Math.sqrt(s*s+1.0); 543 offdiag[i+1]=g*r; 544 c=1/r; 545 s*=c; 546 } 547 g=diag[i+1]-p; 548 r=(diag[i]-g)*s+2.0*c*b; 549 p=s*r; 550 diag[i+1]=g+p; 551 g=c*r-b; 552 } 553 diag[l]=diag[l]-p; 554 offdiag[l]=g; 555 offdiag[m]=0.0; 556 } 557 } while(m!=l); 558 } 559 } 560 564 private static void reduceSymmetric1_SquareToTridiagonal(final double matrix[][],final double diag[],final double offdiag[]) { 565 final int n=diag.length; 566 int i,j,k,l; 567 double f,g,h,hh,scale; 568 for(i=n-1;i>0;i--) { 569 l=i-1; 570 h=scale=0.0; 571 if(l>0) { 572 for(k=0;k<=l;k++) 573 scale+=Math.abs(matrix[i][k]); 574 if(scale==0.0) 575 offdiag[i]=matrix[i][l]; 576 else { 577 for(k=0;k<=l;k++) { 578 matrix[i][k]/=scale; 579 h+=matrix[i][k]*matrix[i][k]; 580 } 581 f=matrix[i][l]; 582 g=(f>=0.0?-Math.sqrt(h):Math.sqrt(h)); 583 offdiag[i]=scale*g; 584 h-=f*g; 585 matrix[i][l]=f-g; 586 f=0.0; 587 for(j=0;j<=l;j++) { 588 g=0.0; 589 for(k=0;k<=j;k++) 590 g+=matrix[j][k]*matrix[i][k]; 591 for(k=j+1;k<=l;k++) 592 g+=matrix[k][j]*matrix[i][k]; 593 offdiag[j]=g/h; 594 f+=offdiag[j]*matrix[i][j]; 595 } 596 hh=f/(h+h); 597 for(j=0;j<=l;j++) { 598 f=matrix[i][j]; 599 offdiag[j]=g=offdiag[j]-hh*f; 600 for(k=0;k<=j;k++) 601 matrix[j][k]-=f*offdiag[k]+g*matrix[i][k]; 602 } 603 } 604 } else 605 offdiag[i]=matrix[i][l]; 606 diag[i]=h; 607 } 608 offdiag[0]=0.0; 609 for(i=0;i<n;i++) 610 diag[i]=matrix[i][i]; 611 } 612 618 private static void eigenSolveSymmetricTridiagonalMatrix(final double diag[],final double offdiag[],final double transf[][]) throws MaximumIterationsExceededException { 619 final int n=diag.length; 620 final int nm1=n-1; 621 int m,l,iteration,i,k; 622 double s,r,p,g,f,dd,c,b; 623 for(l=0;l<n;l++) { 624 iteration=0; 625 do { 626 for(m=l;m<nm1;m++) { 627 dd=Math.abs(diag[m])+Math.abs(diag[m+1]); 628 if(Math.abs(offdiag[m])+dd==dd) 629 break; 630 } 631 if(m!=l) { 632 if(iteration++ == EIGEN_MAX_ITERATIONS) 633 throw new MaximumIterationsExceededException("No convergence after "+EIGEN_MAX_ITERATIONS+" iterations."); 634 g=(diag[l+1]-diag[l])/(2.0*offdiag[l]); 635 r=Math.sqrt(g*g+1.0); 636 g=diag[m]-diag[l]+offdiag[l]/(g+(g<0.0?-Math.abs(r):Math.abs(r))); 637 s=c=1.0; 638 p=0.0; 639 for(i=m-1;i>=l;i--) { 640 f=s*offdiag[i]; 641 b=c*offdiag[i]; 642 if(Math.abs(f)>=Math.abs(g)) { 643 c=g/f; 644 r=Math.sqrt(c*c+1.0); 645 offdiag[i+1]=f*r; 646 s=1/r; 647 c*=s; 648 } else { 649 s=f/g; 650 r=Math.sqrt(s*s+1.0); 651 offdiag[i+1]=g*r; 652 c=1/r; 653 s*=c; 654 } 655 g=diag[i+1]-p; 656 r=(diag[i]-g)*s+2.0*c*b; 657 p=s*r; 658 diag[i+1]=g+p; 659 g=c*r-b; 660 for(k=0;k<n;k++) { 661 f=transf[k][i+1]; 662 transf[k][i+1]=s*transf[k][i]+c*f; 663 transf[k][i]=c*transf[k][i]-s*f; 664 } 665 } 666 diag[l]=diag[l]-p; 667 offdiag[l]=g; 668 offdiag[m]=0.0; 669 } 670 } while(m!=l); 671 } 672 } 673 678 private static void reduceSymmetric2_SquareToTridiagonal(final double matrix[][],final double diag[],final double offdiag[]) { 679 final int n=diag.length; 680 int i,j,k,l; 681 double f,g,h,hh,scale; 682 for(i=n-1;i>0;i--) { 683 l=i-1; 684 h=scale=0.0; 685 if(l>0) { 686 for(k=0;k<=l;k++) 687 scale+=Math.abs(matrix[i][k]); 688 if(scale==0.0) 689 offdiag[i]=matrix[i][l]; 690 else { 691 for(k=0;k<=l;k++) { 692 matrix[i][k]/=scale; 693 h+=matrix[i][k]*matrix[i][k]; 694 } 695 f=matrix[i][l]; 696 g=(f>=0.0?-Math.sqrt(h):Math.sqrt(h)); 697 offdiag[i]=scale*g; 698 h-=f*g; 699 matrix[i][l]=f-g; 700 f=0.0; 701 for(j=0;j<=l;j++) { 702 matrix[j][i]=matrix[i][j]/h; 703 g=0.0; 704 for(k=0;k<=j;k++) 705 g+=matrix[j][k]*matrix[i][k]; 706 for(k=j+1;k<=l;k++) 707 g+=matrix[k][j]*matrix[i][k]; 708 offdiag[j]=g/h; 709 f+=offdiag[j]*matrix[i][j]; 710 } 711 hh=f/(h+h); 712 for(j=0;j<=l;j++) { 713 f=matrix[i][j]; 714 offdiag[j]=g=offdiag[j]-hh*f; 715 for(k=0;k<=j;k++) 716 matrix[j][k]-=f*offdiag[k]+g*matrix[i][k]; 717 } 718 } 719 } else 720 offdiag[i]=matrix[i][l]; 721 diag[i]=h; 722 } 723 diag[0]=offdiag[0]=0.0; 724 for(i=0;i<n;i++) { 725 l=i-1; 726 if(diag[i]!=0.0) { 727 for(j=0;j<=l;j++) { 728 g=0.0; 729 for(k=0;k<=l;k++) 730 g+=matrix[i][k]*matrix[k][j]; 731 for(k=0;k<=l;k++) 732 matrix[k][j]-=g*matrix[k][i]; 733 } 734 } 735 diag[i]=matrix[i][i]; 736 matrix[i][i]=1.0; 737 for(j=0;j<=l;j++) 738 matrix[j][i]=matrix[i][j]=0.0; 739 } 740 } 741 } 742 | Popular Tags |