|
cgma
|
#include <AnalyticGeometryTool.hpp>
Public Member Functions | |
| mgcEigen (int _size) | |
| ~mgcEigen () | |
| mgcEigen & | Matrix (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] |
Definition at line 865 of file AnalyticGeometryTool.hpp.
| mgcEigen::mgcEigen | ( | int | _size | ) |
| 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.
| 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;
}
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] |
{
"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.
1.7.6.1