real THBlas_(dot)(int64_t n, real *x, int64_t incx, real *y, int64_t incy)
{
  if(n == 1)
  {
    incx = 1;
    incy = 1;
  }

#if defined(USE_BLAS) && (defined(TH_REAL_IS_DOUBLE) || defined(TH_REAL_IS_FLOAT))
  if( (n <= INT_MAX) && (incx <= INT_MAX) && (incy <= INT_MAX) )
  {
    int i_n = (int)n;
    int i_incx = (int)incx;
    int i_incy = (int)incy;

#if defined(TH_REAL_IS_DOUBLE)
    return (real) ddot_(&i_n, x, &i_incx, y, &i_incy);
#else
    return (real) sdot_(&i_n, x, &i_incx, y, &i_incy);
#endif
  }
#endif
  {
    int64_t i;
    real sum = 0;
    for(i = 0; i < n; i++)
    sum += x[i*incx]*y[i*incy];
    return sum;
  }
}

int main()
{
  const unsigned int N = 100;
  const unsigned int ONE = 1;
  float a[N];
  float b[N];

  for (unsigned int i = 0; i < N; ++i) {
      a[i] = float(i);
      b[i] = a[i] / 2.;
  }

  float exact_result = 0.;
  for (unsigned int i = 0; i < N; ++i)
    exact_result += a[i] * b[i];

  float result = sdot_(&N, a, &ONE, b, &ONE);
  if (std::abs(result - exact_result) < 1e-6 * fabs(exact_result)) {
      std::cout << "SUCCESS" << std::endl;
      return 0;
  } else {
      std::cout << "FAILED" << std::endl;
      return 1;
  }
}

doublereal 
f2c_sdot(integer* N, 
         real* X, integer* incX, 
         real* Y, integer* incY)
{
    return sdot_(N, X, incX, Y, incY);
}

/* This routine computes the pure pseudoinverse from the svd returned
by svdcmp (U*S*V^T) so Agi = V*S^-1*U^T.  The algorithm used is a
little overly tricky using an internally allocated work vector to 
make the routine nondestructive to the input matrices.  This computes
the "pure" pseudoinverse by setting the svd cutoff value based on
float epsilon (from float.h).

Note input U is mxn, s is an n vector, V is nxn, and the output
Agi is nxm.  

Function returns the number of singular values actually used to 
compute Agi.  

Author:  Gary L. Pavlis
*/
int pseudoinverse(float **U, float *s, float **V, int m, int n, float **Agi)
{
	int i,j, k;  /* counters*/
	float *work;  /* work space */
	float smax;
	float sinv;
	double sv_cutoff;
	int nsv_used;
#ifndef SUNPERF
	int one=1;
#endif


        if((work=(float *)calloc(n,sizeof(float))) == NULL)
                elog_die(1,"Pseudoinverse computation: cannot alloc work array of length %d\n",
                        n);

	/* first find the larges singular value, then just zero
	all those smaller than the cutoff determined as the ratio
	wrt to largest singular value */
	smax = 0.0;
	for(i=0;i<n;++i) 
		if(s[i] > smax) smax = s[i];
	sv_cutoff = (double)smax*FLT_EPSILON;
	/* This is a copy operation */
	for(i=0;i<m;++i)
		for(j=0;j<n;++j) Agi[j][i] = U[i][j];
	/* this works because of C storage order, but is strange.
	It is the multiply by S^-1 */
	for(j=0,nsv_used=0;j<n;++j)
	{
		if( (double)s[j] > sv_cutoff)
		{
			sinv = 1.0/s[j];
			++nsv_used;
		}
		else
			sinv = 0.0;
#ifdef SUNPERF
		sscal(m,sinv,Agi[j],1);
#else
		sscal_(&m,&sinv,Agi[j],&one);
#endif
	}
	/* multiply by V using a column work vector*/
	for(j=0;j<m;++j)
	{
		for(k=0;k<n;++k) work[k] = Agi[k][j];
		for(i=0;i<n;++i)
#ifdef SUNPERF
			Agi[i][j] = sdot(n,work,1,V[i],1);
#else
			Agi[i][j] = sdot_(&n,work,&one,V[i],&one);
#endif
	}
	free(work);
	return(nsv_used);
}

float STARPU_SDOT(const int n, const float *x, const int incx, const float *y, const int incy)
{
	float retVal = 0;

	/* GOTOBLAS will return a FLOATRET which is a double, not a float */
	retVal = (float)sdot_(&n, x, &incx, y, &incy);

	return retVal;
}

GURLS_EXPORT float dot(const gVec<float>& x, const gVec<float>& y)
{
    if ( x.getSize() != y.getSize() )
        throw gException(gurls::Exception_Inconsistent_Size);

    int n = x.getSize();
    int incr = 1;

    return sdot_(&n, const_cast<float*>(x.getData()), &incr, const_cast<float*>(y.getData()), &incr);
}

/* Subroutine */
int slapll_(integer *n, real *x, integer *incx, real *y, integer *incy, real *ssmin)
{
    /* System generated locals */
    integer i__1;
    /* Local variables */
    real c__, a11, a12, a22, tau;
    extern real sdot_(integer *, real *, integer *, real *, integer *);
    extern /* Subroutine */
    int slas2_(real *, real *, real *, real *, real *) ;
    real ssmax;
    extern /* Subroutine */
    int saxpy_(integer *, real *, real *, integer *, real *, integer *), slarfg_(integer *, real *, real *, integer *, real *);
    /* -- LAPACK auxiliary routine (version 3.4.2) -- */
    /* -- LAPACK is a software package provided by Univ. of Tennessee, -- */
    /* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- */
    /* September 2012 */
    /* .. Scalar Arguments .. */
    /* .. */
    /* .. Array Arguments .. */
    /* .. */
    /* ===================================================================== */
    /* .. Parameters .. */
    /* .. */
    /* .. Local Scalars .. */
    /* .. */
    /* .. External Functions .. */
    /* .. */
    /* .. External Subroutines .. */
    /* .. */
    /* .. Executable Statements .. */
    /* Quick return if possible */
    /* Parameter adjustments */
    --y;
    --x;
    /* Function Body */
    if (*n <= 1)
    {
        *ssmin = 0.f;
        return 0;
    }
    /* Compute the QR factorization of the N-by-2 matrix ( X Y ) */
    slarfg_(n, &x[1], &x[*incx + 1], incx, &tau);
    a11 = x[1];
    x[1] = 1.f;
    c__ = -tau * sdot_(n, &x[1], incx, &y[1], incy);
    saxpy_(n, &c__, &x[1], incx, &y[1], incy);
    i__1 = *n - 1;
    slarfg_(&i__1, &y[*incy + 1], &y[(*incy << 1) + 1], incy, &tau);
    a12 = y[1];
    a22 = y[*incy + 1];
    /* Compute the SVD of 2-by-2 Upper triangular matrix. */
    slas2_(&a11, &a12, &a22, ssmin, &ssmax);
    return 0;
    /* End of SLAPLL */
}

/* Subroutine */ int sspgst_(integer *itype, char *uplo, integer *n, real *ap, 
	 real *bp, integer *info)
{
    /* System generated locals */
    integer i__1, i__2;
    real r__1;

    /* Local variables */
    integer j, k, j1, k1, jj, kk;
    real ct, ajj;
    integer j1j1;
    real akk;
    integer k1k1;
    real bjj, bkk;
    extern doublereal sdot_(integer *, real *, integer *, real *, integer *);
    extern /* Subroutine */ int sspr2_(char *, integer *, real *, real *, 
	    integer *, real *, integer *, real *);
    extern logical lsame_(char *, char *);
    extern /* Subroutine */ int sscal_(integer *, real *, real *, integer *);
    logical upper;
    extern /* Subroutine */ int saxpy_(integer *, real *, real *, integer *, 
	    real *, integer *), sspmv_(char *, integer *, real *, real *, 
	    real *, integer *, real *, real *, integer *), stpmv_(
	    char *, char *, char *, integer *, real *, real *, integer *), stpsv_(char *, char *, char *, integer *, 
	     real *, real *, integer *), xerbla_(char 
	    *, integer *);


/*  -- LAPACK routine (version 3.2) -- */
/*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
/*     November 2006 */

/*     .. Scalar Arguments .. */
/*     .. */
/*     .. Array Arguments .. */
/*     .. */

/*  Purpose */
/*  ======= */

/*  SSPGST reduces a real symmetric-definite generalized eigenproblem */
/*  to standard form, using packed storage. */

/*  If ITYPE = 1, the problem is A*x = lambda*B*x, */
/*  and A is overwritten by inv(U**T)*A*inv(U) or inv(L)*A*inv(L**T) */

/*  If ITYPE = 2 or 3, the problem is A*B*x = lambda*x or */
/*  B*A*x = lambda*x, and A is overwritten by U*A*U**T or L**T*A*L. */

/*  B must have been previously factorized as U**T*U or L*L**T by SPPTRF. */

/*  Arguments */
/*  ========= */

/*  ITYPE   (input) INTEGER */
/*          = 1: compute inv(U**T)*A*inv(U) or inv(L)*A*inv(L**T); */
/*          = 2 or 3: compute U*A*U**T or L**T*A*L. */

/*  UPLO    (input) CHARACTER*1 */
/*          = 'U':  Upper triangle of A is stored and B is factored as */
/*                  U**T*U; */
/*          = 'L':  Lower triangle of A is stored and B is factored as */
/*                  L*L**T. */

/*  N       (input) INTEGER */
/*          The order of the matrices A and B.  N >= 0. */

/*  AP      (input/output) REAL array, dimension (N*(N+1)/2) */
/*          On entry, the upper or lower triangle of the symmetric matrix */
/*          A, packed columnwise in a linear array.  The j-th column of A */
/*          is stored in the array AP as follows: */
/*          if UPLO = 'U', AP(i + (j-1)*j/2) = A(i,j) for 1<=i<=j; */
/*          if UPLO = 'L', AP(i + (j-1)*(2n-j)/2) = A(i,j) for j<=i<=n. */

/*          On exit, if INFO = 0, the transformed matrix, stored in the */
/*          same format as A. */

/*  BP      (input) REAL array, dimension (N*(N+1)/2) */
/*          The triangular factor from the Cholesky factorization of B, */
/*          stored in the same format as A, as returned by SPPTRF. */

/*  INFO    (output) INTEGER */
/*          = 0:  successful exit */
/*          < 0:  if INFO = -i, the i-th argument had an illegal value */

/*  ===================================================================== */

/*     .. Parameters .. */
/*     .. */
/*     .. Local Scalars .. */
/*     .. */
/*     .. External Subroutines .. */
/*     .. */
/*     .. External Functions .. */
/*     .. */
/*     .. Executable Statements .. */

/*     Test the input parameters. */

    /* Parameter adjustments */
//.........这里部分代码省略.........

/* Subroutine */ int slarfy_(char *uplo, integer *n, real *v, integer *incv,
                             real *tau, real *c__, integer *ldc, real *work)
{
    /* System generated locals */
    integer c_dim1, c_offset;
    real r__1;

    /* Local variables */
    extern doublereal sdot_(integer *, real *, integer *, real *, integer *);
    extern /* Subroutine */ int ssyr2_(char *, integer *, real *, real *,
                                       integer *, real *, integer *, real *, integer *);
    real alpha;
    extern /* Subroutine */ int saxpy_(integer *, real *, real *, integer *,
                                       real *, integer *), ssymv_(char *, integer *, real *, real *,
                                               integer *, real *, integer *, real *, real *, integer *);


    /*  -- LAPACK auxiliary test routine (version 3.1) -- */
    /*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
    /*     November 2006 */

    /*     .. Scalar Arguments .. */
    /*     .. */
    /*     .. Array Arguments .. */
    /*     .. */

    /*  Purpose */
    /*  ======= */

    /*  SLARFY applies an elementary reflector, or Householder matrix, H, */
    /*  to an n x n symmetric matrix C, from both the left and the right. */

    /*  H is represented in the form */

    /*     H = I - tau * v * v' */

    /*  where  tau  is a scalar and  v  is a vector. */

    /*  If  tau  is  zero, then  H  is taken to be the unit matrix. */

    /*  Arguments */
    /*  ========= */

    /*  UPLO    (input) CHARACTER*1 */
    /*          Specifies whether the upper or lower triangular part of the */
    /*          symmetric matrix C is stored. */
    /*          = 'U':  Upper triangle */
    /*          = 'L':  Lower triangle */

    /*  N       (input) INTEGER */
    /*          The number of rows and columns of the matrix C.  N >= 0. */

    /*  V       (input) REAL array, dimension */
    /*                  (1 + (N-1)*abs(INCV)) */
    /*          The vector v as described above. */

    /*  INCV    (input) INTEGER */
    /*          The increment between successive elements of v.  INCV must */
    /*          not be zero. */

    /*  TAU     (input) REAL */
    /*          The value tau as described above. */

    /*  C       (input/output) REAL array, dimension (LDC, N) */
    /*          On entry, the matrix C. */
    /*          On exit, C is overwritten by H * C * H'. */

    /*  LDC     (input) INTEGER */
    /*          The leading dimension of the array C.  LDC >= max( 1, N ). */

    /*  WORK    (workspace) REAL array, dimension (N) */

    /*  ===================================================================== */

    /*     .. Parameters .. */
    /*     .. */
    /*     .. Local Scalars .. */
    /*     .. */
    /*     .. External Subroutines .. */
    /*     .. */
    /*     .. External Functions .. */
    /*     .. */
    /*     .. Executable Statements .. */

    /* Parameter adjustments */
    --v;
    c_dim1 = *ldc;
    c_offset = 1 + c_dim1;
    c__ -= c_offset;
    --work;

    /* Function Body */
    if (*tau == 0.f) {
        return 0;
    }

    /*     Form  w:= C * v */

    ssymv_(uplo, n, &c_b2, &c__[c_offset], ldc, &v[1], incv, &c_b3, &work[1],
           &c__1);
//.........这里部分代码省略.........

/* Subroutine */ int spptri_(char *uplo, integer *n, real *ap, integer *info)
{
    /* System generated locals */
    integer i__1, i__2;

    /* Local variables */
    integer j, jc, jj;
    real ajj;
    integer jjn;
    extern doublereal sdot_(integer *, real *, integer *, real *, integer *);
    extern /* Subroutine */ int sspr_(char *, integer *, real *, real *, 
	    integer *, real *);
    extern logical lsame_(char *, char *);
    extern /* Subroutine */ int sscal_(integer *, real *, real *, integer *);
    logical upper;
    extern /* Subroutine */ int stpmv_(char *, char *, char *, integer *, 
	    real *, real *, integer *), xerbla_(char *
, integer *), stptri_(char *, char *, integer *, real *, 
	    integer *);


/*  -- LAPACK routine (version 3.2) -- */
/*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
/*     November 2006 */

/*     .. Scalar Arguments .. */
/*     .. */
/*     .. Array Arguments .. */
/*     .. */

/*  Purpose */
/*  ======= */

/*  SPPTRI computes the inverse of a real symmetric positive definite */
/*  matrix A using the Cholesky factorization A = U**T*U or A = L*L**T */
/*  computed by SPPTRF. */

/*  Arguments */
/*  ========= */

/*  UPLO    (input) CHARACTER*1 */
/*          = 'U':  Upper triangular factor is stored in AP; */
/*          = 'L':  Lower triangular factor is stored in AP. */

/*  N       (input) INTEGER */
/*          The order of the matrix A.  N >= 0. */

/*  AP      (input/output) REAL array, dimension (N*(N+1)/2) */
/*          On entry, the triangular factor U or L from the Cholesky */
/*          factorization A = U**T*U or A = L*L**T, packed columnwise as */
/*          a linear array.  The j-th column of U or L is stored in the */
/*          array AP as follows: */
/*          if UPLO = 'U', AP(i + (j-1)*j/2) = U(i,j) for 1<=i<=j; */
/*          if UPLO = 'L', AP(i + (j-1)*(2n-j)/2) = L(i,j) for j<=i<=n. */

/*          On exit, the upper or lower triangle of the (symmetric) */
/*          inverse of A, overwriting the input factor U or L. */

/*  INFO    (output) INTEGER */
/*          = 0:  successful exit */
/*          < 0:  if INFO = -i, the i-th argument had an illegal value */
/*          > 0:  if INFO = i, the (i,i) element of the factor U or L is */
/*                zero, and the inverse could not be computed. */

/*  ===================================================================== */

/*     .. Parameters .. */
/*     .. */
/*     .. Local Scalars .. */
/*     .. */
/*     .. External Functions .. */
/*     .. */
/*     .. External Subroutines .. */
/*     .. */
/*     .. Executable Statements .. */

/*     Test the input parameters. */

    /* Parameter adjustments */
    --ap;

    /* Function Body */
    *info = 0;
    upper = lsame_(uplo, "U");
    if (! upper && ! lsame_(uplo, "L")) {
	*info = -1;
    } else if (*n < 0) {
	*info = -2;
    }
    if (*info != 0) {
	i__1 = -(*info);
	xerbla_("SPPTRI", &i__1);
	return 0;
    }

/*     Quick return if possible */

    if (*n == 0) {
	return 0;
    }
//.........这里部分代码省略.........

//.........这里部分代码省略.........
    The contents of A on exit are illustrated by the following examples   
    with n = 5:   

    if UPLO = 'U':                       if UPLO = 'L':   

      (  d   e   v2  v3  v4 )              (  d                  )   
      (      d   e   v3  v4 )              (  e   d              )   
      (          d   e   v4 )              (  v1  e   d          )   
      (              d   e  )              (  v1  v2  e   d      )   
      (                  d  )              (  v1  v2  v3  e   d  )   

    where d and e denote diagonal and off-diagonal elements of T, and vi 
  
    denotes an element of the vector defining H(i).   

    ===================================================================== 
  


       Test the input parameters   

    
   Parameter adjustments   
       Function Body */
    /* Table of constant values */
    static integer c__1 = 1;
    static real c_b8 = 0.f;
    static real c_b14 = -1.f;
    
    /* System generated locals */
    integer a_dim1, a_offset, i__1, i__2, i__3;
    /* Local variables */
    static real taui;
    extern doublereal sdot_(integer *, real *, integer *, real *, integer *);
    static integer i;
    extern /* Subroutine */ int ssyr2_(char *, integer *, real *, real *, 
	    integer *, real *, integer *, real *, integer *);
    static real alpha;
    extern logical lsame_(char *, char *);
    static logical upper;
    extern /* Subroutine */ int saxpy_(integer *, real *, real *, integer *, 
	    real *, integer *), ssymv_(char *, integer *, real *, real *, 
	    integer *, real *, integer *, real *, real *, integer *), 
	    xerbla_(char *, integer *), slarfg_(integer *, real *, 
	    real *, integer *, real *);



#define D(I) d[(I)-1]
#define E(I) e[(I)-1]
#define TAU(I) tau[(I)-1]

#define A(I,J) a[(I)-1 + ((J)-1)* ( *lda)]

    *info = 0;
    upper = lsame_(uplo, "U");
    if (! upper && ! lsame_(uplo, "L")) {
	*info = -1;
    } else if (*n < 0) {
	*info = -2;
    } else if (*lda < max(1,*n)) {
	*info = -4;
    }
    if (*info != 0) {
	i__1 = -(*info);
	xerbla_("SSYTD2", &i__1);

//.........这里部分代码省略.........

    This routine is a further developed implementation of algorithm   
    BSOLVE in [1] using complete pivoting in the LU factorization.   

    [1] Bo Kagstrom and Lars Westin,   
        Generalized Schur Methods with Condition Estimators for   
        Solving the Generalized Sylvester Equation, IEEE Transactions   
        on Automatic Control, Vol. 34, No. 7, July 1989, pp 745-751.   

    [2] Peter Poromaa,   
        On Efficient and Robust Estimators for the Separation   
        between two Regular Matrix Pairs with Applications in   
        Condition Estimation. Report IMINF-95.05, Departement of   
        Computing Science, Umea University, S-901 87 Umea, Sweden, 1995.   

    =====================================================================   


       Parameter adjustments */
    /* Table of constant values */
    static integer c__1 = 1;
    static integer c_n1 = -1;
    static real c_b23 = 1.f;
    static real c_b37 = -1.f;
    
    /* System generated locals */
    integer z_dim1, z_offset, i__1, i__2;
    real r__1;
    /* Builtin functions */
    double sqrt(doublereal);
    /* Local variables */
    static integer info;
    static real temp;
    extern doublereal sdot_(integer *, real *, integer *, real *, integer *);
    static real work[32];
    static integer i__, j, k;
    extern /* Subroutine */ int sscal_(integer *, real *, real *, integer *);
    static real pmone;
    extern doublereal sasum_(integer *, real *, integer *);
    static real sminu;
    static integer iwork[8];
    extern /* Subroutine */ int scopy_(integer *, real *, integer *, real *, 
	    integer *), saxpy_(integer *, real *, real *, integer *, real *, 
	    integer *);
    static real splus;
    extern /* Subroutine */ int sgesc2_(integer *, real *, integer *, real *, 
	    integer *, integer *, real *);
    static real bm, bp, xm[8], xp[8];
    extern /* Subroutine */ int sgecon_(char *, integer *, real *, integer *, 
	    real *, real *, real *, integer *, integer *), slassq_(
	    integer *, real *, integer *, real *, real *), slaswp_(integer *, 
	    real *, integer *, integer *, integer *, integer *, integer *);
#define z___ref(a_1,a_2) z__[(a_2)*z_dim1 + a_1]


    z_dim1 = *ldz;
    z_offset = 1 + z_dim1 * 1;
    z__ -= z_offset;
    --rhs;
    --ipiv;
    --jpiv;

    /* Function Body */
    if (*ijob != 2) {

/*        Apply permutations IPIV to RHS */

/* Subroutine */ int spst01_(char *uplo, integer *n, real *a, integer *lda, 
	real *afac, integer *ldafac, real *perm, integer *ldperm, integer *
	piv, real *rwork, real *resid, integer *rank)
{
    /* System generated locals */
    integer a_dim1, a_offset, afac_dim1, afac_offset, perm_dim1, perm_offset, 
	    i__1, i__2;

    /* Local variables */
    integer i__, j, k;
    real t, eps;
    extern doublereal sdot_(integer *, real *, integer *, real *, integer *);
    extern /* Subroutine */ int ssyr_(char *, integer *, real *, real *, 
	    integer *, real *, integer *);
    extern logical lsame_(char *, char *);
    extern /* Subroutine */ int sscal_(integer *, real *, real *, integer *);
    real anorm;
    extern /* Subroutine */ int strmv_(char *, char *, char *, integer *, 
	    real *, integer *, real *, integer *);
    extern doublereal slamch_(char *), slansy_(char *, char *, 
	    integer *, real *, integer *, real *);


/*  -- LAPACK test routine (version 3.1) -- */
/*     Craig Lucas, University of Manchester / NAG Ltd. */
/*     October, 2008 */

/*     .. Scalar Arguments .. */
/*     .. */
/*     .. Array Arguments .. */
/*     .. */

/*  Purpose */
/*  ======= */

/*  SPST01 reconstructs a symmetric positive semidefinite matrix A */
/*  from its L or U factors and the permutation matrix P and computes */
/*  the residual */
/*     norm( P*L*L'*P' - A ) / ( N * norm(A) * EPS ) or */
/*     norm( P*U'*U*P' - A ) / ( N * norm(A) * EPS ), */
/*  where EPS is the machine epsilon. */

/*  Arguments */
/*  ========== */

/*  UPLO    (input) CHARACTER*1 */
/*          Specifies whether the upper or lower triangular part of the */
/*          symmetric matrix A is stored: */
/*          = 'U':  Upper triangular */
/*          = 'L':  Lower triangular */

/*  N       (input) INTEGER */
/*          The number of rows and columns of the matrix A.  N >= 0. */

/*  A       (input) REAL array, dimension (LDA,N) */
/*          The original symmetric matrix A. */

/*  LDA     (input) INTEGER */
/*          The leading dimension of the array A.  LDA >= max(1,N) */

/*  AFAC    (input) REAL array, dimension (LDAFAC,N) */
/*          The factor L or U from the L*L' or U'*U */
/*          factorization of A. */

/*  LDAFAC  (input) INTEGER */
/*          The leading dimension of the array AFAC.  LDAFAC >= max(1,N). */

/*  PERM    (output) REAL array, dimension (LDPERM,N) */
/*          Overwritten with the reconstructed matrix, and then with the */
/*          difference P*L*L'*P' - A (or P*U'*U*P' - A) */

/*  LDPERM  (input) INTEGER */
/*          The leading dimension of the array PERM. */
/*          LDAPERM >= max(1,N). */

/*  PIV     (input) INTEGER array, dimension (N) */
/*          PIV is such that the nonzero entries are */
/*          P( PIV( K ), K ) = 1. */

/*  RWORK   (workspace) REAL array, dimension (N) */

/*  RESID   (output) REAL */
/*          If UPLO = 'L', norm(L*L' - A) / ( N * norm(A) * EPS ) */
/*          If UPLO = 'U', norm(U'*U - A) / ( N * norm(A) * EPS ) */

/*  ===================================================================== */

/*     .. Parameters .. */
/*     .. */
/*     .. Local Scalars .. */
/*     .. */
/*     .. External Functions .. */
/*     .. */
/*     .. External Subroutines .. */
/*     .. */
/*     .. Intrinsic Functions .. */
/*     .. */
/*     .. Executable Statements .. */

/*     Quick exit if N = 0. */
//.........这里部分代码省略.........

//.........这里部分代码省略.........
    if UPLO = 'U':                       if UPLO = 'L':   

      (  a   a   a   v4  v5 )              (  d                  )   
      (      a   a   v4  v5 )              (  1   d              )   
      (          a   1   v5 )              (  v1  1   a          )   
      (              d   1  )              (  v1  v2  a   a      )   
      (                  d  )              (  v1  v2  a   a   a  )   

    where d denotes a diagonal element of the reduced matrix, a denotes   
    an element of the original matrix that is unchanged, and vi denotes   
    an element of the vector defining H(i).   

    ===================================================================== 
  


       Quick return if possible   

    
   Parameter adjustments   
       Function Body */
    /* Table of constant values */
    static real c_b5 = -1.f;
    static real c_b6 = 1.f;
    static int c__1 = 1;
    static real c_b16 = 0.f;
    
    /* System generated locals */
/*  Unused variables commented out by MDG on 03-09-05
    int a_dim1, a_offset, w_dim1, w_offset;
*/
    int i__1, i__2, i__3;
    /* Local variables */
    extern doublereal sdot_(int *, real *, int *, real *, int *);
    static int i;
    static real alpha;
    extern logical lsame_(char *, char *);
    extern /* Subroutine */ int sscal_(int *, real *, real *, int *), 
	    sgemv_(char *, int *, int *, real *, real *, int *, 
	    real *, int *, real *, real *, int *), saxpy_(
	    int *, real *, real *, int *, real *, int *), ssymv_(
	    char *, int *, real *, real *, int *, real *, int *, 
	    real *, real *, int *);
    static int iw;
    extern /* Subroutine */ int slarfg_(int *, real *, real *, int *, 
	    real *);



#define E(I) e[(I)-1]
#define TAU(I) tau[(I)-1]

#define A(I,J) a[(I)-1 + ((J)-1)* ( *lda)]
#define W(I,J) w[(I)-1 + ((J)-1)* ( *ldw)]

    if (*n <= 0) {
	return 0;
    }

    if (lsame_(uplo, "U")) {

/*        Reduce last NB columns of upper triangle */

	i__1 = *n - *nb + 1;
	for (i = *n; i >= *n-*nb+1; --i) {
	    iw = i - *n + *nb;

/* Subroutine */
int ssytri_rook_(char *uplo, integer *n, real *a, integer * lda, integer *ipiv, real *work, integer *info)
{
    /* System generated locals */
    integer a_dim1, a_offset, i__1;
    real r__1;
    /* Local variables */
    real d__;
    integer k;
    real t, ak;
    integer kp;
    real akp1, temp;
    extern real sdot_(integer *, real *, integer *, real *, integer *);
    real akkp1;
    extern logical lsame_(char *, char *);
    integer kstep;
    logical upper;
    extern /* Subroutine */
    int scopy_(integer *, real *, integer *, real *, integer *), sswap_(integer *, real *, integer *, real *, integer * ), ssymv_(char *, integer *, real *, real *, integer *, real *, integer *, real *, real *, integer *), xerbla_(char *, integer *);
    /* -- LAPACK computational routine (version 3.4.1) -- */
    /* -- LAPACK is a software package provided by Univ. of Tennessee, -- */
    /* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- */
    /* April 2012 */
    /* .. Scalar Arguments .. */
    /* .. */
    /* .. Array Arguments .. */
    /* .. */
    /* ===================================================================== */
    /* .. Parameters .. */
    /* .. */
    /* .. Local Scalars .. */
    /* .. */
    /* .. External Functions .. */
    /* .. */
    /* .. External Subroutines .. */
    /* .. */
    /* .. Intrinsic Functions .. */
    /* .. */
    /* .. Executable Statements .. */
    /* Test the input parameters. */
    /* Parameter adjustments */
    a_dim1 = *lda;
    a_offset = 1 + a_dim1;
    a -= a_offset;
    --ipiv;
    --work;
    /* Function Body */
    *info = 0;
    upper = lsame_(uplo, "U");
    if (! upper && ! lsame_(uplo, "L"))
    {
        *info = -1;
    }
    else if (*n < 0)
    {
        *info = -2;
    }
    else if (*lda < max(1,*n))
    {
        *info = -4;
    }
    if (*info != 0)
    {
        i__1 = -(*info);
        xerbla_("SSYTRI_ROOK", &i__1);
        return 0;
    }
    /* Quick return if possible */
    if (*n == 0)
    {
        return 0;
    }
    /* Check that the diagonal matrix D is nonsingular. */
    if (upper)
    {
        /* Upper triangular storage: examine D from bottom to top */
        for (*info = *n;
                *info >= 1;
                --(*info))
        {
            if (ipiv[*info] > 0 && a[*info + *info * a_dim1] == 0.f)
            {
                return 0;
            }
            /* L10: */
        }
    }
    else
    {
        /* Lower triangular storage: examine D from top to bottom. */
        i__1 = *n;
        for (*info = 1;
                *info <= i__1;
                ++(*info))
        {
            if (ipiv[*info] > 0 && a[*info + *info * a_dim1] == 0.f)
            {
                return 0;
            }
            /* L20: */
//.........这里部分代码省略.........

/* Subroutine */ int snaitr_(integer *ido, char *bmat, integer *n, integer *k,
	 integer *np, integer *nb, real *resid, real *rnorm, real *v, integer 
	*ldv, real *h__, integer *ldh, integer *ipntr, real *workd, integer *
	info, ftnlen bmat_len)
{
    /* Initialized data */

    static logical first = TRUE_;

    /* System generated locals */
    integer h_dim1, h_offset, v_dim1, v_offset, i__1, i__2;
    real r__1, r__2;

    /* Builtin functions */
    double sqrt(doublereal);

    /* Local variables */
    static integer i__, j;
    static real t0, t1, t2, t3, t4, t5;
    static integer jj, ipj, irj, ivj;
    static real ulp, tst1;
    static integer ierr, iter;
    static real unfl, ovfl;
    extern doublereal sdot_(integer *, real *, integer *, real *, integer *);
    static integer itry;
    static real temp1;
    static logical orth1, orth2, step3, step4;
    extern doublereal snrm2_(integer *, real *, integer *);
    static real betaj;
    extern /* Subroutine */ int sscal_(integer *, real *, real *, integer *);
    static integer infol;
    extern /* Subroutine */ int sgemv_(char *, integer *, integer *, real *, 
	    real *, integer *, real *, integer *, real *, real *, integer *, 
	    ftnlen);
    static real xtemp[2];
    extern /* Subroutine */ int scopy_(integer *, real *, integer *, real *, 
	    integer *);
    static real wnorm;
    extern /* Subroutine */ int saxpy_(integer *, real *, real *, integer *, 
	    real *, integer *), ivout_(integer *, integer *, integer *, 
	    integer *, char *, ftnlen), smout_(integer *, integer *, integer *
	    , real *, integer *, integer *, char *, ftnlen), svout_(integer *,
	     integer *, real *, integer *, char *, ftnlen), sgetv0_(integer *,
	     char *, integer *, logical *, integer *, integer *, real *, 
	    integer *, real *, real *, integer *, real *, integer *, ftnlen);
    static real rnorm1;
    extern /* Subroutine */ int slabad_(real *, real *);
    extern doublereal slamch_(char *, ftnlen);
    extern /* Subroutine */ int second_(real *), slascl_(char *, integer *, 
	    integer *, real *, real *, integer *, integer *, real *, integer *
	    , integer *, ftnlen);
    static logical rstart;
    static integer msglvl;
    static real smlnum;
    extern doublereal slanhs_(char *, integer *, real *, integer *, real *, 
	    ftnlen);


/*     %----------------------------------------------------% */
/*     | Include files for debugging and timing information | */
/*     %----------------------------------------------------% */


/* \SCCS Information: @(#) */
/* FILE: debug.h   SID: 2.3   DATE OF SID: 11/16/95   RELEASE: 2 */

/*     %---------------------------------% */
/*     | See debug.doc for documentation | */
/*     %---------------------------------% */

/*     %------------------% */
/*     | Scalar Arguments | */
/*     %------------------% */

/*     %--------------------------------% */
/*     | See stat.doc for documentation | */
/*     %--------------------------------% */

/* \SCCS Information: @(#) */
/* FILE: stat.h   SID: 2.2   DATE OF SID: 11/16/95   RELEASE: 2 */



/*     %-----------------% */
/*     | Array Arguments | */
/*     %-----------------% */


/*     %------------% */
/*     | Parameters | */
/*     %------------% */


/*     %---------------% */
/*     | Local Scalars | */
/*     %---------------% */


/*     %-----------------------% */
/*     | Local Array Arguments | */
//.........这里部分代码省略.........

GURLS_EXPORT float dot(const int N, const float *X, const int incX, const float *Y, const int incY)
{
    return sdot_(const_cast<int*>(&N), const_cast<float*>(X), const_cast<int*>(&incX), const_cast<float*>(Y), const_cast<int*>(&incY));
}

 int cggbal_(char *job, int *n, complex *a, int *lda, 
	complex *b, int *ldb, int *ilo, int *ihi, float *lscale, 
	float *rscale, float *work, int *info)
{
    /* System generated locals */
    int a_dim1, a_offset, b_dim1, b_offset, i__1, i__2, i__3, i__4;
    float r__1, r__2, r__3;

    /* Builtin functions */
    double r_lg10(float *), r_imag(complex *), c_abs(complex *), r_sign(float *,
	     float *), pow_ri(float *, int *);

    /* Local variables */
    int i__, j, k, l, m;
    float t;
    int jc;
    float ta, tb, tc;
    int ir;
    float ew;
    int it, nr, ip1, jp1, lm1;
    float cab, rab, ewc, cor, sum;
    int nrp2, icab, lcab;
    float beta, coef;
    int irab, lrab;
    float basl, cmax;
    extern double sdot_(int *, float *, int *, float *, int *);
    float coef2, coef5, gamma, alpha;
    extern int lsame_(char *, char *);
    extern  int sscal_(int *, float *, float *, int *);
    float sfmin;
    extern  int cswap_(int *, complex *, int *, 
	    complex *, int *);
    float sfmax;
    int iflow, kount;
    extern  int saxpy_(int *, float *, float *, int *, 
	    float *, int *);
    float pgamma;
    extern int icamax_(int *, complex *, int *);
    extern double slamch_(char *);
    extern  int csscal_(int *, float *, complex *, int 
	    *), xerbla_(char *, int *);
    int lsfmin, lsfmax;


/*  -- LAPACK routine (version 3.2) -- */
/*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
/*     November 2006 */

/*     .. Scalar Arguments .. */
/*     .. */
/*     .. Array Arguments .. */
/*     .. */

/*  Purpose */
/*  ======= */

/*  CGGBAL balances a pair of general complex matrices (A,B).  This */
/*  involves, first, permuting A and B by similarity transformations to */
/*  isolate eigenvalues in the first 1 to ILO$-$1 and last IHI+1 to N */
/*  elements on the diagonal; and second, applying a diagonal similarity */
/*  transformation to rows and columns ILO to IHI to make the rows */
/*  and columns as close in norm as possible. Both steps are optional. */

/*  Balancing may reduce the 1-norm of the matrices, and improve the */
/*  accuracy of the computed eigenvalues and/or eigenvectors in the */
/*  generalized eigenvalue problem A*x = lambda*B*x. */

/*  Arguments */
/*  ========= */

/*  JOB     (input) CHARACTER*1 */
/*          Specifies the operations to be performed on A and B: */
/*          = 'N':  none:  simply set ILO = 1, IHI = N, LSCALE(I) = 1.0 */
/*                  and RSCALE(I) = 1.0 for i=1,...,N; */
/*          = 'P':  permute only; */
/*          = 'S':  scale only; */
/*          = 'B':  both permute and scale. */

/*  N       (input) INTEGER */
/*          The order of the matrices A and B.  N >= 0. */

/*  A       (input/output) COMPLEX array, dimension (LDA,N) */
/*          On entry, the input matrix A. */
/*          On exit, A is overwritten by the balanced matrix. */
/*          If JOB = 'N', A is not referenced. */

/*  LDA     (input) INTEGER */
/*          The leading dimension of the array A. LDA >= MAX(1,N). */

/*  B       (input/output) COMPLEX array, dimension (LDB,N) */
/*          On entry, the input matrix B. */
/*          On exit, B is overwritten by the balanced matrix. */
/*          If JOB = 'N', B is not referenced. */

/*  LDB     (input) INTEGER */
/*          The leading dimension of the array B. LDB >= MAX(1,N). */

/*  ILO     (output) INTEGER */
/*  IHI     (output) INTEGER */
/*          ILO and IHI are set to ints such that on exit */
//.........这里部分代码省略.........

/* Subroutine */
int slatrs_(char *uplo, char *trans, char *diag, char * normin, integer *n, real *a, integer *lda, real *x, real *scale, real *cnorm, integer *info)
{
    /* System generated locals */
    integer a_dim1, a_offset, i__1, i__2, i__3;
    real r__1, r__2, r__3;
    /* Local variables */
    integer i__, j;
    real xj, rec, tjj;
    integer jinc;
    real xbnd;
    integer imax;
    real tmax, tjjs;
    extern real sdot_(integer *, real *, integer *, real *, integer *);
    real xmax, grow, sumj;
    extern logical lsame_(char *, char *);
    extern /* Subroutine */
    int sscal_(integer *, real *, real *, integer *);
    real tscal, uscal;
    integer jlast;
    extern real sasum_(integer *, real *, integer *);
    logical upper;
    extern /* Subroutine */
    int saxpy_(integer *, real *, real *, integer *, real *, integer *), strsv_(char *, char *, char *, integer *, real *, integer *, real *, integer *);
    extern real slamch_(char *);
    extern /* Subroutine */
    int xerbla_(char *, integer *);
    real bignum;
    extern integer isamax_(integer *, real *, integer *);
    logical notran;
    integer jfirst;
    real smlnum;
    logical nounit;
    /* -- LAPACK auxiliary routine (version 3.4.2) -- */
    /* -- LAPACK is a software package provided by Univ. of Tennessee, -- */
    /* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- */
    /* September 2012 */
    /* .. Scalar Arguments .. */
    /* .. */
    /* .. Array Arguments .. */
    /* .. */
    /* ===================================================================== */
    /* .. Parameters .. */
    /* .. */
    /* .. Local Scalars .. */
    /* .. */
    /* .. External Functions .. */
    /* .. */
    /* .. External Subroutines .. */
    /* .. */
    /* .. Intrinsic Functions .. */
    /* .. */
    /* .. Executa