jocobi method recipe not working


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

Saul Teukolsky
07-30-2007, 07:33 PM
Hi,

Your code is not NR code, but your own translation. At least one possible source of error is your use of built-in C vectors instead of the NR vector library. (Among other things, zero-based vs. 1-based.) See Chapter 1 of the book.

Saul Teukolsky

me_here_me
07-31-2007, 03:47 AM
Hi,
I converted the code to zero based from the 1 based given in the book.

I will use the NR vector library to see if things get better :)

can I download the code in the book from some website. It will be much more faster then typing it from the book?

me_here_me
07-31-2007, 06:16 AM
Thanks, u were right :)

the error was coming from the way I converted the code.
Its fine and running now

:)

mahcs2010
07-13-2008, 06:50 AM
Hi all
I use this function, it gives me correct results for symatric matrix and different results than orther methods for unsymatric
could you help me please

void jacobi(double **a, double **v, double *d, int n, int *nrot)
{
int j,iq,ip,i;
double tresh,theta,tau,t,sm,s,h,g,c,*b,*z;
// b=vector(1,n);
// z=vector(1,n);
b=(void*)xalloc(n * sizeof (double));
z=(void*)xalloc(n * sizeof (double));

for (ip=0;ip<n;ip++)
{
//Initialize to the identity matrix.
for (iq=0;iq<n;iq++)
v[ip][iq]=0.0;
v[ip][ip]=1.0;
}

for (ip=0;ip<n;ip++)
{
//Initialize b and d to the diagonal of a. This vector will
//accumulate terms of the form tapq as in equation (11.1.14).
b[ip]=d[ip]=a[ip][ip];
z[ip]=0.0;
}
*nrot=0;

for (i=1;i<=50;i++)
{
sm=0.0;
for (ip=0;ip<n-1;ip++)
{
//Sum off-diagonal elements.
for (iq=ip+1;iq<n;iq++)
sm += (double)fabs(a[ip][iq]);
}
/*printmatrix(a,3,3);*/
// printf("Sum=%f\t%d\n",sm,*nrot);
if (sm == 0.0)
{
//The normal return, which relies on
//quadratic convergence to machine underflow.
// free_vector(z,1,n);
// free_vector(b,1,n);
free(z);
free(b);
return;
}
if (i < 4)
tresh=(double)0.2*sm/(n*n); //...on the first three sweeps.
else
tresh=0.0; //...thereafter.
for (ip=0;ip<n-1;ip++)
{
for (iq=ip+1;iq<n;iq++)
{
g=(double)( 100.0*fabs(a[ip][iq]) );
//After four sweeps, skip the rotation if the o -diagonal element is small..
if (i > 4 && (double)(fabs(d[ip])+g)==(double)fabs(d[ip]) && (double)(fabs(d[iq])+g)==(double)fabs(d[iq]))
a[ip][iq]=0.0;
else if (fabs(a[ip][iq]) > tresh)
{
h=d[iq]-d[ip];
if ((double)(fabs(h)+g) == (double)fabs(h))
t=(a[ip][iq])/h; //t=1/(2*Theta)
else
{
theta=(double)( 0.5*h/(a[ip][iq]) ); //Equation (11.1.10).
t=(double)(1.0/(fabs(theta)+sqrt(1.0+theta*theta)) );
if (theta < 0.0)
t = -t;
}
c=(double)(1.0/sqrt(1+t*t));
s=t*c;
tau=(double)(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;j++)
{
//Case of rotations 1 <= j < p.
ROTATE(a,j,ip,j,iq);
}
for (j=ip+1;j<iq;j++)
{
//Case of rotations p < j < q.
ROTATE(a,ip,j,j,iq);
}
for (j=iq+1;j<n;j++)
{
//Case of rotations q < j <=n.
ROTATE(a,ip,j,iq,j);
}
for (j=0;j<n;j++)
{
ROTATE(v,j,ip,j,iq);
}
++(*nrot);
}
}

}//i=1...50

for (ip=0;ip<n;ip++)
{
//Update d with the sum of tapq, and reinitialize z.
b[ip] += z[ip];
d[ip]=b[ip];
z[ip]=0.0;
}

}
error("Too many iterations in routine jacobi");

}