template<class C>
| void Impala::Core::Matrix::tred2 | ( | MatrixTem< C > & | a, | |
| VectorTem< C > & | d, | |||
| VectorTem< C > & | e | |||
| ) |
Definition at line 636 of file MatrixTem.h.
References Impala::Core::Vector::VectorTem< ElemT >::Size().
00637 { 00638 int i,j,k,l; 00639 double scale,hh,h,g,f; 00640 00641 int n = d.Size() ; 00642 00643 for (i=n-1; i>0; i--) 00644 { 00645 l = i-1; 00646 h = scale = 0.0; 00647 if (l>0) 00648 { 00649 for (k=0; k<=l; k++) scale += fabs(a[i][k]); 00650 if (scale == 0.0) e[i] = a[i][l]; 00651 else 00652 { 00653 for (k=0; k<=l; k++) 00654 { 00655 a[i][k] /= scale; //use scale a's for transformation 00656 h += a[i][k]*a[i][k]; 00657 } 00658 f = a[i][l]; 00659 g = (f >=0.0 ? -sqrt(h) : sqrt(h)); 00660 e[i] = scale*g; 00661 h -= f*g; 00662 a[i][l] = f-g; 00663 f = 0.0; 00664 for (j=0; j<=l; j++) 00665 { 00666 a[j][i] = a[i][j]/h; 00667 g = 0.0; 00668 for (k=0; k<=j; k++) g += a[j][k]*a[i][k]; 00669 for (k=j+1; k<=l; k++) g += a[k][j]*a[i][k]; 00670 e[j] = g/h; //form element of p in temporarily unused element of e 00671 f += e[j]*a[i][j]; 00672 } 00673 hh = f/(h+h); 00674 for (j=0; j<=l; j++) 00675 { 00676 f = a[i][j]; 00677 e[j] = g = e[j]-hh*f; 00678 for (k=0; k<=j; k++) a[j][k] -= (f*e[k]+g*a[i][k]); 00679 } 00680 } 00681 } 00682 else e[i] = a[i][l]; 00683 d[i] = h; 00684 } 00685 00686 d[0] = 0.0; 00687 e[0] = 0.0; 00688 00689 for (i=0; i<n; i++) 00690 { 00691 l = i-1; 00692 if (d[i]) 00693 { 00694 for (j=0; j<=l; j++) 00695 { 00696 g = 0.0; 00697 for (k=0; k<=l; k++) g += a[i][k]*a[k][j]; 00698 for (k=0; k<=l; k++) a[k][j] -= g*a[k][i]; 00699 } 00700 } 00701 d[i] = a[i][i]; 00702 a[i][i] = 1.0; 00703 for (j=0; j<=l; j++) a[j][i] = a[i][j] = 0.0; 00704 } 00705 00706 }
Here is the call graph for this function:
