#include <stdio.h>
#include <stdlib.h>
#include <time.h>
#include <math.h>
/*******************************************************************/
void print_vector(double *c, int n) {
   int i;
   for(i = 0; i < n ; i++)
      printf("% 6.2e  ", c[i]);
   printf("\n");
}
/*******************************************************************/

int main() {
   double **a , **a_tmp, *b, *x, *e, *p, *r, *r_old;
   double alpha, beta, tmp1, tmp2, n_res;
   double eps;
   int    i, j, k, n;

   n   = 11;       /* dimension of the problem */
   eps = 1.e-16;   /* stop for iteration */

   a     = (double **)malloc( n * sizeof(double*) );
   a_tmp = (double **)malloc( n * sizeof(double*) );
   b     = (double *) malloc( n * sizeof(double) );   
   p     = (double *) malloc( n * sizeof(double) );
   x     = (double *) malloc( n * sizeof(double) ); 
   r     = (double *) malloc( n * sizeof(double) );
   r_old = (double *) malloc( n * sizeof(double) );
   e     = (double *) malloc( n * sizeof(double) );

   for(i=0; i < n; i++) {
      a[i]     = (double *)malloc( n * sizeof(double) );
      a_tmp[i] = (double *)malloc( n * sizeof(double) );
   }

   /* initialization */
   /* srand( (unsigned)time( (long *)0 )); */
   for(i = 0; i < n ; i++) {
      for(j = 0; j < n ; j++)
         a_tmp[i][j] = rand() /2147483647.0;
      b[i] = rand() /2147483647.0;              /* right hand side of equation */
      x[i] = 0.0;                               /* start vector */
   }
   for(i = 0; i < n ; i++) 
      for(j = 0; j < n ; j++)
         a[i][j]  = a_tmp[j][i] * a_tmp[i][j];  /* matrix "a" is now symmetric and positiv definit */
 
   for(i = 0; i < n ; i++) {
      r[i] = b[i];                              /* b-Ax = b because of initialization x_0=0 */
      p[i] = r[i];
   }
   k = 1;
   do  {
     tmp1 = 0.0;
     tmp2 = 0.0;
     if( k > 1 ) {
        for(i = 0; i < n ; i++) {
           tmp1 += r[i] * r[i];
           tmp2 += r_old[i] * r_old[i];
        }
        beta = tmp1 / tmp2 ;  
     } else {
        beta = 0.0;
     }
     k++;
     memcpy( r_old , r , n * sizeof(double) ); 

     for(i = 0; i < n ; i++)
        p[i] = r[i] + beta * p[i];

     tmp1 = 0.0;
     tmp2 = 0.0;
     for(i = 0; i < n ; i++) {
        tmp1 += r[i] * r[i];
        for(j = 0; j < n ; j++)
           tmp2 += p[i] * a[i][j] * p[j];
     }
     alpha = tmp1 / tmp2 ;

     for(i = 0; i < n ; i++)
        x[i] += alpha * p[i];

     for(i = 0; i < n ; i++) {
        tmp1 = 0.0;
        for(j = 0; j < n ; j++)
           tmp1 += a[i][j] * p[j];
        r[i] -= alpha * tmp1;
     }
     n_res = 0.0;
     for(i = 0; i < n ; i++)
        n_res += r[i] * r[i];      /* calculate norm of residuum */
     n_res = sqrt(n_res);
     print_vector(x,n); 
   } while ( n_res > eps ) ;

   for(i = 0; i < n; i++) {        /* calculate and print the final error */
      e[i] = 0.0; 
      for(j = 0; j < n; j++) 
         e[i] += a[i][j]*x[j];
      e[i] = b[i]-e[i];
   }
   printf("\n");
   print_vector(e,n);  
}
/*******************************************************************/