Método de Jacobi en lenguaje C

//----------------------------------------------------------------------------//
//Autores.- Sergio Romero Sánchez & Alejandro E. Castro Robredo
//Asignatura.- Análisis Numérico
//Práctica 4.- Implementación del Método de Jacobi
//Especialidad.- I.T.I.G.
//----------------------------------------------------------------------------//

//Librerías necesarias para la implementación//
#include
#include
#include

//Para leer los datos de un fichero//
double** LeerDatos(char *nombrefichero,int *dimension)
{double **matriz;
int i,j;
FILE *f;
if(f=fopen(nombrefichero,"r"),!f)
{
printf("Problema con la lectura del fichero\n");
return matriz;
}
fscanf(f,"%d\n",&(*dimension));

//Asignación de memoria dinamica
matriz=(double**)malloc(*dimension*sizeof(double*));
for(i=0;i<*dimension;i++)
{
matriz[i]=(double*)malloc((*dimension)*sizeof(double));
}

//Lectura de elementos de la matriz
for(i=0;i<*dimension;i++) //Recorremos las Filas
{
for(j=0;j<*dimension;j++) //Recorremos las Columnas
{
fscanf(f,"%lf\n",&(matriz[i][j])); //Almacenamos el elemento en la
} //posición ij.
}
fclose(f);
return matriz;
}

//Procedimiento que muestra una matriz
void muestra(double **matriz,int Tam)
{int i,j;
printf("\n\t");
for(i=0;i {
for(j=0;j {
printf("%lf ",matriz[i][j]);
}
printf("\n\t");
}
printf("\n");
}

//Halla el mayor elemento en valor absoluto de la triangular superior
double mayor(double **matriz,int dimension,int *p, int *q)
{int i,j;
double temp;

temp=fabs(matriz[1][0]);*p=1;*q=0;
for(i=2;i { //de la diagonal principal y por encima de esta.
for(j=0;j { //de la diagonal principal y por encima de esta.
if(fabs(matriz[i][j])>temp)
{
temp=fabs(matriz[i][j]);
*p=i;*q=j;
}
}
}
return temp; //Devolvemos en temp el valor númerico en valor absoluto más
} //alto.
//Ejemplo de lo que hace la Función Mayor, con una matriz de 3x3
//(1,2,3
// 4,5,-6 //Los elementos de la diagonal Superior serían 2,3 y -6
// 7,8,9) //Al tratarse de valor absoluto el valor que nos dará será 6.


//Procedimiento de Jacobi
double ** Jacobi(double **autoval,double **autovec,int dim,double tol,int Nmax)
{int p,q,i,j,k,n=0;
double cot,sen,tan,cos,temp,max,ip,iq;
double zero=0.0000000000000001;
//Primero inicializamos la matriz de rotacion a la identidad
for (i=0;i {
for (j=0;j<=i;j++)
{
if (i==j) autovec[i][j]=1.;
else autovec[i][j]=autovec[j][i]=0.;
}
}

max=mayor(autoval,dim,&p,&q);
while ( (ntol) )
{// Condiciones de parada
cot = (autoval[q][q] - autoval[p][p]) / (2*autoval[p][q]);
if (cot >= zero)
tan = -cot + sqrt(1 + cot * cot);
else
tan = -cot - sqrt(1 + cot *cot);

cos = 1 / sqrt(1+tan*tan);
sen = tan * cos;
for (j=0; j { //Vamos rotando la matriz original, hasta que todos
//sus elementos sean ceros,excepto la diagonal
if ((j!=p) && (j!=q))
{
temp = autoval[p][j];
autoval[p][j]=autoval[p][j]*cos - autoval[q][j] * sen;
autoval[j][p]=autoval[p][j];
autoval[q][j] = temp * sen + autoval[q][j] *cos;
autoval[j][q]=autoval[q][j];
}
}//En la diagonal tendremos los autovalores
autoval[q][q] = autoval[q][q] + tan * autoval[p][q];
autoval[p][p] = autoval[p][p] - tan * autoval[p][q];
autoval[p][q] = 0.;
autoval[q][p] = 0.;


for (i=0; i { //Vamos calculando la nueva matriz de rotación
ip=cos*autovec[i][p]+sen*autovec[i][q];
iq=-sen*autovec[i][p]+cos*autovec[i][q];
autovec[i][p]=ip; //Las columnas seran los autovectores
autovec[i][q]=iq;
}
n++;
max=mayor(autoval,dim,&p,&q);
}
//Paramos porque hemos alcanzado el nº maximo de iteraciones
if (n==Nmax) printf("Hemos superado el Número Máximo\n");
return autoval;
}


int main()
{
double **matautovalores,**matautovectores;
double tol;
int dim,Nmax,i,j;

printf("\n\n ******* METODO DE JACOBI ***********");
printf("\n\n");
printf("Vamos a calcular los autovalores y autovectores de una matriz dada");
printf("\n\n");
printf("Introduzca la tolerancia: \n");scanf("%lf",&tol);
printf("Introduzca el numero maximo de iteraciones (Nmax): \n");scanf("%i",&Nmax);
printf("\n");
matautovalores=LeerDatos("matriz10.dat",&dim);
printf("La dimension de la matriz es %d\n",dim);
printf("\n");

matautovectores=(double**)malloc(dim*sizeof(double*));

printf("La matriz de los autovectores es la primera y la de los autovalores la segunda:\n");
for(i=0;i {
matautovectores[i]=(double*)malloc(dim*sizeof(double));
}
Jacobi(matautovalores,matautovectores,dim,tol,Nmax);

muestra(matautovalores,dim);
muestra(matautovectores,dim);

printf("Los autovalores son: \n\n");
for (i=0;i {
printf(" <%lf>",matautovalores[i][i]);
}

printf("\n");
printf("\n");
printf("Los autovectores son: \n\n");
for (i=0;i {
printf("<");
for (j=0;j {
printf("<%lf>",matautovectores[j][i]);
}
printf("> \n");
}
printf("\n");
printf("\n");
printf(" ************** FIN DEL PROGRAMA ************** ");

free (matautovalores);
free(matautovectores);
return 0;
}