text: #include "blas_extended.h"
  #include "blas_extended_private.h"
  void		BLAS_ddot_x(enum blas_conj_type conj, int n, double alpha,
  	     		const		double *x, int incx, double beta,
  		  		const		double *y, int incy,
  		  		double       *r, enum blas_prec_type prec)
  /*
   * Purpose
   * =======
   *
   * This routine computes the inner product:
   *
   *     r <- beta * r + alpha * SUM_{i=0, n-1} x[i] * y[i].
   *
   * Arguments
   * =========
   *
   * conj   (input) enum blas_conj_type
   *        When x and y are complex vectors, specifies whether vector
   *        components x[i] are used unconjugated or conjugated.
   *
   * n      (input) int
   *        The length of vectors x and y.
   *
   * alpha  (input) double
   *
   * x      (input) const double*
   *        Array of length n.
   *
   * incx   (input) int
   *        The stride used to access components x[i].
   *
   * beta   (input) double
   *
   * y      (input) const double*
   *        Array of length n.
   *
   * incy   (input) int
   *        The stride used to access components y[i].
   *
   * r      (input/output) double*
   *
   * prec   (input) enum blas_prec_type
   *        Specifies the internal precision to be used.
   *        = blas_prec_single: single precision.
   *        = blas_prec_double: double precision.
   *        = blas_prec_extra : anything at least 1.5 times as accurate
   *                            than double, and wider than 80-bits.
   *                            We use double-double in our implementation.
   *
   */
  {
    static const char routine_name[] = "BLAS_ddot_x";
  
    switch (prec) {
    case blas_prec_single:
    case blas_prec_double:
    case blas_prec_indigenous:
      {
        int	      i    , ix = 0, iy = 0;
        double         *r_i = r;
        const double   *x_i = x;
        const double   *y_i = y;
        double	      alpha_i = alpha;
        double	      beta_i = beta;
        double	      x_ii;
        double	      y_ii;
        double	      r_v;
        double	      prod;
        double	      sum;
        double	      tmp1;
        double	      tmp2;
  
  
        /* Test the input parameters. */
        if (n < 0)
  	BLAS_error(routine_name, -2, n, NULL);
        else if (incx == 0)
  	BLAS_error(routine_name, -5, incx, NULL);
        else if (incy == 0)
  	BLAS_error(routine_name, -8, incy, NULL);
  
        /* Immediate return. */
        if ((beta_i == 1.0) && (n == 0 || (alpha_i == 0.0)))
  	return;
  
  
  
        r_v = r_i[0];
        sum = 0.0;
  
  
        if (incx < 0)
  	ix = (-n + 1) * incx;
        if (incy < 0)
  	iy = (-n + 1) * incy;
  
        for (i = 0; i < n; ++i) {
  	x_ii = x_i[ix];
  	y_ii = y_i[iy];
  
  	prod = x_ii * y_ii;	/* prod = x[i]*y[i] */
  	sum = sum + prod;	/* sum = sum+prod */
  	ix += incx;
  	iy += incy;
        }				/* endfor */
  
  
        tmp1 = sum * alpha_i;	/* tmp1 = sum*alpha */
        tmp2 = r_v * beta_i;	/* tmp2 = r*beta */
        tmp1 = tmp1 + tmp2;	/* tmp1 = tmp1+tmp2 */
        *r = tmp1;		/* r = tmp1 */
  
  
      }
      break;
    case blas_prec_extra:
      {
        int	      i    , ix = 0, iy = 0;
        double         *r_i = r;
        const double   *x_i = x;
        const double   *y_i = y;
        double	      alpha_i = alpha;
        double	      beta_i = beta;
        double	      x_ii;
        double	      y_ii;
        double	      r_v;
        double	      head_prod, tail_prod;
        double	      head_sum, tail_sum;
        double	      head_tmp1, tail_tmp1;
        double	      head_tmp2, tail_tmp2;
        FPU_FIX_DECL;
  
        /* Test the input parameters. */
        if (n < 0)
  	BLAS_error(routine_name, -2, n, NULL);
        else if (incx == 0)
  	BLAS_error(routine_name, -5, incx, NULL);
        else if (incy == 0)
  	BLAS_error(routine_name, -8, incy, NULL);
  
        /* Immediate return. */
        if ((beta_i == 1.0) && (n == 0 || (alpha_i == 0.0)))
  	return;
  
        FPU_FIX_START;
  
        r_v = r_i[0];
        head_sum = tail_sum = 0.0;
  
  
        if (incx < 0)
  	ix = (-n + 1) * incx;
        if (incy < 0)
  	iy = (-n + 1) * incy;
  
        for (i = 0; i < n; ++i) {
  	x_ii = x_i[ix];
  	y_ii = y_i[iy];
  
  	{
  	  /* Compute double_double = double * double. */
  	  double	  a1    , a2, b1, b2, con;
  
  	  con = x_ii * split;
  	  a1 = con - x_ii;
  	  a1 = con - a1;
  	  a2 = x_ii - a1;
  	  con = y_ii * split;
  	  b1 = con - y_ii;
  	  b1 = con - b1;
  	  b2 = y_ii - b1;
  
  	  head_prod = x_ii * y_ii;
  	  tail_prod = (((a1 * b1 - head_prod) + a1 * b2) + a2 * b1) + a2 * b2;
  	}			/* prod = x[i]*y[i] */
  	{
  	  /* Compute double-double = double-double + double-double. */
  	  double	  bv;
  	  double	  s1    , s2, t1, t2;
  
  	  /* Add two hi words. */
  	  s1 = head_sum + head_prod;
  	  bv = s1 - head_sum;
  	  s2 = ((head_prod - bv) + (head_sum - (s1 - bv)));
  
  	  /* Add two lo words. */
  	  t1 = tail_sum + tail_prod;
  	  bv = t1 - tail_sum;
  	  t2 = ((tail_prod - bv) + (tail_sum - (t1 - bv)));
  
  	  s2 += t1;
  
  	  /* Renormalize (s1, s2)  to  (t1, s2) */
  	  t1 = s1 + s2;
  	  s2 = s2 - (t1 - s1);
  
  	  t2 += s2;
  
  	  /* Renormalize (t1, t2)  */
  	  head_sum = t1 + t2;
  	  tail_sum = t2 - (head_sum - t1);
  	}			/* sum = sum+prod */
  	ix += incx;
  	iy += incy;
        }				/* endfor */
  
  
        {
  	/* Compute double-double = double-double * double. */
  	double		a11    , a21, b1, b2, c11, c21, c2, con, t1, t2;
  
  	con = head_sum * split;
  	a11 = con - head_sum;
  	a11 = con - a11;
  	a21 = head_sum - a11;
  	con = alpha_i * split;
  	b1 = con - alpha_i;
  	b1 = con - b1;
  	b2 = alpha_i - b1;
  
  	c11 = head_sum * alpha_i;
  	c21 = (((a11 * b1 - c11) + a11 * b2) + a21 * b1) + a21 * b2;
  
  	c2 = tail_sum * alpha_i;
  	t1 = c11 + c2;
  	t2 = (c2 - (t1 - c11)) + c21;
  
  	head_tmp1 = t1 + t2;
  	tail_tmp1 = t2 - (head_tmp1 - t1);
        }				/* tmp1 = sum*alpha */
        {
  	/* Compute double_double = double * double. */
  	double		a1     , a2, b1, b2, con;
  
  	con = r_v * split;
  	a1 = con - r_v;
  	a1 = con - a1;
  	a2 = r_v - a1;
  	con = beta_i * split;
  	b1 = con - beta_i;
  	b1 = con - b1;
  	b2 = beta_i - b1;
  
  	head_tmp2 = r_v * beta_i;
  	tail_tmp2 = (((a1 * b1 - head_tmp2) + a1 * b2) + a2 * b1) + a2 * b2;
        }				/* tmp2 = r*beta */
        {
  	/* Compute double-double = double-double + double-double. */
  	double		bv;
  	double		s1     , s2, t1, t2;
  
  	/* Add two hi words. */
  	s1 = head_tmp1 + head_tmp2;
  	bv = s1 - head_tmp1;
  	s2 = ((head_tmp2 - bv) + (head_tmp1 - (s1 - bv)));
  
  	/* Add two lo words. */
  	t1 = tail_tmp1 + tail_tmp2;
  	bv = t1 - tail_tmp1;
  	t2 = ((tail_tmp2 - bv) + (tail_tmp1 - (t1 - bv)));
  
  	s2 += t1;
  
  	/* Renormalize (s1, s2)  to  (t1, s2) */
  	t1 = s1 + s2;
  	s2 = s2 - (t1 - s1);
  
  	t2 += s2;
  
  	/* Renormalize (t1, t2)  */
  	head_tmp1 = t1 + t2;
  	tail_tmp1 = t2 - (head_tmp1 - t1);
        }				/* tmp1 = tmp1+tmp2 */
        *r = head_tmp1;		/* r = tmp1 */
  
        FPU_FIX_STOP;
      }
      break;
    }
  }
  