#include <stdio.h>
#include <stdlib.h>
#include <time.h>
#include <math.h>


double random_number(double *rmin, double *rmax){
  return (*rmin) + (float)rand()/(float)RAND_MAX*((*rmax)-(*rmin));
}

/* This function prints the matrix A to the screen. */
void print_matrix(double **A, int n) {
  int i, j;
  for(i=0; i<n; i++) {
    for(j=0; j<n; j++){
      printf("% 6.2e  ",A[i][j]);
    }
    printf("\n");
  }
  printf("\n");
}


/* This function prints the vector c to the screen. */
void print_vector(double *c, int n) {
  int i;
  for(i=0; i<n; i++){
    printf("% 6.2e  ", c[i]);
  }
  printf("\n");
}


/* This function computes the scalar product a*b and returns the result. */
double scalar(double *a, double *b, int n){
  int i;
  double sum=0.;
  for(i=0; i<n; i++){
    sum += a[i]*b[i];
  }
  return sum;
}


/* This function computes the matrix-vector product Ax. The result
   is returned as f. */
void matrix_vector(double **A, double *x, double *f, int n){
  int i,j;
  for (i=0; i<n; i++){
    f[i] = 0.;
    for (j=0; j<n; j++){
      f[i] += A[i][j] * x[j];
    }
  }
}


/* This function computes x+al*y and stores the result in u. */
void addmult(double *u, double *x, double *y, double al, int n){
  int i;
  for (i=0; i<n; i++){
    u[i] = x[i] + al*y[i];
  }
}


int main() {
  double **A, *d, *x, *e, *p, *r, *f;
  double *Af;
  double om = 1.7;     /* omega: over-relaxation */
  double eps;          /* stop criterion */
  double rmin,rmax;    /* for generation of randum numbers */
  double tmp,res;      /* res is the resiuum */
  double al,be;
  int i, j, k, n;
  int k_stop;          /* maximal iteration number */

  FILE *fp;

  n = 8;   /* dimension of the problem */
  k_stop = 100;
  eps = 1.0e-20;

  rmin = 0;
  rmax = 1;

  A  = (double**)malloc(n*sizeof(double*));
  Af = (double*)malloc(n*n*sizeof(double));

  d = (double *) malloc(n*sizeof(double) );   
  x = (double *) malloc(n*sizeof(double) ); 
  e = (double *) malloc(n*sizeof(double) );
  r = (double *) malloc(n*sizeof(double) );   
  p = (double *) malloc(n*sizeof(double) );   
  f = (double *) malloc(n*sizeof(double) );   

  for(i=0; i<n; i++){
    A[i] = &Af[i*n];
  }

  /* initialisation */
  srand( (unsigned)time(0));
  for(i=0; i<n; i++) {
    for(j=0; j<=i; j++){
      A[i][j] = 0.;
      while( fabs(A[i][j]) < 1.0e-10 ){
        A[i][j] = random_number(&rmin,&rmax);
        A[j][i] = A[i][j];
      }
    }
    A[i][i] += (float)n+3.;   /* ensure diagonal dominant matrix */ 
    d[i] = random_number(&rmin,&rmax);
    x[i] = 0.0;      /* start vector */
  }

  k=0;
  res = 1.;

  matrix_vector(A,x,f,n);
  for(i=0; i<n; i++){
    r[i] = d[i] - f[i];
    p[i] = r[i];
  }
  tmp = scalar(r,r,n);
  
  fp = fopen("residual.dat","w");
  while( res>eps && k<=k_stop ){
    matrix_vector(A,p,f,n);
    al = -tmp/scalar(p,f,n);
    addmult(x,x,p,-al,n);
    addmult(r,r,f,al,n);
    res = scalar(r,r,n);
    be = res/tmp;
    tmp = res;
    addmult(p,r,p,be,n);
    fprintf(fp,"%d %le\n",k,res);
    k++;
  }
  fclose(fp);

  /* finally check the result */
  for(i=0; i<n; i++){
    e[i] = 0.0; 
    for(j=0; j<n; j++){
      e[i] += A[i][j]*x[j];
    }
    e[i] = d[i]-e[i];
  }
  printf("\n");
  print_vector(e,n);  
}