GenericMathTemplateLibrary: Eigen.h Source File

Go to the documentation of this file.
00001 // GMTL is (C) Copyright 2001-2010 by Allen Bierbaum
00002 // Distributed under the GNU Lesser General Public License 2.1 with an
00003 // addendum covering inlined code. (See accompanying files LICENSE and
00004 // LICENSE.addendum or http://www.gnu.org/copyleft/lesser.txt)
00005 
00006 //
00007 // XXX: Based on source code from: Magic Software, Inc.
00008 //
00009 #ifndef _EIGEN_H
00010 #define _EIGEN_H
00011 
00012 #include <gmtl/gmtlConfig.h>
00013 
00014 namespace gmtl
00015 {
00016 
00017 class Eigen
00018 {
00019 public:
00020     Eigen (int iSize);
00021     ~Eigen ();
00022 
00023     // set the matrix for eigensolving
00024     float& Matrix (int iRow, int iCol);
00025     void SetMatrix (float** aafMat);
00026 
00027     // get the results of eigensolving (eigenvectors are columns of matrix)
00028     float GetEigenvalue (int i) const;
00029     float GetEigenvector (int iRow, int iCol) const;
00030     float* GetEigenvalue ();
00031     float** GetEigenvector ();
00032 
00033     // solve eigensystem
00034     void EigenStuff2 ();
00035     void EigenStuff3 ();
00036     void EigenStuff4 ();
00037     void EigenStuffN ();
00038     void EigenStuff  ();
00039 
00040     // solve eigensystem, use decreasing sort on eigenvalues
00041     void DecrSortEigenStuff2 ();
00042     void DecrSortEigenStuff3 ();
00043     void DecrSortEigenStuff4 ();
00044     void DecrSortEigenStuffN ();
00045     void DecrSortEigenStuff  ();
00046 
00047     // solve eigensystem, use increasing sort on eigenvalues
00048     void IncrSortEigenStuff2 ();
00049     void IncrSortEigenStuff3 ();
00050     void IncrSortEigenStuff4 ();
00051     void IncrSortEigenStuffN ();
00052     void IncrSortEigenStuff  ();
00053 
00054 protected:
00055     int m_iSize;
00056     float** m_aafMat;
00057     float* m_afDiag;
00058     float* m_afSubd;
00059 
00060     // Householder reduction to tridiagonal form
00061     static void Tridiagonal2 (float** aafMat, float* afDiag,
00062         float* afSubd);
00063     static void Tridiagonal3 (float** aafMat, float* afDiag,
00064         float* afSubd);
00065     static void Tridiagonal4 (float** aafMat, float* afDiag,
00066         float* afSubd);
00067     static void TridiagonalN (int iSize, float** aafMat, float* afDiag,
00068         float* afSubd);
00069 
00070     // QL algorithm with implicit shifting, applies to tridiagonal matrices
00071     static bool QLAlgorithm (int iSize, float* afDiag, float* afSubd,
00072         float** aafMat);
00073 
00074     // sort eigenvalues from largest to smallest
00075     static void DecreasingSort (int iSize, float* afEigval,
00076         float** aafEigvec);
00077 
00078     // sort eigenvalues from smallest to largest
00079     static void IncreasingSort (int iSize, float* afEigval,
00080         float** aafEigvec);
00081 };
00082 
00083 //---------------------------------------------------------------------------
00084 inline float& Eigen::Matrix (int iRow, int iCol)
00085 {
00086     return m_aafMat[iRow][iCol];
00087 }
00088 //---------------------------------------------------------------------------
00089 inline float Eigen::GetEigenvalue (int i) const
00090 {
00091     return m_afDiag[i];
00092 }
00093 //---------------------------------------------------------------------------
00094 inline float Eigen::GetEigenvector (int iRow, int iCol) const
00095 {
00096     return m_aafMat[iRow][iCol];
00097 }
00098 //---------------------------------------------------------------------------
00099 inline float* Eigen::GetEigenvalue ()
00100 {
00101     return m_afDiag;
00102 }
00103 //---------------------------------------------------------------------------
00104 inline float** Eigen::GetEigenvector ()
00105 {
00106     return m_aafMat;
00107 }
00108 //---------------------------------------------------------------------------
00109 
00110 
00111 
00112 
00113 
00114 //---------------------------------------------------------------------------
00115 Eigen::Eigen (int iSize)
00116 {
00117     assert( iSize >= 2 );
00118     m_iSize = iSize;
00119 
00120     m_aafMat = new float*[m_iSize];
00121     for (int i = 0; i < m_iSize; i++)
00122         m_aafMat[i] = new float[m_iSize];
00123 
00124     m_afDiag = new float[m_iSize];
00125     m_afSubd = new float[m_iSize];
00126 }
00127 //---------------------------------------------------------------------------
00128 Eigen::~Eigen ()
00129 {
00130     delete[] m_afSubd;
00131     delete[] m_afDiag;
00132     for (int i = 0; i < m_iSize; i++)
00133         delete[] m_aafMat[i];
00134     delete[] m_aafMat;
00135 }
00136 //---------------------------------------------------------------------------
00137 void Eigen::Tridiagonal2 (float** m_aafMat, float* m_afDiag,
00138     float* m_afSubd)
00139 {
00140     // matrix is already tridiagonal
00141     m_afDiag[0] = m_aafMat[0][0];
00142     m_afDiag[1] = m_aafMat[1][1];
00143     m_afSubd[0] = m_aafMat[0][1];
00144     m_afSubd[1] = 0.0;
00145     m_aafMat[0][0] = 1.0;
00146     m_aafMat[0][1] = 0.0;
00147     m_aafMat[1][0] = 0.0;
00148     m_aafMat[1][1] = 1.0;
00149 }
00150 //---------------------------------------------------------------------------
00151 void Eigen::Tridiagonal3 (float** m_aafMat, float* m_afDiag,
00152     float* m_afSubd)
00153 {
00154     float fM00 = m_aafMat[0][0];
00155     float fM01 = m_aafMat[0][1];
00156     float fM02 = m_aafMat[0][2];
00157     float fM11 = m_aafMat[1][1];
00158     float fM12 = m_aafMat[1][2];
00159     float fM22 = m_aafMat[2][2];
00160 
00161     m_afDiag[0] = fM00;
00162     m_afSubd[2] = 0.0;
00163     if ( fM02 != 0.0 )
00164     {
00165         float fLength = Math::sqrt(fM01*fM01+fM02*fM02);
00166         float fInvLength = 1.0/fLength;
00167         fM01 *= fInvLength;
00168         fM02 *= fInvLength;
00169         float fQ = 2.0*fM01*fM12+fM02*(fM22-fM11);
00170         m_afDiag[1] = fM11+fM02*fQ;
00171         m_afDiag[2] = fM22-fM02*fQ;
00172         m_afSubd[0] = fLength;
00173         m_afSubd[1] = fM12-fM01*fQ;
00174         m_aafMat[0][0] = 1.0; m_aafMat[0][1] = 0.0;  m_aafMat[0][2] = 0.0;
00175         m_aafMat[1][0] = 0.0; m_aafMat[1][1] = fM01; m_aafMat[1][2] = fM02;
00176         m_aafMat[2][0] = 0.0; m_aafMat[2][1] = fM02; m_aafMat[2][2] = -fM01;
00177     }
00178     else
00179     {
00180         m_afDiag[1] = fM11;
00181         m_afDiag[2] = fM22;
00182         m_afSubd[0] = fM01;
00183         m_afSubd[1] = fM12;
00184         m_aafMat[0][0] = 1.0; m_aafMat[0][1] = 0.0; m_aafMat[0][2] = 0.0;
00185         m_aafMat[1][0] = 0.0; m_aafMat[1][1] = 1.0; m_aafMat[1][2] = 0.0;
00186         m_aafMat[2][0] = 0.0; m_aafMat[2][1] = 0.0; m_aafMat[2][2] = 1.0;
00187     }
00188 }
00189 //---------------------------------------------------------------------------
00190 void Eigen::Tridiagonal4 (float** m_aafMat, float* m_afDiag,
00191     float* m_afSubd)
00192 {
00193     // save matrix M
00194     float fM00 = m_aafMat[0][0];
00195     float fM01 = m_aafMat[0][1];
00196     float fM02 = m_aafMat[0][2];
00197     float fM03 = m_aafMat[0][3];
00198     float fM11 = m_aafMat[1][1];
00199     float fM12 = m_aafMat[1][2];
00200     float fM13 = m_aafMat[1][3];
00201     float fM22 = m_aafMat[2][2];
00202     float fM23 = m_aafMat[2][3];
00203     float fM33 = m_aafMat[3][3];
00204 
00205     m_afDiag[0] = fM00;
00206     m_afSubd[3] = 0.0;
00207 
00208     m_aafMat[0][0] = 1.0;
00209     m_aafMat[0][1] = 0.0;
00210     m_aafMat[0][2] = 0.0;
00211     m_aafMat[0][3] = 0.0;
00212     m_aafMat[1][0] = 0.0;
00213     m_aafMat[2][0] = 0.0;
00214     m_aafMat[3][0] = 0.0;
00215 
00216     float fLength, fInvLength;
00217 
00218     if ( fM02 != 0.0 || fM03 != 0.0 )
00219     {
00220         float fQ11, fQ12, fQ13;
00221         float fQ21, fQ22, fQ23;
00222         float fQ31, fQ32, fQ33;
00223 
00224         // build column Q1
00225         fLength = Math::sqrt(fM01*fM01 + fM02*fM02 + fM03*fM03);
00226         fInvLength = 1.0/fLength;
00227         fQ11 = fM01*fInvLength;
00228         fQ21 = fM02*fInvLength;
00229         fQ31 = fM03*fInvLength;
00230 
00231         m_afSubd[0] = fLength;
00232 
00233         // compute S*Q1
00234         float fV0 = fM11*fQ11+fM12*fQ21+fM13*fQ31;
00235         float fV1 = fM12*fQ11+fM22*fQ21+fM23*fQ31;
00236         float fV2 = fM13*fQ11+fM23*fQ21+fM33*fQ31;
00237 
00238         m_afDiag[1] = fQ11*fV0+fQ21*fV1+fQ31*fV2;
00239 
00240         // build column Q3 = Q1x(S*Q1)
00241         fQ13 = fQ21*fV2-fQ31*fV1;
00242         fQ23 = fQ31*fV0-fQ11*fV2;
00243         fQ33 = fQ11*fV1-fQ21*fV0;
00244         fLength = Math::sqrt(fQ13*fQ13+fQ23*fQ23+fQ33*fQ33);
00245         if ( fLength > 0.0 )
00246         {
00247             fInvLength = 1.0/fLength;
00248             fQ13 *= fInvLength;
00249             fQ23 *= fInvLength;
00250             fQ33 *= fInvLength;
00251 
00252             // build column Q2 = Q3xQ1
00253             fQ12 = fQ23*fQ31-fQ33*fQ21;
00254             fQ22 = fQ33*fQ11-fQ13*fQ31;
00255             fQ32 = fQ13*fQ21-fQ23*fQ11;
00256 
00257             fV0 = fQ12*fM11+fQ22*fM12+fQ32*fM13;
00258             fV1 = fQ12*fM12+fQ22*fM22+fQ32*fM23;
00259             fV2 = fQ12*fM13+fQ22*fM23+fQ32*fM33;
00260             m_afSubd[1] = fQ11*fV0+fQ21*fV1+fQ31*fV2;
00261             m_afDiag[2] = fQ12*fV0+fQ22*fV1+fQ32*fV2;
00262             m_afSubd[2] = fQ13*fV0+fQ23*fV1+fQ33*fV2;
00263 
00264             fV0 = fQ13*fM11+fQ23*fM12+fQ33*fM13;
00265             fV1 = fQ13*fM12+fQ23*fM22+fQ33*fM23;
00266             fV2 = fQ13*fM13+fQ23*fM23+fQ33*fM33;
00267             m_afDiag[3] = fQ13*fV0+fQ23*fV1+fQ33*fV2;
00268         }
00269         else
00270         {
00271             // S*Q1 parallel to Q1, choose any valid Q2 and Q3
00272             m_afSubd[1] = 0;
00273 
00274             fLength = fQ21*fQ21+fQ31*fQ31;
00275             if ( fLength > 0.0 )
00276             {
00277                 fInvLength = 1.0/fLength;
00278                 float fTmp = fQ11-1.0;
00279                 fQ12 = -fQ21;
00280                 fQ22 = 1.0+fTmp*fQ21*fQ21*fInvLength;
00281                 fQ32 = fTmp*fQ21*fQ31*fInvLength;
00282 
00283                 fQ13 = -fQ31;
00284                 fQ23 = fQ32;
00285                 fQ33 = 1.0+fTmp*fQ31*fQ31*fInvLength;
00286 
00287                 fV0 = fQ12*fM11+fQ22*fM12+fQ32*fM13;
00288                 fV1 = fQ12*fM12+fQ22*fM22+fQ32*fM23;
00289                 fV2 = fQ12*fM13+fQ22*fM23+fQ32*fM33;
00290                 m_afDiag[2] = fQ12*fV0+fQ22*fV1+fQ32*fV2;
00291                 m_afSubd[2] = fQ13*fV0+fQ23*fV1+fQ33*fV2;
00292 
00293                 fV0 = fQ13*fM11+fQ23*fM12+fQ33*fM13;
00294                 fV1 = fQ13*fM12+fQ23*fM22+fQ33*fM23;
00295                 fV2 = fQ13*fM13+fQ23*fM23+fQ33*fM33;
00296                 m_afDiag[3] = fQ13*fV0+fQ23*fV1+fQ33*fV2;
00297             }
00298             else
00299             {
00300                 // Q1 = (+-1,0,0)
00301                 fQ12 = 0.0; fQ22 = 1.0; fQ32 = 0.0;
00302                 fQ13 = 0.0; fQ23 = 0.0; fQ33 = 1.0;
00303 
00304                 m_afDiag[2] = fM22;
00305                 m_afDiag[3] = fM33;
00306                 m_afSubd[2] = fM23;
00307             }
00308         }
00309 
00310         m_aafMat[1][1] = fQ11; m_aafMat[1][2] = fQ12; m_aafMat[1][3] = fQ13;
00311         m_aafMat[2][1] = fQ21; m_aafMat[2][2] = fQ22; m_aafMat[2][3] = fQ23;
00312         m_aafMat[3][1] = fQ31; m_aafMat[3][2] = fQ32; m_aafMat[3][3] = fQ33;
00313     }
00314     else
00315     {
00316         m_afDiag[1] = fM11;
00317         m_afSubd[0] = fM01;
00318         m_aafMat[1][1] = 1.0;
00319         m_aafMat[2][1] = 0.0;
00320         m_aafMat[3][1] = 0.0;
00321 
00322         if ( fM13 != 0.0 )
00323         {
00324             fLength = Math::sqrt(fM12*fM12+fM13*fM13);
00325             fInvLength = 1.0/fLength;
00326             fM12 *= fInvLength;
00327             fM13 *= fInvLength;
00328             float fQ = 2.0*fM12*fM23+fM13*(fM33-fM22);
00329 
00330             m_afDiag[2] = fM22+fM13*fQ;
00331             m_afDiag[3] = fM33-fM13*fQ;
00332             m_afSubd[1] = fLength;
00333             m_afSubd[2] = fM23-fM12*fQ;
00334             m_aafMat[1][2] = 0.0;
00335             m_aafMat[1][3] = 0.0;
00336             m_aafMat[2][2] = fM12;
00337             m_aafMat[2][3] = fM13;
00338             m_aafMat[3][2] = fM13;
00339             m_aafMat[3][3] = -fM12;
00340         }
00341         else
00342         {
00343             m_afDiag[2] = fM22;
00344             m_afDiag[3] = fM33;
00345             m_afSubd[1] = fM12;
00346             m_afSubd[2] = fM23;
00347             m_aafMat[1][2] = 0.0;
00348             m_aafMat[1][3] = 0.0;
00349             m_aafMat[2][2] = 1.0;
00350             m_aafMat[2][3] = 0.0;
00351             m_aafMat[3][2] = 0.0;
00352             m_aafMat[3][3] = 1.0;
00353         }
00354     }
00355 }
00356 //---------------------------------------------------------------------------
00357 void Eigen::TridiagonalN (int iSize, float** m_aafMat,
00358     float* m_afDiag, float* m_afSubd)
00359 {
00360     int i0, i1, i2, i3;
00361 
00362     for (i0 = iSize-1, i3 = iSize-2; i0 >= 1; i0--, i3--)
00363     {
00364         float fH = 0.0, fScale = 0.0;
00365 
00366         if ( i3 > 0 )
00367         {
00368             for (i2 = 0; i2 <= i3; i2++)
00369                 fScale += Math::abs(m_aafMat[i0][i2]);
00370             if ( fScale == 0 )
00371             {
00372                 m_afSubd[i0] = m_aafMat[i0][i3];
00373             }
00374             else
00375             {
00376                 float fInvScale = 1.0/fScale;
00377                 for (i2 = 0; i2 <= i3; i2++)
00378                 {
00379                     m_aafMat[i0][i2] *= fInvScale;
00380                     fH += m_aafMat[i0][i2]*m_aafMat[i0][i2];
00381                 }
00382                 float fF = m_aafMat[i0][i3];
00383                 float fG = Math::sqrt(fH);
00384                 if ( fF > 0.0 )
00385                     fG = -fG;
00386                 m_afSubd[i0] = fScale*fG;
00387                 fH -= fF*fG;
00388                 m_aafMat[i0][i3] = fF-fG;
00389                 fF = 0.0;
00390                 float fInvH = 1.0/fH;
00391                 for (i1 = 0; i1 <= i3; i1++)
00392                 {
00393                     m_aafMat[i1][i0] = m_aafMat[i0][i1]*fInvH;
00394                     fG = 0.0;
00395                     for (i2 = 0; i2 <= i1; i2++)
00396                         fG += m_aafMat[i1][i2]*m_aafMat[i0][i2];
00397                     for (i2 = i1+1; i2 <= i3; i2++)
00398                         fG += m_aafMat[i2][i1]*m_aafMat[i0][i2];
00399                     m_afSubd[i1] = fG*fInvH;
00400                     fF += m_afSubd[i1]*m_aafMat[i0][i1];
00401                 }
00402                 float fHalfFdivH = 0.5*fF*fInvH;
00403                 for (i1 = 0; i1 <= i3; i1++)
00404                 {
00405                     fF = m_aafMat[i0][i1];
00406                     fG = m_afSubd[i1] - fHalfFdivH*fF;
00407                     m_afSubd[i1] = fG;
00408                     for (i2 = 0; i2 <= i1; i2++)
00409                     {
00410                         m_aafMat[i1][i2] -= fF*m_afSubd[i2] +
00411                             fG*m_aafMat[i0][i2];
00412                     }
00413                 }
00414             }
00415         }
00416         else
00417         {
00418             m_afSubd[i0] = m_aafMat[i0][i3];
00419         }
00420 
00421         m_afDiag[i0] = fH;
00422     }
00423 
00424     m_afDiag[0] = m_afSubd[0] = 0;
00425     for (i0 = 0, i3 = -1; i0 <= iSize-1; i0++, i3++)
00426     {
00427         if ( m_afDiag[i0] )
00428         {
00429             for (i1 = 0; i1 <= i3; i1++)
00430             {
00431                 float fSum = 0;
00432                 for (i2 = 0; i2 <= i3; i2++)
00433                     fSum += m_aafMat[i0][i2]*m_aafMat[i2][i1];
00434                 for (i2 = 0; i2 <= i3; i2++)
00435                     m_aafMat[i2][i1] -= fSum*m_aafMat[i2][i0];
00436             }
00437         }
00438         m_afDiag[i0] = m_aafMat[i0][i0];
00439         m_aafMat[i0][i0] = 1;
00440         for (i1 = 0; i1 <= i3; i1++)
00441             m_aafMat[i1][i0] = m_aafMat[i0][i1] = 0;
00442     }
00443 
00444     // re-ordering if Eigen::QLAlgorithm is used subsequently
00445     for (i0 = 1, i3 = 0; i0 < iSize; i0++, i3++)
00446         m_afSubd[i3] = m_afSubd[i0];
00447     m_afSubd[iSize-1] = 0;
00448 }
00449 //---------------------------------------------------------------------------
00450 bool Eigen::QLAlgorithm (int iSize, float* m_afDiag, float* m_afSubd,
00451     float** m_aafMat)
00452 {
00453     const int iMaxIter = 32;
00454 
00455     for (int i0 = 0; i0 < iSize; i0++)
00456     {
00457         int i1;
00458         for (i1 = 0; i1 < iMaxIter; i1++)
00459         {
00460             int i2;
00461             for (i2 = i0; i2 <= iSize-2; i2++)
00462             {
00463                 float fTmp =
00464                     Math::abs(m_afDiag[i2])+ Math::abs(m_afDiag[i2+1]);
00465                 if ( Math::abs(m_afSubd[i2]) + fTmp == fTmp )
00466                     break;
00467             }
00468             if ( i2 == i0 )
00469                 break;
00470 
00471             float fG = (m_afDiag[i0+1]-m_afDiag[i0])/(2.0*m_afSubd[i0]);
00472             float fR = Math::sqrt(fG*fG+1.0);
00473             if ( fG < 0.0 )
00474                 fG = m_afDiag[i2]-m_afDiag[i0]+m_afSubd[i0]/(fG-fR);
00475             else
00476                 fG = m_afDiag[i2]-m_afDiag[i0]+m_afSubd[i0]/(fG+fR);
00477             float fSin = 1.0, fCos = 1.0, fP = 0.0;
00478             for (int i3 = i2-1; i3 >= i0; i3--)
00479             {
00480                 float fF = fSin*m_afSubd[i3];
00481                 float fB = fCos*m_afSubd[i3];
00482                 if ( Math::abs(fF) >= Math::abs(fG) )
00483                 {
00484                     fCos = fG/fF;
00485                     fR = sqrt(fCos*fCos+1.0);
00486                     m_afSubd[i3+1] = fF*fR;
00487                     fSin = 1.0/fR;
00488                     fCos *= fSin;
00489                 }
00490                 else
00491                 {
00492                     fSin = fF/fG;
00493                     fR = Math::sqrt(fSin*fSin+1.0);
00494                     m_afSubd[i3+1] = fG*fR;
00495                     fCos = 1.0/fR;
00496                     fSin *= fCos;
00497                 }
00498                 fG = m_afDiag[i3+1]-fP;
00499                 fR = (m_afDiag[i3]-fG)*fSin+2.0*fB*fCos;
00500                 fP = fSin*fR;
00501                 m_afDiag[i3+1] = fG+fP;
00502                 fG = fCos*fR-fB;
00503 
00504                 for (int i4 = 0; i4 < iSize; i4++)
00505                 {
00506                     fF = m_aafMat[i4][i3+1];
00507                     m_aafMat[i4][i3+1] = fSin*m_aafMat[i4][i3]+fCos*fF;
00508                     m_aafMat[i4][i3] = fCos*m_aafMat[i4][i3]-fSin*fF;
00509                 }
00510             }
00511             m_afDiag[i0] -= fP;
00512             m_afSubd[i0] = fG;
00513             m_afSubd[i2] = 0.0;
00514         }
00515         if ( i1 == iMaxIter )
00516             return false;
00517     }
00518 
00519     return true;
00520 }
00521 //---------------------------------------------------------------------------
00522 void Eigen::DecreasingSort (int iSize, float* afEigval,
00523     float** aafEigvec)
00524 {
00525     // sort eigenvalues in decreasing order, e[0] >= ... >= e[iSize-1]
00526     for (int i0 = 0, i1; i0 <= iSize-2; i0++)
00527     {
00528         // locate maximum eigenvalue
00529         i1 = i0;
00530         float fMax = afEigval[i1];
00531         int i2;
00532         for (i2 = i0+1; i2 < iSize; i2++)
00533         {
00534             if ( afEigval[i2] > fMax )
00535             {
00536                 i1 = i2;
00537                 fMax = afEigval[i1];
00538             }
00539         }
00540 
00541         if ( i1 != i0 )
00542         {
00543             // swap eigenvalues
00544             afEigval[i1] = afEigval[i0];
00545             afEigval[i0] = fMax;
00546 
00547             // swap eigenvectors
00548             for (i2 = 0; i2 < iSize; i2++)
00549             {
00550                 float fTmp = aafEigvec[i2][i0];
00551                 aafEigvec[i2][i0] = aafEigvec[i2][i1];
00552                 aafEigvec[i2][i1] = fTmp;
00553             }
00554         }
00555     }
00556 }
00557 //---------------------------------------------------------------------------
00558 void Eigen::IncreasingSort (int iSize, float* afEigval,
00559     float** aafEigvec)
00560 {
00561     // sort eigenvalues in increasing order, e[0] <= ... <= e[iSize-1]
00562     for (int i0 = 0, i1; i0 <= iSize-2; i0++)
00563     {
00564         // locate minimum eigenvalue
00565         i1 = i0;
00566         float fMin = afEigval[i1];
00567         int i2;
00568         for (i2 = i0+1; i2 < iSize; i2++)
00569         {
00570             if ( afEigval[i2] < fMin )
00571             {
00572                 i1 = i2;
00573                 fMin = afEigval[i1];
00574             }
00575         }
00576 
00577         if ( i1 != i0 )
00578         {
00579             // swap eigenvalues
00580             afEigval[i1] = afEigval[i0];
00581             afEigval[i0] = fMin;
00582 
00583             // swap eigenvectors
00584             for (i2 = 0; i2 < iSize; i2++)
00585             {
00586                 float fTmp = aafEigvec[i2][i0];
00587                 aafEigvec[i2][i0] = aafEigvec[i2][i1];
00588                 aafEigvec[i2][i1] = fTmp;
00589             }
00590         }
00591     }
00592 }
00593 //---------------------------------------------------------------------------
00594 void Eigen::SetMatrix (float** aafMat)
00595 {
00596     for (int iRow = 0; iRow < m_iSize; iRow++)
00597     {
00598         for (int iCol = 0; iCol < m_iSize; iCol++)
00599             m_aafMat[iRow][iCol] = aafMat[iRow][iCol];
00600     }
00601 }
00602 //---------------------------------------------------------------------------
00603 void Eigen::EigenStuff2 ()
00604 {
00605     Tridiagonal2(m_aafMat,m_afDiag,m_afSubd);
00606     QLAlgorithm(m_iSize,m_afDiag,m_afSubd,m_aafMat);
00607 }
00608 //---------------------------------------------------------------------------
00609 void Eigen::EigenStuff3 ()
00610 {
00611     Tridiagonal3(m_aafMat,m_afDiag,m_afSubd);
00612     QLAlgorithm(m_iSize,m_afDiag,m_afSubd,m_aafMat);
00613 }
00614 //---------------------------------------------------------------------------
00615 void Eigen::EigenStuff4 ()
00616 {
00617     Tridiagonal4(m_aafMat,m_afDiag,m_afSubd);
00618     QLAlgorithm(m_iSize,m_afDiag,m_afSubd,m_aafMat);
00619 }
00620 //---------------------------------------------------------------------------
00621 void Eigen::EigenStuffN ()
00622 {
00623     TridiagonalN(m_iSize,m_aafMat,m_afDiag,m_afSubd);
00624     QLAlgorithm(m_iSize,m_afDiag,m_afSubd,m_aafMat);
00625 }
00626 //---------------------------------------------------------------------------
00627 void Eigen::EigenStuff ()
00628 {
00629     switch ( m_iSize )
00630     {
00631         case 2:
00632             Tridiagonal2(m_aafMat,m_afDiag,m_afSubd);
00633             break;
00634         case 3:
00635             Tridiagonal3(m_aafMat,m_afDiag,m_afSubd);
00636             break;
00637         case 4:
00638             Tridiagonal4(m_aafMat,m_afDiag,m_afSubd);
00639             break;
00640         default:
00641             TridiagonalN(m_iSize,m_aafMat,m_afDiag,m_afSubd);
00642             break;
00643     }
00644     QLAlgorithm(m_iSize,m_afDiag,m_afSubd,m_aafMat);
00645 }
00646 //---------------------------------------------------------------------------
00647 void Eigen::DecrSortEigenStuff2 ()
00648 {
00649     Tridiagonal2(m_aafMat,m_afDiag,m_afSubd);
00650     QLAlgorithm(m_iSize,m_afDiag,m_afSubd,m_aafMat);
00651     DecreasingSort(m_iSize,m_afDiag,m_aafMat);
00652 }
00653 //---------------------------------------------------------------------------
00654 void Eigen::DecrSortEigenStuff3 ()
00655 {
00656     Tridiagonal3(m_aafMat,m_afDiag,m_afSubd);
00657     QLAlgorithm(m_iSize,m_afDiag,m_afSubd,m_aafMat);
00658     DecreasingSort(m_iSize,m_afDiag,m_aafMat);
00659 }
00660 //---------------------------------------------------------------------------
00661 void Eigen::DecrSortEigenStuff4 ()
00662 {
00663     Tridiagonal4(m_aafMat,m_afDiag,m_afSubd);
00664     QLAlgorithm(m_iSize,m_afDiag,m_afSubd,m_aafMat);
00665     DecreasingSort(m_iSize,m_afDiag,m_aafMat);
00666 }
00667 //---------------------------------------------------------------------------
00668 void Eigen::DecrSortEigenStuffN ()
00669 {
00670     TridiagonalN(m_iSize,m_aafMat,m_afDiag,m_afSubd);
00671     QLAlgorithm(m_iSize,m_afDiag,m_afSubd,m_aafMat);
00672     DecreasingSort(m_iSize,m_afDiag,m_aafMat);
00673 }
00674 //---------------------------------------------------------------------------
00675 void Eigen::DecrSortEigenStuff ()
00676 {
00677     switch ( m_iSize )
00678     {
00679         case 2:
00680             Tridiagonal2(m_aafMat,m_afDiag,m_afSubd);
00681             break;
00682         case 3:
00683             Tridiagonal3(m_aafMat,m_afDiag,m_afSubd);
00684             break;
00685         case 4:
00686             Tridiagonal4(m_aafMat,m_afDiag,m_afSubd);
00687             break;
00688         default:
00689             TridiagonalN(m_iSize,m_aafMat,m_afDiag,m_afSubd);
00690             break;
00691     }
00692     QLAlgorithm(m_iSize,m_afDiag,m_afSubd,m_aafMat);
00693     DecreasingSort(m_iSize,m_afDiag,m_aafMat);
00694 }
00695 //---------------------------------------------------------------------------
00696 void Eigen::IncrSortEigenStuff2 ()
00697 {
00698     Tridiagonal2(m_aafMat,m_afDiag,m_afSubd);
00699     QLAlgorithm(m_iSize,m_afDiag,m_afSubd,m_aafMat);
00700     IncreasingSort(m_iSize,m_afDiag,m_aafMat);
00701 }
00702 //---------------------------------------------------------------------------
00703 void Eigen::IncrSortEigenStuff3 ()
00704 {
00705     Tridiagonal3(m_aafMat,m_afDiag,m_afSubd);
00706     QLAlgorithm(m_iSize,m_afDiag,m_afSubd,m_aafMat);
00707     IncreasingSort(m_iSize,m_afDiag,m_aafMat);
00708 }
00709 //---------------------------------------------------------------------------
00710 void Eigen::IncrSortEigenStuff4 ()
00711 {
00712     Tridiagonal4(m_aafMat,m_afDiag,m_afSubd);
00713     QLAlgorithm(m_iSize,m_afDiag,m_afSubd,m_aafMat);
00714     IncreasingSort(m_iSize,m_afDiag,m_aafMat);
00715 }
00716 //---------------------------------------------------------------------------
00717 void Eigen::IncrSortEigenStuffN ()
00718 {
00719     TridiagonalN(m_iSize,m_aafMat,m_afDiag,m_afSubd);
00720     QLAlgorithm(m_iSize,m_afDiag,m_afSubd,m_aafMat);
00721     IncreasingSort(m_iSize,m_afDiag,m_aafMat);
00722 }
00723 //---------------------------------------------------------------------------
00724 void Eigen::IncrSortEigenStuff ()
00725 {
00726     switch ( m_iSize )
00727     {
00728         case 2:
00729             Tridiagonal2(m_aafMat,m_afDiag,m_afSubd);
00730             break;
00731         case 3:
00732             Tridiagonal3(m_aafMat,m_afDiag,m_afSubd);
00733             break;
00734         case 4:
00735             Tridiagonal4(m_aafMat,m_afDiag,m_afSubd);
00736             break;
00737         default:
00738             TridiagonalN(m_iSize,m_aafMat,m_afDiag,m_afSubd);
00739             break;
00740     }
00741     QLAlgorithm(m_iSize,m_afDiag,m_afSubd,m_aafMat);
00742     IncreasingSort(m_iSize,m_afDiag,m_aafMat);
00743 }
00744 //---------------------------------------------------------------------------
00745 
00746 };
00747 
00748 
00749 /*
00750 #ifdef EIGEN_TEST
00751 
00752 int main ()
00753 {
00754     Eigen kES(3);
00755 
00756     kES.Matrix(0,0) = 2.0;  kES.Matrix(0,1) = 1.0;  kES.Matrix(0,2) = 1.0;
00757     kES.Matrix(1,0) = 1.0;  kES.Matrix(1,1) = 2.0;  kES.Matrix(1,2) = 1.0;
00758     kES.Matrix(2,0) = 1.0;  kES.Matrix(2,1) = 1.0;  kES.Matrix(2,2) = 2.0;
00759 
00760     kES.IncrSortEigenStuff3();
00761 
00762     cout.setf(ios::fixed);
00763 
00764     cout << "eigenvalues = " << endl;
00765     int iRow;
00766     for (iRow = 0; iRow < 3; iRow++)
00767         cout << kES.GetEigenvalue(iRow) << ' ';
00768     cout << endl;
00769 
00770     cout << "eigenvectors = " << endl;
00771     for (iRow = 0; iRow < 3; iRow++)
00772     {
00773         for (int iCol = 0; iCol < 3; iCol++)
00774             cout << kES.GetEigenvector(iRow,iCol) << ' ';
00775         cout << endl;
00776     }
00777 
00778     // eigenvalues =
00779     //    1.000000 1.000000 4.000000
00780     // eigenvectors =
00781     //    0.411953  0.704955 0.577350
00782     //    0.404533 -0.709239 0.577350
00783     //   -0.816485  0.004284 0.577350
00784 
00785     return 0;
00786 }
00787 
00788 #endif
00789 */
00790 
00791 #endif