Método de Gauss-Seidel en Lenguaje C (Mj_CJ)

#include
#include
#include
#define TOL 0.00001


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));


for(i=0;i<*dimension;i++)
for(j=0;j<*dimension;j++)
fscanf(f,"%lf\n",&(matriz[i][j]));


for(i=0;i<*dimension;i++)
fscanf(f,"%lf\n",&(matriz[i][*dimension]));

fclose(f);
return matriz;
}


void liberaMatriz(double **matriz,int dimension){
int i;

for(i=0;i free(matriz[i]);
free(matriz);
}


double **multiplica(double **A,double **B,int dimension){
int i,j,k;
double **C;
C=(double**)malloc(dimension*sizeof(double*));
for(i=0;i C[i]=(double*)malloc((dimension+1)*sizeof(double));

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;
}

double **matrizAmpliada(double **matriz,double *vector,int dimension){
int i,j;
double **C;
C=(double**)malloc(dimension*sizeof(double*));
for(i=0;i C[i]=(double*)malloc((dimension+1)*sizeof(double));

for (i=0;i for (j=0;j C[i][j]=matriz[i][j];

for (i=0;i C[i][dimension]=vector[i];
return C;
}
double *matrizPorVector(double **matriz,double *vector,int dimension){
double *resultado,a,b,mult;
int i,j;

resultado=(double*)malloc(dimension*sizeof(double));

for(i=0;i mult=0;
for(j=0;j mult+=(matriz[i][j]*vector[j]);
resultado[i]=mult;
}

return resultado;
}


double* resolverSistemaGauss(double **A,double *b,int dimension){ /*funcion q la usa inversa*/
double M,t,*X,suma,**aux;
int i,j,p,k,c,*Fila;

Fila=(int*)malloc(dimension*sizeof(int*));
X=(double*)malloc(dimension*sizeof(double));


aux=(double**)malloc(dimension*sizeof(double*));
for(i=0;i aux[i]=(double*)malloc((dimension+1)*sizeof(double));

for(i=0;i for(j=0;j aux[i][j]=A[i][j];

aux=matrizAmpliada(aux,b,dimension);

for(j=0;j Fila[j]=j;

for(p=0;p for(k=p+1;k if(fabs(aux[Fila[k]][p])>(fabs(aux[Fila[p]][p]))){
t=Fila[p];
Fila[p]=Fila[k];
Fila[k]=t;
}

if(fabs(aux[Fila[p]][p]) < TOL){
printf("\n\n\"resolverSistemaGauss\" informa:\n");
printf("La matriz es singular\n");
return NULL;
}

for(k=p+1;k M=(aux[Fila[k]][p])/(aux[Fila[p]][p]);
for(c=p;c<=dimension;c++)
aux[Fila[k]][c]=(aux[Fila[k]][c])-(M*aux[Fila[p]][c]);
}
}

if(fabs(aux[Fila[dimension-1]][dimension-1]) < TOL){
printf("\n\n\"resolverSistemaGauss\" informa:\n");
printf("La matriz es singular\n");
return NULL;
}

X[dimension-1]=aux[Fila[dimension-1]][dimension]/(aux[Fila[dimension-1]][dimension-1]);

for(k=dimension-2;k>=0;k--){
suma=0.;
for(c=dimension-1;c>k;c--)
suma+=X[c]*aux[Fila[k]][c];
X[k]=(aux[Fila[k]][dimension]-suma)/aux[Fila[k]][k];
}
return X;
}


double **inversa(double **A,int dimension){
double **invA,*C,**aux,*b;
int i,j;

b=(double*)malloc((dimension)*sizeof(double));

invA=(double**)malloc(dimension*sizeof(double*));
for(i=0;i invA[i]=(double*)malloc((dimension)*sizeof(double));

aux=(double**)malloc(dimension*sizeof(double*));
for(i=0;i aux[i]=(double*)malloc((dimension+1)*sizeof(double));
for(i=0;i for(j=0;j aux[i][j]=A[i][j];

for(i=0;i for(j=0;j if (j==i)
b[j]=1;
else
b[j]=0;
}

if (C=resolverSistemaGauss(aux,b,dimension))
for(j=0;j invA[j][i]=C[j];
else
return NULL;
}
liberaMatriz(aux,dimension);
return invA;
}


double *suma_vector(double *A,double *B,int dimension){

int i;
double *C;

C=(double*)malloc(dimension*sizeof(double));

for (i=0;i C[i] = A[i]+B[i];

return C;
}

void Gauss_Seidel_2 (double **m, double *v, double *aprox, int dimension, int iter){
double **DL, **U, *C;
int i,j;


DL=(double**)malloc(dimension*sizeof(double*));
for(i=0;i DL[i]=(double*)malloc((dimension)*sizeof(double));

U=(double**)malloc(dimension*sizeof(double*));
for(i=0;i U[i]=(double*)malloc((dimension)*sizeof(double));


/*Hallamos la matriz D + L */

for (i=0;i for(j=0;j if (j>i)
DL[i][j] = 0.;
else
DL[i][j] = m[i][j];


/*Hallamos la inversa de D+L*/

DL = inversa(DL, dimension);

/*Obtenemos la matriz -U*/

for (i=0;i for(j=0;j if (j>i)
U[i][j] = -m[i][j];
else
U[i][j] = 0.;


/*OBTENEMOS LA MATRIZ DE GAUSS SEIDEL Mj*/

liberaMatriz(m,dimension);
m = multiplica(DL, U, dimension);
printf("\La matriz Mgs da: \n\n");

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



/*CALCULAMOS EL VECTOR Cj*/
C=matrizPorVector(DL,v,dimension);
printf("\nEl vector C da: \n");
for (j=0;j {
printf(" %4.3f ",C[j]);
}


for(i=0;i aprox = matrizPorVector(m,aprox,dimension);
aprox = suma_vector(aprox, C, dimension);}

printf("\nEl vector solucion para ese numero de iteraciones es:\n");
for(i=0;i printf("\nvector Un[%d] = %lf",i,aprox[i]);
liberaMatriz(m,dimension);
liberaMatriz(U,dimension);
liberaMatriz(DL,dimension);

}


int main(){



double **matriz;
double *v, *aprox, aux;
int dimension, iter, i, j;

matriz=LeerDatos("matriz4.dat",&dimension);

printf("\t---------------------------------------------\n");
printf("\t| METODO DE GAUSS-SEIDEL_MC |\n");
printf("\t---------------------------------------------\n");

v = (double*)malloc(dimension*sizeof(double));
aprox = (double*)malloc(dimension*sizeof(double));

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");
}

for (i=0;i v[i] = matriz[i][dimension];

printf("\nIntroduzca la aproximacion inicial:\n");
for(i=0;i scanf("\n%f", &aux);
aprox[i] = aux;}


printf("\nIntroduzca el numero de interaciones deseado:\n");
scanf("%d",&iter);


Gauss_Seidel_2 (matriz, v ,aprox, dimension, iter);


/* destruimos los vectores de punteros */
/*la matriz ya la libero arriba en gauss*/
free(v);
free(aprox);

return 0;

}