Método de Relajación en lenguaje C

#include
#include
#include

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));
/* Asignacion de memoria dinamica*/
matriz=(double**)malloc(*dimension*sizeof(double*));
for(i=0;i<*dimension;i++)
{
matriz[i]=(double*)malloc((*dimension+1)*sizeof(double));
}
/* Lectura de elementos de la matriz*/
for(i=0;i<*dimension;i++)
{
for(j=0;j<*dimension;j++)
{
fscanf(f,"%lf\n",&(matriz[i][j]));
}
}
/* Leemos el vector de terminos independientes */
for(i=0;i<*dimension;i++)
{
fscanf(f,"%lf\n",&(matriz[i][*dimension]));
}
fclose(f);
return matriz;
}

/* Procedimiento para mostrar matrices */
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");
}
/* Procedimiento que multiplica matrices */
double** mulmat(double **A, double **B,int dimension)
{
int i,j,k;
double **C=(double**)malloc(dimension*sizeof(double*));

for(i=0;i {
C[i]=(double*)malloc((dimension)*sizeof(double));
}
for(i=0;i {
for(j=0;j {
C[i][j]=0;
}
}
for (i=0;i {
for (j=0;j {
for (k=0;k {
C[i][j]=C[i][j]+(A[i][k]*B[k][j]);
}
}
}
return (C);
}
/* funcion que construye la matriz de Jacobi Mj*/
double** MJ(double **matriz,int dimension)

{
int i,j;
double **mj,**L_U,**Dinv;
/* Tomamos memoria */
L_U=(double**)malloc(dimension*sizeof(double*));
Dinv=(double**)malloc(dimension*sizeof(double*));
for(i=0;i {
L_U[i]=(double*)malloc((dimension)*sizeof(double));
Dinv[i]=(double*)malloc((dimension)*sizeof(double));
}
/* calculamos inv(D), y L_U=-L-U*/
for(i=0;i {
for(j=0;j {
if (j==i)
{
Dinv[i][j]=1/matriz[i][j];
L_U[i][j]=0;
}
else
{
Dinv[i][j]=0;
L_U[i][j]=(-1)*matriz[i][j];
}
}
}
/* mj=inv(d)*(-L-U)*/
mj=mulmat(Dinv,L_U,dimension);
free(L_U);
free(Dinv);
return (mj);
}
double norma(double *vec,int dimension,int Tnorma)
{
int i,cero=0.00000000001;
double temp=0;
if (Tnorma==1)
{/*Norma ecuclidea*/
for(i=0;i {
temp+=vec[i]*vec[i];
}
return sqrt(temp);
}
else
{/*Norma infinita*/
temp=fabs(vec[0]);
for (i=1;i {
if ((fabs(vec[i])-fabs(temp)) > cero)
temp=fabs(vec[i]);
}
return temp;
}
}

/* Devuelve el signo del autovalor*/
double sign(double *v1, double *v2, int dimension)
{
int i;
double temp=0;

for (i=0; i {
temp = temp + v1[i] * v2[i];
}
if (temp>=0)
return (1);
else
return (-1);
}

/* Multiplica una matriz por un vector*/
double* MatxVec(double** matriz,double* vec,int dimension)
{
int i,j;
double *sol=(double*)malloc(dimension*sizeof(double));

for (i=0;i {
sol[i]=0;
for (j=0;j {
sol[i]= sol[i] + matriz[i][j]*vec[j];
}
}
return sol;
}

double Potencia(double **matriz,int dimension,double *u_n,double tolerancia,int Nmax)
{
double* temp=(double*)malloc(dimension*sizeof(double));
double* u_n_1=(double*)malloc(dimension*sizeof(double));
double Norma,Norma_Anterior=0,solucion;
int n=1,j;

Norma=norma(u_n,dimension,1);
while((n<=Nmax)&&(fabs(Norma-Norma_Anterior)>tolerancia))
{
/*Guardamos la iteracion anterior*/
for (j=0;j {
u_n_1[j]=u_n[j];
}
/*Dividimos entre la norma la iteracion actual*/
for (j=0;j {
u_n[j]=u_n[j]/Norma;
}
/* Y resolvemos la iteracion siguiente*/
temp = MatxVec(matriz,u_n,dimension);
free(u_n);
u_n=temp;
Norma_Anterior=Norma;
/* Hallamos el signo del autovalor maximo*/
solucion=Norma*sign(u_n, u_n_1, dimension);
n=n+1;
Norma=norma(u_n,dimension,1);
}
if (n>Nmax)
printf("Se ha superado el numero de iteraciones maximo.\n");
free(u_n_1);
free(temp);
return (solucion);
}

void relajacion(double **matriz,double *b,double *u_n,int dimension,int Nmax,double w,double Tol)
{
double dif,temp,maxdif=1.;
double *anterior;
int i,j,k,n=1;

anterior=(double*)malloc(dimension*sizeof(double));
while((n<=Nmax)&&(maxdif>Tol)) /*Bucle para calcular las aproximaciones hasta el valor numero de iteraciones*/
{
printf(" u%d[ ",n);
for (i=0;i {
printf(" %lf",u_n[i]);
anterior[i]=u_n[i]; /*guardamos en el vector anterior la posicion anterior del vector de aproximacion*/
}
printf("]\n");

/*Calculamos la siguiente iteracion*/
for (i=0;i {
temp=0;
for (j=0;j {
if (i!=j)
{
temp=temp-(matriz[i][j]*u_n[j]);
}
}
u_n[i]=w*( (temp + b[i]) / matriz[i][i])+((1-w)*u_n[i]); /*Aqui se almacena la siguiente aprox obtenida*/
}
/*Calculamos la mayor diferencia para el calculo del error*/
maxdif=fabs(u_n[0]-anterior[0]);
for (k=1;k {
dif=fabs(u_n[k]-anterior[k]);
if (dif>maxdif)
{
maxdif=dif;
}
}
n++;
}
if (n>Nmax)
printf("Se ha superado el numero de iteraciones maximo.\n");
free(anterior);
}


int main()
{
int dimension,Nmax,i,j;
double Tol,wopt;
double **matriz,**Mj,*u,*b;

printf("\t---------------------------------------------\n");
printf("\t| METODO DE RELAJACION |\n");
printf("\t---------------------------------------------\n");
printf("\n Introduzca tolerancia: ");
scanf("%lf",&Tol);
printf(" Introduzca el numero maximo de iteraciones: ");
scanf("%i",&Nmax);

matriz=LeerDatos("matriz4.dat",&dimension);
u=(double*)malloc(dimension*sizeof(double));/*Vector de aproximaciones iniciales*/
b=(double*)malloc(dimension*sizeof(double));/*Vector de terminos independientes*/
printf(" Terminos independientes: ");
printf("\n b[");
for(i=0;i {
b[i]=matriz[i][dimension];
printf(" %lf",b[i]);
}
printf(" ]\n");
printf("\n La matriz de coeficientes es: \n\n");

for (i=0;i {
printf("{ ");
for (j=0;j {
printf("%4.3f ",matriz[i][j]);
}
printf("}\n");
}

/*Aplicamos el Metodo de la Potencia para despues calcular el Wopt.\n");*/
printf("\n Calculamos Mj=ind(D)(-L-U): \n");
Mj=MJ(matriz,dimension); /* Calculamos Mj=inv(D)(-L-U)*/
printf(" La matriz Mj da: \n");
muestra(Mj,dimension);
printf(" Introduce el vector de aproximacion inicial para el metodo de la potencia:\n");
for (i=0;i {
printf(" U[%d]=",i);
scanf("%lf",&(u[i]));
}

wopt=Potencia(Mj,dimension,u,Tol,Nmax); /* Calculamos el autovalor maximo*/

wopt=(2/(1+sqrt(1-(wopt*wopt))));/* Calculamos el w optimo*/

printf(" El Wopt es: %lf \n",wopt);
printf(" Aplicamos Relajacion.\n");
printf(" Introduce el vector de aproximacion inicial para el metodo de Relajacion:\n");

u=(double*)malloc(dimension*sizeof(double));
for (i=0;i {
printf(" U[%d]=",i);
scanf("%lf",&(u[i]));
}
/* Aplicamos relajacion */
relajacion(matriz,b,u,dimension,Nmax,wopt,Tol); /*procedure q nos devuelve de salida u*/
printf("\n El resultado[");
for(i=0;i {
printf(" %lf",u[i]);
}
printf(" ]\n");
return 0;
}