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");

}

The Jacobi method is only for symmetric matrices. Read the description in the book for more info.

Saul Teukolsky

Thank you so much
your comment help me so much