cgma
mgcEigen Class Reference

#include <AnalyticGeometryTool.hpp>

List of all members.

Public Member Functions

 mgcEigen (int _size)
 ~mgcEigen ()
mgcEigenMatrix (float **inmat)

Static Public Member Functions

static void Report (std::ostream &ostr)

Static Public Attributes

static int verbose1 = 0
static unsigned error = 0

Private Member Functions

void Tridiagonal2 (float **pmat, float *pdiag, float *psubd)
void Tridiagonal3 (float **pmat, float *pdiag, float *psubd)
void Tridiagonal4 (float **pmat, float *pdiag, float *psubd)
void TridiagonalN (int n, float **mat, float *diag, float *subd)
void QLAlgorithm (int n, float *pdiag, float *psubd, float **pmat)
void DecreasingSort (int n, float *eigval, float **eigvec)
void IncreasingSort (int n, float *eigval, float **eigvec)

Static Private Member Functions

static int Number (unsigned single_error)
static void Report (unsigned single_error)

Private Attributes

int size
float ** mat
float * diag
float * subd

Static Private Attributes

static const unsigned invalid_size = 0x00000001
static const unsigned allocation_failed = 0x00000002
static const unsigned ql_exceeded = 0x00000004
static const char * message [3]

Detailed Description

Definition at line 865 of file AnalyticGeometryTool.hpp.


Constructor & Destructor Documentation

mgcEigen::mgcEigen ( int  _size)

Member Function Documentation

void mgcEigen::DecreasingSort ( int  n,
float *  eigval,
float **  eigvec 
) [private]

Definition at line 3810 of file AnalyticGeometryTool.cpp.

{
        // sort eigenvalues in decreasing order, e[0] >= ... >= e[n-1]
        for (int i = 0, k; i <= n-2; i++) {
                // locate maximum eigenvalue
                float max = eigval[k=i];
                int j;
                for (j = i+1; j < n; j++)
                        if ( eigval[j] > max )
                                max = eigval[k=j];

                if ( k != i ) {
                        // swap eigenvalues
                        eigval[k] = eigval[i];
                        eigval[i] = max;

                        // swap eigenvectors
                        for (j = 0; j < n; j++) {
                                float tmp = eigvec[j][i];
                                eigvec[j][i] = eigvec[j][k];
                                eigvec[j][k] = tmp;
                        }
                }
        }
}
void mgcEigen::IncreasingSort ( int  n,
float *  eigval,
float **  eigvec 
) [private]

Definition at line 3837 of file AnalyticGeometryTool.cpp.

{
        // sort eigenvalues in increasing order, e[0] <= ... <= e[n-1]
        for (int i = 0, k; i <= n-2; i++) {
                // locate minimum eigenvalue
                float min = eigval[k=i];
        int j;
                for (j = i+1; j < n; j++)
                        if ( eigval[j] < min )
                                min = eigval[k=j];

                if ( k != i ) {
                        // swap eigenvalues
                        eigval[k] = eigval[i];
                        eigval[i] = min;

                        // swap eigenvectors
                        for (j = 0; j < n; j++) {
                                float tmp = eigvec[j][i];
                                eigvec[j][i] = eigvec[j][k];
                                eigvec[j][k] = tmp;
                        }
                }
        }
}
mgcEigen& mgcEigen::Matrix ( float **  inmat)
int mgcEigen::Number ( unsigned  single_error) [static, private]

Definition at line 3864 of file AnalyticGeometryTool.cpp.

{
        int result;
        for (result = -1; single_error; single_error >>= 1)
                result++;
        return result;
}
void mgcEigen::QLAlgorithm ( int  n,
float *  pdiag,
float *  psubd,
float **  pmat 
) [private]

Definition at line 3749 of file AnalyticGeometryTool.cpp.

{
        const int eigen_maxiter = 30;

        for (int ell = 0; ell < n; ell++) {
                int iter;
                for (iter = 0; iter < eigen_maxiter; iter++) {
                        int m;
                        for (m = ell; m <= n-2; m++) {
                                float dd = float(fabs(pdiag[m])+fabs(pdiag[m+1]));
                                if ( (float)(fabs(psubd[m])+dd) == dd )
                                        break;
                        }
                        if ( m == ell )
                                break;

                        float g = (pdiag[ell+1]-pdiag[ell])/(2*psubd[ell]);
                        float r = float(sqrt(g*g+1));
                        if ( g < 0 )
                                g = pdiag[m]-pdiag[ell]+psubd[ell]/(g-r);
                        else
                                g = pdiag[m]-pdiag[ell]+psubd[ell]/(g+r);
                        float s = 1, c = 1, p = 0;
                        for (int i = m-1; i >= ell; i--) {
                                float f = s*psubd[i], b = c*psubd[i];
                                if ( fabs(f) >= fabs(g) ) {
                                        c = g/f;
                                        r = float(sqrt(c*c+1));
                                        psubd[i+1] = f*r;
                                        c *= (s = 1/r);
                                }
                                else {
                                        s = f/g;
                                        r = float(sqrt(s*s+1));
                                        psubd[i+1] = g*r;
                                        s *= (c = 1/r);
                                }
                                g = pdiag[i+1]-p;
                                r = (pdiag[i]-g)*s+2*b*c;
                                p = s*r;
                                pdiag[i+1] = g+p;
                                g = c*r-b;

                                for (int k = 0; k < n; k++) {
                                        f = pmat[k][i+1];
                                        pmat[k][i+1] = s*pmat[k][i]+c*f;
                                        pmat[k][i] = c*pmat[k][i]-s*f;
                                }
                        }
                        pdiag[ell] -= p;
                        psubd[ell] = g;
                        psubd[m] = 0;
                }
                if ( iter == eigen_maxiter ) {
                        Report(ql_exceeded);
                        return;
                }
        }
}
static void mgcEigen::Report ( std::ostream &  ostr) [static]
void mgcEigen::Report ( unsigned  single_error) [static, private]

Definition at line 3873 of file AnalyticGeometryTool.cpp.

{
        if ( mgcEigen::verbose1 )
                cout << "mgcEigen: " << message[Number(single_error)] << endl;
        else  {
                ofstream ostr("eigen.err",ios::out|ios::app);
                ostr << "mgcEigen: " << message[Number(single_error)] << endl;
        }
        error |= single_error;
}
void mgcEigen::Tridiagonal2 ( float **  pmat,
float *  pdiag,
float *  psubd 
) [private]

Definition at line 3494 of file AnalyticGeometryTool.cpp.

{
        // matrix is already tridiagonal

        pdiag[0] = pmat[0][0];
        pdiag[1] = pmat[1][1];
        psubd[0] = pmat[0][1];
        psubd[1] = 0;
        pmat[0][0] = 1;  pmat[0][1] = 0;
        pmat[1][0] = 0;  pmat[1][1] = 1;
}
void mgcEigen::Tridiagonal3 ( float **  pmat,
float *  pdiag,
float *  psubd 
) [private]

Definition at line 3507 of file AnalyticGeometryTool.cpp.

{
        float a = pmat[0][0], b = pmat[0][1], c = pmat[0][2],
                                                 d = pmat[1][1], e = pmat[1][2],
                                                                                f = pmat[2][2];

        pdiag[0] = a;
        psubd[2] = 0;
        if ( c != 0 ) {
                float ell = float(sqrt(b*b+c*c));
                b /= ell;
                c /= ell;
                float q = 2*b*e+c*(f-d);
                pdiag[1] = d+c*q;
                pdiag[2] = f-c*q;
                psubd[0] = ell;
                psubd[1] = e-b*q;
                pmat[0][0] = 1; pmat[0][1] = 0; pmat[0][2] = 0;
                pmat[1][0] = 0; pmat[1][1] = b; pmat[1][2] = c;
                pmat[2][0] = 0; pmat[2][1] = c; pmat[2][2] = -b;
        }
        else {
                pdiag[1] = d;
                pdiag[2] = f;
                psubd[0] = b;
                psubd[1] = e;
                pmat[0][0] = 1; pmat[0][1] = 0; pmat[0][2] = 0;
                pmat[1][0] = 0; pmat[1][1] = 1; pmat[1][2] = 0;
                pmat[2][0] = 0; pmat[2][1] = 0; pmat[2][2] = 1;
        }
}
void mgcEigen::Tridiagonal4 ( float **  pmat,
float *  pdiag,
float *  psubd 
) [private]

Definition at line 3540 of file AnalyticGeometryTool.cpp.

{
        // save pmatrix M
        float
        a = pmat[0][0], b = pmat[0][1], c = pmat[0][2], d = pmat[0][3],
                                   e = pmat[1][1], f = pmat[1][2], g = pmat[1][3],
                                                                  h = pmat[2][2], i = pmat[2][3],
                                                                                                 j = pmat[3][3];

        pdiag[0] = a;
        psubd[3] = 0;

        pmat[0][0] = 1; pmat[0][1] = 0; pmat[0][2] = 0; pmat[0][3] = 0;
        pmat[1][0] = 0;
        pmat[2][0] = 0;
        pmat[3][0] = 0;

        if ( c != 0 || d != 0 ) {
                float q11, q12, q13;
                float q21, q22, q23;
                float q31, q32, q33;

                // build column Q1
                float len = float(sqrt(b*b+c*c+d*d));
                q11 = b/len;
                q21 = c/len;
                q31 = d/len;

                psubd[0] = len;

                // compute S*Q1
                float v0 = e*q11+f*q21+g*q31;
                float v1 = f*q11+h*q21+i*q31;
                float v2 = g*q11+i*q21+j*q31;

                pdiag[1] = q11*v0+q21*v1+q31*v2;

        // build column Q3 = Q1x(S*Q1)
                q13 = q21*v2-q31*v1;
                q23 = q31*v0-q11*v2;
                q33 = q11*v1-q21*v0;
                len = float(sqrt(q13*q13+q23*q23+q33*q33));
                if ( len > 0 ) {
                        q13 /= len;
                        q23 /= len;
                        q33 /= len;

                        // build column Q2 = Q3xQ1
                        q12 = q23*q31-q33*q21; 
                        q22 = q33*q11-q13*q31;
                        q32 = q13*q21-q23*q11;

                        v0 = q12*e+q22*f+q32*g;
                        v1 = q12*f+q22*h+q32*i;
                        v2 = q12*g+q22*i+q32*j;
                        psubd[1] = q11*v0+q21*v1+q31*v2;
                        pdiag[2] = q12*v0+q22*v1+q32*v2;
                        psubd[2] = q13*v0+q23*v1+q33*v2;

                        v0 = q13*e+q23*f+q33*g;
                        v1 = q13*f+q23*h+q33*i;
                        v2 = q13*g+q23*i+q33*j;
                        pdiag[3] = q13*v0+q23*v1+q33*v2;
                }
                else {  // S*Q1 parallel to Q1, choose any valid Q2 and Q3
                        psubd[1] = 0;

                        len = q21*q21+q31*q31;
                        if ( len > 0 ) {
                                float tmp = q11-1;
                                q12 = -q21;
                                q22 = 1+tmp*q21*q21/len;
                q32 = tmp*q21*q31/len;

                                q13 = -q31;
                                q23 = q32;
                                q33 = 1+tmp*q31*q31/len;

                                v0 = q12*e+q22*f+q32*g;
                                v1 = q12*f+q22*h+q32*i;
                                v2 = q12*g+q22*i+q32*j;
                                pdiag[2] = q12*v0+q22*v1+q32*v2;
                                psubd[2] = q13*v0+q23*v1+q33*v2;

                                v0 = q13*e+q23*f+q33*g;
                                v1 = q13*f+q23*h+q33*i;
                                v2 = q13*g+q23*i+q33*j;
                                pdiag[3] = q13*v0+q23*v1+q33*v2;
                        }
                        else {  // Q1 = (+-1,0,0)
                                q12 = 0; q22 = 1; q32 = 0;
                                q13 = 0; q23 = 0; q33 = 1;

                                pdiag[2] = h;
                                pdiag[3] = j;
                                psubd[2] = i;
                        }
                }

                pmat[1][1] = q11; pmat[1][2] = q12; pmat[1][3] = q13;
                pmat[2][1] = q21; pmat[2][2] = q22; pmat[2][3] = q23;
                pmat[3][1] = q31; pmat[3][2] = q32; pmat[3][3] = q33;
        }
        else {
                pdiag[1] = e;
                psubd[0] = b;
                pmat[1][1] = 1;
                pmat[2][1] = 0;
                pmat[3][1] = 0; 

                if ( g != 0 ) {
                        float ell = float(sqrt(f*f+g*g));
                        f /= ell;
                        g /= ell;
                        float Q = 2*f*i+g*(j-h);

                        pdiag[2] = h+g*Q;
                        pdiag[3] = j-g*Q;
                        psubd[1] = ell;
                        psubd[2] = i-f*Q;
                        pmat[1][2] = 0;  pmat[1][3] = 0;
                        pmat[2][2] = f;  pmat[2][3] = g;
                        pmat[3][2] = g;  pmat[3][3] = -f;
                }
                else {
                        pdiag[2] = h;
                        pdiag[3] = j;
                        psubd[1] = f;
                        psubd[2] = i;
                        pmat[1][2] = 0;  pmat[1][3] = 0;
                        pmat[2][2] = 1;  pmat[2][3] = 0;
                        pmat[3][2] = 0;  pmat[3][3] = 1;
                }
        }
}
void mgcEigen::TridiagonalN ( int  n,
float **  mat,
float *  diag,
float *  subd 
) [private]

Definition at line 3677 of file AnalyticGeometryTool.cpp.

{
        int i, j, k, ell;

        for (i = n-1, ell = n-2; i >= 1; i--, ell--) {
                float h = 0, scale = 0;

                if ( ell > 0 ) {
                        for (k = 0; k <= ell; k++)
                                scale += float(fabs(pmat[i][k]));
                        if ( scale == 0 )
                                psubd[i] = pmat[i][ell];
                        else {
                                for (k = 0; k <= ell; k++) {
                                        pmat[i][k] /= scale;
                                        h += pmat[i][k]*pmat[i][k];
                                }
                                float f = pmat[i][ell];
                                float g = ( f > 0 ? -float(sqrt(h)) : float(sqrt(h)) );
                                psubd[i] = scale*g;
                                h -= f*g;
                                pmat[i][ell] = f-g;
                                f = 0;
                                for (j = 0; j <= ell; j++) {
                                        pmat[j][i] = pmat[i][j]/h;
                                        g = 0;
                                        for (k = 0; k <= j; k++)
                                                g += pmat[j][k]*pmat[i][k];
                                        for (k = j+1; k <= ell; k++)
                                                g += pmat[k][j]*pmat[i][k];
                                        psubd[j] = g/h;
                                        f += psubd[j]*pmat[i][j];
                                }
                                float hh = f/(h+h);
                                for (j = 0; j <= ell; j++) {
                                        f = pmat[i][j];
                                        psubd[j] = g = psubd[j] - hh*f;
                                        for (k = 0; k <= j; k++)
                                                pmat[j][k] -= f*psubd[k]+g*pmat[i][k];
                                }
            }
                }
                else
                        psubd[i] = pmat[i][ell];

                pdiag[i] = h;
        }

        pdiag[0] = psubd[0] = 0;
        for (i = 0, ell = -1; i <= n-1; i++, ell++) {
                if ( pdiag[i] ) {
                        for (j = 0; j <= ell; j++) {
                                float sum = 0;
                                for (k = 0; k <= ell; k++)
                                        sum += pmat[i][k]*pmat[k][j];
                                for (k = 0; k <= ell; k++)
                                        pmat[k][j] -= sum*pmat[k][i];
                        }
                }
                pdiag[i] = pmat[i][i];
                pmat[i][i] = 1;
                for (j = 0; j <= ell; j++)
                        pmat[j][i] = pmat[i][j] = 0;
        }

        // re-ordering if mgcEigen::QLAlgorithm is used subsequently
        for (i = 1, ell = 0; i < n; i++, ell++)
                psubd[ell] = psubd[i];
        psubd[n-1] = 0;
}

Member Data Documentation

const unsigned mgcEigen::allocation_failed = 0x00000002 [static, private]

Definition at line 901 of file AnalyticGeometryTool.hpp.

float* mgcEigen::diag [private]

Definition at line 876 of file AnalyticGeometryTool.hpp.

unsigned mgcEigen::error = 0 [static]

Definition at line 897 of file AnalyticGeometryTool.hpp.

const unsigned mgcEigen::invalid_size = 0x00000001 [static, private]

Definition at line 900 of file AnalyticGeometryTool.hpp.

float** mgcEigen::mat [private]

Definition at line 875 of file AnalyticGeometryTool.hpp.

const char * mgcEigen::message [static, private]
Initial value:
 {
        "invalid matrix size",
        "allocation failed",
        "QL algorithm - exceeded maximum iterations"
}

Definition at line 903 of file AnalyticGeometryTool.hpp.

const unsigned mgcEigen::ql_exceeded = 0x00000004 [static, private]

Definition at line 902 of file AnalyticGeometryTool.hpp.

int mgcEigen::size [private]

Definition at line 874 of file AnalyticGeometryTool.hpp.

float* mgcEigen::subd [private]

Definition at line 877 of file AnalyticGeometryTool.hpp.

int mgcEigen::verbose1 = 0 [static]

Definition at line 896 of file AnalyticGeometryTool.hpp.


The documentation for this class was generated from the following files:
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines