C/C++ Matrix/Vector Calculation

[arperidot]

Hey, I have a program which is done but it gives me the wrong output. It would be greatly appreciated if somebody could help me get the correct output. The question asks:

Develop a program which computes the current value of the vector {x} based on the following forward iteration: {x}(n+1) = [K] {x}(n)+ {z}, n = 0,1,2,3,4.

In other words, the next vector {x} is equal to the product of [K] and the current vector {x} + constant vector {z}.

Perform the matrix multiplication by using the function:

void matrix_mult(dou ble a[3][3], double *b, double *c, int rows);

You should pass the matrix [K], and a pointer array {x} to matrix mult(), which will pass back a vector {y}=[K]{x}.

[K] is a constant 3x3 matrix defined by:
double k[3][3] = { { 1.2, -3.0, 1.1},
{-4.1, 6.2, -2.2},
{ 3.4, -2.5, -3.3} };

The initial components of the vector {x}(0) are:
double x[3] = { 1./sqrt(2.), 0., -1./sqrt(2.) },

while the constant vector {z} is:
double z[3] = { 1./sqrt(3.), -1./sqrt(3.), 1./sqrt(3.) }.

Your output should look like:

Initial vector:
x[0] = [ 0.707107 0.000000 -0.707107 ]

New vector:
y[1] = [ 0.648061 -1.920853 5.314966 ]
Normalized new vector:
x[1] = [ 0.113926 -0.337676 0.934343 ]

............... ............... ........…
............... ............... ........…
............... ............... ........…

My code:

#include <stdio.h>
#include <math.h>
void unit_norm(double *vec, int cols);
void matrix_mult(double a[3][3], double *b, double *c, int rows);
main()
{
 int row, i, j;
 double k[3][3] = {{1.2, -3.0,  1.1},
                  {-4.1, 6.2, -2.2},
                  {3.4, -2.5, -3.3}};
 double y[3], x[3] = { 1./sqrt(2.), 0., -1./sqrt(2.) };

 double w[3], z[3] = { 1./sqrt(3.), -1/sqrt(3.), 1./sqrt(3.) };

 printf("\nInitial vector:\n\n");
 printf("x[0] = [");
 for(row=0; row<3; row++)
        printf(" %.6f ", x[row]);
 printf("]\n\n");

  for(i=0; i<=4; i++)
  {
   matrix_mult(k,x,y,3);
   printf("New Vector:\ny[%d] = [", i+1);
   for(j=0; j<=2; j++)
        printf(" %.6f ", y[j]);
   printf("]\n");

   unit_norm(y,3);
   printf("Normalized new vector:\nx[%d] = [", i+1);
   for(j=0; j<=2; j++)
        printf(" %.6f ", y[j]);
   printf("]\n");
   for(j=0; j<=2; j++)
        x[j] = y[j];
   printf("\n");
  }
  return(0);
}

void matrix_mult(double a[3][3], double *b, double *c, int rows)
{
        int k, s, r;
        double dot_p, res_vec[50];
        double w[3], z[3] = { 1./sqrt(3.), -1/sqrt(3.), 1./sqrt(3.) };
        for(r=0; r<rows; r++)
        {
                for(s=0; s<3; s++)
                {
                        for(k=0, dot_p=0; k<3; k++)
                        { dot_p += a[s][k]*b[k]; }
                }
                { c[s] = dot_p; }
        { c[r] += z[r]; }
        }
        return;
}

void unit_norm(double *vec, int cols)
{
        double norm;
        double sum=0.;
        int n;
        for(n=0; n<cols; n++)
         sum += vec[n] * vec[n];
        sum = sqrt(sum);
        for(n=0; n<cols; n++)
            vec[n] = vec[n]/sum;
return;
}

[jkmyoung]

Interchange lines 52 and 53!
}
{ c[s] = dot_p; }

[arperidot]

thanks for the quick reply.
after doing that, it gives a different value but it is still wrong. any other suggestions??

[jkmyoung]

There's something wrong with your for loops.
You're have 3: -> 3 * 3 * 3 = 27 multiplications .
A [3X3] array * [3X1] array should only have 9 multiplications .

The r for loop is misplaced somehow. If you got rid of that and used
c[s] += z[s];
instead, it'll probably work out.

[arperidot]

thanks alot!!!! that along with rearranging my loop produced the desired output. this is what the final matric_mult looked like:

void matrix_mult(double a[3][3], double *b, double *c, int rows)
{
int k, s, r;
double dot_p, res_vec[5];
double w[3], z[3] = { 1./sqrt(3.), -1/sqrt(3.), 1./sqrt(3.) };
        for(s=0; s<3; s++)
        {
                for(k=0, dot_p=0; k<3; k++) 
                {
                dot_p += a[s][k]*b[k];
                }
        c[s] = dot_p;
        }
        for(r=0; r<rows; r++)
        {
                c[r] += z[r];
        }
return;
}