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?
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.