Método de Cholesky en Lenguaje C

//Principal.c
//Autores.- Sergio Romero Sánchez & Alejandro E. Castro Robredo
//Práctica Método de Cholesky.-
//Especialidad.- I.T.I.G.
//Fecha.- 18/01/2005

#include
#include
#include "cholesky.h"
#include "Leer.h"
int main()
{
double **A,*x;
int n,i,j;

A=LeerDatos("Dto2.dat",&n);
if ( (x=Cholesky(A,n))!=NULL)
{
printf ("\nEstos son los resultados: \n");
for(i=0;i printf("\n X%d = %.14lf",i,x[i]);
free(x);
printf ("\n");
}
else
{
printf ("No se ha podido resolver el sistema \n");
}
for (i=0;i {
free(A[i]);
}
free(A);
return 1;
}




//Leer.h
double** LeerDatos(char *nombrefichero,int *dimension);

//Leer.c

#include
#include "Leer.h"
#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;
}



//Cholesky.h

int Diagonalizar(double **A,double **B,int dimension);
void MatrizTraspuesta(double **A,double **At,int dimension);
void DescensoCholesky(double **B,double *y,double *b,int dimension);
void AscensoCholesky(double **Bt,double *x,double *y,int dimension);
int Simetrica (double **A,int n);
double** LeerDatos(char *nombrefichero,int *dimension);
double* Cholesky (double **A, int dimension);

//Cholesky.c

#include "cholesky.h"
#include
#include
#include

int Diagonalizar(double **A,double **B,int dimension)
/* Ambas matrices deben ya disponer de espacio de memoria dinamica
// A = matriz original de coeficientes
// B = triangular inferior (los de arriba todos a cero) despues de
aplicar Cholesky
// dimension = dimension de la matriz cuadrada*/
{
int i,j,k;
double suma,producto,z;

for(j=0;j {
for (i=j;i {
if(i==j)
{
suma = 0.;
for(k=0;k {
producto = B[j][k]*B[j][k];
suma = suma + producto;
}
if ((z=A[j][j]-suma)<=0)
{
return 0;
}
else
{
B[j][j] = sqrt(z);
}
}
else
{
suma = 0.;
for(k=0;k {
producto = B[i][k]*B[j][k];
suma = suma + producto;
}
B[i][j] = (double)(A[i][j]-suma)/B[j][j];
}
printf("\nB[%d][%d] = %.14lf",i,j,B[i][j]);
}
}
return 1;
}

void MatrizTraspuesta(double **A,double **At,int dimension)
/* Ambas matrices deben ya disponer de espacio de memoria dinamica
// A = matriz original
// At = matriz traspuesta de la matriz A
// dimension = dimension de la matriz cuadrada*/
{
int i,j;

for (j=0;j {
for(i=0;i At[i][j] = A[j][i];
}
return;
}


void DescensoCholesky(double **B,double *y,double *b,int dimension)
/*Ambas matrices deben ya disponer de espacio de memoria dinamica
// B = matriz B
// y = soluciones del descenso
// b = vector de terminos independientes
// dimension = dimension de la matriz cuadrada*/
{
int i,j;
double suma;

for(i=0;i {
suma = 0.;
for(j=0;j {
suma = suma + (B[i][j]*y[j]);
}
y[i] = (double)(b[i]-suma)/B[i][i];
}

for(i=0;i<4;i++)
printf("\ny%d = %.14lf",i,y[i]);

return;
}

void AscensoCholesky(double **Bt,double *x,double *y,int dimension)
/*Ambas matrices deben ya disponer de espacio de memoria dinamica
// Bt = matriz traspuesta de B
// x = soluciones del sistemas
// y = soluciones del descenso
// dimension = dimension de la matriz cuadrada */
{
int i,j;
double suma;

for(i=(dimension-1);i>=0;i--)
{
suma = 0.;
for(j=i;j {
suma = suma + (Bt[i][j]*x[j]);
}
x[i] = (double)(y[i]-suma)/Bt[i][i];
}
return;
}
int Simetrica (double **A,int n)

/* A=Matriz de coeficientes si es simétrica devuelve 1, en caso contrario
devuelve 0 */
{
int i,j;

for (i=0;i {
for (j=i+1;j {
if (A[i][j]!=A[j][i])
return 0;
}
}
return 1;
}

double* Cholesky (double **A, int dimension)
{
/* Matriz A contiene coefcientes y tréminos independientes
dimension de la matriz
Si A cumple condiciones de Cholesky entonces
Devuelve X con la solucion del sistema
en caso contrario
Devuelve NULL */

double *b,**B,**Bt,*x,*y;
int h,i;
if (Simetrica(A,dimension)==1)
{
b=(double*)malloc(dimension*sizeof(double));
B = (double**)malloc(dimension*sizeof(double*));
for(h=0;h B[h] = (double*)malloc(dimension*sizeof(double));

Bt = (double**)malloc(dimension*sizeof(double*));
for(h=0;h Bt[h] = (double*)malloc(dimension*sizeof(double));
x = (double*)malloc(dimension*sizeof(double));
y = (double*)malloc(dimension*sizeof(double));
for(i=0;i b[i]=A[i][dimension];

if (Diagonalizar(A,B,dimension)==1)
{
MatrizTraspuesta(B,Bt,4);
DescensoCholesky(B,y,b,dimension);
AscensoCholesky(Bt,x,y,dimension);

for (i=0;i {
free(B[i]);free(Bt[i]);
}
free(y);
free(B);
free(Bt);
free(b);
return x;
}
else
{
printf ("La matriz de coeficientes no es definida positiva,y por
tanto no se puede resolver a través de Cholesky \n");
for (i=0;i {
free(B[i]);free(Bt[i]);
}
free(y);free(B);
free(Bt);free(b);
return NULL;
}}
else
{
printf ("La matriz de coeficientes no es simétrica,y por tanto no
se puede resolver a través de Cholesky \n");
return NULL;
}}