Commit 3fd22904 authored by Will Schroeder's avatar Will Schroeder
Browse files

ENH: Added jacobi method.

parent 3a908aed
......@@ -357,3 +357,146 @@ void vtkMath::SingularValueBackSubstitution(double **u, double *w, double **v,
delete [] tmp;
}
// Description:
// Jacobi iteration for the solution of eigenvectors/eigenvalues of a real
// symmetric matrix. Square matrix a; system size n; output eigenvalues in d;
// eigenvectors in v; and number of jacobi rotations (output) nrot.
#define ROTATE(a,i,j,k,l) g=a[i][j];h=a[k][l];a[i][j]=g-s*(h+g*tau);\
a[k][l]=h+s*(g-h*tau);
int vtkMath::Jacobi(float **a, int n, float *d, float **v, int *nrot)
{
int j,iq,ip,i;
float tresh,theta,tau,t,sm,s,h,g,c;
static float *b = {NULL};
static float *z = {NULL};
static int nMax = {0};
//
// Check on allocation of working vectors
//
if ( n > nMax )
{
delete [] b; b = NULL;
delete [] z; z = NULL;
}
if ( b == NULL )
{
b = new float[n];
z = new float[n];
nMax = n;
}
for (ip=0;ip<n;ip++)
{
for (iq=0;iq<n;iq++) v[ip][iq]=0.0;
v[ip][ip]=1.0;
}
for (ip=0;ip<n;ip++)
{
b[ip]=d[ip]=a[ip][ip];
z[ip]=0.0;
}
*nrot=0;
for (i=0;i<50;i++)
{
sm=0.0;
for (ip=0;ip<n-1;ip++)
{
for (iq=ip+1;iq<n;iq++)
sm += fabs(a[ip][iq]);
}
if (sm == 0.0) return 1;
if (i < 4)
tresh=0.2*sm/(n*n);
else
tresh=0.0;
for (ip=0;ip<n-1;ip++)
{
for (iq=ip+1;iq<n;iq++)
{
g=100.0*fabs(a[ip][iq]);
if (i > 4 && (float)(fabs(d[ip])+g) == (float)fabs(d[ip])
&& (float)(fabs(d[iq])+g) == (float)fabs(d[iq]))
a[ip][iq]=0.0;
else if (fabs(a[ip][iq]) > tresh)
{
h=d[iq]-d[ip];
if ((float)(fabs(h)+g) == (float)fabs(h))
t=(a[ip][iq])/h;
else
{
theta=0.5*h/(a[ip][iq]);
t=1.0/(fabs(theta)+sqrt(1.0+theta*theta));
if (theta < 0.0) t = -t;
}
c=1.0/sqrt(1+t*t);
s=t*c;
tau=s/(1.0+c);
h=t*a[ip][iq];
z[ip] -= h;
z[iq] += h;
d[ip] -= h;
d[iq] += h;
a[ip][iq]=0.0;
for (j=0;j<ip-1;j++)
{
ROTATE(a,j,ip,j,iq)
}
for (j=ip+1;j<iq-1;j++)
{
ROTATE(a,ip,j,j,iq)
}
for (j=iq+1;j<n;j++)
{
ROTATE(a,ip,j,iq,j)
}
for (j=0;j<n;j++)
{
ROTATE(v,j,ip,j,iq)
}
++(*nrot);
}
}
}
for (ip=0;ip<n;ip++)
{
b[ip] += z[ip];
d[ip]=b[ip];
z[ip]=0.0;
}
}
// vtkErrorMacro(<<"Error extracting eigenfunctions")'
return 0;
}
#undef ROTATE
// Description:
// Sort eigenvalues/eigenvectors into descending order. Call after Jacobi.
// Vector d are the eigenvalues; matrix v are eigenvectors (in columns; n
// is the size of the matrix).
void vtkMath::Eigsrt(float *d, float **v, int n)
{
int k,j,i;
float p;
for (i=0;i<n;i++)
{
p=d[k=i];
for (j=i; j<n; j++)
if (d[j] >= p) p=d[k=j];
if (k != i)
{
d[k]=d[i];
d[i]=p;
for (j=0;j<n;j++)
{
p=v[j][i];
v[j][i]=v[j][k];
v[j][k]=p;
}
}
}
}
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment