7.4 Jacobi’s method 219
{
// Setting up the eigenvector matrix
for(inti = 0; i < n; i++ ){
for(intj = 0; j < n; j++ ){
if( i == j ){
R[i][j] = 1.0;
}else{
R[i][j] = 0.0;
}
}
}
intk, l;
doubleepsilon = 1.0e-8;
doublemax_number_terations = (double) n(double) n(double) n;
intiterations = 0;
doublemax_offdiag = maxoffdiag ( A, &k, &l, n );
while( fabs(max_offdiag) > epsilon && (double) iterations < max_number_iterations ){
max:offdiag = maxoffdiag ( A, &k, &l, n );
rotate ( A, R, k, l, n );
iterations++;
}
std::cout <<"Number of iterations: "<< iterations <<"\n";
return;
}
// Function to find the maximum matrix element. Can you figure out a more
// elegant algorithm?
doublemaxoffdiag (doubleA,intk,intl,intn )
{
doublemax = 0.0;
for(inti = 0; i < n; i++ ){
for(intj = i + 1; j < n; j++ ){
if( fabs(A[i][j]) > max ){
max = fabs(A[i][j]);
l = i;
k = j;
}
}
}
returnmax;
}
// Function to find the values of cos and sin
voidrotate (doubleA,double*R,intk,intl,intn )
{
doubles, c;
if( A[k][l] != 0.0 ){
doublet, tau;
tau = (A[l][l] - A[k][k])/(2A[k][l]);
if( tau > 0 ) {
t = 1.0/(tau + sqrt(1.0 + tautau);
}else{
t = -1.0/( -tau + sqrt(1.0 + tautau);
}
c = 1/sqrt(1+tt);
s = ct;
}else{
c = 1.0;
s = 0.0;