Thread: Cholesky decomposition example

Hello guys,

i have the cholesky decomposition to solve an equation Ax=b
I dont really understand it so can someone show me on an example
like A=
1,1,0,0
1,1,1,0
0,1,1,1
0,0,1,1
and
(1,1,0,0)=b
how it works? i cant get the main program for it so i guess a simple example will help me understand a lot.

here are the subroutines

Code:

void nrerror(char error_text[])
/* Numerical Recipes standard error handler */
{
    fprintf(stderr,"Numerical Recipes run-time error...\n");
    fprintf(stderr,"%s\n",error_text);
    fprintf(stderr,"...now exiting to system...\n");
    exit(1);
}

void choldc(float **a, int n ,float p[])
{
    void nrerror(char error_text[]);
    int i,j,k;
    float sum;

    for (i=1;i<=n;i++) {
        for (j=i;j<=n;j++) {
            for (sum=a[i][j],k=i-1;k>=1;k--) sum -= a[i][k]*a[j][k];
            if (i == j) {
                if (sum <= 0.0)
                    nrerror("choldc failed");
                p[i]=sqrt(sum);
            } else a[j][i]=sum/p[i];
        }
    }
}

Code:

void cholsl(float **a, int n, float p[], float b[], float x[])
{
    int i,k;
    float sum;

    for (i=1;i<=n;i++) {
        for (sum=b[i],k=i-1;k>=1;k--) sum -= a[i][k]*x[k];
        x[i]=sum/p[i];
    }
    for (i=n;i>=1;i--) {
        for (sum=x[i],k=i+1;k<=n;k++) sum -= a[k][i]*x[k];
        x[i]=sum/p[i];
    }
}

This are the functions i got from numerical recipes and should be correct. I try to understand them by example, like a simple main programm which takes these 2 functions to solve the equation.
Take a matrix and its dimension into the choldc and print result, then take the results and a vector b and use it in cholsl and print the final result.

there are functions from nemerical recipes and should be correct. i basically want to understand them by looking at a testing programm. it shouldnt be hard but i somehow dont get a correct main

Code:

int main{

int n = 4 int a[]  = { 1, 1, 0, 0, 
             1, 1, 1, 0,    
             0, 1, 1, 1,    
             0, 0, 1, 1,   
            };

choldc(&a, n ,p[]);
printing(&a,n)


}

Code:

int printing(int *a, int n)
{
    int i, j;
    for (i = 0; i < n; i++)
    {
        for(j = 0; j < n; j++)
        {
            printf("%d\t", mat[j*n+i]);
        }
        printf("\n");
    }
}

with a print function

sry there was some major nonsense.

Code:

    int n=4;
    float L,p;
    int line;
    int row;
    float a[line][row];

a[1][1]=3;   a[1][2]= -2; a[1][3]=0;     a[1][4]=0;
a[2][1]= -2; a[2][2]=3;   a[2][3]= 6;    a[2][4]=0;
a[3][1]=0;   a[3][2]=6;   a[3][3]=3;     a[3][4]= -2;
a[4][1]=0;   a[4][2]=0;   a[4][3]= -2;   a[4][4]=3;



    L = choldc(&a, n, p[]);
    printf("  L %\n ", L);

    return 0;

Later i need to use this result and a vector for the second function and get my solution

Code:

int n=4;
    float L,p;
    int line =4;

    int row =4;        

    float a[line][row];
 
a[1][1]=3;   a[1][2]= -2; a[1][3]=0;     a[1][4]=0;
a[2][1]= -2; a[2][2]=3;   a[2][3]= 6;    a[2][4]=0;
a[3][1]=0;   a[3][2]=6;   a[3][3]=3;     a[3][4]= -2;
a[4][1]=0;   a[4][2]=0;   a[4][3]= -2;   a[4][4]=3;
 
 
 
    L = choldc(a, n, p[]);

    printf("  L %\n ", L);
 
    return 0;

So line and row are constants now, i killed the "&" so its copy and paste;
problem is at the p[] cause it expects something in front of ']'
can i simply change the int **a to int a[][] in the choldc function so it matches?

next problem that arises is printing it. i need to create a print function for that i guess?

Im basically just a starter and want to see how these functions can be used in main and how to apply.
L = ... is actually not necessary and its just ... i guess.

Code:

choldc(float **a, int n ,float p[]);

in main
choldc(&a, n , p[]);

I came to the following main.
Problem, it runs but it shows error even though its not positive definite

Code:

void nrerror(char error_text[])
/* Numerical Recipes standard error handler */
{
    fprintf(stderr,"Numerical Recipes run-time error...\n");
    fprintf(stderr,"%s\n",error_text);
    fprintf(stderr,"...now exiting to system...\n");
    exit(1);
}

float choldc(float **a, int n ,float p[])
{
    void nrerror(char error_text[]);
    int i,j,k;
    float sum;

    for (i=1;i<=n;i++) {
        for (j=i;j<=n;j++) {
            for (sum=a[i][j],k=i-1;k>=1;k--) sum -= a[i][k]*a[j][k];
            if (i == j) {
                if (sum <= 0.0)
                    nrerror("choldc failed");
                p[i]=sqrt(sum);
            } else a[j][i]=sum/p[i];
        }
    }
    return 1;
}



float cholsl(float **a, int n, float p[], float b[], float x[])
{
    int i,k;
    float sum;

    for (i=1;i<=n;i++) {
        for (sum=b[i],k=i-1;k>=1;k--) sum -= a[i][k]*x[k];
        x[i]=sum/p[i];
    }
    for (i=n;i>=1;i--) {
        for (sum=x[i],k=i+1;k<=n;k++) sum -= a[k][i]*x[k];
        x[i]=sum/p[i];
    }
    return 1;
}


int main() {



    int n = 4;
    float L, T;
    float p[n], x[n];
    int zeile =4;
    int spalte =4;
    float a[zeile][spalte];
    float b[n] = {4, -1,0,0};

a[1][1]=1;     a[1][2]= -0.2; a[1][3]=0;     a[1][4]=0;
a[2][1]= -0.2; a[2][2]=2;     a[2][3]= -0.3; a[2][4]=0;
a[3][1]=0;     a[3][2]= -0.3; a[3][3]=3;     a[3][4]= -0.4;
a[4][1]=0;     a[4][2]=0;     a[4][3]= -0.4; a[4][4]=4;



    L = choldc(a, n, p);
    printf(" Untere Dreiecksmatrix L %\n ", L);

    T = cholsl( a, n, p, b, x)
    printf(" Lösungsvektor T%\n ", T)



    return 0;


}