Método de Crout en lenguaje C

//principal.c

#include
#include
#include
#include "matrices.h"
#include "crout.h"

void captura_matriz(double **mat,int nfil,int ncol)
{
int i,j;
double valor;

for (i=0;i {
for (j=0;j {
printf("introduzca el valor[%d,%d]=",i,j);
scanf("%f",&valor);
mat[i][j]=valor;

}
}}
void muestra_matriz (double **mat,int nfil,int ncol)
{
int i,j;

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

{
double **matriz,**matriz_u,**matriz_l,*vect,*termindp;
FILE *fich;
int i,j,tam;
/*apertura de fichero, devolviendo el tamaño de la matriz*/
tam = abrirFicheroMatriz("matriz4.dat",&fich);

/*se reserva para las matrices vectores de punteros a double que serían las filas (i)*/
matriz_u = (double**)malloc(tam*sizeof(double*));
matriz_l = (double**)malloc(tam*sizeof(double*));
vect = (double*)malloc(tam*sizeof(double));
termindp = (double*)malloc(tam*sizeof(double));
/*una vez reservadas las filas se reservan para cada puntero del vector un vector que seran
las columnas*/
for(i=0;i matriz_u[i]=(double*)malloc(tam*sizeof(double));
matriz_l[i]=(double*)malloc(tam*sizeof(double));
}
printf("\t---------------------------------------------\n");
printf("\t| METODO DE CROUT |\n");
printf("\t---------------------------------------------\n");

/*Se procedera en primer lugar a descomponer la matriz A, en matriz B y matriz B'\n");
usando para ello el metodo de Crout.\n");
Se realizara el descenso y el remonte de las misma, obtiendo asi las soluciones del sistema\n\n");*/


matriz=leerDatos(fich,tam);
/*se lee la matriz A y se muestra*/
printf("La matriz de coeficientes es: \n ");
muestra_matriz(matriz,tam,tam);

crout(matriz,matriz_u,matriz_l,tam);
/*una vez hallado B y B' o L y U son mostradas por pantalla*/

printf("\n\n El metodo de Crout concluido con exito \n");
printf("\nEsta son las dos matrices obtenidas : ");

printf("\nMatriz B' : \n ");
muestra_matriz(matriz_l,tam,tam);
printf("\nmatriz B: \n ");
muestra_matriz(matriz_u,tam,tam);
/*vector de terminos independientes*/

printf("Vector de terminos independientes :\n");
for(i=0;i {
printf("\n");
termindp[i]=matriz[i][tam];
printf("( %4.3lf )", termindp[i]);
}
printf("\n");
vect =descenso(matriz_l,termindp,tam);
vect =remonte(matriz_u,vect,tam);
printf("La resolucion del sistema de ecuaciones a concluido con exito, las soluciones son: ");
printf("\n");
for(i=0;i printf("( %4.3lf )", vect[i]);
printf("\n");
for(i=0; i {
free(matriz[i]);
free(matriz_l[i]);
free(matriz_u[i]);
}

/* destruimos losvectorres de punteros */
free(matriz);
free(matriz_l);
free(matriz_u);
free(vect);
free(termindp);
return 0;
}



//matrices.h
#ifndef MATRICES_DEF
#define MATRICES_DEF

/*--------------------------------------------------------------------------------*/
/* devuelve la matriz resultante de A*B */
double **multiplica(double **A,double **B,int dimension);
/*--------------------------------------------------------------------------------*/


/*--------------------------------------------------------------------------------*/
/* devuelve la matriz resultante de A*x, siendo x un vector */
double *matrizPorVector(double **matriz,double *vector,int dimension);
/*--------------------------------------------------------------------------------*/


/*--------------------------------------------------------------------------------*/
/* devuelve la matriz resultante de A+B */
double **suma(double **A,double **B,int dimension);
/*--------------------------------------------------------------------------------*/


/*--------------------------------------------------------------------------------*/
/* devuelve la matriz resultante de A-B */
double **resta(double **A,double **B,int dimension);
/*--------------------------------------------------------------------------------*/


/*--------------------------------------------------------------------------------*/
/* Le resta a la diagonal el autovalor pasado */
/* si se le pasa la matriz ampliada solo devulve la matriz */
/* sin la columna de terminos independientes */
double **matrizMenosAutovalor(double **matriz,double autovalor,int dimension);
/*--------------------------------------------------------------------------------*/


/*--------------------------------------------------------------------------------*/
/* devuelve el vector de terminos independientes de la matriz ampliada */
double *terminosIndependientes(double **matriz,int dimension);
/*--------------------------------------------------------------------------------*/


/*--------------------------------------------------------------------------------*/
/* crea una matriz ampliada pasandole una matriz y */
/* el vector de terminos independientes */
double **matrizAmpliada(double **matriz,double *vector,int dimension);
/*--------------------------------------------------------------------------------*/


/*--------------------------------------------------------------------------------*/
/* devuelve la matriz traspuesta de una matriz de entrada */
double **traspuesta(double **A,int dimension);
/*--------------------------------------------------------------------------------*/


/*--------------------------------------------------------------------------------*/
/* devuelve el vector de soluciones obtenidas al realizar
el descenso en una matriz triangular superior */
double *remonte(double **B,double *b,int dimension);
/*--------------------------------------------------------------------------------*/


/*--------------------------------------------------------------------------------*/
/* devuelve el vector de soluciones obtenidas al realizar
el descenso en una matriz triangular inferior */
double *descenso(double **B,double *b,int dimension);
/*--------------------------------------------------------------------------------*/


/*--------------------------------------------------------------------------------*/
/* comprueba si una matriz es simetrica */
/* 1: SI es simetrica */
/* 0: NO es simetrica */
int simetrica(double **A,int dimension);
/*--------------------------------------------------------------------------------*/


/*--------------------------------------------------------------------------------*/
/* comprueba si dos matrices son iguales */
/* 1: SI son iguales */
/* 0: NO son iguales */
int iguales(double **A,double **B,int dimension);
/*--------------------------------------------------------------------------------*/


/*--------------------------------------------------------------------------------*/
/* comprueba si una matriz es tridiagonal */
/* 1: SI es tridiagonal */
/* 0: NO es tridiagonal */
int tridiagonal(double **A,int dimension);
/*--------------------------------------------------------------------------------*/


/*--------------------------------------------------------------------------------*/
/* imprime por pantalla la matriz A */
/* La dimension maxima para poder imprirse es 6 */
void imprimirMatrizPantalla(double **A,int dimension);
/*--------------------------------------------------------------------------------*/


/*--------------------------------------------------------------------------------*/
/* imprime a un fichero la matriz C (maximo dimension 10) */
/* requiere que el fichero haya sido abierto previamente */
void imprimirMatrizFichero(FILE *fichero,double **C,int dimension);
/*--------------------------------------------------------------------------------*/

/*--------------------------------------------------------------------------------*/
/* Analogo que imprimirMatrizPantalla pero con un vector */
void imprimirVectorPantalla(double *vector,int dimension);
/*--------------------------------------------------------------------------------*/


/*--------------------------------------------------------------------------------*/
/* Analogo que imprimirMatrizFichero pero con un vector */
void imprimirVectorFichero(FILE *fichero,double *vector,int dimension);
/*--------------------------------------------------------------------------------*/


/*--------------------------------------------------------------------------------*/
/* devuelve la matriz B de la descomposicion A=B*Bt */
/* comprueba que la matriz A sea simetrica, devolviendo
un NULL en caso de no ser simetrica */
double **descomposicionCholesky (double **A,int dimension);
/*--------------------------------------------------------------------------------*/


/*--------------------------------------------------------------------------------*/
/* lee de un fichero de texto una matriz */
/* la matriz tiene que venir de la siguiente manera */
/* en la primera linea la dimension, despues en cada linea
se van introduciendo los elementos de las distintas filas
ordenadamente. Al final y tambien en cada linea los terminos
independientes; Ejemplo:
2
1
0
0
1
1
1
Esta seria la matriz identidad de orden dos y con el vector (1,1) como
terminos independientes. */
double** leerDatos(FILE *f,int dimension);
/*--------------------------------------------------------------------------------*/


/*--------------------------------------------------------------------------------*/
/* abre un fichero que contiene una matriz con la representacion
anterior devolviendo:
la dimension si pudo abrir el fichero satisfactoriamente
0 en caso contrario */
int abrirFicheroMatriz(char *nombreFichero,FILE **fentrada);
/*--------------------------------------------------------------------------------*/


/*--------------------------------------------------------------------------------*/
/* crea un fichero con una matriz con la representacion comentada en leerDatos
para luego poder se leida por leerDatos. Como parametros requiere la matriz,
el vector de terminos independientes y la dimension de la matriz.
Devuelve:
1: pudo crear el fichero correctamente
0: no pudo crear el fichero */
int crearFicheroMatriz(double **matriz,double *vector,int dimension);
/*--------------------------------------------------------------------------------*/


/*--------------------------------------------------------------------------------*/
/* elimina el espacio reservado para una matriz */
void eliminaMatriz(double **matriz,int dimension);
/*--------------------------------------------------------------------------------*/


/*--------------------------------------------------------------------------------*/
/* Esta funcion se apoya en la funcion descomposicionCholesky para obtener la
/* matriz B y luego resolver el sistema B*z=b mediante un remonte. La solucion
/* de este sistema la utilizamos para resolver el sistema Bt*u=z mediante un
/* descenso. Requiere que la matriz sea definida positiva y simetrica */
double *resolverSistemaCholesky(double **A,double *vector,int dimension);
/*--------------------------------------------------------------------------------*/


double* resolverSistemaGauss(double **A,double *b,int dimension);

double* resolverSistemaGaussSeidel(double **A,double *b,int dimension,int *Nmax);

double *autovaloresJacobi(double **matriz,int dimension,int *Nmax);

double determinante(double **matriz,int dimension);


/*--------------------------------------------------------------------------------*/
/* Calcula la inversa de una matriz. Para ello resuelve tantos sistemas como sea
/* la dimension de entrada. La peculiaridad de estos sistemas es que tiene como
/* terminos independientes las columnas de la matriz identidad. Los vectores obtenidos
/* como soluciones de estos sistemas seran las columnas de la matriz inversa. Para
/* resolver los distintos sistemas utilizamos la funcion resolverSistemaGauss. Si
/* Gauss no es capaz de resolver el sistema devuelve NULL */
double **inversa(double **A,int dimension);
/*--------------------------------------------------------------------------------*/


/*--------------------------------------------------------------------------------*/
/* calcula el autovalor maximo de la matriz de entrada que en valor absoluto es */
/* el radio espectral de la matriz de entrada */
double potencia(double **matriz, double *x,double *autovector,int dimension,int *Nmax);
/*--------------------------------------------------------------------------------*/


/*--------------------------------------------------------------------------------*/
/* resuelve un sistema de ecuaciones por el metodo de la relajacion, que es muy
/* parecido al metodo de Gauss-Seidel pero introduciendo un parametro de
/* relajacion w. Para el calculo de este w nos apoyamos en el metodo de la
/* potencia para el calculo del autovalor maximo. Para ello es imprescindible
/* que la matriz de entrada sea tridiagonal. En caso de no ser tridiagonal
/* resolvemos el sistema con una w predefinida (ver implementacion) */
double* resolverSistemaRelajacion(double **A,double *b,int dimension,int *Nmax);
/*--------------------------------------------------------------------------------*/

double potenciaInversa(double **matriz, double *uInicial,double *autovector,int dimension,int *Nmax);

#endif

//matrices.c

#include
#include
#include
#include
#include "matrices.h"
#define TOL 0.00000001

int iguales(double **A,double **B,int dimension){
int i,j;
for(i=0;i for(j=0;j if(fabs((A[i][j])-(B[i][j])) > TOL)
return 0;
return 1;
printf("\nLas matrices estan bien");
}

int abrirFicheroMatriz(char *nombreFichero,FILE **fentrada){
int dimension;
if(*fentrada=fopen(nombreFichero,"rb"),!*fentrada){
printf("\n\n\"abrirFicheroMatriz\" informa:\n");
printf("ERROR: No se ha podido abrir el fichero \"%s\"",nombreFichero);
return 0;
}
fscanf(*fentrada,"%d\n",&dimension);
return dimension;
}

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** leerDatos(FILE *f,int dimension){
double **matriz;
int i,j;

/* Asignacion de memoria dinamica*/
matriz=(double**)malloc(dimension*sizeof(double*));
for(i=0;i matriz[i]=(double*)malloc((dimension+1)*sizeof(double));

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

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

fclose(f);
return matriz;
}

int tridiagonal(double **A,int dimension){
int i,j;
for(i=2;i for(j=0;j if((fabs(A[j][i])+fabs(A[i][j])) > TOL)
return 0;
return 1;
}

double *descenso(double **B,double *b,int dimension){
int i,j;
double suma,*vec;
vec=(double*)malloc(dimension*sizeof(double));
vec[0]=b[0]/(B[0][0]);
for(i=1;i suma=0.;
for(j=0;j suma+=vec[j]*B[i][j];
vec[i]=(b[i]-suma)/B[i][i];
}
return vec;
}

double *remonte(double **B,double *b,int dimension){
int i,j;
double suma,*vec;
vec=(double*)malloc(dimension*sizeof(double));
vec[dimension-1]=b[dimension-1]/(B[dimension-1][dimension-1]);
for(i=dimension-2;i>=0;i--){
suma=0.;
for(j=dimension-1;j>i;j--)
suma+=vec[j]*B[i][j];
vec[i]=(b[i]-suma)/B[i][i];
}
return vec;
}

//crout.h

#ifndef CROUT
#define CROUT

/**************METODO DE FACTORIZACION DE CROUT*****************/

void crout(double **mat,double **matriz_u,double **matriz_l,int tam);
/*procedimiento que dada una matriz y un tamaño devuelve dos matrices triangulares que su multiplicacion dan mat, dado que crout es un metodo de factorizacion*/
#endif


//crout.c

void crout(double **mat,double **mat_u,double **mat_l,int tam)
{

int i,j,k;


for (i=0; i {
for(j=0;j {
if (i==j)
mat_u[i][j]=1.;
else
mat_u[i][j]=0.;

mat_l[i][j]=0.;
}
}

mat_l[0][0]=mat[0][0];
mat_u[0][1]=(mat[0][1]/mat_l[0][0]);
mat_l[1][0]=mat[1][0];

/*Se rellenan las matrices*/

for(i=1;i {
mat_l[i][i]=mat[i][i]-(mat_l[i][i-1]*mat_u[i-1][i]);
mat_u[i][i+1]=(mat[i][i+1]/mat_l[i][i]);
mat_l[i+1][i]=mat[i+1][i];
}
mat_l[tam-1][tam-2]=mat[tam-1][tam-2];
mat_l[tam-1][tam-1]=mat[tam-1][tam-1]-mat_l[tam-1][tam-2]*mat_u[tam-2][tam-1];
}