source: projectionDesigner/trunk/projdesigner/include/gmtl/Numerics/Eigen.h @ 217

Last change on this file since 217 was 4, checked in by Torben Dannhauer, 15 years ago
File size: 24.8 KB
Line 
1/************************************************************** ggt-head beg
2 *
3 * GGT: Generic Graphics Toolkit
4 *
5 * Original Authors:
6 *   Allen Bierbaum
7 *
8 * -----------------------------------------------------------------
9 * File:          Eigen.h,v
10 * Date modified: 2002/02/10 04:38:07
11 * Version:       1.4
12 * -----------------------------------------------------------------
13 *
14 *********************************************************** ggt-head end */
15/*************************************************************** ggt-cpr beg
16*
17* GGT: The Generic Graphics Toolkit
18* Copyright (C) 2001,2002 Allen Bierbaum
19*
20* This library is free software; you can redistribute it and/or
21* modify it under the terms of the GNU Lesser General Public
22* License as published by the Free Software Foundation; either
23* version 2.1 of the License, or (at your option) any later version.
24*
25* This library is distributed in the hope that it will be useful,
26* but WITHOUT ANY WARRANTY; without even the implied warranty of
27* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
28* Lesser General Public License for more details.
29*
30* You should have received a copy of the GNU Lesser General Public
31* License along with this library; if not, write to the Free Software
32* Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
33*
34 ************************************************************ ggt-cpr end */
35//
36// XXX: Based on source code from: Magic Software, Inc.
37//
38#ifndef _EIGEN_H
39#define _EIGEN_H
40
41#include <gmtl/gmtlConfig.h>
42
43namespace gmtl
44{
45
46class Eigen
47{
48public:
49    Eigen (int iSize);
50    ~Eigen ();
51
52    // set the matrix for eigensolving
53    float& Matrix (int iRow, int iCol);
54    void SetMatrix (float** aafMat);
55
56    // get the results of eigensolving (eigenvectors are columns of matrix)
57    float GetEigenvalue (int i) const;
58    float GetEigenvector (int iRow, int iCol) const;
59    float* GetEigenvalue ();
60    float** GetEigenvector ();
61
62    // solve eigensystem
63    void EigenStuff2 ();
64    void EigenStuff3 ();
65    void EigenStuff4 ();
66    void EigenStuffN ();
67    void EigenStuff  ();
68
69    // solve eigensystem, use decreasing sort on eigenvalues
70    void DecrSortEigenStuff2 ();
71    void DecrSortEigenStuff3 ();
72    void DecrSortEigenStuff4 ();
73    void DecrSortEigenStuffN ();
74    void DecrSortEigenStuff  ();
75
76    // solve eigensystem, use increasing sort on eigenvalues
77    void IncrSortEigenStuff2 ();
78    void IncrSortEigenStuff3 ();
79    void IncrSortEigenStuff4 ();
80    void IncrSortEigenStuffN ();
81    void IncrSortEigenStuff  ();
82
83protected:
84    int m_iSize;
85    float** m_aafMat;
86    float* m_afDiag;
87    float* m_afSubd;
88
89    // Householder reduction to tridiagonal form
90    static void Tridiagonal2 (float** aafMat, float* afDiag,
91        float* afSubd);
92    static void Tridiagonal3 (float** aafMat, float* afDiag,
93        float* afSubd);
94    static void Tridiagonal4 (float** aafMat, float* afDiag,
95        float* afSubd);
96    static void TridiagonalN (int iSize, float** aafMat, float* afDiag,
97        float* afSubd);
98
99    // QL algorithm with implicit shifting, applies to tridiagonal matrices
100    static bool QLAlgorithm (int iSize, float* afDiag, float* afSubd,
101        float** aafMat);
102
103    // sort eigenvalues from largest to smallest
104    static void DecreasingSort (int iSize, float* afEigval,
105        float** aafEigvec);
106
107    // sort eigenvalues from smallest to largest
108    static void IncreasingSort (int iSize, float* afEigval,
109        float** aafEigvec);
110};
111
112//---------------------------------------------------------------------------
113inline float& Eigen::Matrix (int iRow, int iCol)
114{
115    return m_aafMat[iRow][iCol];
116}
117//---------------------------------------------------------------------------
118inline float Eigen::GetEigenvalue (int i) const
119{
120    return m_afDiag[i];
121}
122//---------------------------------------------------------------------------
123inline float Eigen::GetEigenvector (int iRow, int iCol) const
124{
125    return m_aafMat[iRow][iCol];
126}
127//---------------------------------------------------------------------------
128inline float* Eigen::GetEigenvalue ()
129{
130    return m_afDiag;
131}
132//---------------------------------------------------------------------------
133inline float** Eigen::GetEigenvector ()
134{
135    return m_aafMat;
136}
137//---------------------------------------------------------------------------
138
139
140
141
142
143//---------------------------------------------------------------------------
144Eigen::Eigen (int iSize)
145{
146    assert( iSize >= 2 );
147    m_iSize = iSize;
148
149    m_aafMat = new float*[m_iSize];
150    for (int i = 0; i < m_iSize; i++)
151        m_aafMat[i] = new float[m_iSize];
152
153    m_afDiag = new float[m_iSize];
154    m_afSubd = new float[m_iSize];
155}
156//---------------------------------------------------------------------------
157Eigen::~Eigen ()
158{
159    delete[] m_afSubd;
160    delete[] m_afDiag;
161    for (int i = 0; i < m_iSize; i++)
162        delete[] m_aafMat[i];
163    delete[] m_aafMat;
164}
165//---------------------------------------------------------------------------
166void Eigen::Tridiagonal2 (float** m_aafMat, float* m_afDiag,
167    float* m_afSubd)
168{
169    // matrix is already tridiagonal
170    m_afDiag[0] = m_aafMat[0][0];
171    m_afDiag[1] = m_aafMat[1][1];
172    m_afSubd[0] = m_aafMat[0][1];
173    m_afSubd[1] = 0.0;
174    m_aafMat[0][0] = 1.0;
175    m_aafMat[0][1] = 0.0;
176    m_aafMat[1][0] = 0.0;
177    m_aafMat[1][1] = 1.0;
178}
179//---------------------------------------------------------------------------
180void Eigen::Tridiagonal3 (float** m_aafMat, float* m_afDiag,
181    float* m_afSubd)
182{
183    float fM00 = m_aafMat[0][0];
184    float fM01 = m_aafMat[0][1];
185    float fM02 = m_aafMat[0][2];
186    float fM11 = m_aafMat[1][1];
187    float fM12 = m_aafMat[1][2];
188    float fM22 = m_aafMat[2][2];
189
190    m_afDiag[0] = fM00;
191    m_afSubd[2] = 0.0;
192    if ( fM02 != 0.0 )
193    {
194        float fLength = Math::sqrt(fM01*fM01+fM02*fM02);
195        float fInvLength = 1.0/fLength;
196        fM01 *= fInvLength;
197        fM02 *= fInvLength;
198        float fQ = 2.0*fM01*fM12+fM02*(fM22-fM11);
199        m_afDiag[1] = fM11+fM02*fQ;
200        m_afDiag[2] = fM22-fM02*fQ;
201        m_afSubd[0] = fLength;
202        m_afSubd[1] = fM12-fM01*fQ;
203        m_aafMat[0][0] = 1.0; m_aafMat[0][1] = 0.0;  m_aafMat[0][2] = 0.0;
204        m_aafMat[1][0] = 0.0; m_aafMat[1][1] = fM01; m_aafMat[1][2] = fM02;
205        m_aafMat[2][0] = 0.0; m_aafMat[2][1] = fM02; m_aafMat[2][2] = -fM01;
206    }
207    else
208    {
209        m_afDiag[1] = fM11;
210        m_afDiag[2] = fM22;
211        m_afSubd[0] = fM01;
212        m_afSubd[1] = fM12;
213        m_aafMat[0][0] = 1.0; m_aafMat[0][1] = 0.0; m_aafMat[0][2] = 0.0;
214        m_aafMat[1][0] = 0.0; m_aafMat[1][1] = 1.0; m_aafMat[1][2] = 0.0;
215        m_aafMat[2][0] = 0.0; m_aafMat[2][1] = 0.0; m_aafMat[2][2] = 1.0;
216    }
217}
218//---------------------------------------------------------------------------
219void Eigen::Tridiagonal4 (float** m_aafMat, float* m_afDiag,
220    float* m_afSubd)
221{
222    // save matrix M
223    float fM00 = m_aafMat[0][0];
224    float fM01 = m_aafMat[0][1];
225    float fM02 = m_aafMat[0][2];
226    float fM03 = m_aafMat[0][3];
227    float fM11 = m_aafMat[1][1];
228    float fM12 = m_aafMat[1][2];
229    float fM13 = m_aafMat[1][3];
230    float fM22 = m_aafMat[2][2];
231    float fM23 = m_aafMat[2][3];
232    float fM33 = m_aafMat[3][3];
233
234    m_afDiag[0] = fM00;
235    m_afSubd[3] = 0.0;
236
237    m_aafMat[0][0] = 1.0;
238    m_aafMat[0][1] = 0.0;
239    m_aafMat[0][2] = 0.0;
240    m_aafMat[0][3] = 0.0;
241    m_aafMat[1][0] = 0.0;
242    m_aafMat[2][0] = 0.0;
243    m_aafMat[3][0] = 0.0;
244
245    float fLength, fInvLength;
246
247    if ( fM02 != 0.0 || fM03 != 0.0 )
248    {
249        float fQ11, fQ12, fQ13;
250        float fQ21, fQ22, fQ23;
251        float fQ31, fQ32, fQ33;
252
253        // build column Q1
254        fLength = Math::sqrt(fM01*fM01 + fM02*fM02 + fM03*fM03);
255        fInvLength = 1.0/fLength;
256        fQ11 = fM01*fInvLength;
257        fQ21 = fM02*fInvLength;
258        fQ31 = fM03*fInvLength;
259
260        m_afSubd[0] = fLength;
261
262        // compute S*Q1
263        float fV0 = fM11*fQ11+fM12*fQ21+fM13*fQ31;
264        float fV1 = fM12*fQ11+fM22*fQ21+fM23*fQ31;
265        float fV2 = fM13*fQ11+fM23*fQ21+fM33*fQ31;
266
267        m_afDiag[1] = fQ11*fV0+fQ21*fV1+fQ31*fV2;
268
269        // build column Q3 = Q1x(S*Q1)
270        fQ13 = fQ21*fV2-fQ31*fV1;
271        fQ23 = fQ31*fV0-fQ11*fV2;
272        fQ33 = fQ11*fV1-fQ21*fV0;
273        fLength = Math::sqrt(fQ13*fQ13+fQ23*fQ23+fQ33*fQ33);
274        if ( fLength > 0.0 )
275        {
276            fInvLength = 1.0/fLength;
277            fQ13 *= fInvLength;
278            fQ23 *= fInvLength;
279            fQ33 *= fInvLength;
280
281            // build column Q2 = Q3xQ1
282            fQ12 = fQ23*fQ31-fQ33*fQ21;
283            fQ22 = fQ33*fQ11-fQ13*fQ31;
284            fQ32 = fQ13*fQ21-fQ23*fQ11;
285
286            fV0 = fQ12*fM11+fQ22*fM12+fQ32*fM13;
287            fV1 = fQ12*fM12+fQ22*fM22+fQ32*fM23;
288            fV2 = fQ12*fM13+fQ22*fM23+fQ32*fM33;
289            m_afSubd[1] = fQ11*fV0+fQ21*fV1+fQ31*fV2;
290            m_afDiag[2] = fQ12*fV0+fQ22*fV1+fQ32*fV2;
291            m_afSubd[2] = fQ13*fV0+fQ23*fV1+fQ33*fV2;
292
293            fV0 = fQ13*fM11+fQ23*fM12+fQ33*fM13;
294            fV1 = fQ13*fM12+fQ23*fM22+fQ33*fM23;
295            fV2 = fQ13*fM13+fQ23*fM23+fQ33*fM33;
296            m_afDiag[3] = fQ13*fV0+fQ23*fV1+fQ33*fV2;
297        }
298        else
299        {
300            // S*Q1 parallel to Q1, choose any valid Q2 and Q3
301            m_afSubd[1] = 0;
302
303            fLength = fQ21*fQ21+fQ31*fQ31;
304            if ( fLength > 0.0 )
305            {
306                fInvLength = 1.0/fLength;
307                float fTmp = fQ11-1.0;
308                fQ12 = -fQ21;
309                fQ22 = 1.0+fTmp*fQ21*fQ21*fInvLength;
310                fQ32 = fTmp*fQ21*fQ31*fInvLength;
311
312                fQ13 = -fQ31;
313                fQ23 = fQ32;
314                fQ33 = 1.0+fTmp*fQ31*fQ31*fInvLength;
315
316                fV0 = fQ12*fM11+fQ22*fM12+fQ32*fM13;
317                fV1 = fQ12*fM12+fQ22*fM22+fQ32*fM23;
318                fV2 = fQ12*fM13+fQ22*fM23+fQ32*fM33;
319                m_afDiag[2] = fQ12*fV0+fQ22*fV1+fQ32*fV2;
320                m_afSubd[2] = fQ13*fV0+fQ23*fV1+fQ33*fV2;
321
322                fV0 = fQ13*fM11+fQ23*fM12+fQ33*fM13;
323                fV1 = fQ13*fM12+fQ23*fM22+fQ33*fM23;
324                fV2 = fQ13*fM13+fQ23*fM23+fQ33*fM33;
325                m_afDiag[3] = fQ13*fV0+fQ23*fV1+fQ33*fV2;
326            }
327            else
328            {
329                // Q1 = (+-1,0,0)
330                fQ12 = 0.0; fQ22 = 1.0; fQ32 = 0.0;
331                fQ13 = 0.0; fQ23 = 0.0; fQ33 = 1.0;
332
333                m_afDiag[2] = fM22;
334                m_afDiag[3] = fM33;
335                m_afSubd[2] = fM23;
336            }
337        }
338
339        m_aafMat[1][1] = fQ11; m_aafMat[1][2] = fQ12; m_aafMat[1][3] = fQ13;
340        m_aafMat[2][1] = fQ21; m_aafMat[2][2] = fQ22; m_aafMat[2][3] = fQ23;
341        m_aafMat[3][1] = fQ31; m_aafMat[3][2] = fQ32; m_aafMat[3][3] = fQ33;
342    }
343    else
344    {
345        m_afDiag[1] = fM11;
346        m_afSubd[0] = fM01;
347        m_aafMat[1][1] = 1.0;
348        m_aafMat[2][1] = 0.0;
349        m_aafMat[3][1] = 0.0;
350
351        if ( fM13 != 0.0 )
352        {
353            fLength = Math::sqrt(fM12*fM12+fM13*fM13);
354            fInvLength = 1.0/fLength;
355            fM12 *= fInvLength;
356            fM13 *= fInvLength;
357            float fQ = 2.0*fM12*fM23+fM13*(fM33-fM22);
358
359            m_afDiag[2] = fM22+fM13*fQ;
360            m_afDiag[3] = fM33-fM13*fQ;
361            m_afSubd[1] = fLength;
362            m_afSubd[2] = fM23-fM12*fQ;
363            m_aafMat[1][2] = 0.0;
364            m_aafMat[1][3] = 0.0;
365            m_aafMat[2][2] = fM12;
366            m_aafMat[2][3] = fM13;
367            m_aafMat[3][2] = fM13;
368            m_aafMat[3][3] = -fM12;
369        }
370        else
371        {
372            m_afDiag[2] = fM22;
373            m_afDiag[3] = fM33;
374            m_afSubd[1] = fM12;
375            m_afSubd[2] = fM23;
376            m_aafMat[1][2] = 0.0;
377            m_aafMat[1][3] = 0.0;
378            m_aafMat[2][2] = 1.0;
379            m_aafMat[2][3] = 0.0;
380            m_aafMat[3][2] = 0.0;
381            m_aafMat[3][3] = 1.0;
382        }
383    }
384}
385//---------------------------------------------------------------------------
386void Eigen::TridiagonalN (int iSize, float** m_aafMat,
387    float* m_afDiag, float* m_afSubd)
388{
389    int i0, i1, i2, i3;
390
391    for (i0 = iSize-1, i3 = iSize-2; i0 >= 1; i0--, i3--)
392    {
393        float fH = 0.0, fScale = 0.0;
394
395        if ( i3 > 0 )
396        {
397            for (i2 = 0; i2 <= i3; i2++)
398                fScale += Math::abs(m_aafMat[i0][i2]);
399            if ( fScale == 0 )
400            {
401                m_afSubd[i0] = m_aafMat[i0][i3];
402            }
403            else
404            {
405                float fInvScale = 1.0/fScale;
406                for (i2 = 0; i2 <= i3; i2++)
407                {
408                    m_aafMat[i0][i2] *= fInvScale;
409                    fH += m_aafMat[i0][i2]*m_aafMat[i0][i2];
410                }
411                float fF = m_aafMat[i0][i3];
412                float fG = Math::sqrt(fH);
413                if ( fF > 0.0 )
414                    fG = -fG;
415                m_afSubd[i0] = fScale*fG;
416                fH -= fF*fG;
417                m_aafMat[i0][i3] = fF-fG;
418                fF = 0.0;
419                float fInvH = 1.0/fH;
420                for (i1 = 0; i1 <= i3; i1++)
421                {
422                    m_aafMat[i1][i0] = m_aafMat[i0][i1]*fInvH;
423                    fG = 0.0;
424                    for (i2 = 0; i2 <= i1; i2++)
425                        fG += m_aafMat[i1][i2]*m_aafMat[i0][i2];
426                    for (i2 = i1+1; i2 <= i3; i2++)
427                        fG += m_aafMat[i2][i1]*m_aafMat[i0][i2];
428                    m_afSubd[i1] = fG*fInvH;
429                    fF += m_afSubd[i1]*m_aafMat[i0][i1];
430                }
431                float fHalfFdivH = 0.5*fF*fInvH;
432                for (i1 = 0; i1 <= i3; i1++)
433                {
434                    fF = m_aafMat[i0][i1];
435                    fG = m_afSubd[i1] - fHalfFdivH*fF;
436                    m_afSubd[i1] = fG;
437                    for (i2 = 0; i2 <= i1; i2++)
438                    {
439                        m_aafMat[i1][i2] -= fF*m_afSubd[i2] +
440                            fG*m_aafMat[i0][i2];
441                    }
442                }
443            }
444        }
445        else
446        {
447            m_afSubd[i0] = m_aafMat[i0][i3];
448        }
449
450        m_afDiag[i0] = fH;
451    }
452
453    m_afDiag[0] = m_afSubd[0] = 0;
454    for (i0 = 0, i3 = -1; i0 <= iSize-1; i0++, i3++)
455    {
456        if ( m_afDiag[i0] )
457        {
458            for (i1 = 0; i1 <= i3; i1++)
459            {
460                float fSum = 0;
461                for (i2 = 0; i2 <= i3; i2++)
462                    fSum += m_aafMat[i0][i2]*m_aafMat[i2][i1];
463                for (i2 = 0; i2 <= i3; i2++)
464                    m_aafMat[i2][i1] -= fSum*m_aafMat[i2][i0];
465            }
466        }
467        m_afDiag[i0] = m_aafMat[i0][i0];
468        m_aafMat[i0][i0] = 1;
469        for (i1 = 0; i1 <= i3; i1++)
470            m_aafMat[i1][i0] = m_aafMat[i0][i1] = 0;
471    }
472
473    // re-ordering if Eigen::QLAlgorithm is used subsequently
474    for (i0 = 1, i3 = 0; i0 < iSize; i0++, i3++)
475        m_afSubd[i3] = m_afSubd[i0];
476    m_afSubd[iSize-1] = 0;
477}
478//---------------------------------------------------------------------------
479bool Eigen::QLAlgorithm (int iSize, float* m_afDiag, float* m_afSubd,
480    float** m_aafMat)
481{
482    const int iMaxIter = 32;
483
484    for (int i0 = 0; i0 < iSize; i0++)
485    {
486        int i1;
487        for (i1 = 0; i1 < iMaxIter; i1++)
488        {
489            int i2;
490            for (i2 = i0; i2 <= iSize-2; i2++)
491            {
492                float fTmp =
493                    Math::abs(m_afDiag[i2])+ Math::abs(m_afDiag[i2+1]);
494                if ( Math::abs(m_afSubd[i2]) + fTmp == fTmp )
495                    break;
496            }
497            if ( i2 == i0 )
498                break;
499
500            float fG = (m_afDiag[i0+1]-m_afDiag[i0])/(2.0*m_afSubd[i0]);
501            float fR = Math::sqrt(fG*fG+1.0);
502            if ( fG < 0.0 )
503                fG = m_afDiag[i2]-m_afDiag[i0]+m_afSubd[i0]/(fG-fR);
504            else
505                fG = m_afDiag[i2]-m_afDiag[i0]+m_afSubd[i0]/(fG+fR);
506            float fSin = 1.0, fCos = 1.0, fP = 0.0;
507            for (int i3 = i2-1; i3 >= i0; i3--)
508            {
509                float fF = fSin*m_afSubd[i3];
510                float fB = fCos*m_afSubd[i3];
511                if ( Math::abs(fF) >= Math::abs(fG) )
512                {
513                    fCos = fG/fF;
514                    fR = sqrt(fCos*fCos+1.0);
515                    m_afSubd[i3+1] = fF*fR;
516                    fSin = 1.0/fR;
517                    fCos *= fSin;
518                }
519                else
520                {
521                    fSin = fF/fG;
522                    fR = Math::sqrt(fSin*fSin+1.0);
523                    m_afSubd[i3+1] = fG*fR;
524                    fCos = 1.0/fR;
525                    fSin *= fCos;
526                }
527                fG = m_afDiag[i3+1]-fP;
528                fR = (m_afDiag[i3]-fG)*fSin+2.0*fB*fCos;
529                fP = fSin*fR;
530                m_afDiag[i3+1] = fG+fP;
531                fG = fCos*fR-fB;
532
533                for (int i4 = 0; i4 < iSize; i4++)
534                {
535                    fF = m_aafMat[i4][i3+1];
536                    m_aafMat[i4][i3+1] = fSin*m_aafMat[i4][i3]+fCos*fF;
537                    m_aafMat[i4][i3] = fCos*m_aafMat[i4][i3]-fSin*fF;
538                }
539            }
540            m_afDiag[i0] -= fP;
541            m_afSubd[i0] = fG;
542            m_afSubd[i2] = 0.0;
543        }
544        if ( i1 == iMaxIter )
545            return false;
546    }
547
548    return true;
549}
550//---------------------------------------------------------------------------
551void Eigen::DecreasingSort (int iSize, float* afEigval,
552    float** aafEigvec)
553{
554    // sort eigenvalues in decreasing order, e[0] >= ... >= e[iSize-1]
555    for (int i0 = 0, i1; i0 <= iSize-2; i0++)
556    {
557        // locate maximum eigenvalue
558        i1 = i0;
559        float fMax = afEigval[i1];
560        int i2;
561        for (i2 = i0+1; i2 < iSize; i2++)
562        {
563            if ( afEigval[i2] > fMax )
564            {
565                i1 = i2;
566                fMax = afEigval[i1];
567            }
568        }
569
570        if ( i1 != i0 )
571        {
572            // swap eigenvalues
573            afEigval[i1] = afEigval[i0];
574            afEigval[i0] = fMax;
575
576            // swap eigenvectors
577            for (i2 = 0; i2 < iSize; i2++)
578            {
579                float fTmp = aafEigvec[i2][i0];
580                aafEigvec[i2][i0] = aafEigvec[i2][i1];
581                aafEigvec[i2][i1] = fTmp;
582            }
583        }
584    }
585}
586//---------------------------------------------------------------------------
587void Eigen::IncreasingSort (int iSize, float* afEigval,
588    float** aafEigvec)
589{
590    // sort eigenvalues in increasing order, e[0] <= ... <= e[iSize-1]
591    for (int i0 = 0, i1; i0 <= iSize-2; i0++)
592    {
593        // locate minimum eigenvalue
594        i1 = i0;
595        float fMin = afEigval[i1];
596        int i2;
597        for (i2 = i0+1; i2 < iSize; i2++)
598        {
599            if ( afEigval[i2] < fMin )
600            {
601                i1 = i2;
602                fMin = afEigval[i1];
603            }
604        }
605
606        if ( i1 != i0 )
607        {
608            // swap eigenvalues
609            afEigval[i1] = afEigval[i0];
610            afEigval[i0] = fMin;
611
612            // swap eigenvectors
613            for (i2 = 0; i2 < iSize; i2++)
614            {
615                float fTmp = aafEigvec[i2][i0];
616                aafEigvec[i2][i0] = aafEigvec[i2][i1];
617                aafEigvec[i2][i1] = fTmp;
618            }
619        }
620    }
621}
622//---------------------------------------------------------------------------
623void Eigen::SetMatrix (float** aafMat)
624{
625    for (int iRow = 0; iRow < m_iSize; iRow++)
626    {
627        for (int iCol = 0; iCol < m_iSize; iCol++)
628            m_aafMat[iRow][iCol] = aafMat[iRow][iCol];
629    }
630}
631//---------------------------------------------------------------------------
632void Eigen::EigenStuff2 ()
633{
634    Tridiagonal2(m_aafMat,m_afDiag,m_afSubd);
635    QLAlgorithm(m_iSize,m_afDiag,m_afSubd,m_aafMat);
636}
637//---------------------------------------------------------------------------
638void Eigen::EigenStuff3 ()
639{
640    Tridiagonal3(m_aafMat,m_afDiag,m_afSubd);
641    QLAlgorithm(m_iSize,m_afDiag,m_afSubd,m_aafMat);
642}
643//---------------------------------------------------------------------------
644void Eigen::EigenStuff4 ()
645{
646    Tridiagonal4(m_aafMat,m_afDiag,m_afSubd);
647    QLAlgorithm(m_iSize,m_afDiag,m_afSubd,m_aafMat);
648}
649//---------------------------------------------------------------------------
650void Eigen::EigenStuffN ()
651{
652    TridiagonalN(m_iSize,m_aafMat,m_afDiag,m_afSubd);
653    QLAlgorithm(m_iSize,m_afDiag,m_afSubd,m_aafMat);
654}
655//---------------------------------------------------------------------------
656void Eigen::EigenStuff ()
657{
658    switch ( m_iSize )
659    {
660        case 2:
661            Tridiagonal2(m_aafMat,m_afDiag,m_afSubd);
662            break;
663        case 3:
664            Tridiagonal3(m_aafMat,m_afDiag,m_afSubd);
665            break;
666        case 4:
667            Tridiagonal4(m_aafMat,m_afDiag,m_afSubd);
668            break;
669        default:
670            TridiagonalN(m_iSize,m_aafMat,m_afDiag,m_afSubd);
671            break;
672    }
673    QLAlgorithm(m_iSize,m_afDiag,m_afSubd,m_aafMat);
674}
675//---------------------------------------------------------------------------
676void Eigen::DecrSortEigenStuff2 ()
677{
678    Tridiagonal2(m_aafMat,m_afDiag,m_afSubd);
679    QLAlgorithm(m_iSize,m_afDiag,m_afSubd,m_aafMat);
680    DecreasingSort(m_iSize,m_afDiag,m_aafMat);
681}
682//---------------------------------------------------------------------------
683void Eigen::DecrSortEigenStuff3 ()
684{
685    Tridiagonal3(m_aafMat,m_afDiag,m_afSubd);
686    QLAlgorithm(m_iSize,m_afDiag,m_afSubd,m_aafMat);
687    DecreasingSort(m_iSize,m_afDiag,m_aafMat);
688}
689//---------------------------------------------------------------------------
690void Eigen::DecrSortEigenStuff4 ()
691{
692    Tridiagonal4(m_aafMat,m_afDiag,m_afSubd);
693    QLAlgorithm(m_iSize,m_afDiag,m_afSubd,m_aafMat);
694    DecreasingSort(m_iSize,m_afDiag,m_aafMat);
695}
696//---------------------------------------------------------------------------
697void Eigen::DecrSortEigenStuffN ()
698{
699    TridiagonalN(m_iSize,m_aafMat,m_afDiag,m_afSubd);
700    QLAlgorithm(m_iSize,m_afDiag,m_afSubd,m_aafMat);
701    DecreasingSort(m_iSize,m_afDiag,m_aafMat);
702}
703//---------------------------------------------------------------------------
704void Eigen::DecrSortEigenStuff ()
705{
706    switch ( m_iSize )
707    {
708        case 2:
709            Tridiagonal2(m_aafMat,m_afDiag,m_afSubd);
710            break;
711        case 3:
712            Tridiagonal3(m_aafMat,m_afDiag,m_afSubd);
713            break;
714        case 4:
715            Tridiagonal4(m_aafMat,m_afDiag,m_afSubd);
716            break;
717        default:
718            TridiagonalN(m_iSize,m_aafMat,m_afDiag,m_afSubd);
719            break;
720    }
721    QLAlgorithm(m_iSize,m_afDiag,m_afSubd,m_aafMat);
722    DecreasingSort(m_iSize,m_afDiag,m_aafMat);
723}
724//---------------------------------------------------------------------------
725void Eigen::IncrSortEigenStuff2 ()
726{
727    Tridiagonal2(m_aafMat,m_afDiag,m_afSubd);
728    QLAlgorithm(m_iSize,m_afDiag,m_afSubd,m_aafMat);
729    IncreasingSort(m_iSize,m_afDiag,m_aafMat);
730}
731//---------------------------------------------------------------------------
732void Eigen::IncrSortEigenStuff3 ()
733{
734    Tridiagonal3(m_aafMat,m_afDiag,m_afSubd);
735    QLAlgorithm(m_iSize,m_afDiag,m_afSubd,m_aafMat);
736    IncreasingSort(m_iSize,m_afDiag,m_aafMat);
737}
738//---------------------------------------------------------------------------
739void Eigen::IncrSortEigenStuff4 ()
740{
741    Tridiagonal4(m_aafMat,m_afDiag,m_afSubd);
742    QLAlgorithm(m_iSize,m_afDiag,m_afSubd,m_aafMat);
743    IncreasingSort(m_iSize,m_afDiag,m_aafMat);
744}
745//---------------------------------------------------------------------------
746void Eigen::IncrSortEigenStuffN ()
747{
748    TridiagonalN(m_iSize,m_aafMat,m_afDiag,m_afSubd);
749    QLAlgorithm(m_iSize,m_afDiag,m_afSubd,m_aafMat);
750    IncreasingSort(m_iSize,m_afDiag,m_aafMat);
751}
752//---------------------------------------------------------------------------
753void Eigen::IncrSortEigenStuff ()
754{
755    switch ( m_iSize )
756    {
757        case 2:
758            Tridiagonal2(m_aafMat,m_afDiag,m_afSubd);
759            break;
760        case 3:
761            Tridiagonal3(m_aafMat,m_afDiag,m_afSubd);
762            break;
763        case 4:
764            Tridiagonal4(m_aafMat,m_afDiag,m_afSubd);
765            break;
766        default:
767            TridiagonalN(m_iSize,m_aafMat,m_afDiag,m_afSubd);
768            break;
769    }
770    QLAlgorithm(m_iSize,m_afDiag,m_afSubd,m_aafMat);
771    IncreasingSort(m_iSize,m_afDiag,m_aafMat);
772}
773//---------------------------------------------------------------------------
774
775};
776
777
778/*
779#ifdef EIGEN_TEST
780
781int main ()
782{
783    Eigen kES(3);
784
785    kES.Matrix(0,0) = 2.0;  kES.Matrix(0,1) = 1.0;  kES.Matrix(0,2) = 1.0;
786    kES.Matrix(1,0) = 1.0;  kES.Matrix(1,1) = 2.0;  kES.Matrix(1,2) = 1.0;
787    kES.Matrix(2,0) = 1.0;  kES.Matrix(2,1) = 1.0;  kES.Matrix(2,2) = 2.0;
788
789    kES.IncrSortEigenStuff3();
790
791    cout.setf(ios::fixed);
792
793    cout << "eigenvalues = " << endl;
794    int iRow;
795    for (iRow = 0; iRow < 3; iRow++)
796        cout << kES.GetEigenvalue(iRow) << ' ';
797    cout << endl;
798
799    cout << "eigenvectors = " << endl;
800    for (iRow = 0; iRow < 3; iRow++)
801    {
802        for (int iCol = 0; iCol < 3; iCol++)
803            cout << kES.GetEigenvector(iRow,iCol) << ' ';
804        cout << endl;
805    }
806
807    // eigenvalues =
808    //    1.000000 1.000000 4.000000
809    // eigenvectors =
810    //    0.411953  0.704955 0.577350
811    //    0.404533 -0.709239 0.577350
812    //   -0.816485  0.004284 0.577350
813
814    return 0;
815}
816
817#endif
818*/
819
820#endif
Note: See TracBrowser for help on using the repository browser.