me_here_me
07-30-2007, 02:57 PM
HI
The jacobi method reipe does not seem to work here with me. I took the code from the book but the results that I get seem to differ. Below is the code that I am using:
#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);
void jacobi (float a[][4], int n, float d[], float v[][4], int *nrot){
int j, iq, ip, i;
float tresh, theta, tau, t, sm, s, h, g, c, *b, *z;
b = new float[n];
z = new float[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;
}
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;
}
}
}
I tried running the following example matrix:
/*
4 -30 60 -35
-30 300 -675 420
60 -675 1620 -1050
-35 420 -1050 700
*/
The answer to this matrix should be:
eigen values:
(2585.25 )
( 37.10 )
( 0.17 )
( 1.48 )
and eigen vectors:
( 0.03 ) (-0.18 ) ( 0.79 ) (-0.58 )
(-0.33 ) ( 0.74 ) ( 0.45 ) ( 0.37 )
( 0.79 ) (-0.10 ) ( 0.32 ) ( 0.51 )
(-0.51 ) (-0.64 ) ( 0.25 ) ( 0.51 )
as given in: http://en.wikipedia.org/wiki/Jacobi_eigenvalue_algorithm
but the output that I got is:
eigen values:
(2585.25 )
( 39.99 )
( 1.99 )
( -2.8 )
and eigen vectors:
( 0.03 ) (-0.33 ) ( 0.79 ) (-0.51 )
( 0.15 ) (-0.74 ) ( 0.09 ) ( 0.63 )
( 0.85 ) ( 0.38 ) ( 0.26 ) ( 0.21 )
(-0.51 ) ( 0.42 ) ( 0.54 ) ( 0.51 )
I tried a couple of more matrices and cross checked by the online eigen vector calculator given at:
http://www.bluebit.gr/matrix-calculator/
and there was always big diffference, can some one help me with this matter. I need the solution kindda urgent
thanks in advance
The jacobi method reipe does not seem to work here with me. I took the code from the book but the results that I get seem to differ. Below is the code that I am using:
#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);
void jacobi (float a[][4], int n, float d[], float v[][4], int *nrot){
int j, iq, ip, i;
float tresh, theta, tau, t, sm, s, h, g, c, *b, *z;
b = new float[n];
z = new float[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;
}
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;
}
}
}
I tried running the following example matrix:
/*
4 -30 60 -35
-30 300 -675 420
60 -675 1620 -1050
-35 420 -1050 700
*/
The answer to this matrix should be:
eigen values:
(2585.25 )
( 37.10 )
( 0.17 )
( 1.48 )
and eigen vectors:
( 0.03 ) (-0.18 ) ( 0.79 ) (-0.58 )
(-0.33 ) ( 0.74 ) ( 0.45 ) ( 0.37 )
( 0.79 ) (-0.10 ) ( 0.32 ) ( 0.51 )
(-0.51 ) (-0.64 ) ( 0.25 ) ( 0.51 )
as given in: http://en.wikipedia.org/wiki/Jacobi_eigenvalue_algorithm
but the output that I got is:
eigen values:
(2585.25 )
( 39.99 )
( 1.99 )
( -2.8 )
and eigen vectors:
( 0.03 ) (-0.33 ) ( 0.79 ) (-0.51 )
( 0.15 ) (-0.74 ) ( 0.09 ) ( 0.63 )
( 0.85 ) ( 0.38 ) ( 0.26 ) ( 0.21 )
(-0.51 ) ( 0.42 ) ( 0.54 ) ( 0.51 )
I tried a couple of more matrices and cross checked by the online eigen vector calculator given at:
http://www.bluebit.gr/matrix-calculator/
and there was always big diffference, can some one help me with this matter. I need the solution kindda urgent
thanks in advance