matrix arrays 1

[aaitch]

i wrote this program that inputs a square matrix dimensions n x n. it takes it and makes it into a lower triangular matrix and an upper triangular matrix. if done properly, when the lower and upper triangular matrices are multiplied, the result should be the original matrix. however, when i multiply, all the values are right except for the ones in the middle of the diagonal, for example:
a a a a a
a b a a a
a a b a a
a a a b a
a a a a a the a's match the values in the original matrix but the b's are the ones that are different. i don't know whether i calculated the upper and lower matrices wrong, or i multiplied wrong.

#include <stdio.h>
void mat_in(double **a, int n);
void print_mat(double **a, int n);
void mat_mult(double **c, double **a, double **b, int n);

int
main(void)
{
double **num, **U, **L, **UL, sigma;
int i, j, k, n;
char name;

scanf("%d", &n);

num = (double **) calloc (n, sizeof(double*));
for(i=0; i<n; i++)
num = (double*)calloc(n, sizeof(double));

U = (double **) calloc (n, sizeof(double*));
for(i=0; i<n; i++)
U = (double*)calloc(n, sizeof(double));

L = (double **) calloc (n, sizeof(double*));
for(i=0; i<n; i++)
L = (double*)calloc(n, sizeof(double));

UL = (double **) calloc (n, sizeof(double*));
for(i=0; i<n; i++)
UL = (double*)calloc(n, sizeof(double));

mat_in(num, n);

/* main diagonal of U to 1 */
j = 0;
for(i=0; i<n; i++)
{
U[j] = 1;
j = j+1;
}

/* first L entry */
L[0][0] = num[0][0];

/* first row U, first column L */
for(i=1; i<n; i++)
{
for(j=1; j<n; j++)
{
U[0][j] = num[0][j]/L[0][0];
}

L[0] = num[0]/U[0][0];
}

/* main diagonal L */
for(i=1; i<(n-1); i++)
{
for(k=0; k<(i-1); k++)
{
sigma = sigma + (L[k]*U[k]);
printf("%lf\n", sigma);
}
L = num - sigma;
sigma = 0;

}

/* ith row of U */
for(i=1; i<(n-1); i++)
{
sigma = 0;
for(j=i; j<n; j++)
{
for(k=0; k<(n-1); k++)
{
sigma = sigma + (L[k] * U[k][j]);
}
U[j] = (1/L) * (num[j] - sigma);
sigma = 0;

if (j!=i)
{
for(k=0; k<(n-1); k++)
{
sigma = sigma + (L[j][k] * U[k]);
}
L[j] = (1/U) * (num[j] - sigma);
sigma = 0;
}
}
}

/* last value of n */
for(k=0; k<(n-1); k++)
sigma = sigma + (L[n-1][k] * U[k][n-1]);
L[n-1][n-1] = num[n-1][n-1] - sigma;

mat_mult(UL, L, U, n);

/* printout of U L */
printf("U*L:\n");
for(i=0; i<n; i++)
for(j=0; j<n; j++)
{
printf("%8.4lf", UL[j]);
if (j == (n-1))
printf("\n");
}
}

void
mat_in(double **a, int n)
{
int i, j;
double x;
for(i=0; i<n; i++)
for(j=0; j<n; j++)
{
scanf("%lf", &x);
a[j] = x;
}
}


void
print_mat(double **a, int n)
{
int i, j;

for(i=0; i<n; i++)
for(j=0; j<n; j++)
{
printf("%8.4lf", a[j]);
if (j == (n-1))
printf("\n");
}
}

void
mat_mult(double **c, double **a, double **b, int n)
{
int i, j, k;

for(i=0; i<n; i++)
for(j=0; j<n; j++)
{
c[j] = 0;
for(k=0; k<n; k++)
{
c[j] = c[j] + (a[k] * b[k][j]);
}
}

for(i=0; i<n; i++)
{
c = 0;
for(k=0; k<n; k++)
{
c = c + (a[k] * b[k]);
}
}
}

 

[Salem]

Please use the [tt][ignore]

[/ignore][/tt] tags when posting code.

> num = (double **) calloc (n, sizeof(double*));
Remove the cast - that's the (double**) and add
#include <stdlib.h>
to the start of the program.
If you're still getting complaints from the compiler (like cannot convert void* or something), then try using a C compiler instead of a C++ compiler.

Also FYI, calloc does not guarantee that the floating point values stored in your array actually equal 0.0.

http://www.eskimo.com/~scs/C-faq/q7.31.html

> all the values are right except for the ones in the middle of the diagonal, for example
How are they wrong?
Two things spring to mind
1. You have an "off-by-1" error in one of your loops. That is, you're looping one too few, or one too many times.
2. The nature of your computation is such that what appears to be mathematically correct is in fact computationally in error. For example, ((1.0/3.0)*3.0) != 1.0 (computationally). A few of those, and you can be a long way off pretty quickly.

--