cgma
mgcEigenD Class Reference

#include <AnalyticGeometryTool.hpp>

List of all members.

Public Member Functions

 mgcEigenD (int _size)
 ~mgcEigenD ()
double & Matrix (int row, int col)
mgcEigenDMatrix (double **inmat)
double Eigenvector (int row, int col)
const double ** Eigenvector ()
void EigenStuff3 ()

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 (double **mat, double *diag, double *subd)
void Tridiagonal3 (double **mat, double *diag, double *subd)
void Tridiagonal4 (double **mat, double *diag, double *subd)
void TridiagonalN (int n, double **mat, double *diag, double *subd)
void QLAlgorithm (int n, double *diag, double *subd, double **mat)
void DecreasingSort (int n, double *eigval, double **eigvec)
void IncreasingSort (int n, double *eigval, double **eigvec)

Static Private Member Functions

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

Private Attributes

int size
double ** mat
double * diag
double * 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 909 of file AnalyticGeometryTool.hpp.


Constructor & Destructor Documentation

mgcEigenD::mgcEigenD ( int  _size)

Definition at line 3907 of file AnalyticGeometryTool.cpp.

{
        if ( (size = _size) <= 1 ) {
                Report(invalid_size);
                return;
        }
        if ( (mat = new double*[size]) == 0 ) {
                Report(allocation_failed);
                return;
        }
        for (int d = 0; d < size; d++)
                if ( (mat[d] = new double[size]) == 0 ) {
                        Report(allocation_failed);
                        return;
                }
        if ( (diag = new double[size]) == 0 ) {
                Report(allocation_failed);
                return;
        }
        if ( (subd = new double[size]) == 0 ) {
                Report(allocation_failed);
                return;
        }
}

Definition at line 3933 of file AnalyticGeometryTool.cpp.

{
        delete[] subd;
        delete[] diag;
        for (int d = 0; d < size; d++)
                delete[] mat[d];
        delete[] mat;
}

Member Function Documentation

void mgcEigenD::DecreasingSort ( int  n,
double *  eigval,
double **  eigvec 
) [private]

Definition at line 4257 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
                double 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++) {
                                double tmp = eigvec[j][i];
                                eigvec[j][i] = eigvec[j][k];
                                eigvec[j][k] = tmp;
                        }
                }
        }
}
double mgcEigenD::Eigenvector ( int  row,
int  col 
) [inline]

Definition at line 920 of file AnalyticGeometryTool.hpp.

{ return mat[row][col]; }
const double** mgcEigenD::Eigenvector ( ) [inline]

Definition at line 921 of file AnalyticGeometryTool.hpp.

{ return (const double**) mat; }
void mgcEigenD::IncreasingSort ( int  n,
double *  eigval,
double **  eigvec 
) [private]

Definition at line 4284 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
                double 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++) {
                                double tmp = eigvec[j][i];
                                eigvec[j][i] = eigvec[j][k];
                                eigvec[j][k] = tmp;
                        }
                }
        }
}
double& mgcEigenD::Matrix ( int  row,
int  col 
) [inline]

Definition at line 916 of file AnalyticGeometryTool.hpp.

{ return mat[row][col]; }
mgcEigenD & mgcEigenD::Matrix ( double **  inmat)

Definition at line 4311 of file AnalyticGeometryTool.cpp.

{
        for (int row = 0; row < size; row++)
                for (int col = 0; col < size; col++)
                        mat[row][col] = inmat[row][col];
        return *this;
}
int mgcEigenD::Number ( unsigned  single_error) [static, private]

Definition at line 4327 of file AnalyticGeometryTool.cpp.

{
        int result;
        for (result = -1; single_error; single_error >>= 1)
                result++;
        return result;
}
void mgcEigenD::QLAlgorithm ( int  n,
double *  diag,
double *  subd,
double **  mat 
) [private]

Definition at line 4196 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++) {
                                double dd = fabs(pdiag[m])+fabs(pdiag[m+1]);
                                if ( (double)(fabs(psubd[m])+dd) == dd )
                                        break;
                        }
                        if ( m == ell )
                                break;

                        double g = (pdiag[ell+1]-pdiag[ell])/(2*psubd[ell]);
                        double r = 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);
                        double s = 1, c = 1, p = 0;
                        for (int i = m-1; i >= ell; i--) {
                                double f = s*psubd[i], b = c*psubd[i];
                                if ( fabs(f) >= fabs(g) ) {
                                        c = g/f;
                                        r = sqrt(c*c+1);
                                        psubd[i+1] = f*r;
                                        c *= (s = 1/r);
                                }
                                else {
                                        s = f/g;
                                        r = 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 mgcEigenD::Report ( std::ostream &  ostr) [static]
void mgcEigenD::Report ( unsigned  single_error) [static, private]

Definition at line 4336 of file AnalyticGeometryTool.cpp.

{
        if ( mgcEigenD::verbose1 )
                cout << "mgcEigenD: " << message[Number(single_error)] << endl;
        else {
          ofstream ostr("eigen.err",ios::out|ios::app);
          ostr << "mgcEigenD: " << message[Number(single_error)] << endl;
        }
        error |= single_error;
}
void mgcEigenD::Tridiagonal2 ( double **  mat,
double *  diag,
double *  subd 
) [private]

Definition at line 3943 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 mgcEigenD::Tridiagonal3 ( double **  mat,
double *  diag,
double *  subd 
) [private]

Definition at line 3956 of file AnalyticGeometryTool.cpp.

{
        double 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 ) {
                double ell = sqrt(b*b+c*c);
                b /= ell;
                c /= ell;
                double 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 mgcEigenD::Tridiagonal4 ( double **  mat,
double *  diag,
double *  subd 
) [private]

Definition at line 3988 of file AnalyticGeometryTool.cpp.

{
        // save pmatrix M
        double
        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 ) {
                double q11, q12, q13;
                double q21, q22, q23;
                double q31, q32, q33;

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

                psubd[0] = len;

                // compute S*Q1
                double v0 = e*q11+f*q21+g*q31;
                double v1 = f*q11+h*q21+i*q31;
                double 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 = 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 ) {
                                double 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 ) {
                        double ell = sqrt(f*f+g*g);
                        f /= ell;
                        g /= ell;
                        double 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 mgcEigenD::TridiagonalN ( int  n,
double **  mat,
double *  diag,
double *  subd 
) [private]

Definition at line 4124 of file AnalyticGeometryTool.cpp.

{
        int i, j, k, ell;

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

                if ( ell > 0 ) {
                        for (k = 0; k <= ell; k++)
                                scale += 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];
                                }
                                double f = pmat[i][ell];
                                double g = ( f > 0 ? -sqrt(h) : 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];
                                }
                                double 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++) {
                                double 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 mgcEigenD::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 mgcEigenD::allocation_failed = 0x00000002 [static, private]

Definition at line 953 of file AnalyticGeometryTool.hpp.

double* mgcEigenD::diag [private]

Definition at line 928 of file AnalyticGeometryTool.hpp.

unsigned mgcEigenD::error = 0 [static]

Definition at line 949 of file AnalyticGeometryTool.hpp.

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

Definition at line 952 of file AnalyticGeometryTool.hpp.

double** mgcEigenD::mat [private]

Definition at line 927 of file AnalyticGeometryTool.hpp.

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

Definition at line 955 of file AnalyticGeometryTool.hpp.

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

Definition at line 954 of file AnalyticGeometryTool.hpp.

int mgcEigenD::size [private]

Definition at line 926 of file AnalyticGeometryTool.hpp.

double* mgcEigenD::subd [private]

Definition at line 929 of file AnalyticGeometryTool.hpp.

int mgcEigenD::verbose1 = 0 [static]

Definition at line 948 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