|
template<class C>
Definition at line 628 of file Matrix.h. Referenced by CxMatrixTem< C >::klm(). 00629 { 00630 int i,j,k,l; 00631 double scale,hh,h,g,f; 00632 00633 int n = d.Size() ; 00634 00635 for (i=n-1; i>0; i--) 00636 { 00637 l = i-1; 00638 h = scale = 0.0; 00639 if (l>0) 00640 { 00641 for (k=0; k<=l; k++) scale += fabs(a[i][k]); 00642 if (scale == 0.0) e[i] = a[i][l]; 00643 else 00644 { 00645 for (k=0; k<=l; k++) 00646 { 00647 a[i][k] /= scale; //use scale a's for transformation 00648 h += a[i][k]*a[i][k]; 00649 } 00650 f = a[i][l]; 00651 g = (f >=0.0 ? -sqrt(h) : sqrt(h)); 00652 e[i] = scale*g; 00653 h -= f*g; 00654 a[i][l] = f-g; 00655 f = 0.0; 00656 for (j=0; j<=l; j++) 00657 { 00658 a[j][i] = a[i][j]/h; 00659 g = 0.0; 00660 for (k=0; k<=j; k++) g += a[j][k]*a[i][k]; 00661 for (k=j+1; k<=l; k++) g += a[k][j]*a[i][k]; 00662 e[j] = g/h; //form element of p in temporarily unused element of e 00663 f += e[j]*a[i][j]; 00664 } 00665 hh = f/(h+h); 00666 for (j=0; j<=l; j++) 00667 { 00668 f = a[i][j]; 00669 e[j] = g = e[j]-hh*f; 00670 for (k=0; k<=j; k++) a[j][k] -= (f*e[k]+g*a[i][k]); 00671 } 00672 } 00673 } 00674 else e[i] = a[i][l]; 00675 d[i] = h; 00676 } 00677 00678 d[0] = 0.0; 00679 e[0] = 0.0; 00680 00681 for (i=0; i<n; i++) 00682 { 00683 l = i-1; 00684 if (d[i]) 00685 { 00686 for (j=0; j<=l; j++) 00687 { 00688 g = 0.0; 00689 for (k=0; k<=l; k++) g += a[i][k]*a[k][j]; 00690 for (k=0; k<=l; k++) a[k][j] -= g*a[k][i]; 00691 } 00692 } 00693 d[i] = a[i][i]; 00694 a[i][i] = 1.0; 00695 for (j=0; j<=l; j++) a[j][i] = a[i][j] = 0.0; 00696 } 00697 00698 }
|
1.5.1