Main Page   Modules   Namespace List   Class Hierarchy   Alphabetical List   Compound List   File List   Namespace Members   Compound Members   File Members   Related Pages  

Eigen.h

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

Generated on Mon Apr 7 15:28:54 2003 for GenericMathTemplateLibrary by doxygen1.2.14 written by Dimitri van Heesch, © 1997-2002