Hello,

I took the Jacobi routine from the book "Numerical Recipes in C (page 467-468)" to obatin eigenvalues and corresponding eigenvectors from a real symmetric matrix.

There are no compilation errors, but the program does not produce correct eigenvalues and eigenvectors. It's completely wrong.

Given a real symmetric matrix is
5.0 -2.0
-2.0 2.0

The eigenvalues and eigenvectors generated by the program are
0.000000
2.000000

and

0.000000 0.000000
0.000000 1.000000

Actually, the program produces same eigenvectors regardless of what a given symmetric matrix is. Also, it always produces the last element of a given symmetric matrix instead of eigenvalues.
I already checked several times.

That means, the program really does not work.
Would you tell me what I need to fix the routine?

<IMPORTANT>

Followings are main program constructed by myself and a subroutine taken from the book "Numerical Recipes" respectively.

************* main program *************

#include <stdio.h>

/*declare global variables */

float a[2][2]={5.0,-2.0,-2.0,2.0};

float d[2];
float v[2][2];
int n=2;
int nrot=0;

void jacobi(float (*a)[2],int n, float* d,
float (*v)[2], int *nrot);

int main(void)
{
int x=0;
int y=0;

jacobi(a,n,d,v,&nrot);

for(x=0;x<2;x++)
{
printf("%lf \n",d[x]);
//returns the eigenvalues of a
}
printf("\n");

for(x=0;x<2;x++)
{
for(y=0;y<2;y++)
{

printf("%lf",v[x][y]);
/ / returns the eigenvectors of a
}

return 0;
}


***** subroutine (refer to page 467-468)*****

#include <math.h>
#include "nrutil.h"
#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)[2], int n, float* d, float (*v)[2], int *nrot)
{
int j,iq,ip,i;
float tresh,theta,tau,t,sm,s,h,g,c,*b,*z;

b=vector(1,n);
z=vector(1,n);
for (ip=1;ip<=n;ip++) {
for (iq=1;iq<=n;iq++) v[ip][iq]=0.0;
v[ip][ip]=1.0;
}



for (ip=1;ip<=n;ip++)
{
b[ip]=d[ip]=a[ip][ip];
z[ip]=0.0;
}
*nrot=0;
for (i=1;i<=50;i++) {
sm=0.0;
for (ip=1;ip<=n-1;ip++) {
for (iq=ip+1;iq<=n;iq++)
sm += fabs(a[ip][iq]);
}
if (sm == 0.0) {
free_vector(z,1,n);
free_vector(b,1,n);
return;
}
if (i < 4)
tresh=0.2*sm/(n*n);
else
tresh=0.0;
for (ip=1;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=1;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=1;j<=n;j++) {
ROTATE(v,j,ip,j,iq)
}
++(*nrot);
}
}
}
for (ip=1;ip<=n;ip++) {
b[ip] += z[ip];
d[ip]=b[ip];
z[ip]=0.0;
}
}
nrerror("Too many iterations in routine jacobi");


}

I took the Jacobi routine from the book ...It's completely wrong.
If you are going to use the old version, you have two choices:


Play by the rules. Use the functions as listed and call them as the book tells you. See Footnote.


Make your own rules: convert all code from the index-1 based matrix and vector functions to index-0 base. Change "pointer to pointer" arguments to 2-D arrays if you want to.

However...

You have to do these things for all of the functions that are called by the functions that you call.


Here's an example showing how you can use the function from the book:

/*
Driver for Numerical Recipes Version 2 jacobi() function

davekw7x
*/
#include <stdio.h>
#include "nr.h"
#include "nrutil.h"

int main(void)
{
int i, j, nrot;

int n = 2;
/*
We need a Numerical Recipes-style "matrix" to use
as an argument to call the jacobi function.
Here are the matrix values stored in row-major
order in a 1-D array. This array name will be used
as an argument for the convert_matrix function.
*/
float a_array[] = {
5.0, -2.0, /* First row */
-2.0, 2.0 /* Second row */
};

float *d = vector(1, n);
float *r = vector(1, n);
float **v = matrix(1, n, 1, n);
float **a = convert_matrix(&a_array[0], 1, n, 1, n);

printf("Matrix:\n");
for (i = 1; i <= n; i++) {
for (j = 1; j <= n; j++) {
printf(" %+13.6e", a[i][j]);
}
printf("\n");
}
printf("\n");

jacobi(a, n, d, v, &nrot);

printf("Eigenvalues:\n");
for (i = 1; i <= n; i++) {
printf(" Number %d:%+14.6e\n", i, d[i]);
}
printf("\n");

printf("Eigenvectors:\n");
for (i = 1; i <= n; i++) {
printf(" Number %d:", i);
for (j = 1; j <= n; j++) {
printf("%+14.6e", v[j][i]);
}
printf("\n");
}
printf("\n");

printf("Number of Jacobi rotations = %d\n", nrot);

free_convert_matrix(a, 1, n, 1, n);
free_matrix(v, 1, n, 1, n);
free_vector(r, 1, n);
free_vector(d, 1, n);

return 0;
}


Output:

Matrix:
+5.000000e+00 -2.000000e+00
-2.000000e+00 +2.000000e+00

Eigenvalues:
Number 1: +6.000000e+00
Number 2: +1.000000e+00

Eigenvectors:
Number 1: +8.944272e-01 -4.472136e-01
Number 2: +4.472136e-01 +8.944272e-01

Number of Jacobi rotations = 1



Regards,

Dave

Footnote:
The old stuff in the "Numerical Recipes for C" consists mainly of Fortran functions that were transmogrified into something that looks like C to a C compiler but doesn't appeal to C programmers. The C++ versions of the Numerical Recipes books and code have eliminated the ugly 1-based-index-for-arrays stuff.

[/begin Opinion]
The old code may be usable (if you use it in the way that it was designed), but upgrading to Numerical Recipes Version 3 is recommended (if you can afford it). It's C++. In addition to new stuff, there are lots of changes (improvements) in the code and in many of the algorithms themselves way above and beyond vector and matrix notation.
[/end Opinion]


Statement of full disclosure:
I have no connection with the authors of the books or of the code, so there's nothing in it for me whether you upgrade or not.

I appreciate your help. The program works now. I will refer to the book "Numerical Recipes in C++" for future works.

Thank you very much again.
Best regards


If you are going to use the old version, you have two choices:


Play by the rules. Use the functions as listed and call them as the book tells you. See Footnote.


Make your own rules: convert all code from the index-1 based matrix and vector functions to index-0 base. Change "pointer to pointer" arguments to 2-D arrays if you want to.

However...

You have to do these things for all of the functions that are called by the functions that you call.


Here's an example showing how you can use the function from the book:

/*
Driver for Numerical Recipes Version 2 jacobi() function

davekw7x
*/
#include <stdio.h>
#include "nr.h"
#include "nrutil.h"

int main(void)
{
int i, j, nrot;

int n = 2;
/*
We need a Numerical Recipes-style "matrix" to use
as an argument to call the jacobi function.
Here are the matrix values stored in row-major
order in a 1-D array. This array name will be used
as an argument for the convert_matrix function.
*/
float a_array[] = {
5.0, -2.0, /* First row */
-2.0, 2.0 /* Second row */
};

float *d = vector(1, n);
float *r = vector(1, n);
float **v = matrix(1, n, 1, n);
float **a = convert_matrix(&a_array[0], 1, n, 1, n);

printf("Matrix:\n");
for (i = 1; i <= n; i++) {
for (j = 1; j <= n; j++) {
printf(" %+13.6e", a[i][j]);
}
printf("\n");
}
printf("\n");

jacobi(a, n, d, v, &nrot);

printf("Eigenvalues:\n");
for (i = 1; i <= n; i++) {
printf(" Number %d:%+14.6e\n", i, d[i]);
}
printf("\n");

printf("Eigenvectors:\n");
for (i = 1; i <= n; i++) {
printf(" Number %d:", i);
for (j = 1; j <= n; j++) {
printf("%+14.6e", v[j][i]);
}
printf("\n");
}
printf("\n");

printf("Number of Jacobi rotations = %d\n", nrot);

free_convert_matrix(a, 1, n, 1, n);
free_matrix(v, 1, n, 1, n);
free_vector(r, 1, n);
free_vector(d, 1, n);

return 0;
}


Output:

Matrix:
+5.000000e+00 -2.000000e+00
-2.000000e+00 +2.000000e+00

Eigenvalues:
Number 1: +6.000000e+00
Number 2: +1.000000e+00

Eigenvectors:
Number 1: +8.944272e-01 -4.472136e-01
Number 2: +4.472136e-01 +8.944272e-01

Number of Jacobi rotations = 1



Regards,

Dave

Footnote:
The old stuff in the "Numerical Recipes for C" consists mainly of Fortran functions that were transmogrified into something that looks like C to a C compiler but doesn't appeal to C programmers. The C++ versions of the Numerical Recipes books and code have eliminated the ugly 1-based-index-for-arrays stuff.

[/begin Opinion]
The old code may be usable (if you use it in the way that it was designed), but upgrading to Numerical Recipes Version 3 is recommended (if you can afford it). It's C++. In addition to new stuff, there are lots of changes (improvements) in the code and in many of the algorithms themselves way above and beyond vector and matrix notation.
[/end Opinion]


Statement of full disclosure:
I have no connection with the authors of the books or of the code, so there's nothing in it for me whether you upgrade or not.

Hello,

I am using NR in C++ second edition with printing corrected till software version 2.11. So maybe this issue is already corrected in one of the new versions.

In formula 11.1.18, tau is defined as s/(1+c). However when I do the derivation myself, I find tau=(1-c)/s. This can be easily found by substracting equation 11.1.4 from equation 11.1.16

Strangely enough, tau=s/(1+c) is adopted in the code and nobody complains about the outcome of the code.

Any thoughts/corrections?


Edit: It follows from s²+c²=1