-1

I am implementing the Cholesky Method in C but the program quits when it arrives at this point.

After the answers : Now it works thanks to the answers of (devnull & piotruś) but it doens't give me the right answer

/* Ax=b
 *This algorithm does:
 *   A = U * U' 
 * with
 *   U  := lower left triangle matrix
 *   U' := the transposed form of U.
 */

double** cholesky(double **A, int N) //this gives me the U (edited)
{
    int i, j, k;
    double sum, **p, **U;
    U=(double**)malloc(N*sizeof(double*));
    for(p=U; p<U+N; p++)
        *p=(double*)malloc(N*sizeof(double));
    for (j=0; j<N; j++) {
        sum = A[j][j];
        for (k=0; k<(j-1); k++) sum -= U[k][j]*U[k][j];
            U[j][j] = sqrt(sum);
            for (i=j; i<N; i++) {
                sum = A[j][i];
                for (k=0; k<(j-1); k++)
                    sum -= U[k][j]*U[k][i];
                U[j][i] = sum/U[j][j];
            }
    }
return U;
}

am I doing something wrong here?

user3316322
  • 3
  • 1
  • 3
  • Did you try a debugger? Using a debugger, you could at least pin down the offending line and that would likely make it clear to yourself what the problem is. That could also tell you what is the error: does it simply skip execution of the function and hence terminate? Does it get killed with a segmentation fault? A floating point exception? – Shahbaz Feb 16 '14 at 17:02
  • One definite problem is that you are assuming that arrays in C are indexed from 1. – devnull Feb 16 '14 at 17:02
  • what it means it quits ? at what point? what is called? – 4pie0 Feb 16 '14 at 17:02
  • 1
    When I used to write code like this in C I extensively used Valgrind to check for memory related bugs and errors. – shuttle87 Feb 16 '14 at 17:03
  • You are doing a lot "wrong", but the main thing is that for modern C there is no reason at all to emulate a 2D matrix with tables of pointers. It only brings you trouble and has worse performance. – Jens Gustedt Feb 16 '14 at 17:19

1 Answers1

0
double** cholesky(double **A, int N)

in this function you assume array length is N. This means the last index of array is at N-1 not at N. Change the code into:

for ( j = 0; j < N; ++j)

and the rest similarly.

4pie0
  • 29,204
  • 9
  • 82
  • 118