Método de Jacobi Mj_Cj 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));


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 *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 *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 jacobi2 (double **m, double *v, double *aprox, int dimension, int iter){
double **D1,**LR, *C;
int i,j;


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

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


/*HALLAMOS LA MATRIZ D INVERSA*/

for (i=0;i for(j=0;j if (j==i)
D1[i][j] = 1./m[i][j];
else
D1[i][j] = 0.;



/*HALLAMOS LA MATRIZ -L-U*/

for (i=0;i for(j=0;j if (j!=i)
LR[i][j] = -m[i][j];
else
LR[i][j] = 0.;


/*OBTENEMOS LA MATRIZ DE JACOBI Mj*/
liberaMatriz(m,dimension);
m = multiplica(D1, LR, dimension);
printf("\n La matriz Mj da como resultado :\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(D1,v,dimension);

printf("\n El vector C da como resultado :\n\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("\n\nEl vector solucion es:\n");
for(i=0;i printf("\nvector Un[%d] = %lf",i,aprox[i]);


}

int main(){
double **matriz;
double *v, *aprox, aux;
int dimension, iter, i, j;

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

printf("\t---------------------------------------------\n");
printf("\t| METODO DE GAUSS-JACOBI_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);
jacobi2 (matriz, v ,aprox, dimension, iter);
/* destruimos los vectores de punteros */
/*la matriz ya la libero arriba en gauss*/
free(v);
free(aprox);
system ("PAUSE");
return 0;
}