text: #include <math.h>
  #include "blas_extended.h"
  #include "blas_extended_private.h"
  
  void		BLAS_ctrsv_s_x(enum blas_order_type order, enum blas_uplo_type uplo,
  		     		enum		blas_trans_type trans, enum blas_diag_type diag,
  		     		int		n       , const void *alpha, const float *T, int ldt,
  	     		void         *x, int incx, enum blas_prec_type prec)
  /*
   * Purpose
   * =======
   *
   * This routine solve :
   *
   *     x <- alpha * inverse(T) * x
   *
   * Arguments
   * =========
   *
   * order  (input) enum blas_order_type
   *        column major, row major
   *
   * uplo   (input) enum blas_uplo_type
   *        upper, lower
   *
   * trans  (input) enum blas_trans_type
   *        no trans, trans, conj trans
   *
   * diag   (input) enum blas_diag_type
   *        unit, non unit
   *
   * n      (input) int
   *        the dimension of T
   *
   * alpha  (input) const void*
   *
   * T      (input) float*
   *        Triangular matrix
   *
   * x      (input) const void*
   *           Array of length n.
   *
   * incx   (input) int
   *           The stride used to access components x[i].
   *
   * 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.
   *
   */
  {
    char           *routine_name = "BLAS_ctrsv_s";
  
    int		  i       , j;	/* used to idx matrix */
    int		  ix      , jx;	/* used to idx vector x */
    int		  start_x;	/* used as the starting idx to vector x */
    const float    *T_i = T;	/* internal matrix T */
    float          *x_i = (float *)x;	/* internal x */
    float          *alpha_i = (float *)alpha;	/* internal alpha */
    float		  T_element;	/* temporary variable for an element of
  				 * matrix A */
    int		  incT = 1;	/* internal ldt */
  
    if ((order != blas_rowmajor && order != blas_colmajor) ||
        (uplo != blas_upper && uplo != blas_lower) ||
        (trans != blas_trans && trans !=
         blas_no_trans && trans != blas_conj_trans) ||
        (diag != blas_non_unit_diag && diag != blas_unit_diag) ||
        (ldt < n) ||
        (incx == 0)) {
      BLAS_error(routine_name, 0, 0, NULL);
    }
    if (n <= 0)
      return;
  
  
    incx *= 2;
    /* configuring the vector starting idx */
    if (incx <= 0) {
      start_x = -(n - 1) * incx;
    } else {
      start_x = 0;
    }
  
    /* if alpha is zero, then return x as a zero vector */
    if (alpha_i[0] == 0.0 && alpha_i[1] == 0.0) {
      ix = start_x;
      for (i = 0; i < n; i++) {
        x_i[ix] = 0.0;
        x_i[ix + 1] = 0.0;
        ix += incx;
      }
      return;
    }
    switch (prec) {
    case blas_prec_single:
  
      {
        float	      temp1[2];	/* temporary variable for calculations */
        float	      temp2[2];	/* temporary variable for calculations */
        float	      temp3[2];	/* temporary variable for calculations */
  
        if ((order == blas_rowmajor &&
  	   trans == blas_no_trans && uplo == blas_upper) ||
  	  (order == blas_colmajor &&
  	   trans != blas_no_trans && uplo == blas_lower)) {
  
  	jx = start_x + (n - 1) * incx;
  	for (j = n - 1; j >= 0; j--) {
  
  	  /*
  	   * compute Xj = alpha*Xj - SUM Tij(or Tji) * Xi i=j+1 to n-1
  	   */
  	  temp3[0] = x_i[jx];
  	  temp3[1] = x_i[1 + jx];
  	  {
  	    temp1[0] = temp3[0] * alpha_i[0] - temp3[1] * alpha_i[1];
  	    temp1[1] = temp3[0] * alpha_i[1] + temp3[1] * alpha_i[0];
  	  }
  
  
  	  ix = start_x + (n - 1) * incx;
  	  for (i = n - 1; i >= j + 1; i--) {
  	    T_element = T_i[i * incT + j * ldt * incT];
  
  	    temp3[0] = x_i[ix];
  	    temp3[1] = x_i[1 + ix];
  	    {
  	      temp2[0] = temp3[0] * T_element;
  	      temp2[1] = temp3[1] * T_element;
  	    }
  	    temp1[0] = temp1[0] - temp2[0];
  	    temp1[1] = temp1[1] - temp2[1];
  	    ix -= incx;
  	  }			/* for j<n */
  
  	  /*
  	   * if the diagonal entry is not equal to one, then divide Xj by the
  	   * entry
  	   */
  	  if (diag == blas_non_unit_diag) {
  	    T_element = T_i[j * incT + j * ldt * incT];
  
  
  	    temp1[0] = temp1[0] / T_element;
  	    temp1[1] = temp1[1] / T_element;
  
  	  }			/* if (diag == blas_non_unit_diag) */
  	  x_i[jx] = temp1[0];
  	  x_i[jx + 1] = temp1[1];
  
  	  jx -= incx;
  	}			/* for j>=0 */
        } else if ((order == blas_rowmajor &&
  		  trans == blas_no_trans && uplo == blas_lower) ||
  		 (order == blas_colmajor &&
  		  trans != blas_no_trans && uplo == blas_upper)) {
  
  	jx = start_x;
  	for (j = 0; j < n; j++) {
  
  	  /*
  	   * compute Xj = alpha*Xj - SUM Aij(or Aji) * Xi i=j+1 to n-1
  	   */
  	  temp3[0] = x_i[jx];
  	  temp3[1] = x_i[1 + jx];
  	  /* multiply by alpha */
  	  {
  	    temp1[0] = temp3[0] * alpha_i[0] - temp3[1] * alpha_i[1];
  	    temp1[1] = temp3[0] * alpha_i[1] + temp3[1] * alpha_i[0];
  	  }
  
  
  	  ix = start_x;
  	  for (i = 0; i < j; i++) {
  	    T_element = T_i[i * incT + j * ldt * incT];
  
  	    temp3[0] = x_i[ix];
  	    temp3[1] = x_i[1 + ix];
  	    {
  	      temp2[0] = temp3[0] * T_element;
  	      temp2[1] = temp3[1] * T_element;
  	    }
  	    temp1[0] = temp1[0] - temp2[0];
  	    temp1[1] = temp1[1] - temp2[1];
  	    ix += incx;
  	  }			/* for i<j */
  
  	  /*
  	   * if the diagonal entry is not equal to one, then divide Xj by the
  	   * entry
  	   */
  	  if (diag == blas_non_unit_diag) {
  	    T_element = T_i[j * incT + j * ldt * incT];
  
  
  	    temp1[0] = temp1[0] / T_element;
  	    temp1[1] = temp1[1] / T_element;
  
  	  }			/* if (diag == blas_non_unit_diag) */
  	  x_i[jx] = temp1[0];
  	  x_i[jx + 1] = temp1[1];
  	  jx += incx;
  	}			/* for j<n */
        } else if ((order == blas_rowmajor &&
  		  trans != blas_no_trans && uplo == blas_lower) ||
  		 (order == blas_colmajor &&
  		  trans == blas_no_trans && uplo == blas_upper)) {
  
  	jx = start_x + (n - 1) * incx;
  	for (j = n - 1; j >= 0; j--) {
  
  	  /*
  	   * compute Xj = alpha*Xj - SUM Tij(or Tji) * Xi i=j+1 to n-1
  	   */
  	  temp3[0] = x_i[jx];
  	  temp3[1] = x_i[1 + jx];
  	  {
  	    temp1[0] = temp3[0] * alpha_i[0] - temp3[1] * alpha_i[1];
  	    temp1[1] = temp3[0] * alpha_i[1] + temp3[1] * alpha_i[0];
  	  }
  
  
  	  ix = start_x + (n - 1) * incx;
  	  for (i = n - 1; i >= j + 1; i--) {
  	    T_element = T_i[j * incT + i * ldt * incT];
  
  	    temp3[0] = x_i[ix];
  	    temp3[1] = x_i[1 + ix];
  	    {
  	      temp2[0] = temp3[0] * T_element;
  	      temp2[1] = temp3[1] * T_element;
  	    }
  	    temp1[0] = temp1[0] - temp2[0];
  	    temp1[1] = temp1[1] - temp2[1];
  	    ix -= incx;
  	  }			/* for j<n */
  
  	  /*
  	   * if the diagonal entry is not equal to one, then divide Xj by the
  	   * entry
  	   */
  	  if (diag == blas_non_unit_diag) {
  	    T_element = T_i[j * incT + j * ldt * incT];
  
  
  	    temp1[0] = temp1[0] / T_element;
  	    temp1[1] = temp1[1] / T_element;
  
  	  }			/* if (diag == blas_non_unit_diag) */
  	  x_i[jx] = temp1[0];
  	  x_i[jx + 1] = temp1[1];
  
  	  jx -= incx;
  	}			/* for j>=0 */
        } else if ((order == blas_rowmajor &&
  		  trans != blas_no_trans && uplo == blas_upper) ||
  		 (order == blas_colmajor &&
  		  trans == blas_no_trans && uplo == blas_lower)) {
  
  	jx = start_x;
  	for (j = 0; j < n; j++) {
  
  	  /*
  	   * compute Xj = alpha*Xj - SUM Aij(or Aji) * Xi i=j+1 to n-1
  	   */
  	  temp3[0] = x_i[jx];
  	  temp3[1] = x_i[1 + jx];
  	  /* multiply by alpha */
  	  {
  	    temp1[0] = temp3[0] * alpha_i[0] - temp3[1] * alpha_i[1];
  	    temp1[1] = temp3[0] * alpha_i[1] + temp3[1] * alpha_i[0];
  	  }
  
  
  	  ix = start_x;
  	  for (i = 0; i < j; i++) {
  	    T_element = T_i[j * incT + i * ldt * incT];
  
  	    temp3[0] = x_i[ix];
  	    temp3[1] = x_i[1 + ix];
  	    {
  	      temp2[0] = temp3[0] * T_element;
  	      temp2[1] = temp3[1] * T_element;
  	    }
  	    temp1[0] = temp1[0] - temp2[0];
  	    temp1[1] = temp1[1] - temp2[1];
  	    ix += incx;
  	  }			/* for i<j */
  
  	  /*
  	   * if the diagonal entry is not equal to one, then divide Xj by the
  	   * entry
  	   */
  	  if (diag == blas_non_unit_diag) {
  	    T_element = T_i[j * incT + j * ldt * incT];
  
  
  	    temp1[0] = temp1[0] / T_element;
  	    temp1[1] = temp1[1] / T_element;
  
  	  }			/* if (diag == blas_non_unit_diag) */
  	  x_i[jx] = temp1[0];
  	  x_i[jx + 1] = temp1[1];
  	  jx += incx;
  	}			/* for j<n */
        }
      }
      break;
    case blas_prec_double:
    case blas_prec_indigenous:
  
      {
        int	      inc_intx;	/* inc for intx */
        double	      temp1[2];	/* temporary variable for calculations */
        double	      temp2[2];	/* temporary variable for calculations */
        double	      temp3[2];	/* temporary variable for calculations */
        double         *intx;	/* copy of x used for calculations */
  
        /* allocate space for intx */
        intx = (double *)blas_malloc(n * sizeof(double) * 2);
        if (n > 0 && intx == NULL) {
  	BLAS_error("blas_malloc", 0, 0, "malloc failed.\n");
        }
        /*
         * since intx is for internal usage, set it to 1 and then adjust it if
         * necessary
         */
        inc_intx = 1;
        inc_intx *= 2;
  
        /* copy x to intx */
        ix = start_x;
        jx = 0;
        for (i = 0; i < n; i++) {
  	temp1[0] = x_i[ix];
  	temp1[1] = x_i[1 + ix];
  	intx[jx] = temp1[0];
  	intx[1 + jx] = temp1[1];
  	ix += incx;
  	jx += inc_intx;
        }
  
        if ((order == blas_rowmajor &&
  	   trans == blas_no_trans && uplo == blas_upper) ||
  	  (order == blas_colmajor &&
  	   trans != blas_no_trans && uplo == blas_lower)) {
  
  	jx = (n - 1) * inc_intx;
  	for (j = n - 1; j >= 0; j--) {
  
  	  /*
  	   * compute Xj = alpha*Xj - SUM Aij(or Aji) * Xi i=j+1 to n-1
  	   */
  	  temp3[0] = intx[jx];
  	  temp3[1] = intx[1 + jx];
  	  /* multiply by alpha */
  	  {
  	    temp1[0] = (double)temp3[0] * alpha_i[0] - (double)temp3[1] * alpha_i[1];
  	    temp1[1] = (double)temp3[0] * alpha_i[1] + (double)temp3[1] * alpha_i[0];
  	  }
  
  	  ix = (n - 1) * inc_intx;
  	  for (i = n - 1; i >= j + 1; i--) {
  	    T_element = T_i[i * incT + j * ldt * incT];
  
  	    temp3[0] = intx[ix];
  	    temp3[1] = intx[1 + ix];
  	    {
  	      temp2[0] = temp3[0] * T_element;
  	      temp2[1] = temp3[1] * T_element;
  	    }
  	    temp1[0] = temp1[0] - temp2[0];
  	    temp1[1] = temp1[1] - temp2[1];
  	    ix -= inc_intx;
  	  }			/* for j<n */
  
  	  /*
  	   * if the diagonal entry is not equal to one, then divide Xj by the
  	   * entry
  	   */
  	  if (diag == blas_non_unit_diag) {
  	    T_element = T_i[j * incT + j * ldt * incT];
  
  
  	    {
  	      double	      dt = T_element;
  	      temp1[0] = temp1[0] / dt;
  	      temp1[1] = temp1[1] / dt;
  	    }
  
  	  }			/* if (diag == blas_non_unit_diag) */
  	  intx[jx] = temp1[0];
  	  intx[1 + jx] = temp1[1];
  
  	  jx -= inc_intx;
  	}			/* for j>=0 */
        } else if ((order == blas_rowmajor &&
  		  trans == blas_no_trans && uplo == blas_lower) ||
  		 (order == blas_colmajor &&
  		  trans != blas_no_trans && uplo == blas_upper)) {
  
  	jx = 0;
  	for (j = 0; j < n; j++) {
  
  	  /*
  	   * compute Xj = Xj - SUM Aij(or Aji) * Xi i=j+1 to n-1
  	   */
  	  temp3[0] = intx[jx];
  	  temp3[1] = intx[1 + jx];
  	  /* multiply by alpha */
  	  {
  	    temp1[0] = (double)temp3[0] * alpha_i[0] - (double)temp3[1] * alpha_i[1];
  	    temp1[1] = (double)temp3[0] * alpha_i[1] + (double)temp3[1] * alpha_i[0];
  	  }
  
  	  ix = 0;
  	  for (i = 0; i < j; i++) {
  	    T_element = T_i[i * incT + j * ldt * incT];
  
  	    temp3[0] = intx[ix];
  	    temp3[1] = intx[1 + ix];
  	    {
  	      temp2[0] = temp3[0] * T_element;
  	      temp2[1] = temp3[1] * T_element;
  	    }
  	    temp1[0] = temp1[0] - temp2[0];
  	    temp1[1] = temp1[1] - temp2[1];
  	    ix += inc_intx;
  	  }			/* for i<j */
  
  	  /*
  	   * if the diagonal entry is not equal to one, then divide Xj by the
  	   * entry
  	   */
  	  if (diag == blas_non_unit_diag) {
  	    T_element = T_i[j * incT + j * ldt * incT];
  
  
  	    {
  	      double	      dt = T_element;
  	      temp1[0] = temp1[0] / dt;
  	      temp1[1] = temp1[1] / dt;
  	    }
  
  	  }			/* if (diag == blas_non_unit_diag) */
  	  intx[jx] = temp1[0];
  	  intx[1 + jx] = temp1[1];
  	  jx += inc_intx;
  	}			/* for j<n */
        } else if ((order == blas_rowmajor &&
  		  trans != blas_no_trans && uplo == blas_lower) ||
  		 (order == blas_colmajor &&
  		  trans == blas_no_trans && uplo == blas_upper)) {
  
  	jx = (n - 1) * inc_intx;
  	for (j = n - 1; j >= 0; j--) {
  
  	  /*
  	   * compute Xj = alpha*Xj - SUM Aij(or Aji) * Xi i=j+1 to n-1
  	   */
  	  temp3[0] = intx[jx];
  	  temp3[1] = intx[1 + jx];
  	  /* multiply by alpha */
  	  {
  	    temp1[0] = (double)temp3[0] * alpha_i[0] - (double)temp3[1] * alpha_i[1];
  	    temp1[1] = (double)temp3[0] * alpha_i[1] + (double)temp3[1] * alpha_i[0];
  	  }
  
  	  ix = (n - 1) * inc_intx;
  	  for (i = n - 1; i >= j + 1; i--) {
  	    T_element = T_i[j * incT + i * ldt * incT];
  
  	    temp3[0] = intx[ix];
  	    temp3[1] = intx[1 + ix];
  	    {
  	      temp2[0] = temp3[0] * T_element;
  	      temp2[1] = temp3[1] * T_element;
  	    }
  	    temp1[0] = temp1[0] - temp2[0];
  	    temp1[1] = temp1[1] - temp2[1];
  	    ix -= inc_intx;
  	  }			/* for j<n */
  
  	  /*
  	   * if the diagonal entry is not equal to one, then divide Xj by the
  	   * entry
  	   */
  	  if (diag == blas_non_unit_diag) {
  	    T_element = T_i[j * incT + j * ldt * incT];
  
  
  	    {
  	      double	      dt = T_element;
  	      temp1[0] = temp1[0] / dt;
  	      temp1[1] = temp1[1] / dt;
  	    }
  
  	  }			/* if (diag == blas_non_unit_diag) */
  	  intx[jx] = temp1[0];
  	  intx[1 + jx] = temp1[1];
  
  	  jx -= inc_intx;
  	}			/* for j>=0 */
        } else if ((order == blas_rowmajor &&
  		  trans != blas_no_trans && uplo == blas_upper) ||
  		 (order == blas_colmajor &&
  		  trans == blas_no_trans && uplo == blas_lower)) {
  
  	jx = 0;
  	for (j = 0; j < n; j++) {
  
  	  /*
  	   * compute Xj = Xj - SUM Aij(or Aji) * Xi i=j+1 to n-1
  	   */
  	  temp3[0] = intx[jx];
  	  temp3[1] = intx[1 + jx];
  	  /* multiply by alpha */
  	  {
  	    temp1[0] = (double)temp3[0] * alpha_i[0] - (double)temp3[1] * alpha_i[1];
  	    temp1[1] = (double)temp3[0] * alpha_i[1] + (double)temp3[1] * alpha_i[0];
  	  }
  
  	  ix = 0;
  	  for (i = 0; i < j; i++) {
  	    T_element = T_i[j * incT + i * ldt * incT];
  
  	    temp3[0] = intx[ix];
  	    temp3[1] = intx[1 + ix];
  	    {
  	      temp2[0] = temp3[0] * T_element;
  	      temp2[1] = temp3[1] * T_element;
  	    }
  	    temp1[0] = temp1[0] - temp2[0];
  	    temp1[1] = temp1[1] - temp2[1];
  	    ix += inc_intx;
  	  }			/* for i<j */
  
  	  /*
  	   * if the diagonal entry is not equal to one, then divide Xj by the
  	   * entry
  	   */
  	  if (diag == blas_non_unit_diag) {
  	    T_element = T_i[j * incT + j * ldt * incT];
  
  
  	    {
  	      double	      dt = T_element;
  	      temp1[0] = temp1[0] / dt;
  	      temp1[1] = temp1[1] / dt;
  	    }
  
  	  }			/* if (diag == blas_non_unit_diag) */
  	  intx[jx] = temp1[0];
  	  intx[1 + jx] = temp1[1];
  	  jx += inc_intx;
  	}			/* for j<n */
        }
        /* copy the final results from intx to x */
        ix = start_x;
        jx = 0;
        for (i = 0; i < n; i++) {
  	temp1[0] = intx[jx];
  	temp1[1] = intx[1 + jx];
  	x_i[ix] = temp1[0];
  	x_i[ix + 1] = temp1[1];
  	ix += incx;
  	jx += inc_intx;
        }
  
        blas_free(intx);
      }
      break;
    case blas_prec_extra:
      {
        FPU_FIX_DECL;
        FPU_FIX_START;
        {
  	{
  	  int		  inc_intx;	/* inc for intx */
  	  double	  head_temp1[2], tail_temp1[2];	/* temporary variable
  							 * for calculations */
  	  double	  head_temp2[2], tail_temp2[2];	/* temporary variable
  							 * for calculations */
  	  double	  head_temp3[2], tail_temp3[2];	/* temporary variable
  							 * for calculations */
  	  double         *head_intx, *tail_intx;
  	  /* copy of x used for calculations */
  
  	  /* allocate space for intx */
  	  head_intx = (double *)blas_malloc(n * sizeof(double) * 2);
  	  tail_intx = (double *)blas_malloc(n * sizeof(double) * 2);
  	  if (n > 0 && (head_intx == NULL || tail_intx == NULL)) {
  	    BLAS_error("blas_malloc", 0, 0, "malloc failed.\n");
  	  }
  	  /*
  	   * since intx is for internal usage, set it to 1 and then adjust it
  	   * if necessary
  	   */
  	  inc_intx = 1;
  	  inc_intx *= 2;
  
  	  /* copy x to intx */
  	  ix = start_x;
  	  jx = 0;
  	  for (i = 0; i < n; i++) {
  	    head_temp1[0] = x_i[ix];
  	    tail_temp1[0] = 0.0;
  	    head_temp1[1] = x_i[1 + ix];
  	    tail_temp1[1] = 0.0;
  	    head_intx[jx] = head_temp1[0];
  	    tail_intx[jx] = tail_temp1[0];
  	    head_intx[1 + jx] = head_temp1[1];
  	    tail_intx[1 + jx] = tail_temp1[1];
  	    ix += incx;
  	    jx += inc_intx;
  	  }
  
  	  if ((order == blas_rowmajor &&
  	       trans == blas_no_trans && uplo == blas_upper) ||
  	      (order == blas_colmajor &&
  	       trans != blas_no_trans && uplo == blas_lower)) {
  
  	    jx = (n - 1) * inc_intx;
  	    for (j = n - 1; j >= 0; j--) {
  
  	      /*
  	       * compute Xj = alpha*Xj - SUM Aij(or Aji) * Xi i=j+1 to n-1
  	       */
  	      head_temp3[0] = head_intx[jx];
  	      head_temp3[1] = head_intx[1 + jx];
  	      tail_temp3[0] = tail_intx[jx];
  	      tail_temp3[1] = tail_intx[1 + jx];
  	      /* multiply by alpha */
  	      {
  		double		cd      [2];
  		cd[0] = (double)alpha_i[0];
  		cd[1] = (double)alpha_i[1];
  		{
  		  /* Compute complex-extra = complex-extra * complex-double. */
  		  double	  head_a0, tail_a0;
  		  double	  head_a1, tail_a1;
  		  double	  head_t1, tail_t1;
  		  double	  head_t2, tail_t2;
  		  head_a0 = head_temp3[0];
  		  tail_a0 = tail_temp3[0];
  		  head_a1 = head_temp3[1];
  		  tail_a1 = tail_temp3[1];
  		  /* real part */
  		  {
  		    /* Compute double-double = double-double * double. */
  		    double	    a11 , a21, b1, b2, c11, c21, c2, con, t1,
  		    		    t2;
  
  		    con = head_a0 * split;
  		    a11 = con - head_a0;
  		    a11 = con - a11;
  		    a21 = head_a0 - a11;
  		    con = cd[0] * split;
  		    b1 = con - cd[0];
  		    b1 = con - b1;
  		    b2 = cd[0] - b1;
  
  		    c11 = head_a0 * cd[0];
  		    c21 = (((a11 * b1 - c11) + a11 * b2) + a21 * b1) + a21 * b2;
  
  		    c2 = tail_a0 * cd[0];
  		    t1 = c11 + c2;
  		    t2 = (c2 - (t1 - c11)) + c21;
  
  		    head_t1 = t1 + t2;
  		    tail_t1 = t2 - (head_t1 - t1);
  		  }
  		  {
  		    /* Compute double-double = double-double * double. */
  		    double	    a11 , a21, b1, b2, c11, c21, c2, con, t1,
  		    		    t2;
  
  		    con = head_a1 * split;
  		    a11 = con - head_a1;
  		    a11 = con - a11;
  		    a21 = head_a1 - a11;
  		    con = cd[1] * split;
  		    b1 = con - cd[1];
  		    b1 = con - b1;
  		    b2 = cd[1] - b1;
  
  		    c11 = head_a1 * cd[1];
  		    c21 = (((a11 * b1 - c11) + a11 * b2) + a21 * b1) + a21 * b2;
  
  		    c2 = tail_a1 * cd[1];
  		    t1 = c11 + c2;
  		    t2 = (c2 - (t1 - c11)) + c21;
  
  		    head_t2 = t1 + t2;
  		    tail_t2 = t2 - (head_t2 - t1);
  		  }
  		  head_t2 = -head_t2;
  		  tail_t2 = -tail_t2;
  		  {
  		    /* Compute double-double = double-double + double-double. */
  		    double	    bv;
  		    double	    s1  , s2, t1, t2;
  
  		    /* Add two hi words. */
  		    s1 = head_t1 + head_t2;
  		    bv = s1 - head_t1;
  		    s2 = ((head_t2 - bv) + (head_t1 - (s1 - bv)));
  
  		    /* Add two lo words. */
  		    t1 = tail_t1 + tail_t2;
  		    bv = t1 - tail_t1;
  		    t2 = ((tail_t2 - bv) + (tail_t1 - (t1 - bv)));
  
  		    s2 += t1;
  
  		    /* Renormalize (s1, s2)  to  (t1, s2) */
  		    t1 = s1 + s2;
  		    s2 = s2 - (t1 - s1);
  
  		    t2 += s2;
  
  		    /* Renormalize (t1, t2)  */
  		    head_t1 = t1 + t2;
  		    tail_t1 = t2 - (head_t1 - t1);
  		  }
  		  head_temp1[0] = head_t1;
  		  tail_temp1[0] = tail_t1;
  		  /* imaginary part */
  		  {
  		    /* Compute double-double = double-double * double. */
  		    double	    a11 , a21, b1, b2, c11, c21, c2, con, t1,
  		    		    t2;
  
  		    con = head_a1 * split;
  		    a11 = con - head_a1;
  		    a11 = con - a11;
  		    a21 = head_a1 - a11;
  		    con = cd[0] * split;
  		    b1 = con - cd[0];
  		    b1 = con - b1;
  		    b2 = cd[0] - b1;
  
  		    c11 = head_a1 * cd[0];
  		    c21 = (((a11 * b1 - c11) + a11 * b2) + a21 * b1) + a21 * b2;
  
  		    c2 = tail_a1 * cd[0];
  		    t1 = c11 + c2;
  		    t2 = (c2 - (t1 - c11)) + c21;
  
  		    head_t1 = t1 + t2;
  		    tail_t1 = t2 - (head_t1 - t1);
  		  }
  		  {
  		    /* Compute double-double = double-double * double. */
  		    double	    a11 , a21, b1, b2, c11, c21, c2, con, t1,
  		    		    t2;
  
  		    con = head_a0 * split;
  		    a11 = con - head_a0;
  		    a11 = con - a11;
  		    a21 = head_a0 - a11;
  		    con = cd[1] * split;
  		    b1 = con - cd[1];
  		    b1 = con - b1;
  		    b2 = cd[1] - b1;
  
  		    c11 = head_a0 * cd[1];
  		    c21 = (((a11 * b1 - c11) + a11 * b2) + a21 * b1) + a21 * b2;
  
  		    c2 = tail_a0 * cd[1];
  		    t1 = c11 + c2;
  		    t2 = (c2 - (t1 - c11)) + c21;
  
  		    head_t2 = t1 + t2;
  		    tail_t2 = t2 - (head_t2 - t1);
  		  }
  		  {
  		    /* Compute double-double = double-double + double-double. */
  		    double	    bv;
  		    double	    s1  , s2, t1, t2;
  
  		    /* Add two hi words. */
  		    s1 = head_t1 + head_t2;
  		    bv = s1 - head_t1;
  		    s2 = ((head_t2 - bv) + (head_t1 - (s1 - bv)));
  
  		    /* Add two lo words. */
  		    t1 = tail_t1 + tail_t2;
  		    bv = t1 - tail_t1;
  		    t2 = ((tail_t2 - bv) + (tail_t1 - (t1 - bv)));
  
  		    s2 += t1;
  
  		    /* Renormalize (s1, s2)  to  (t1, s2) */
  		    t1 = s1 + s2;
  		    s2 = s2 - (t1 - s1);
  
  		    t2 += s2;
  
  		    /* Renormalize (t1, t2)  */
  		    head_t1 = t1 + t2;
  		    tail_t1 = t2 - (head_t1 - t1);
  		  }
  		  head_temp1[1] = head_t1;
  		  tail_temp1[1] = tail_t1;
  		}
  
  	      }
  
  	      ix = (n - 1) * inc_intx;
  	      for (i = n - 1; i >= j + 1; i--) {
  		T_element = T_i[i * incT + j * ldt * incT];
  
  		head_temp3[0] = head_intx[ix];
  		head_temp3[1] = head_intx[1 + ix];
  		tail_temp3[0] = tail_intx[ix];
  		tail_temp3[1] = tail_intx[1 + ix];
  		{
  		  double	  dt = (double)T_element;
  		  {
  		    /* Compute complex-extra = complex-extra * real. */
  		    double	    head_a0, tail_a0;
  		    double	    head_a1, tail_a1;
  		    double	    head_t, tail_t;
  		    head_a0 = head_temp3[0];
  		    tail_a0 = tail_temp3[0];
  		    head_a1 = head_temp3[1];
  		    tail_a1 = tail_temp3[1];
  		    {
  		      /* Compute double-double = double-double * double. */
  		      double	      a11, a21, b1, b2, c11, c21, c2, con,
  		      		      t1     , t2;
  
  		      con = head_a0 * split;
  		      a11 = con - head_a0;
  		      a11 = con - a11;
  		      a21 = head_a0 - a11;
  		      con = dt * split;
  		      b1 = con - dt;
  		      b1 = con - b1;
  		      b2 = dt - b1;
  
  		      c11 = head_a0 * dt;
  		      c21 = (((a11 * b1 - c11) + a11 * b2) + a21 * b1) + a21 * b2;
  
  		      c2 = tail_a0 * dt;
  		      t1 = c11 + c2;
  		      t2 = (c2 - (t1 - c11)) + c21;
  
  		      head_t = t1 + t2;
  		      tail_t = t2 - (head_t - t1);
  		    }
  		    head_temp2[0] = head_t;
  		    tail_temp2[0] = tail_t;
  		    {
  		      /* Compute double-double = double-double * double. */
  		      double	      a11, a21, b1, b2, c11, c21, c2, con,
  		      		      t1     , t2;
  
  		      con = head_a1 * split;
  		      a11 = con - head_a1;
  		      a11 = con - a11;
  		      a21 = head_a1 - a11;
  		      con = dt * split;
  		      b1 = con - dt;
  		      b1 = con - b1;
  		      b2 = dt - b1;
  
  		      c11 = head_a1 * dt;
  		      c21 = (((a11 * b1 - c11) + a11 * b2) + a21 * b1) + a21 * b2;
  
  		      c2 = tail_a1 * dt;
  		      t1 = c11 + c2;
  		      t2 = (c2 - (t1 - c11)) + c21;
  
  		      head_t = t1 + t2;
  		      tail_t = t2 - (head_t - t1);
  		    }
  		    head_temp2[1] = head_t;
  		    tail_temp2[1] = tail_t;
  		  }
  
  		}
  		{
  		  double	  head_at, tail_at;
  		  double	  head_bt, tail_bt;
  		  double	  head_ct, tail_ct;
  
  		  /* Real part */
  		  head_at = head_temp1[0];
  		  tail_at = tail_temp1[0];
  		  head_bt = -head_temp2[0];
  		  tail_bt = -tail_temp2[0];
  		  {
  		    /* Compute double-double = double-double + double-double. */
  		    double	    bv;
  		    double	    s1  , s2, t1, t2;
  
  		    /* Add two hi words. */
  		    s1 = head_at + head_bt;
  		    bv = s1 - head_at;
  		    s2 = ((head_bt - bv) + (head_at - (s1 - bv)));
  
  		    /* Add two lo words. */
  		    t1 = tail_at + tail_bt;
  		    bv = t1 - tail_at;
  		    t2 = ((tail_bt - bv) + (tail_at - (t1 - bv)));
  
  		    s2 += t1;
  
  		    /* Renormalize (s1, s2)  to  (t1, s2) */
  		    t1 = s1 + s2;
  		    s2 = s2 - (t1 - s1);
  
  		    t2 += s2;
  
  		    /* Renormalize (t1, t2)  */
  		    head_ct = t1 + t2;
  		    tail_ct = t2 - (head_ct - t1);
  		  }
  		  head_temp1[0] = head_ct;
  		  tail_temp1[0] = tail_ct;
  		  /* Imaginary part */
  		  head_at = head_temp1[1];
  		  tail_at = tail_temp1[1];
  		  head_bt = -head_temp2[1];
  		  tail_bt = -tail_temp2[1];
  		  {
  		    /* Compute double-double = double-double + double-double. */
  		    double	    bv;
  		    double	    s1  , s2, t1, t2;
  
  		    /* Add two hi words. */
  		    s1 = head_at + head_bt;
  		    bv = s1 - head_at;
  		    s2 = ((head_bt - bv) + (head_at - (s1 - bv)));
  
  		    /* Add two lo words. */
  		    t1 = tail_at + tail_bt;
  		    bv = t1 - tail_at;
  		    t2 = ((tail_bt - bv) + (tail_at - (t1 - bv)));
  
  		    s2 += t1;
  
  		    /* Renormalize (s1, s2)  to  (t1, s2) */
  		    t1 = s1 + s2;
  		    s2 = s2 - (t1 - s1);
  
  		    t2 += s2;
  
  		    /* Renormalize (t1, t2)  */
  		    head_ct = t1 + t2;
  		    tail_ct = t2 - (head_ct - t1);
  		  }
  		  head_temp1[1] = head_ct;
  		  tail_temp1[1] = tail_ct;
  		}
  		ix -= inc_intx;
  	      }			/* for j<n */
  
  	      /*
  	       * if the diagonal entry is not equal to one, then divide Xj by
  	       * the entry
  	       */
  	      if (diag == blas_non_unit_diag) {
  		T_element = T_i[j * incT + j * ldt * incT];
  
  
  		{
  		  double	  head_a, tail_a;
  		  double	  head_b, tail_b;
  		  head_a = head_temp1[0];
  		  tail_a = tail_temp1[0];
  		  {
  		    /*
  		     * Compute double-double = double-double / double, using
  		     * a Newton iteration scheme.
  		     */
  		    double	    b1  , b2, con, e, t1, t2, t11, t21, t12,
  		    		    t22;
  
  		    /* Compute a DP approximation to the quotient. */
  		    t1 = head_a / T_element;
  
  		    /*
  		     * Split t1 and b into two parts with at most 26 bits
  		     * each, using the Dekker-Veltkamp method.
  		     */
  		    con = t1 * split;
  		    t11 = con - (con - t1);
  		    t21 = t1 - t11;
  		    con = T_element * split;
  		    b1 = con - (con - T_element);
  		    b2 = T_element - b1;
  
  		    /* Compute t1 * b using Dekker method. */
  		    t12 = t1 * T_element;
  		    t22 = (((t11 * b1 - t12) + t11 * b2) + t21 * b1) + t21 * b2;
  
  		    /* Compute dda - (t12, t22) using Knuth trick. */
  		    t11 = head_a - t12;
  		    e = t11 - head_a;
  		    t21 = ((-t12 - e) + (head_a - (t11 - e))) + tail_a - t22;
  
  		    /* Compute high-order word of (t11, t21) and divide by b. */
  		    t2 = (t11 + t21) / T_element;
  
  		    /* The result is t1 + t2, after normalization. */
  		    head_b = t1 + t2;
  		    tail_b = t2 - (head_b - t1);
  		  }
  		  head_temp1[0] = head_b;
  		  tail_temp1[0] = tail_b;
  		  head_a = head_temp1[1];
  		  tail_a = tail_temp1[1];
  		  {
  		    /*
  		     * Compute double-double = double-double / double, using
  		     * a Newton iteration scheme.
  		     */
  		    double	    b1  , b2, con, e, t1, t2, t11, t21, t12,
  		    		    t22;
  
  		    /* Compute a DP approximation to the quotient. */
  		    t1 = head_a / T_element;
  
  		    /*
  		     * Split t1 and b into two parts with at most 26 bits
  		     * each, using the Dekker-Veltkamp method.
  		     */
  		    con = t1 * split;
  		    t11 = con - (con - t1);
  		    t21 = t1 - t11;
  		    con = T_element * split;
  		    b1 = con - (con - T_element);
  		    b2 = T_element - b1;
  
  		    /* Compute t1 * b using Dekker method. */
  		    t12 = t1 * T_element;
  		    t22 = (((t11 * b1 - t12) + t11 * b2) + t21 * b1) + t21 * b2;
  
  		    /* Compute dda - (t12, t22) using Knuth trick. */
  		    t11 = head_a - t12;
  		    e = t11 - head_a;
  		    t21 = ((-t12 - e) + (head_a - (t11 - e))) + tail_a - t22;
  
  		    /* Compute high-order word of (t11, t21) and divide by b. */
  		    t2 = (t11 + t21) / T_element;
  
  		    /* The result is t1 + t2, after normalization. */
  		    head_b = t1 + t2;
  		    tail_b = t2 - (head_b - t1);
  		  }
  		  head_temp1[1] = head_b;
  		  tail_temp1[1] = tail_b;
  		}
  
  	      }			/* if (diag == blas_non_unit_diag) */
  	      head_intx[jx] = head_temp1[0];
  	      tail_intx[jx] = tail_temp1[0];
  	      head_intx[1 + jx] = head_temp1[1];
  	      tail_intx[1 + jx] = tail_temp1[1];
  
  	      jx -= inc_intx;
  	    }			/* for j>=0 */
  	  } else if ((order == blas_rowmajor &&
  		      trans == blas_no_trans && uplo == blas_lower) ||
  		     (order == blas_colmajor &&
  		      trans != blas_no_trans && uplo == blas_upper)) {
  
  	    jx = 0;
  	    for (j = 0; j < n; j++) {
  
  	      /*
  	       * compute Xj = Xj - SUM Aij(or Aji) * Xi i=j+1 to n-1
  	       */
  	      head_temp3[0] = head_intx[jx];
  	      head_temp3[1] = head_intx[1 + jx];
  	      tail_temp3[0] = tail_intx[jx];
  	      tail_temp3[1] = tail_intx[1 + jx];
  	      /* multiply by alpha */
  	      {
  		double		cd      [2];
  		cd[0] = (double)alpha_i[0];
  		cd[1] = (double)alpha_i[1];
  		{
  		  /* Compute complex-extra = complex-extra * complex-double. */
  		  double	  head_a0, tail_a0;
  		  double	  head_a1, tail_a1;
  		  double	  head_t1, tail_t1;
  		  double	  head_t2, tail_t2;
  		  head_a0 = head_temp3[0];
  		  tail_a0 = tail_temp3[0];
  		  head_a1 = head_temp3[1];
  		  tail_a1 = tail_temp3[1];
  		  /* real part */
  		  {
  		    /* Compute double-double = double-double * double. */
  		    double	    a11 , a21, b1, b2, c11, c21, c2, con, t1,
  		    		    t2;
  
  		    con = head_a0 * split;
  		    a11 = con - head_a0;
  		    a11 = con - a11;
  		    a21 = head_a0 - a11;
  		    con = cd[0] * split;
  		    b1 = con - cd[0];
  		    b1 = con - b1;
  		    b2 = cd[0] - b1;
  
  		    c11 = head_a0 * cd[0];
  		    c21 = (((a11 * b1 - c11) + a11 * b2) + a21 * b1) + a21 * b2;
  
  		    c2 = tail_a0 * cd[0];
  		    t1 = c11 + c2;
  		    t2 = (c2 - (t1 - c11)) + c21;
  
  		    head_t1 = t1 + t2;
  		    tail_t1 = t2 - (head_t1 - t1);
  		  }
  		  {
  		    /* Compute double-double = double-double * double. */
  		    double	    a11 , a21, b1, b2, c11, c21, c2, con, t1,
  		    		    t2;
  
  		    con = head_a1 * split;
  		    a11 = con - head_a1;
  		    a11 = con - a11;
  		    a21 = head_a1 - a11;
  		    con = cd[1] * split;
  		    b1 = con - cd[1];
  		    b1 = con - b1;
  		    b2 = cd[1] - b1;
  
  		    c11 = head_a1 * cd[1];
  		    c21 = (((a11 * b1 - c11) + a11 * b2) + a21 * b1) + a21 * b2;
  
  		    c2 = tail_a1 * cd[1];
  		    t1 = c11 + c2;
  		    t2 = (c2 - (t1 - c11)) + c21;
  
  		    head_t2 = t1 + t2;
  		    tail_t2 = t2 - (head_t2 - t1);
  		  }
  		  head_t2 = -head_t2;
  		  tail_t2 = -tail_t2;
  		  {
  		    /* Compute double-double = double-double + double-double. */
  		    double	    bv;
  		    double	    s1  , s2, t1, t2;
  
  		    /* Add two hi words. */
  		    s1 = head_t1 + head_t2;
  		    bv = s1 - head_t1;
  		    s2 = ((head_t2 - bv) + (head_t1 - (s1 - bv)));
  
  		    /* Add two lo words. */
  		    t1 = tail_t1 + tail_t2;
  		    bv = t1 - tail_t1;
  		    t2 = ((tail_t2 - bv) + (tail_t1 - (t1 - bv)));
  
  		    s2 += t1;
  
  		    /* Renormalize (s1, s2)  to  (t1, s2) */
  		    t1 = s1 + s2;
  		    s2 = s2 - (t1 - s1);
  
  		    t2 += s2;
  
  		    /* Renormalize (t1, t2)  */
  		    head_t1 = t1 + t2;
  		    tail_t1 = t2 - (head_t1 - t1);
  		  }
  		  head_temp1[0] = head_t1;
  		  tail_temp1[0] = tail_t1;
  		  /* imaginary part */
  		  {
  		    /* Compute double-double = double-double * double. */
  		    double	    a11 , a21, b1, b2, c11, c21, c2, con, t1,
  		    		    t2;
  
  		    con = head_a1 * split;
  		    a11 = con - head_a1;
  		    a11 = con - a11;
  		    a21 = head_a1 - a11;
  		    con = cd[0] * split;
  		    b1 = con - cd[0];
  		    b1 = con - b1;
  		    b2 = cd[0] - b1;
  
  		    c11 = head_a1 * cd[0];
  		    c21 = (((a11 * b1 - c11) + a11 * b2) + a21 * b1) + a21 * b2;
  
  		    c2 = tail_a1 * cd[0];
  		    t1 = c11 + c2;
  		    t2 = (c2 - (t1 - c11)) + c21;
  
  		    head_t1 = t1 + t2;
  		    tail_t1 = t2 - (head_t1 - t1);
  		  }
  		  {
  		    /* Compute double-double = double-double * double. */
  		    double	    a11 , a21, b1, b2, c11, c21, c2, con, t1,
  		    		    t2;
  
  		    con = head_a0 * split;
  		    a11 = con - head_a0;
  		    a11 = con - a11;
  		    a21 = head_a0 - a11;
  		    con = cd[1] * split;
  		    b1 = con - cd[1];
  		    b1 = con - b1;
  		    b2 = cd[1] - b1;
  
  		    c11 = head_a0 * cd[1];
  		    c21 = (((a11 * b1 - c11) + a11 * b2) + a21 * b1) + a21 * b2;
  
  		    c2 = tail_a0 * cd[1];
  		    t1 = c11 + c2;
  		    t2 = (c2 - (t1 - c11)) + c21;
  
  		    head_t2 = t1 + t2;
  		    tail_t2 = t2 - (head_t2 - t1);
  		  }
  		  {
  		    /* Compute double-double = double-double + double-double. */
  		    double	    bv;
  		    double	    s1  , s2, t1, t2;
  
  		    /* Add two hi words. */
  		    s1 = head_t1 + head_t2;
  		    bv = s1 - head_t1;
  		    s2 = ((head_t2 - bv) + (head_t1 - (s1 - bv)));
  
  		    /* Add two lo words. */
  		    t1 = tail_t1 + tail_t2;
  		    bv = t1 - tail_t1;
  		    t2 = ((tail_t2 - bv) + (tail_t1 - (t1 - bv)));
  
  		    s2 += t1;
  
  		    /* Renormalize (s1, s2)  to  (t1, s2) */
  		    t1 = s1 + s2;
  		    s2 = s2 - (t1 - s1);
  
  		    t2 += s2;
  
  		    /* Renormalize (t1, t2)  */
  		    head_t1 = t1 + t2;
  		    tail_t1 = t2 - (head_t1 - t1);
  		  }
  		  head_temp1[1] = head_t1;
  		  tail_temp1[1] = tail_t1;
  		}
  
  	      }
  
  	      ix = 0;
  	      for (i = 0; i < j; i++) {
  		T_element = T_i[i * incT + j * ldt * incT];
  
  		head_temp3[0] = head_intx[ix];
  		head_temp3[1] = head_intx[1 + ix];
  		tail_temp3[0] = tail_intx[ix];
  		tail_temp3[1] = tail_intx[1 + ix];
  		{
  		  double	  dt = (double)T_element;
  		  {
  		    /* Compute complex-extra = complex-extra * real. */
  		    double	    head_a0, tail_a0;
  		    double	    head_a1, tail_a1;
  		    double	    head_t, tail_t;
  		    head_a0 = head_temp3[0];
  		    tail_a0 = tail_temp3[0];
  		    head_a1 = head_temp3[1];
  		    tail_a1 = tail_temp3[1];
  		    {
  		      /* Compute double-double = double-double * double. */
  		      double	      a11, a21, b1, b2, c11, c21, c2, con,
  		      		      t1     , t2;
  
  		      con = head_a0 * split;
  		      a11 = con - head_a0;
  		      a11 = con - a11;
  		      a21 = head_a0 - a11;
  		      con = dt * split;
  		      b1 = con - dt;
  		      b1 = con - b1;
  		      b2 = dt - b1;
  
  		      c11 = head_a0 * dt;
  		      c21 = (((a11 * b1 - c11) + a11 * b2) + a21 * b1) + a21 * b2;
  
  		      c2 = tail_a0 * dt;
  		      t1 = c11 + c2;
  		      t2 = (c2 - (t1 - c11)) + c21;
  
  		      head_t = t1 + t2;
  		      tail_t = t2 - (head_t - t1);
  		    }
  		    head_temp2[0] = head_t;
  		    tail_temp2[0] = tail_t;
  		    {
  		      /* Compute double-double = double-double * double. */
  		      double	      a11, a21, b1, b2, c11, c21, c2, con,
  		      		      t1     , t2;
  
  		      con = head_a1 * split;
  		      a11 = con - head_a1;
  		      a11 = con - a11;
  		      a21 = head_a1 - a11;
  		      con = dt * split;
  		      b1 = con - dt;
  		      b1 = con - b1;
  		      b2 = dt - b1;
  
  		      c11 = head_a1 * dt;
  		      c21 = (((a11 * b1 - c11) + a11 * b2) + a21 * b1) + a21 * b2;
  
  		      c2 = tail_a1 * dt;
  		      t1 = c11 + c2;
  		      t2 = (c2 - (t1 - c11)) + c21;
  
  		      head_t = t1 + t2;
  		      tail_t = t2 - (head_t - t1);
  		    }
  		    head_temp2[1] = head_t;
  		    tail_temp2[1] = tail_t;
  		  }
  
  		}
  		{
  		  double	  head_at, tail_at;
  		  double	  head_bt, tail_bt;
  		  double	  head_ct, tail_ct;
  
  		  /* Real part */
  		  head_at = head_temp1[0];
  		  tail_at = tail_temp1[0];
  		  head_bt = -head_temp2[0];
  		  tail_bt = -tail_temp2[0];
  		  {
  		    /* Compute double-double = double-double + double-double. */
  		    double	    bv;
  		    double	    s1  , s2, t1, t2;
  
  		    /* Add two hi words. */
  		    s1 = head_at + head_bt;
  		    bv = s1 - head_at;
  		    s2 = ((head_bt - bv) + (head_at - (s1 - bv)));
  
  		    /* Add two lo words. */
  		    t1 = tail_at + tail_bt;
  		    bv = t1 - tail_at;
  		    t2 = ((tail_bt - bv) + (tail_at - (t1 - bv)));
  
  		    s2 += t1;
  
  		    /* Renormalize (s1, s2)  to  (t1, s2) */
  		    t1 = s1 + s2;
  		    s2 = s2 - (t1 - s1);
  
  		    t2 += s2;
  
  		    /* Renormalize (t1, t2)  */
  		    head_ct = t1 + t2;
  		    tail_ct = t2 - (head_ct - t1);
  		  }
  		  head_temp1[0] = head_ct;
  		  tail_temp1[0] = tail_ct;
  		  /* Imaginary part */
  		  head_at = head_temp1[1];
  		  tail_at = tail_temp1[1];
  		  head_bt = -head_temp2[1];
  		  tail_bt = -tail_temp2[1];
  		  {
  		    /* Compute double-double = double-double + double-double. */
  		    double	    bv;
  		    double	    s1  , s2, t1, t2;
  
  		    /* Add two hi words. */
  		    s1 = head_at + head_bt;
  		    bv = s1 - head_at;
  		    s2 = ((head_bt - bv) + (head_at - (s1 - bv)));
  
  		    /* Add two lo words. */
  		    t1 = tail_at + tail_bt;
  		    bv = t1 - tail_at;
  		    t2 = ((tail_bt - bv) + (tail_at - (t1 - bv)));
  
  		    s2 += t1;
  
  		    /* Renormalize (s1, s2)  to  (t1, s2) */
  		    t1 = s1 + s2;
  		    s2 = s2 - (t1 - s1);
  
  		    t2 += s2;
  
  		    /* Renormalize (t1, t2)  */
  		    head_ct = t1 + t2;
  		    tail_ct = t2 - (head_ct - t1);
  		  }
  		  head_temp1[1] = head_ct;
  		  tail_temp1[1] = tail_ct;
  		}
  		ix += inc_intx;
  	      }			/* for i<j */
  
  	      /*
  	       * if the diagonal entry is not equal to one, then divide Xj by
  	       * the entry
  	       */
  	      if (diag == blas_non_unit_diag) {
  		T_element = T_i[j * incT + j * ldt * incT];
  
  
  		{
  		  double	  head_a, tail_a;
  		  double	  head_b, tail_b;
  		  head_a = head_temp1[0];
  		  tail_a = tail_temp1[0];
  		  {
  		    /*
  		     * Compute double-double = double-double / double, using
  		     * a Newton iteration scheme.
  		     */
  		    double	    b1  , b2, con, e, t1, t2, t11, t21, t12,
  		    		    t22;
  
  		    /* Compute a DP approximation to the quotient. */
  		    t1 = head_a / T_element;
  
  		    /*
  		     * Split t1 and b into two parts with at most 26 bits
  		     * each, using the Dekker-Veltkamp method.
  		     */
  		    con = t1 * split;
  		    t11 = con - (con - t1);
  		    t21 = t1 - t11;
  		    con = T_element * split;
  		    b1 = con - (con - T_element);
  		    b2 = T_element - b1;
  
  		    /* Compute t1 * b using Dekker method. */
  		    t12 = t1 * T_element;
  		    t22 = (((t11 * b1 - t12) + t11 * b2) + t21 * b1) + t21 * b2;
  
  		    /* Compute dda - (t12, t22) using Knuth trick. */
  		    t11 = head_a - t12;
  		    e = t11 - head_a;
  		    t21 = ((-t12 - e) + (head_a - (t11 - e))) + tail_a - t22;
  
  		    /* Compute high-order word of (t11, t21) and divide by b. */
  		    t2 = (t11 + t21) / T_element;
  
  		    /* The result is t1 + t2, after normalization. */
  		    head_b = t1 + t2;
  		    tail_b = t2 - (head_b - t1);
  		  }
  		  head_temp1[0] = head_b;
  		  tail_temp1[0] = tail_b;
  		  head_a = head_temp1[1];
  		  tail_a = tail_temp1[1];
  		  {
  		    /*
  		     * Compute double-double = double-double / double, using
  		     * a Newton iteration scheme.
  		     */
  		    double	    b1  , b2, con, e, t1, t2, t11, t21, t12,
  		    		    t22;
  
  		    /* Compute a DP approximation to the quotient. */
  		    t1 = head_a / T_element;
  
  		    /*
  		     * Split t1 and b into two parts with at most 26 bits
  		     * each, using the Dekker-Veltkamp method.
  		     */
  		    con = t1 * split;
  		    t11 = con - (con - t1);
  		    t21 = t1 - t11;
  		    con = T_element * split;
  		    b1 = con - (con - T_element);
  		    b2 = T_element - b1;
  
  		    /* Compute t1 * b using Dekker method. */
  		    t12 = t1 * T_element;
  		    t22 = (((t11 * b1 - t12) + t11 * b2) + t21 * b1) + t21 * b2;
  
  		    /* Compute dda - (t12, t22) using Knuth trick. */
  		    t11 = head_a - t12;
  		    e = t11 - head_a;
  		    t21 = ((-t12 - e) + (head_a - (t11 - e))) + tail_a - t22;
  
  		    /* Compute high-order word of (t11, t21) and divide by b. */
  		    t2 = (t11 + t21) / T_element;
  
  		    /* The result is t1 + t2, after normalization. */
  		    head_b = t1 + t2;
  		    tail_b = t2 - (head_b - t1);
  		  }
  		  head_temp1[1] = head_b;
  		  tail_temp1[1] = tail_b;
  		}
  
  	      }			/* if (diag == blas_non_unit_diag) */
  	      head_intx[jx] = head_temp1[0];
  	      tail_intx[jx] = tail_temp1[0];
  	      head_intx[1 + jx] = head_temp1[1];
  	      tail_intx[1 + jx] = tail_temp1[1];
  	      jx += inc_intx;
  	    }			/* for j<n */
  	  } else if ((order == blas_rowmajor &&
  		      trans != blas_no_trans && uplo == blas_lower) ||
  		     (order == blas_colmajor &&
  		      trans == blas_no_trans && uplo == blas_upper)) {
  
  	    jx = (n - 1) * inc_intx;
  	    for (j = n - 1; j >= 0; j--) {
  
  	      /*
  	       * compute Xj = alpha*Xj - SUM Aij(or Aji) * Xi i=j+1 to n-1
  	       */
  	      head_temp3[0] = head_intx[jx];
  	      head_temp3[1] = head_intx[1 + jx];
  	      tail_temp3[0] = tail_intx[jx];
  	      tail_temp3[1] = tail_intx[1 + jx];
  	      /* multiply by alpha */
  	      {
  		double		cd      [2];
  		cd[0] = (double)alpha_i[0];
  		cd[1] = (double)alpha_i[1];
  		{
  		  /* Compute complex-extra = complex-extra * complex-double. */
  		  double	  head_a0, tail_a0;
  		  double	  head_a1, tail_a1;
  		  double	  head_t1, tail_t1;
  		  double	  head_t2, tail_t2;
  		  head_a0 = head_temp3[0];
  		  tail_a0 = tail_temp3[0];
  		  head_a1 = head_temp3[1];
  		  tail_a1 = tail_temp3[1];
  		  /* real part */
  		  {
  		    /* Compute double-double = double-double * double. */
  		    double	    a11 , a21, b1, b2, c11, c21, c2, con, t1,
  		    		    t2;
  
  		    con = head_a0 * split;
  		    a11 = con - head_a0;
  		    a11 = con - a11;
  		    a21 = head_a0 - a11;
  		    con = cd[0] * split;
  		    b1 = con - cd[0];
  		    b1 = con - b1;
  		    b2 = cd[0] - b1;
  
  		    c11 = head_a0 * cd[0];
  		    c21 = (((a11 * b1 - c11) + a11 * b2) + a21 * b1) + a21 * b2;
  
  		    c2 = tail_a0 * cd[0];
  		    t1 = c11 + c2;
  		    t2 = (c2 - (t1 - c11)) + c21;
  
  		    head_t1 = t1 + t2;
  		    tail_t1 = t2 - (head_t1 - t1);
  		  }
  		  {
  		    /* Compute double-double = double-double * double. */
  		    double	    a11 , a21, b1, b2, c11, c21, c2, con, t1,
  		    		    t2;
  
  		    con = head_a1 * split;
  		    a11 = con - head_a1;
  		    a11 = con - a11;
  		    a21 = head_a1 - a11;
  		    con = cd[1] * split;
  		    b1 = con - cd[1];
  		    b1 = con - b1;
  		    b2 = cd[1] - b1;
  
  		    c11 = head_a1 * cd[1];
  		    c21 = (((a11 * b1 - c11) + a11 * b2) + a21 * b1) + a21 * b2;
  
  		    c2 = tail_a1 * cd[1];
  		    t1 = c11 + c2;
  		    t2 = (c2 - (t1 - c11)) + c21;
  
  		    head_t2 = t1 + t2;
  		    tail_t2 = t2 - (head_t2 - t1);
  		  }
  		  head_t2 = -head_t2;
  		  tail_t2 = -tail_t2;
  		  {
  		    /* Compute double-double = double-double + double-double. */
  		    double	    bv;
  		    double	    s1  , s2, t1, t2;
  
  		    /* Add two hi words. */
  		    s1 = head_t1 + head_t2;
  		    bv = s1 - head_t1;
  		    s2 = ((head_t2 - bv) + (head_t1 - (s1 - bv)));
  
  		    /* Add two lo words. */
  		    t1 = tail_t1 + tail_t2;
  		    bv = t1 - tail_t1;
  		    t2 = ((tail_t2 - bv) + (tail_t1 - (t1 - bv)));
  
  		    s2 += t1;
  
  		    /* Renormalize (s1, s2)  to  (t1, s2) */
  		    t1 = s1 + s2;
  		    s2 = s2 - (t1 - s1);
  
  		    t2 += s2;
  
  		    /* Renormalize (t1, t2)  */
  		    head_t1 = t1 + t2;
  		    tail_t1 = t2 - (head_t1 - t1);
  		  }
  		  head_temp1[0] = head_t1;
  		  tail_temp1[0] = tail_t1;
  		  /* imaginary part */
  		  {
  		    /* Compute double-double = double-double * double. */
  		    double	    a11 , a21, b1, b2, c11, c21, c2, con, t1,
  		    		    t2;
  
  		    con = head_a1 * split;
  		    a11 = con - head_a1;
  		    a11 = con - a11;
  		    a21 = head_a1 - a11;
  		    con = cd[0] * split;
  		    b1 = con - cd[0];
  		    b1 = con - b1;
  		    b2 = cd[0] - b1;
  
  		    c11 = head_a1 * cd[0];
  		    c21 = (((a11 * b1 - c11) + a11 * b2) + a21 * b1) + a21 * b2;
  
  		    c2 = tail_a1 * cd[0];
  		    t1 = c11 + c2;
  		    t2 = (c2 - (t1 - c11)) + c21;
  
  		    head_t1 = t1 + t2;
  		    tail_t1 = t2 - (head_t1 - t1);
  		  }
  		  {
  		    /* Compute double-double = double-double * double. */
  		    double	    a11 , a21, b1, b2, c11, c21, c2, con, t1,
  		    		    t2;
  
  		    con = head_a0 * split;
  		    a11 = con - head_a0;
  		    a11 = con - a11;
  		    a21 = head_a0 - a11;
  		    con = cd[1] * split;
  		    b1 = con - cd[1];
  		    b1 = con - b1;
  		    b2 = cd[1] - b1;
  
  		    c11 = head_a0 * cd[1];
  		    c21 = (((a11 * b1 - c11) + a11 * b2) + a21 * b1) + a21 * b2;
  
  		    c2 = tail_a0 * cd[1];
  		    t1 = c11 + c2;
  		    t2 = (c2 - (t1 - c11)) + c21;
  
  		    head_t2 = t1 + t2;
  		    tail_t2 = t2 - (head_t2 - t1);
  		  }
  		  {
  		    /* Compute double-double = double-double + double-double. */
  		    double	    bv;
  		    double	    s1  , s2, t1, t2;
  
  		    /* Add two hi words. */
  		    s1 = head_t1 + head_t2;
  		    bv = s1 - head_t1;
  		    s2 = ((head_t2 - bv) + (head_t1 - (s1 - bv)));
  
  		    /* Add two lo words. */
  		    t1 = tail_t1 + tail_t2;
  		    bv = t1 - tail_t1;
  		    t2 = ((tail_t2 - bv) + (tail_t1 - (t1 - bv)));
  
  		    s2 += t1;
  
  		    /* Renormalize (s1, s2)  to  (t1, s2) */
  		    t1 = s1 + s2;
  		    s2 = s2 - (t1 - s1);
  
  		    t2 += s2;
  
  		    /* Renormalize (t1, t2)  */
  		    head_t1 = t1 + t2;
  		    tail_t1 = t2 - (head_t1 - t1);
  		  }
  		  head_temp1[1] = head_t1;
  		  tail_temp1[1] = tail_t1;
  		}
  
  	      }
  
  	      ix = (n - 1) * inc_intx;
  	      for (i = n - 1; i >= j + 1; i--) {
  		T_element = T_i[j * incT + i * ldt * incT];
  
  		head_temp3[0] = head_intx[ix];
  		head_temp3[1] = head_intx[1 + ix];
  		tail_temp3[0] = tail_intx[ix];
  		tail_temp3[1] = tail_intx[1 + ix];
  		{
  		  double	  dt = (double)T_element;
  		  {
  		    /* Compute complex-extra = complex-extra * real. */
  		    double	    head_a0, tail_a0;
  		    double	    head_a1, tail_a1;
  		    double	    head_t, tail_t;
  		    head_a0 = head_temp3[0];
  		    tail_a0 = tail_temp3[0];
  		    head_a1 = head_temp3[1];
  		    tail_a1 = tail_temp3[1];
  		    {
  		      /* Compute double-double = double-double * double. */
  		      double	      a11, a21, b1, b2, c11, c21, c2, con,
  		      		      t1     , t2;
  
  		      con = head_a0 * split;
  		      a11 = con - head_a0;
  		      a11 = con - a11;
  		      a21 = head_a0 - a11;
  		      con = dt * split;
  		      b1 = con - dt;
  		      b1 = con - b1;
  		      b2 = dt - b1;
  
  		      c11 = head_a0 * dt;
  		      c21 = (((a11 * b1 - c11) + a11 * b2) + a21 * b1) + a21 * b2;
  
  		      c2 = tail_a0 * dt;
  		      t1 = c11 + c2;
  		      t2 = (c2 - (t1 - c11)) + c21;
  
  		      head_t = t1 + t2;
  		      tail_t = t2 - (head_t - t1);
  		    }
  		    head_temp2[0] = head_t;
  		    tail_temp2[0] = tail_t;
  		    {
  		      /* Compute double-double = double-double * double. */
  		      double	      a11, a21, b1, b2, c11, c21, c2, con,
  		      		      t1     , t2;
  
  		      con = head_a1 * split;
  		      a11 = con - head_a1;
  		      a11 = con - a11;
  		      a21 = head_a1 - a11;
  		      con = dt * split;
  		      b1 = con - dt;
  		      b1 = con - b1;
  		      b2 = dt - b1;
  
  		      c11 = head_a1 * dt;
  		      c21 = (((a11 * b1 - c11) + a11 * b2) + a21 * b1) + a21 * b2;
  
  		      c2 = tail_a1 * dt;
  		      t1 = c11 + c2;
  		      t2 = (c2 - (t1 - c11)) + c21;
  
  		      head_t = t1 + t2;
  		      tail_t = t2 - (head_t - t1);
  		    }
  		    head_temp2[1] = head_t;
  		    tail_temp2[1] = tail_t;
  		  }
  
  		}
  		{
  		  double	  head_at, tail_at;
  		  double	  head_bt, tail_bt;
  		  double	  head_ct, tail_ct;
  
  		  /* Real part */
  		  head_at = head_temp1[0];
  		  tail_at = tail_temp1[0];
  		  head_bt = -head_temp2[0];
  		  tail_bt = -tail_temp2[0];
  		  {
  		    /* Compute double-double = double-double + double-double. */
  		    double	    bv;
  		    double	    s1  , s2, t1, t2;
  
  		    /* Add two hi words. */
  		    s1 = head_at + head_bt;
  		    bv = s1 - head_at;
  		    s2 = ((head_bt - bv) + (head_at - (s1 - bv)));
  
  		    /* Add two lo words. */
  		    t1 = tail_at + tail_bt;
  		    bv = t1 - tail_at;
  		    t2 = ((tail_bt - bv) + (tail_at - (t1 - bv)));
  
  		    s2 += t1;
  
  		    /* Renormalize (s1, s2)  to  (t1, s2) */
  		    t1 = s1 + s2;
  		    s2 = s2 - (t1 - s1);
  
  		    t2 += s2;
  
  		    /* Renormalize (t1, t2)  */
  		    head_ct = t1 + t2;
  		    tail_ct = t2 - (head_ct - t1);
  		  }
  		  head_temp1[0] = head_ct;
  		  tail_temp1[0] = tail_ct;
  		  /* Imaginary part */
  		  head_at = head_temp1[1];
  		  tail_at = tail_temp1[1];
  		  head_bt = -head_temp2[1];
  		  tail_bt = -tail_temp2[1];
  		  {
  		    /* Compute double-double = double-double + double-double. */
  		    double	    bv;
  		    double	    s1  , s2, t1, t2;
  
  		    /* Add two hi words. */
  		    s1 = head_at + head_bt;
  		    bv = s1 - head_at;
  		    s2 = ((head_bt - bv) + (head_at - (s1 - bv)));
  
  		    /* Add two lo words. */
  		    t1 = tail_at + tail_bt;
  		    bv = t1 - tail_at;
  		    t2 = ((tail_bt - bv) + (tail_at - (t1 - bv)));
  
  		    s2 += t1;
  
  		    /* Renormalize (s1, s2)  to  (t1, s2) */
  		    t1 = s1 + s2;
  		    s2 = s2 - (t1 - s1);
  
  		    t2 += s2;
  
  		    /* Renormalize (t1, t2)  */
  		    head_ct = t1 + t2;
  		    tail_ct = t2 - (head_ct - t1);
  		  }
  		  head_temp1[1] = head_ct;
  		  tail_temp1[1] = tail_ct;
  		}
  		ix -= inc_intx;
  	      }			/* for j<n */
  
  	      /*
  	       * if the diagonal entry is not equal to one, then divide Xj by
  	       * the entry
  	       */
  	      if (diag == blas_non_unit_diag) {
  		T_element = T_i[j * incT + j * ldt * incT];
  
  
  		{
  		  double	  head_a, tail_a;
  		  double	  head_b, tail_b;
  		  head_a = head_temp1[0];
  		  tail_a = tail_temp1[0];
  		  {
  		    /*
  		     * Compute double-double = double-double / double, using
  		     * a Newton iteration scheme.
  		     */
  		    double	    b1  , b2, con, e, t1, t2, t11, t21, t12,
  		    		    t22;
  
  		    /* Compute a DP approximation to the quotient. */
  		    t1 = head_a / T_element;
  
  		    /*
  		     * Split t1 and b into two parts with at most 26 bits
  		     * each, using the Dekker-Veltkamp method.
  		     */
  		    con = t1 * split;
  		    t11 = con - (con - t1);
  		    t21 = t1 - t11;
  		    con = T_element * split;
  		    b1 = con - (con - T_element);
  		    b2 = T_element - b1;
  
  		    /* Compute t1 * b using Dekker method. */
  		    t12 = t1 * T_element;
  		    t22 = (((t11 * b1 - t12) + t11 * b2) + t21 * b1) + t21 * b2;
  
  		    /* Compute dda - (t12, t22) using Knuth trick. */
  		    t11 = head_a - t12;
  		    e = t11 - head_a;
  		    t21 = ((-t12 - e) + (head_a - (t11 - e))) + tail_a - t22;
  
  		    /* Compute high-order word of (t11, t21) and divide by b. */
  		    t2 = (t11 + t21) / T_element;
  
  		    /* The result is t1 + t2, after normalization. */
  		    head_b = t1 + t2;
  		    tail_b = t2 - (head_b - t1);
  		  }
  		  head_temp1[0] = head_b;
  		  tail_temp1[0] = tail_b;
  		  head_a = head_temp1[1];
  		  tail_a = tail_temp1[1];
  		  {
  		    /*
  		     * Compute double-double = double-double / double, using
  		     * a Newton iteration scheme.
  		     */
  		    double	    b1  , b2, con, e, t1, t2, t11, t21, t12,
  		    		    t22;
  
  		    /* Compute a DP approximation to the quotient. */
  		    t1 = head_a / T_element;
  
  		    /*
  		     * Split t1 and b into two parts with at most 26 bits
  		     * each, using the Dekker-Veltkamp method.
  		     */
  		    con = t1 * split;
  		    t11 = con - (con - t1);
  		    t21 = t1 - t11;
  		    con = T_element * split;
  		    b1 = con - (con - T_element);
  		    b2 = T_element - b1;
  
  		    /* Compute t1 * b using Dekker method. */
  		    t12 = t1 * T_element;
  		    t22 = (((t11 * b1 - t12) + t11 * b2) + t21 * b1) + t21 * b2;
  
  		    /* Compute dda - (t12, t22) using Knuth trick. */
  		    t11 = head_a - t12;
  		    e = t11 - head_a;
  		    t21 = ((-t12 - e) + (head_a - (t11 - e))) + tail_a - t22;
  
  		    /* Compute high-order word of (t11, t21) and divide by b. */
  		    t2 = (t11 + t21) / T_element;
  
  		    /* The result is t1 + t2, after normalization. */
  		    head_b = t1 + t2;
  		    tail_b = t2 - (head_b - t1);
  		  }
  		  head_temp1[1] = head_b;
  		  tail_temp1[1] = tail_b;
  		}
  
  	      }			/* if (diag == blas_non_unit_diag) */
  	      head_intx[jx] = head_temp1[0];
  	      tail_intx[jx] = tail_temp1[0];
  	      head_intx[1 + jx] = head_temp1[1];
  	      tail_intx[1 + jx] = tail_temp1[1];
  
  	      jx -= inc_intx;
  	    }			/* for j>=0 */
  	  } else if ((order == blas_rowmajor &&
  		      trans != blas_no_trans && uplo == blas_upper) ||
  		     (order == blas_colmajor &&
  		      trans == blas_no_trans && uplo == blas_lower)) {
  
  	    jx = 0;
  	    for (j = 0; j < n; j++) {
  
  	      /*
  	       * compute Xj = Xj - SUM Aij(or Aji) * Xi i=j+1 to n-1
  	       */
  	      head_temp3[0] = head_intx[jx];
  	      head_temp3[1] = head_intx[1 + jx];
  	      tail_temp3[0] = tail_intx[jx];
  	      tail_temp3[1] = tail_intx[1 + jx];
  	      /* multiply by alpha */
  	      {
  		double		cd      [2];
  		cd[0] = (double)alpha_i[0];
  		cd[1] = (double)alpha_i[1];
  		{
  		  /* Compute complex-extra = complex-extra * complex-double. */
  		  double	  head_a0, tail_a0;
  		  double	  head_a1, tail_a1;
  		  double	  head_t1, tail_t1;
  		  double	  head_t2, tail_t2;
  		  head_a0 = head_temp3[0];
  		  tail_a0 = tail_temp3[0];
  		  head_a1 = head_temp3[1];
  		  tail_a1 = tail_temp3[1];
  		  /* real part */
  		  {
  		    /* Compute double-double = double-double * double. */
  		    double	    a11 , a21, b1, b2, c11, c21, c2, con, t1,
  		    		    t2;
  
  		    con = head_a0 * split;
  		    a11 = con - head_a0;
  		    a11 = con - a11;
  		    a21 = head_a0 - a11;
  		    con = cd[0] * split;
  		    b1 = con - cd[0];
  		    b1 = con - b1;
  		    b2 = cd[0] - b1;
  
  		    c11 = head_a0 * cd[0];
  		    c21 = (((a11 * b1 - c11) + a11 * b2) + a21 * b1) + a21 * b2;
  
  		    c2 = tail_a0 * cd[0];
  		    t1 = c11 + c2;
  		    t2 = (c2 - (t1 - c11)) + c21;
  
  		    head_t1 = t1 + t2;
  		    tail_t1 = t2 - (head_t1 - t1);
  		  }
  		  {
  		    /* Compute double-double = double-double * double. */
  		    double	    a11 , a21, b1, b2, c11, c21, c2, con, t1,
  		    		    t2;
  
  		    con = head_a1 * split;
  		    a11 = con - head_a1;
  		    a11 = con - a11;
  		    a21 = head_a1 - a11;
  		    con = cd[1] * split;
  		    b1 = con - cd[1];
  		    b1 = con - b1;
  		    b2 = cd[1] - b1;
  
  		    c11 = head_a1 * cd[1];
  		    c21 = (((a11 * b1 - c11) + a11 * b2) + a21 * b1) + a21 * b2;
  
  		    c2 = tail_a1 * cd[1];
  		    t1 = c11 + c2;
  		    t2 = (c2 - (t1 - c11)) + c21;
  
  		    head_t2 = t1 + t2;
  		    tail_t2 = t2 - (head_t2 - t1);
  		  }
  		  head_t2 = -head_t2;
  		  tail_t2 = -tail_t2;
  		  {
  		    /* Compute double-double = double-double + double-double. */
  		    double	    bv;
  		    double	    s1  , s2, t1, t2;
  
  		    /* Add two hi words. */
  		    s1 = head_t1 + head_t2;
  		    bv = s1 - head_t1;
  		    s2 = ((head_t2 - bv) + (head_t1 - (s1 - bv)));
  
  		    /* Add two lo words. */
  		    t1 = tail_t1 + tail_t2;
  		    bv = t1 - tail_t1;
  		    t2 = ((tail_t2 - bv) + (tail_t1 - (t1 - bv)));
  
  		    s2 += t1;
  
  		    /* Renormalize (s1, s2)  to  (t1, s2) */
  		    t1 = s1 + s2;
  		    s2 = s2 - (t1 - s1);
  
  		    t2 += s2;
  
  		    /* Renormalize (t1, t2)  */
  		    head_t1 = t1 + t2;
  		    tail_t1 = t2 - (head_t1 - t1);
  		  }
  		  head_temp1[0] = head_t1;
  		  tail_temp1[0] = tail_t1;
  		  /* imaginary part */
  		  {
  		    /* Compute double-double = double-double * double. */
  		    double	    a11 , a21, b1, b2, c11, c21, c2, con, t1,
  		    		    t2;
  
  		    con = head_a1 * split;
  		    a11 = con - head_a1;
  		    a11 = con - a11;
  		    a21 = head_a1 - a11;
  		    con = cd[0] * split;
  		    b1 = con - cd[0];
  		    b1 = con - b1;
  		    b2 = cd[0] - b1;
  
  		    c11 = head_a1 * cd[0];
  		    c21 = (((a11 * b1 - c11) + a11 * b2) + a21 * b1) + a21 * b2;
  
  		    c2 = tail_a1 * cd[0];
  		    t1 = c11 + c2;
  		    t2 = (c2 - (t1 - c11)) + c21;
  
  		    head_t1 = t1 + t2;
  		    tail_t1 = t2 - (head_t1 - t1);
  		  }
  		  {
  		    /* Compute double-double = double-double * double. */
  		    double	    a11 , a21, b1, b2, c11, c21, c2, con, t1,
  		    		    t2;
  
  		    con = head_a0 * split;
  		    a11 = con - head_a0;
  		    a11 = con - a11;
  		    a21 = head_a0 - a11;
  		    con = cd[1] * split;
  		    b1 = con - cd[1];
  		    b1 = con - b1;
  		    b2 = cd[1] - b1;
  
  		    c11 = head_a0 * cd[1];
  		    c21 = (((a11 * b1 - c11) + a11 * b2) + a21 * b1) + a21 * b2;
  
  		    c2 = tail_a0 * cd[1];
  		    t1 = c11 + c2;
  		    t2 = (c2 - (t1 - c11)) + c21;
  
  		    head_t2 = t1 + t2;
  		    tail_t2 = t2 - (head_t2 - t1);
  		  }
  		  {
  		    /* Compute double-double = double-double + double-double. */
  		    double	    bv;
  		    double	    s1  , s2, t1, t2;
  
  		    /* Add two hi words. */
  		    s1 = head_t1 + head_t2;
  		    bv = s1 - head_t1;
  		    s2 = ((head_t2 - bv) + (head_t1 - (s1 - bv)));
  
  		    /* Add two lo words. */
  		    t1 = tail_t1 + tail_t2;
  		    bv = t1 - tail_t1;
  		    t2 = ((tail_t2 - bv) + (tail_t1 - (t1 - bv)));
  
  		    s2 += t1;
  
  		    /* Renormalize (s1, s2)  to  (t1, s2) */
  		    t1 = s1 + s2;
  		    s2 = s2 - (t1 - s1);
  
  		    t2 += s2;
  
  		    /* Renormalize (t1, t2)  */
  		    head_t1 = t1 + t2;
  		    tail_t1 = t2 - (head_t1 - t1);
  		  }
  		  head_temp1[1] = head_t1;
  		  tail_temp1[1] = tail_t1;
  		}
  
  	      }
  
  	      ix = 0;
  	      for (i = 0; i < j; i++) {
  		T_element = T_i[j * incT + i * ldt * incT];
  
  		head_temp3[0] = head_intx[ix];
  		head_temp3[1] = head_intx[1 + ix];
  		tail_temp3[0] = tail_intx[ix];
  		tail_temp3[1] = tail_intx[1 + ix];
  		{
  		  double	  dt = (double)T_element;
  		  {
  		    /* Compute complex-extra = complex-extra * real. */
  		    double	    head_a0, tail_a0;
  		    double	    head_a1, tail_a1;
  		    double	    head_t, tail_t;
  		    head_a0 = head_temp3[0];
  		    tail_a0 = tail_temp3[0];
  		    head_a1 = head_temp3[1];
  		    tail_a1 = tail_temp3[1];
  		    {
  		      /* Compute double-double = double-double * double. */
  		      double	      a11, a21, b1, b2, c11, c21, c2, con,
  		      		      t1     , t2;
  
  		      con = head_a0 * split;
  		      a11 = con - head_a0;
  		      a11 = con - a11;
  		      a21 = head_a0 - a11;
  		      con = dt * split;
  		      b1 = con - dt;
  		      b1 = con - b1;
  		      b2 = dt - b1;
  
  		      c11 = head_a0 * dt;
  		      c21 = (((a11 * b1 - c11) + a11 * b2) + a21 * b1) + a21 * b2;
  
  		      c2 = tail_a0 * dt;
  		      t1 = c11 + c2;
  		      t2 = (c2 - (t1 - c11)) + c21;
  
  		      head_t = t1 + t2;
  		      tail_t = t2 - (head_t - t1);
  		    }
  		    head_temp2[0] = head_t;
  		    tail_temp2[0] = tail_t;
  		    {
  		      /* Compute double-double = double-double * double. */
  		      double	      a11, a21, b1, b2, c11, c21, c2, con,
  		      		      t1     , t2;
  
  		      con = head_a1 * split;
  		      a11 = con - head_a1;
  		      a11 = con - a11;
  		      a21 = head_a1 - a11;
  		      con = dt * split;
  		      b1 = con - dt;
  		      b1 = con - b1;
  		      b2 = dt - b1;
  
  		      c11 = head_a1 * dt;
  		      c21 = (((a11 * b1 - c11) + a11 * b2) + a21 * b1) + a21 * b2;
  
  		      c2 = tail_a1 * dt;
  		      t1 = c11 + c2;
  		      t2 = (c2 - (t1 - c11)) + c21;
  
  		      head_t = t1 + t2;
  		      tail_t = t2 - (head_t - t1);
  		    }
  		    head_temp2[1] = head_t;
  		    tail_temp2[1] = tail_t;
  		  }
  
  		}
  		{
  		  double	  head_at, tail_at;
  		  double	  head_bt, tail_bt;
  		  double	  head_ct, tail_ct;
  
  		  /* Real part */
  		  head_at = head_temp1[0];
  		  tail_at = tail_temp1[0];
  		  head_bt = -head_temp2[0];
  		  tail_bt = -tail_temp2[0];
  		  {
  		    /* Compute double-double = double-double + double-double. */
  		    double	    bv;
  		    double	    s1  , s2, t1, t2;
  
  		    /* Add two hi words. */
  		    s1 = head_at + head_bt;
  		    bv = s1 - head_at;
  		    s2 = ((head_bt - bv) + (head_at - (s1 - bv)));
  
  		    /* Add two lo words. */
  		    t1 = tail_at + tail_bt;
  		    bv = t1 - tail_at;
  		    t2 = ((tail_bt - bv) + (tail_at - (t1 - bv)));
  
  		    s2 += t1;
  
  		    /* Renormalize (s1, s2)  to  (t1, s2) */
  		    t1 = s1 + s2;
  		    s2 = s2 - (t1 - s1);
  
  		    t2 += s2;
  
  		    /* Renormalize (t1, t2)  */
  		    head_ct = t1 + t2;
  		    tail_ct = t2 - (head_ct - t1);
  		  }
  		  head_temp1[0] = head_ct;
  		  tail_temp1[0] = tail_ct;
  		  /* Imaginary part */
  		  head_at = head_temp1[1];
  		  tail_at = tail_temp1[1];
  		  head_bt = -head_temp2[1];
  		  tail_bt = -tail_temp2[1];
  		  {
  		    /* Compute double-double = double-double + double-double. */
  		    double	    bv;
  		    double	    s1  , s2, t1, t2;
  
  		    /* Add two hi words. */
  		    s1 = head_at + head_bt;
  		    bv = s1 - head_at;
  		    s2 = ((head_bt - bv) + (head_at - (s1 - bv)));
  
  		    /* Add two lo words. */
  		    t1 = tail_at + tail_bt;
  		    bv = t1 - tail_at;
  		    t2 = ((tail_bt - bv) + (tail_at - (t1 - bv)));
  
  		    s2 += t1;
  
  		    /* Renormalize (s1, s2)  to  (t1, s2) */
  		    t1 = s1 + s2;
  		    s2 = s2 - (t1 - s1);
  
  		    t2 += s2;
  
  		    /* Renormalize (t1, t2)  */
  		    head_ct = t1 + t2;
  		    tail_ct = t2 - (head_ct - t1);
  		  }
  		  head_temp1[1] = head_ct;
  		  tail_temp1[1] = tail_ct;
  		}
  		ix += inc_intx;
  	      }			/* for i<j */
  
  	      /*
  	       * if the diagonal entry is not equal to one, then divide Xj by
  	       * the entry
  	       */
  	      if (diag == blas_non_unit_diag) {
  		T_element = T_i[j * incT + j * ldt * incT];
  
  
  		{
  		  double	  head_a, tail_a;
  		  double	  head_b, tail_b;
  		  head_a = head_temp1[0];
  		  tail_a = tail_temp1[0];
  		  {
  		    /*
  		     * Compute double-double = double-double / double, using
  		     * a Newton iteration scheme.
  		     */
  		    double	    b1  , b2, con, e, t1, t2, t11, t21, t12,
  		    		    t22;
  
  		    /* Compute a DP approximation to the quotient. */
  		    t1 = head_a / T_element;
  
  		    /*
  		     * Split t1 and b into two parts with at most 26 bits
  		     * each, using the Dekker-Veltkamp method.
  		     */
  		    con = t1 * split;
  		    t11 = con - (con - t1);
  		    t21 = t1 - t11;
  		    con = T_element * split;
  		    b1 = con - (con - T_element);
  		    b2 = T_element - b1;
  
  		    /* Compute t1 * b using Dekker method. */
  		    t12 = t1 * T_element;
  		    t22 = (((t11 * b1 - t12) + t11 * b2) + t21 * b1) + t21 * b2;
  
  		    /* Compute dda - (t12, t22) using Knuth trick. */
  		    t11 = head_a - t12;
  		    e = t11 - head_a;
  		    t21 = ((-t12 - e) + (head_a - (t11 - e))) + tail_a - t22;
  
  		    /* Compute high-order word of (t11, t21) and divide by b. */
  		    t2 = (t11 + t21) / T_element;
  
  		    /* The result is t1 + t2, after normalization. */
  		    head_b = t1 + t2;
  		    tail_b = t2 - (head_b - t1);
  		  }
  		  head_temp1[0] = head_b;
  		  tail_temp1[0] = tail_b;
  		  head_a = head_temp1[1];
  		  tail_a = tail_temp1[1];
  		  {
  		    /*
  		     * Compute double-double = double-double / double, using
  		     * a Newton iteration scheme.
  		     */
  		    double	    b1  , b2, con, e, t1, t2, t11, t21, t12,
  		    		    t22;
  
  		    /* Compute a DP approximation to the quotient. */
  		    t1 = head_a / T_element;
  
  		    /*
  		     * Split t1 and b into two parts with at most 26 bits
  		     * each, using the Dekker-Veltkamp method.
  		     */
  		    con = t1 * split;
  		    t11 = con - (con - t1);
  		    t21 = t1 - t11;
  		    con = T_element * split;
  		    b1 = con - (con - T_element);
  		    b2 = T_element - b1;
  
  		    /* Compute t1 * b using Dekker method. */
  		    t12 = t1 * T_element;
  		    t22 = (((t11 * b1 - t12) + t11 * b2) + t21 * b1) + t21 * b2;
  
  		    /* Compute dda - (t12, t22) using Knuth trick. */
  		    t11 = head_a - t12;
  		    e = t11 - head_a;
  		    t21 = ((-t12 - e) + (head_a - (t11 - e))) + tail_a - t22;
  
  		    /* Compute high-order word of (t11, t21) and divide by b. */
  		    t2 = (t11 + t21) / T_element;
  
  		    /* The result is t1 + t2, after normalization. */
  		    head_b = t1 + t2;
  		    tail_b = t2 - (head_b - t1);
  		  }
  		  head_temp1[1] = head_b;
  		  tail_temp1[1] = tail_b;
  		}
  
  	      }			/* if (diag == blas_non_unit_diag) */
  	      head_intx[jx] = head_temp1[0];
  	      tail_intx[jx] = tail_temp1[0];
  	      head_intx[1 + jx] = head_temp1[1];
  	      tail_intx[1 + jx] = tail_temp1[1];
  	      jx += inc_intx;
  	    }			/* for j<n */
  	  }
  	  /* copy the final results from intx to x */
  	  ix = start_x;
  	  jx = 0;
  	  for (i = 0; i < n; i++) {
  	    head_temp1[0] = head_intx[jx];
  	    head_temp1[1] = head_intx[1 + jx];
  	    tail_temp1[0] = tail_intx[jx];
  	    tail_temp1[1] = tail_intx[1 + jx];
  	    x_i[ix] = head_temp1[0];
  	    x_i[ix + 1] = head_temp1[1];
  	    ix += incx;
  	    jx += inc_intx;
  	  }
  
  	  blas_free(head_intx);
  	  blas_free(tail_intx);
  	}
        }
        FPU_FIX_STOP;
      }
      break;
    }
  }
  