Impala Documentation

template<class MatrixT, class VectorT>
void Impala::Core::Matrix::Tred2	(	MatrixT *	a,
		VectorT *	d,
		VectorT *	e
	)
Definition at line 24 of file MatKLM.h.
References MatE(), VecE(), and VecNrElem().
Referenced by MatKLM().
00025 {
00026     int i,j,k,l;
00027     double scale,hh,h,g,f;
00028 
00029     int n = VecNrElem(d);
00030 
00031     for (i=n-1; i>0; i--)
00032     {
00033         l = i-1;
00034         h = scale = 0.0;
00035         if (l>0)
00036         {
00037             for (k=0; k<=l; k++)
00038                 scale += fabs(*MatE(a, i, k));
00039             if (scale == 0.0)
00040                 *VecE(e, i) = *MatE(a, i, l);
00041             else
00042             {
00043                 for (k=0; k<=l; k++)
00044                 {
00045                     *MatE(a, i, k) /= scale; //use scale a's for transformation
00046                     h += *MatE(a, i, k) * *MatE(a, i, k);
00047                 }
00048                 f = *MatE(a, i, l);
00049                 g = (f >=0.0 ? -sqrt(h) : sqrt(h));
00050                 *VecE(e, i) = scale*g;
00051                 h -= f*g;
00052                 *MatE(a, i, l) = f-g;
00053                 f = 0.0;
00054                 for (j=0; j<=l; j++)
00055                 {
00056                     *MatE(a, j, i) = *MatE(a, i, j)/h;
00057                     g = 0.0;
00058                     for (k=0; k<=j; k++)
00059                         g += *MatE(a, j, k) * *MatE(a, i, k);
00060                     for (k=j+1; k<=l; k++)
00061                         g += *MatE(a, k, j) * *MatE(a, i, k);
00062                     *VecE(e, j) = g/h; //form element of p in temporarily unused element of e
00063                     f += *VecE(e, j) * *MatE(a, i, j);
00064                 }
00065                 hh = f/(h+h);
00066                 for (j=0; j<=l; j++)
00067                 {
00068                     f = *MatE(a, i, j);
00069                     *VecE(e, j) = g = *VecE(e, j)-hh*f;
00070                     for (k=0; k<=j; k++)
00071                         *MatE(a, j, k) -= (f * *VecE(e, k) 
00072                                              + g * *MatE(a, i, k));
00073                 }
00074             }
00075         }
00076         else
00077             *VecE(e, i) = *MatE(a, i, l);
00078         *VecE(d, i) = h;
00079     }
00080 
00081     *VecE(d, 0) = 0.0;
00082     *VecE(e, 0) = 0.0;
00083 
00084     for (i=0; i<n; i++)
00085     {
00086         l = i-1;
00087         if (*VecE(d, i))
00088         {
00089             for (j=0; j<=l; j++)
00090             {
00091                 g = 0.0;
00092                 for (k=0; k<=l; k++)
00093                     g += *MatE(a, i, k) * *MatE(a, k, j);
00094                 for (k=0; k<=l; k++)
00095                     *MatE(a, k, j) -= g * *MatE(a, k, i);
00096             }
00097         }
00098         *VecE(d, i) = *MatE(a, i, i);
00099         *MatE(a, i, i) = 1.0;
00100         for (j=0; j<=l; j++)
00101             *MatE(a, j, i) = *MatE(a, i, j) = 0.0;
00102     }
00103 }
Here is the call graph for this function: