Thread: Cholesky decomposition example

  1. #1
    Registered User
    Join Date
    Nov 2017
    Posts
    8

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

  2. #2
    Programming Wraith GReaper's Avatar
    Join Date
    Apr 2009
    Location
    Greece
    Posts
    2,739
    I know nothing about the cholesky decomposition, but I detect many loops counting from 1 to n, while array indices run from 0 to n-1. Is this code yours, have you tested it to see if it works?

    Have you read the wikipedia page on this?
    Cholesky decomposition - Wikipedia
    Devoted my life to programming...

  3. #3
    Registered User
    Join Date
    Nov 2017
    Posts
    8
    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.

  4. #4
    Registered User
    Join Date
    Nov 2017
    Posts
    8
    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

  5. #5
    and the hat of int overfl Salem's Avatar
    Join Date
    Aug 2001
    Location
    The edge of the known universe
    Posts
    38,620
    Perhaps you should post your main(), if that's what you're having trouble with.
    If you dance barefoot on the broken glass of undefined behaviour, you've got to expect the occasional cut.
    If at first you don't succeed, try writing your phone number on the exam paper.

  6. #6
    Registered User
    Join Date
    Nov 2017
    Posts
    8
    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

  7. #7
    Registered User
    Join Date
    Nov 2017
    Posts
    8
    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

  8. #8
    and the hat of int overfl Salem's Avatar
    Join Date
    Aug 2001
    Location
    The edge of the known universe
    Posts
    38,620
    Well you need to initialise line and row.
    And they need to be literal constants unless you have a C99 compiler.

    Consider this an example of how to pass an array to a function.
    Code:
    void foo ( int arr[3][4] ) {
    }
    
    int main ( ) {
      int arr[3][4];
      foo(arr);
    }
    Notice that the parameter declaration is a copy/paste of the array you want to pass.

    Now declare suitable arrays in your main which match what your functions expect.

    Note that
    void foo ( int **a )
    is VERY different from
    void foo ( int a[][N] )
    They are not interchangeable.
    If you dance barefoot on the broken glass of undefined behaviour, you've got to expect the occasional cut.
    If at first you don't succeed, try writing your phone number on the exam paper.

  9. #9
    Registered User
    Join Date
    Nov 2017
    Posts
    8
    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?

  10. #10
    and the hat of int overfl Salem's Avatar
    Join Date
    Aug 2001
    Location
    The edge of the known universe
    Posts
    38,620
    Well it's a side effect of having "found" code without having a clue how to use it.

    Where did you figure/learn that
    L = choldc(a, n, p[]);
    was in any way valid syntax.

    Just how much actual C have you done prior to today?
    If you dance barefoot on the broken glass of undefined behaviour, you've got to expect the occasional cut.
    If at first you don't succeed, try writing your phone number on the exam paper.

  11. #11
    Registered User
    Join Date
    Nov 2017
    Posts
    8
    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[]);

  12. #12
    Registered User
    Join Date
    Nov 2017
    Posts
    8
    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;
    
    
    }

Popular pages Recent additions subscribe to a feed

Similar Threads

  1. Integers Decomposition
    By raven2313 in forum C Programming
    Replies: 4
    Last Post: 07-19-2015, 12:02 PM
  2. LU decomposition
    By khdani in forum C Programming
    Replies: 3
    Last Post: 10-15-2010, 08:08 PM
  3. Complex Cholesky Function
    By magda3227 in forum C Programming
    Replies: 3
    Last Post: 07-18-2008, 08:55 PM
  4. Lu Decomposition Help Me Please
    By scottyg in forum C Programming
    Replies: 3
    Last Post: 09-29-2006, 08:33 AM
  5. LU decomposition
    By scottmanc in forum C++ Programming
    Replies: 1
    Last Post: 10-15-2003, 08:25 AM

Tags for this Thread