/*calculates C.x/
/* C     : symetric positive-definite matrix recorded */
/*         (per raws) as a vector with only components */
/*         Cij so that j <= i, i = [0...n-1]*/
/* x     : vector, xi avec i = [0...n-1]*/
/* b     : vector, bi avec i = [0...n-1]*/
/* n     : dimension of matrix Cij*/
/* */
/* The solution vector is returned in b*/

void mat_vec(double *C, double *x, double *b, int n)

{
  int i,k,l;
  double temp;

  for (i = 0; i < n; i++) {
    temp = 0.;

    for (k = 0; k < n; k++) {
      if ( k <= i ) {
	l = k+i*(i+1)/2;
      } else {
	l = i+k*(k+1)/2;
      }
      temp += C[l]*x[k];
    }

    b[i] = temp;
  }

  return;
}