C/C++ - METODO DI JACOBI

Ciao a tutti! ho scritto questo codice del metodo di Jacobi per risolvere i sistemi lineari con metodo iterativo, però non riesco ad avere le soluzioni corrette. Mi potreste aiutare a capire dov'è l'errore? il criterio di fermata è corretto?
questo è il codice:

/*METODO DI JACOBI*/

#include<stdio.h>
#include<stdlib.h>
#include<math.h>

const int size=50;
typedef double matrice[size][size];
typedef double vettore[size];


void leggimatrice(matrice,int);
void leggivettore(vettore,int);
void jacobi(matrice,vettore,vettore,vettore,int);
double diagonaledominante(matrice,int);
double normamax(vettore, vettore, int);
void stampa(vettore,int);

main()
{
      matrice A; vettore b, vettk, vettkp1; double delta;
      
      int n, k=0; double toll, lambda, rho;
      
      printf("\n dammi la dimensione della matrice  n=");scanf("%d",&n);
      
      leggimatrice(A,n);
      
      lambda=diagonaledominante(A,n);//lamda=max-i(somm a[i][j])/a[i][i]
      
      if(lambda>=1)
                   {
                      printf("\n la matrice non e'a diagonale dominante!non puoi usare jacobi\n\n");
                      exit(1);
                   }
      
      leggivettore(b,n);
      
      printf("\n\n Inserire il limite di precisione delle soluzioni.        tolleranza=");
      scanf("%lf", &toll);
      
      for(int i=1; i<=n; i++)         vettk[i]=0;
      
      for(int i=1; i<=n; i++)    vettkp1[i]=0;
      
      delta=normamax(vettkp1, vettk, n);
      
      rho=(log(toll*((1-lambda)/delta))/log(lambda));
      
      do(k++);
      while(k<=rho);        
      
      while(delta>=toll)
             {
                        for(int j=1; j<n; j++)
                                {
                                        for(int i=1;i<=n;i++)   vettk[i]=vettkp1[i];
                   
                                        jacobi(A,b,vettk,vettkp1,n);
                                        
                                        delta=fabs(vettkp1[j]-vettk[j]);
                                 }
             }
     
     stampa(vettkp1,n);
     
     system("PAUSE");
     return 0;
}


void leggimatrice(matrice A,int n)
{
     printf("\n La matrice A e' costituita dai seguenti elementi:");
     
     for(int i=1;i<=n;i++)
             for(int j=1;j<=n;j++)
                     {
                         printf("\n Inserire l'elemento A[%d][%d]=",i,j);
                         scanf("%lf",&A[i][j]);
                     }
     return;
}


void leggivettore(vettore b,int n)
{
     for(int i=1;i<=n;i++)
             {
                  printf("\n Inserire il termine noto b[%d]=",i);
                  scanf("%lf",&b[i]);
             }
     return;
}


void jacobi(matrice A, vettore b,vettore x,vettore y,int n)
{
     for(int i=1;i<=n;i++)
             {
                  y[i]=b[i];
                  
                  for(int j=1;j<=n;j++)
                          if(j!=i)          y[i]=y[i]-A[i][j]*x[j];
                  
                  y[i]=y[i]/A[i][i];
             }
     return;
}



double diagonaledominante(matrice A,int n)
{
       double sum, lambda=0;
       
       for(int i=1;i<=n;i++)
               {
                   sum=0;
                   
                   for(int j=1;j<=n;j++)
                           if(i!=j)         sum=sum+(fabs(A[i][j]));
                           
                   sum=sum/fabs(A[i][i]);
                   
                   if(lambda<sum)           lambda=sum;
               }
  return lambda;
}


double normamax(vettore vettkp1, vettore vettk, int n)
{
       
       double max=(fabs(vettkp1[1]-vettk[1]));
       
       for(int i=2;i<=n;i++)
               if(max<(fabs(vettkp1[i]-vettk[i])))     max=(fabs(vettkp1[i]-vettk[i]));
               
       return max;
}

        
void stampa(vettore y,int n)
{
     printf("\n Le soluzioni del sistema sono:");
     
     for(int i=1;i<=n;i++)   printf("\t%lf",y[i]);
     
     return;
}

grazie mille