|
cgma
|
#include <AnalyticGeometryTool.hpp>
Public Member Functions | |
| mgcEigenD (int _size) | |
| ~mgcEigenD () | |
| double & | Matrix (int row, int col) |
| mgcEigenD & | Matrix (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] |
Definition at line 909 of file AnalyticGeometryTool.hpp.
| 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;
}
}
| 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;
}
}
}
}
| void mgcEigenD::EigenStuff3 | ( | ) |
Definition at line 4320 of file AnalyticGeometryTool.cpp.
{
Tridiagonal3(mat,diag,subd);
QLAlgorithm(size,diag,subd,mat);
}
| 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.
| 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.
| 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;
}
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] |
{
"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.
1.7.6.1