#include <stdio.h>
#include <math.h>
/* jacobie.c : 1999 3 8 by Oguni */
#define N 3
#define pi 3.14159265358979
main(){
  double A[N][N]={3.,2.,1.,2.,3.,2.,1.,2.,3.};
  double tol=0.00001;
  int n=3, itr[2]={20, 0}, i;
  jacobie(A,tol,n,itr);
  printf("iteration=%d\n", itr[1]);
  for (i=0; i<n; i++){
    printf("%f\n", A[i][i]); }
}
/* function jacobie */
jacobie(A,tol,n,itr)
double A[N][N], tol;
int n, itr[2];
{
   double v=1000000.0, vmax, t, c, s, akk, all, akj, alj;
   int i, j, k, l, m;
   for (m=1; m<=itr[0]; m++){
     vmax=0.0;
     for (j=0; j<n; j++){
       for (i=0; i<n; i++){
         if (vmax<fabs(A[i][j]) && i!=j){
           vmax=fabs(A[i][j]); k=i; l=j; }
       }
     }
     v=vmax;
     if (v<tol) return;
     if (A[k][k]==A[l][l]) t=pi/4.0;
     else t=atan(2.*A[k][l]/(A[k][k]-A[l][l]))/2.0;
     c=cos(t); s=sin(t);
     akk=A[k][k]*c*c+A[l][l]*s*s+2.0*A[k][l]*c*s;
     all=A[k][k]*s*s+A[l][l]*c*c-2.0*A[k][l]*c*s;
     A[k][k]=akk; A[k][l]=0.; A[l][k]=0.; A[l][l]=all;
     for (j=0; j<n; j++){
       if (j!=k && j!=l){
         akj=A[k][j]*c+A[l][j]*s;
         alj=-s*A[k][j]+c*A[l][j];
         A[j][l]=alj; A[l][j]=alj; A[k][j]=akj; A[j][k]=akj; }
     }
     itr[1]=m;
   }
   if (m>itr[0]) printf("non-convergence\n");
   return;
 }