C++ saxpy_函数代码示例

int main(void)
{
	
     int i, N=4, inca=1, incb=1;
     float *a, *b, *bb, alpha=0.5;
     a=sf_floatalloc(N);
     b=sf_floatalloc(N);
     bb=sf_floatalloc(N);

     a[0]=1; a[1]=2; a[2]=3; a[3]=4; 
     b[0]=4; b[1]=3; b[2]=2; b[3]=1; 
     bb[0]=4; bb[1]=3; bb[2]=2; bb[3]=1; 
  
     for(i=0;i<N;i++)
     	sf_warning("a[%d]=%f",i,a[i]);
     for(i=0;i<N;i++)
     	sf_warning("b[%d]=%f",i,b[i]);

     cblas_saxpy(N, 0.5 , a, 1, b, 1);   

     for(i=0;i<N;i++)
     	sf_warning("(0.5*a+b)[%d]=%f",i,b[i]);

     saxpy_(&N, &alpha, a, &inca, bb, &incb);
 
     for(i=0;i<N;i++)
     	sf_warning("(0.5*a+bb)[%d]=%f",i,bb[i]);

    exit(0);
}

void THBlas_(axpy)(long n, real a, real *x, long incx, real *y, long 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)
    daxpy_(&i_n, &a, x, &i_incx, y, &i_incy);
#else
    saxpy_(&i_n, &a, x, &i_incx, y, &i_incy);
#endif
    return;
  }
#endif
  {
    long i;
    for(i = 0; i < n; i++)
      y[i*incy] += a*x[i*incx];
  }
}

int
f2c_saxpy(integer* N,
          real* alpha,
          real* X, integer* incX,
          real* Y, integer* incY)
{
    saxpy_(N, alpha, X, incX, Y, incY);
    return 0;
}

INLINE void do_saxpy(float *x, float *y, int iterations, int *limit)
{
  REGISTER int i = 0;
  extern int saxpy_();

  for (;i<iterations;i++)
    {
      saxpy_(limit,&falpha,x,&stride,y,&stride);
    }
}

/* 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 */
}

//.........这里部分代码省略.........
         M(j) <= M(j-1) * ( 1 + CNORM(j) ) / | A(j,j) |   

              <= M(0) * product ( ( 1 + CNORM(i) ) / |A(i,i)| )   
                        1<=i<=j   

    and we can safely call STPSV if 1/M(n) and 1/G(n) are both greater   
    than max(underflow, 1/overflow).   

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


       Parameter adjustments */
    /* Table of constant values */
    static integer c__1 = 1;
    static real c_b36 = .5f;
    
    /* System generated locals */
    integer i__1, i__2, i__3;
    real r__1, r__2, r__3;
    /* Local variables */
    static integer jinc, jlen;
    static real xbnd;
    static integer imax;
    static real tmax, tjjs;
    extern doublereal sdot_(integer *, real *, integer *, real *, integer *);
    static real xmax, grow, sumj;
    static integer i__, j;
    extern logical lsame_(char *, char *);
    extern /* Subroutine */ int sscal_(integer *, real *, real *, integer *);
    static real tscal, uscal;
    static integer jlast;
    extern doublereal sasum_(integer *, real *, integer *);
    static logical upper;
    extern /* Subroutine */ int saxpy_(integer *, real *, real *, integer *, 
	    real *, integer *), stpsv_(char *, char *, char *, integer *, 
	    real *, real *, integer *);
    static integer ip;
    static real xj;
    extern doublereal slamch_(char *);
    extern /* Subroutine */ int xerbla_(char *, integer *);
    static real bignum;
    extern integer isamax_(integer *, real *, integer *);
    static logical notran;
    static integer jfirst;
    static real smlnum;
    static logical nounit;
    static real rec, tjj;


    --cnorm;
    --x;
    --ap;

    /* Function Body */
    *info = 0;
    upper = lsame_(uplo, "U");
    notran = lsame_(trans, "N");
    nounit = lsame_(diag, "N");

/*     Test the input parameters. */

    if (! upper && ! lsame_(uplo, "L")) {
	*info = -1;
    } else if (! notran && ! lsame_(trans, "T") && ! 
	    lsame_(trans, "C")) {
	*info = -2;

 int slarz_(char *side, int *m, int *n, int *l, 
	float *v, int *incv, float *tau, float *c__, int *ldc, float *
	work)
{
    /* System generated locals */
    int c_dim1, c_offset;
    float r__1;

    /* Local variables */
    extern  int sger_(int *, int *, float *, float *, 
	    int *, float *, int *, float *, int *);
    extern int lsame_(char *, char *);
    extern  int sgemv_(char *, int *, int *, float *, 
	    float *, int *, float *, int *, float *, float *, int *), scopy_(int *, float *, int *, float *, int *), 
	    saxpy_(int *, float *, float *, int *, float *, int *);


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

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

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

/*  SLARZ applies a float elementary reflector H to a float M-by-N */
/*  matrix C, from either the left or the right. H is represented in the */
/*  form */

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

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

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


/*  H is a product of k elementary reflectors as returned by STZRZF. */

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

/*  SIDE    (input) CHARACTER*1 */
/*          = 'L': form  H * C */
/*          = 'R': form  C * H */

/*  M       (input) INTEGER */
/*          The number of rows of the matrix C. */

/*  N       (input) INTEGER */
/*          The number of columns of the matrix C. */

/*  L       (input) INTEGER */
/*          The number of entries of the vector V containing */
/*          the meaningful part of the Householder vectors. */
/*          If SIDE = 'L', M >= L >= 0, if SIDE = 'R', N >= L >= 0. */

/*  V       (input) REAL array, dimension (1+(L-1)*ABS(INCV)) */
/*          The vector v in the representation of H as returned by */
/*          STZRZF. V is not used if TAU = 0. */

/*  INCV    (input) INTEGER */
/*          The increment between elements of v. INCV <> 0. */

/*  TAU     (input) REAL */
/*          The value tau in the representation of H. */

/*  C       (input/output) REAL array, dimension (LDC,N) */
/*          On entry, the M-by-N matrix C. */
/*          On exit, C is overwritten by the matrix H * C if SIDE = 'L', */
/*          or C * H if SIDE = 'R'. */

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

/*  WORK    (workspace) REAL array, dimension */
/*                         (N) if SIDE = 'L' */
/*                      or (M) if SIDE = 'R' */

/*  Further Details */
/*  =============== */

/*  Based on contributions by */
/*    A. Petitet, Computer Science Dept., Univ. of Tenn., Knoxville, USA */

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

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

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

//.........这里部分代码省略.........
      (          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);
	return 0;
    }

/*     Quick return if possible */

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

//.........这里部分代码省略.........
            criterion is satisfied, should be at least 1.   

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


       Test the input parameters.   

    
   Parameter adjustments   
       Function Body */
    /* Table of constant values */
    static integer c__2 = 2;
    static integer c__1 = 1;
    static integer c_n1 = -1;
    
    /* System generated locals */
    integer z_dim1, z_offset, i__1, i__2, i__3;
    real r__1, r__2, r__3, r__4, r__5;
    /* Builtin functions */
    double sqrt(doublereal);
    /* Local variables */
    static integer jblk, nblk, jmax;
    extern doublereal sdot_(integer *, real *, integer *, real *, integer *), 
	    snrm2_(integer *, real *, integer *);
    static integer i, j, iseed[4], gpind, iinfo;
    extern /* Subroutine */ int sscal_(integer *, real *, real *, integer *);
    static integer b1;
    extern doublereal sasum_(integer *, real *, integer *);
    static integer j1;
    extern /* Subroutine */ int scopy_(integer *, real *, integer *, real *, 
	    integer *);
    static real ortol;
    extern /* Subroutine */ int saxpy_(integer *, real *, real *, integer *, 
	    real *, integer *);
    static integer indrv1, indrv2, indrv3, indrv4, indrv5, bn;
    static real xj;
    extern doublereal slamch_(char *);
    extern /* Subroutine */ int xerbla_(char *, integer *), slagtf_(
	    integer *, real *, real *, real *, real *, real *, real *, 
	    integer *, integer *);
    static integer nrmchk;
    extern integer isamax_(integer *, real *, integer *);
    extern /* Subroutine */ int slagts_(integer *, integer *, real *, real *, 
	    real *, real *, integer *, real *, real *, integer *);
    static integer blksiz;
    static real onenrm, pertol;
    extern /* Subroutine */ int slarnv_(integer *, integer *, integer *, real 
	    *);
    static real stpcrt, scl, eps, ctr, sep, nrm, tol;
    static integer its;
    static real xjm, eps1;



#define ISEED(I) iseed[(I)]
#define D(I) d[(I)-1]
#define E(I) e[(I)-1]
#define W(I) w[(I)-1]
#define IBLOCK(I) iblock[(I)-1]
#define ISPLIT(I) isplit[(I)-1]
#define WORK(I) work[(I)-1]
#define IWORK(I) iwork[(I)-1]
#define IFAIL(I) ifail[(I)-1]

#define Z(I,J) z[(I)-1 + ((J)-1)* ( *ldz)]

/* Subroutine */ int stbt03_(char *uplo, char *trans, char *diag, integer *n,
                             integer *kd, integer *nrhs, real *ab, integer *ldab, real *scale,
                             real *cnorm, real *tscal, real *x, integer *ldx, real *b, integer *
                             ldb, real *work, real *resid)
{
    /* System generated locals */
    integer ab_dim1, ab_offset, b_dim1, b_offset, x_dim1, x_offset, i__1;
    real r__1, r__2, r__3;

    /* Local variables */
    static integer j;
    extern logical lsame_(char *, char *);
    extern /* Subroutine */ int sscal_(integer *, real *, real *, integer *);
    static real xscal;
    extern /* Subroutine */ int stbmv_(char *, char *, char *, integer *,
                                       integer *, real *, integer *, real *, integer *), scopy_(integer *, real *, integer *, real *, integer *);
    static real tnorm, xnorm;
    extern /* Subroutine */ int saxpy_(integer *, real *, real *, integer *,
                                       real *, integer *), slabad_(real *, real *);
    static integer ix;
    extern doublereal slamch_(char *);
    static real bignum;
    extern integer isamax_(integer *, real *, integer *);
    static real smlnum, eps, err;


#define b_ref(a_1,a_2) b[(a_2)*b_dim1 + a_1]
#define x_ref(a_1,a_2) x[(a_2)*x_dim1 + a_1]
#define ab_ref(a_1,a_2) ab[(a_2)*ab_dim1 + a_1]


    /*  -- LAPACK test routine (version 3.0) --
           Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
           Courant Institute, Argonne National Lab, and Rice University
           February 29, 1992


        Purpose
        =======

        STBT03 computes the residual for the solution to a scaled triangular
        system of equations  A*x = s*b  or  A'*x = s*b  when A is a
        triangular band matrix. Here A' is the transpose of A, s is a scalar,
        and x and b are N by NRHS matrices.  The test ratio is the maximum
        over the number of right hand sides of
           norm(s*b - op(A)*x) / ( norm(op(A)) * norm(x) * EPS ),
        where op(A) denotes A or A' and EPS is the machine epsilon.

        Arguments
        =========

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

        TRANS   (input) CHARACTER*1
                Specifies the operation applied to A.
                = 'N':  A *x = b  (No transpose)
                = 'T':  A'*x = b  (Transpose)
                = 'C':  A'*x = b  (Conjugate transpose = Transpose)

        DIAG    (input) CHARACTER*1
                Specifies whether or not the matrix A is unit triangular.
                = 'N':  Non-unit triangular
                = 'U':  Unit triangular

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

        KD      (input) INTEGER
                The number of superdiagonals or subdiagonals of the
                triangular band matrix A.  KD >= 0.

        NRHS    (input) INTEGER
                The number of right hand sides, i.e., the number of columns
                of the matrices X and B.  NRHS >= 0.

        AB      (input) REAL array, dimension (LDAB,N)
                The upper or lower triangular band matrix A, stored in the
                first kd+1 rows of the array. The j-th column of A is stored
                in the j-th column of the array AB as follows:
                if UPLO = 'U', AB(kd+1+i-j,j) = A(i,j) for max(1,j-kd)<=i<=j;
                if UPLO = 'L', AB(1+i-j,j)    = A(i,j) for j<=i<=min(n,j+kd).

        LDAB    (input) INTEGER
                The leading dimension of the array AB.  LDAB >= KD+1.

        SCALE   (input) REAL
                The scaling factor s used in solving the triangular system.

        CNORM   (input) REAL array, dimension (N)
                The 1-norms of the columns of A, not counting the diagonal.

        TSCAL   (input) REAL
                The scaling factor used in computing the 1-norms in CNORM.
                CNORM actually contains the column norms of TSCAL*A.

        X       (input) REAL array, dimension (LDX,NRHS)
                The computed solution vectors for the system of linear
//.........这里部分代码省略.........

/* 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 .. */
    /* .. */
    /* .. Executable Statements .. */
    /* Parameter adjustments */
    a_dim1 = *lda;
    a_offset = 1 + a_dim1;
    a -= a_offset;
    --x;
    --cnorm;
    /* Function Body */
    *info = 0;
    upper = lsame_(uplo, "U");
    notran = lsame_(trans, "N");
    nounit = lsame_(diag, "N");
    /* Test the input parameters. */
    if (! upper && ! lsame_(uplo, "L"))
    {
        *info = -1;
    }
    else if (! notran && ! lsame_(trans, "T") && ! lsame_(trans, "C"))
    {
        *info = -2;
    }
    else if (! nounit && ! lsame_(diag, "U"))
    {
        *info = -3;
    }
    else if (! lsame_(normin, "Y") && ! lsame_(normin, "N"))
    {
        *info = -4;
    }
    else if (*n < 0)
    {
        *info = -5;
    }
    else if (*lda < max(1,*n))
    {
        *info = -7;
    }
    if (*info != 0)
    {
        i__1 = -(*info);
        xerbla_("SLATRS", &i__1);
        return 0;
    }
    /* Quick return if possible */
    if (*n == 0)
    {
        return 0;
    }
//.........这里部分代码省略.........

void saxpy( int n, float alpha, float *x, int incx,  float *y, int incy)
{
    saxpy_(&n, &alpha, x, &incx, y, &incy);
}

/* Subroutine */ int sptrfs_(integer *n, integer *nrhs, real *d__, real *e, 
	real *df, real *ef, real *b, integer *ldb, real *x, integer *ldx, 
	real *ferr, real *berr, real *work, integer *info)
{
    /* System generated locals */
    integer b_dim1, b_offset, x_dim1, x_offset, i__1, i__2;
    real r__1, r__2, r__3;

    /* Local variables */
    integer i__, j;
    real s, bi, cx, dx, ex;
    integer ix, nz;
    real eps, safe1, safe2;
    integer count;
    extern /* Subroutine */ int saxpy_(integer *, real *, real *, integer *, 
	    real *, integer *);
    extern doublereal slamch_(char *);
    real safmin;
    extern /* Subroutine */ int xerbla_(char *, integer *);
    extern integer isamax_(integer *, real *, integer *);
    real lstres;
    extern /* Subroutine */ int spttrs_(integer *, integer *, real *, real *, 
	    real *, integer *, integer *);


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

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

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

/*  SPTRFS improves the computed solution to a system of linear */
/*  equations when the coefficient matrix is symmetric positive definite */
/*  and tridiagonal, and provides error bounds and backward error */
/*  estimates for the solution. */

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

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

/*  NRHS    (input) INTEGER */
/*          The number of right hand sides, i.e., the number of columns */
/*          of the matrix B.  NRHS >= 0. */

/*  D       (input) REAL array, dimension (N) */
/*          The n diagonal elements of the tridiagonal matrix A. */

/*  E       (input) REAL array, dimension (N-1) */
/*          The (n-1) subdiagonal elements of the tridiagonal matrix A. */

/*  DF      (input) REAL array, dimension (N) */
/*          The n diagonal elements of the diagonal matrix D from the */
/*          factorization computed by SPTTRF. */

/*  EF      (input) REAL array, dimension (N-1) */
/*          The (n-1) subdiagonal elements of the unit bidiagonal factor */
/*          L from the factorization computed by SPTTRF. */

/*  B       (input) REAL array, dimension (LDB,NRHS) */
/*          The right hand side matrix B. */

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

/*  X       (input/output) REAL array, dimension (LDX,NRHS) */
/*          On entry, the solution matrix X, as computed by SPTTRS. */
/*          On exit, the improved solution matrix X. */

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

/*  FERR    (output) REAL array, dimension (NRHS) */
/*          The forward error bound for each solution vector */
/*          X(j) (the j-th column of the solution matrix X). */
/*          If XTRUE is the true solution corresponding to X(j), FERR(j) */
/*          is an estimated upper bound for the magnitude of the largest */
/*          element in (X(j) - XTRUE) divided by the magnitude of the */
/*          largest element in X(j). */

/*  BERR    (output) REAL array, dimension (NRHS) */
/*          The componentwise relative backward error of each solution */
/*          vector X(j) (i.e., the smallest relative change in */
/*          any element of A or B that makes X(j) an exact solution). */

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

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

/*  Internal Parameters */
/*  =================== */
//.........这里部分代码省略.........

//.........这里部分代码省略.........
	return 0;
    }

    eps = slamch_("Epsilon");

/*     Copy the matrix U to WORK. */

    i__1 = *n;
    for (j = 1; j <= i__1; ++j) {
	i__2 = *n;
	for (i__ = 1; i__ <= i__2; ++i__) {
	    work[i__ + j * work_dim1] = 0.f;
/* L10: */
	}
/* L20: */
    }
    i__1 = *n;
    for (i__ = 1; i__ <= i__1; ++i__) {
	if (i__ == 1) {
	    work[i__ + i__ * work_dim1] = df[i__];
	    if (*n >= 2) {
		work[i__ + (i__ + 1) * work_dim1] = duf[i__];
	    }
	    if (*n >= 3) {
		work[i__ + (i__ + 2) * work_dim1] = du2[i__];
	    }
	} else if (i__ == *n) {
	    work[i__ + i__ * work_dim1] = df[i__];
	} else {
	    work[i__ + i__ * work_dim1] = df[i__];
	    work[i__ + (i__ + 1) * work_dim1] = duf[i__];
	    if (i__ < *n - 1) {
		work[i__ + (i__ + 2) * work_dim1] = du2[i__];
	    }
	}
/* L30: */
    }

/*     Multiply on the left by L. */

    lastj = *n;
    for (i__ = *n - 1; i__ >= 1; --i__) {
	li = dlf[i__];
	i__1 = lastj - i__ + 1;
	saxpy_(&i__1, &li, &work[i__ + i__ * work_dim1], ldwork, &work[i__ + 
		1 + i__ * work_dim1], ldwork);
	ip = ipiv[i__];
	if (ip == i__) {
/* Computing MIN */
	    i__1 = i__ + 2;
	    lastj = min(i__1,*n);
	} else {
	    i__1 = lastj - i__ + 1;
	    sswap_(&i__1, &work[i__ + i__ * work_dim1], ldwork, &work[i__ + 1 
		    + i__ * work_dim1], ldwork);
	}
/* L40: */
    }

/*     Subtract the matrix A. */

    work[work_dim1 + 1] -= d__[1];
    if (*n > 1) {
	work[(work_dim1 << 1) + 1] -= du[1];
	work[*n + (*n - 1) * work_dim1] -= dl[*n - 1];
	work[*n + *n * work_dim1] -= d__[*n];
	i__1 = *n - 1;
	for (i__ = 2; i__ <= i__1; ++i__) {
	    work[i__ + (i__ - 1) * work_dim1] -= dl[i__ - 1];
	    work[i__ + i__ * work_dim1] -= d__[i__];
	    work[i__ + (i__ + 1) * work_dim1] -= du[i__];
/* L50: */
	}
    }

/*     Compute the 1-norm of the tridiagonal matrix A. */

    anorm = slangt_("1", n, &dl[1], &d__[1], &du[1]);

/*     Compute the 1-norm of WORK, which is only guaranteed to be */
/*     upper Hessenberg. */

    *resid = slanhs_("1", n, &work[work_offset], ldwork, &rwork[1])
	    ;

/*     Compute norm(L*U - A) / (norm(A) * EPS) */

    if (anorm <= 0.f) {
	if (*resid != 0.f) {
	    *resid = 1.f / eps;
	}
    } else {
	*resid = *resid / anorm / eps;
    }

    return 0;

/*     End of SGTT01 */

} /* sgtt01_ */

/* Subroutine */ int ssapps_(integer *n, integer *kev, integer *np, real *
                             shift, real *v, integer *ldv, real *h__, integer *ldh, real *resid,
                             real *q, integer *ldq, real *workd)
{
    /* Initialized data */

    static logical first = TRUE_;

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

    /* Local variables */
    static real c__, f, g;
    static integer i__, j;
    static real r__, s, a1, a2, a3, a4, t0, t1;
    static integer jj;
    static real big;
    static integer iend, itop;
    extern /* Subroutine */ int sscal_(integer *, real *, real *, integer *),
           sgemv_(char *, integer *, integer *, real *, real *, integer *,
                  real *, integer *, real *, real *, integer *, ftnlen), scopy_(
               integer *, real *, integer *, real *, integer *), saxpy_(integer *
                       , real *, real *, integer *, real *, integer *), ivout_(integer *,
                               integer *, integer *, integer *, char *, ftnlen), svout_(integer
                                       *, integer *, real *, integer *, char *, ftnlen);
    extern doublereal slamch_(char *, ftnlen);
    extern /* Subroutine */ int arscnd_(real *);
    static real epsmch;
    static integer istart, kplusp, msglvl;
    extern /* Subroutine */ int slacpy_(char *, integer *, integer *, real *,
                                        integer *, real *, integer *, ftnlen), slartg_(real *, real *,
                                                real *, real *, real *), slaset_(char *, integer *, integer *,
                                                        real *, real *, real *, integer *, 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 | */
    /*     %---------------% */



    /*     %----------------------% */
    /*     | External Subroutines | */
    /*     %----------------------% */


    /*     %--------------------% */
    /*     | External Functions | */
    /*     %--------------------% */


    /*     %----------------------% */
    /*     | Intrinsics Functions | */
    /*     %----------------------% */


    /*     %----------------% */
    /*     | Data statments | */
    /*     %----------------% */

    /* Parameter adjustments */
    --workd;
    --resid;
//.........这里部分代码省略.........