自学教程：C++ snrm2_函数代码示例

int toScalarF(int code, KFVEC(x), FVEC(r)) {    REQUIRES(rn==1,BAD_SIZE);    DEBUGMSG("toScalarF");    float res;    integer one = 1;    integer n = xn;    switch(code) {        case 0: { res = snrm2_(&n,xp,&one); break; }        case 1: { res = sasum_(&n,xp,&one);  break; }        case 2: { res = vector_max_index_f(V(x));  break; }        case 3: { res = vector_max_f(V(x));  break; }        case 4: { res = vector_min_index_f(V(x)); break; }        case 5: { res = vector_min_f(V(x)); break; }        default: ERROR(BAD_CODE);    }    rp[0] = res;    OK}

/* Subroutine */ int slaed3_(integer *k, integer *n, integer *n1, real *d__, 	real *q, integer *ldq, real *rho, real *dlamda, real *q2, integer *	indx, integer *ctot, real *w, real *s, integer *info){    /* System generated locals */    integer q_dim1, q_offset, i__1, i__2;    real r__1;    /* Builtin functions */    double sqrt(doublereal), r_sign(real *, real *);    /* Local variables */    integer i__, j, n2, n12, ii, n23, iq2;    real temp;    extern doublereal snrm2_(integer *, real *, integer *);    extern /* Subroutine */ int sgemm_(char *, char *, integer *, integer *, 	    integer *, real *, real *, integer *, real *, integer *, real *, 	    real *, integer *), scopy_(integer *, real *, 	    integer *, real *, integer *), slaed4_(integer *, integer *, real 	    *, real *, real *, real *, real *, integer *);    extern doublereal slamc3_(real *, real *);    extern /* Subroutine */ int xerbla_(char *, integer *), slacpy_(	    char *, integer *, integer *, real *, integer *, real *, integer *), slaset_(char *, integer *, integer *, real *, real *, 	    real *, integer *);/*  -- LAPACK routine (version 3.2) -- *//*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. *//*     November 2006 *//*     .. Scalar Arguments .. *//*     .. *//*     .. Array Arguments .. *//*     .. *//*  Purpose *//*  ======= *//*  SLAED3 finds the roots of the secular equation, as defined by the *//*  values in D, W, and RHO, between 1 and K.  It makes the *//*  appropriate calls to SLAED4 and then updates the eigenvectors by *//*  multiplying the matrix of eigenvectors of the pair of eigensystems *//*  being combined by the matrix of eigenvectors of the K-by-K system *//*  which is solved here. *//*  This code makes very mild assumptions about floating point *//*  arithmetic. It will work on machines with a guard digit in *//*  add/subtract, or on those binary machines without guard digits *//*  which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or Cray-2. *//*  It could conceivably fail on hexadecimal or decimal machines *//*  without guard digits, but we know of none. *//*  Arguments *//*  ========= *//*  K       (input) INTEGER *//*          The number of terms in the rational function to be solved by *//*          SLAED4.  K >= 0. *//*  N       (input) INTEGER *//*          The number of rows and columns in the Q matrix. *//*          N >= K (deflation may result in N>K). *//*  N1      (input) INTEGER *//*          The location of the last eigenvalue in the leading submatrix. *//*          min(1,N) <= N1 <= N/2. *//*  D       (output) REAL array, dimension (N) *//*          D(I) contains the updated eigenvalues for *//*          1 <= I <= K. *//*  Q       (output) REAL array, dimension (LDQ,N) *//*          Initially the first K columns are used as workspace. *//*          On output the columns 1 to K contain *//*          the updated eigenvectors. *//*  LDQ     (input) INTEGER *//*          The leading dimension of the array Q.  LDQ >= max(1,N). *//*  RHO     (input) REAL *//*          The value of the parameter in the rank one update equation. *//*          RHO >= 0 required. *//*  DLAMDA  (input/output) REAL array, dimension (K) *//*          The first K elements of this array contain the old roots *//*          of the deflated updating problem.  These are the poles *//*          of the secular equation. May be changed on output by *//*          having lowest order bit set to zero on Cray X-MP, Cray Y-MP, *//*          Cray-2, or Cray C-90, as described above. *//*  Q2      (input) REAL array, dimension (LDQ2, N) *//*          The first K columns of this matrix contain the non-deflated *//*          eigenvectors for the split problem. *//*  INDX    (input) INTEGER array, dimension (N) *//*          The permutation used to arrange the columns of the deflated *//*          Q matrix into three groups (see SLAED2). *//*          The rows of the eigenvectors found by SLAED4 must be likewise *//*          permuted before the matrix multiply can take place. *///.........这里部分代码省略.........

doublereal sqrt12_(integer *m, integer *n, real *a, integer *lda, real *s, 	real *work, integer *lwork){    /* System generated locals */    integer a_dim1, a_offset, i__1, i__2;    real ret_val;    /* Local variables */    integer i__, j, mn, iscl, info;    real anrm;    extern doublereal snrm2_(integer *, real *, integer *), sasum_(integer *, 	    real *, integer *);    real dummy[1];    extern /* Subroutine */ int saxpy_(integer *, real *, real *, integer *, 	    real *, integer *), sgebd2_(integer *, integer *, real *, integer 	    *, real *, real *, real *, real *, real *, integer *), slabad_(	    real *, real *);    extern doublereal slamch_(char *), slange_(char *, integer *, 	    integer *, real *, integer *, real *);    extern /* Subroutine */ int xerbla_(char *, integer *);    real bignum;    extern /* Subroutine */ int slascl_(char *, integer *, integer *, real *, 	    real *, integer *, integer *, real *, integer *, integer *), slaset_(char *, integer *, integer *, real *, real *, 	    real *, integer *), sbdsqr_(char *, integer *, integer *, 	    integer *, integer *, real *, real *, real *, integer *, real *, 	    integer *, real *, integer *, real *, integer *);    real smlnum, nrmsvl;/*  -- LAPACK test routine (version 3.1.1) -- *//*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. *//*     January 2007 *//*     .. Scalar Arguments .. *//*     .. *//*     .. Array Arguments .. *//*     .. *//*  Purpose *//*  ======= *//*  SQRT12 computes the singular values `svlues' of the upper trapezoid *//*  of A(1:M,1:N) and returns the ratio *//*       || s - svlues||/(||svlues||*eps*max(M,N)) *//*  Arguments *//*  ========= *//*  M       (input) INTEGER *//*          The number of rows of the matrix A. *//*  N       (input) INTEGER *//*          The number of columns of the matrix A. *//*  A       (input) REAL array, dimension (LDA,N) *//*          The M-by-N matrix A. Only the upper trapezoid is referenced. *//*  LDA     (input) INTEGER *//*          The leading dimension of the array A. *//*  S       (input) REAL array, dimension (min(M,N)) *//*          The singular values of the matrix A. *//*  WORK    (workspace) REAL array, dimension (LWORK) *//*  LWORK   (input) INTEGER *//*          The length of the array WORK. LWORK >= max(M*N + 4*min(M,N) + *//*          max(M,N), M*N+2*MIN( M, N )+4*N). *//*  ===================================================================== *//*     .. Parameters .. *//*     .. *//*     .. Local Scalars .. *//*     .. *//*     .. External Functions .. *//*     .. *//*     .. External Subroutines .. *//*     .. *//*     .. Intrinsic Functions .. *//*     .. *//*     .. Local Arrays .. *//*     .. *//*     .. Executable Statements .. */    /* Parameter adjustments */    a_dim1 = *lda;    a_offset = 1 + a_dim1;    a -= a_offset;    --s;    --work;    /* Function Body */    ret_val = 0.f;/*     Test that enough workspace is supplied *//* Computing MAX */    i__1 = *m * *n + (min(*m,*n) << 2) + max(*m,*n), i__2 = *m * *n + (min(*m,//.........这里部分代码省略.........

/* Subroutine */ int slarfp_(integer *n, real *alpha, real *x, integer *incx, 	real *tau){    /* System generated locals */    integer i__1;    real r__1;    /* Local variables */    integer j, knt;    real beta;    real xnorm;    real safmin, rsafmn;/*  -- LAPACK auxiliary routine (version 3.2) -- *//*     November 2006 *//*  Purpose *//*  ======= *//*  SLARFP generates a real elementary reflector H of order n, such *//*  that *//*        H * ( alpha ) = ( beta ),   H' * H = I. *//*            (   x   )   (   0  ) *//*  where alpha and beta are scalars, beta is non-negative, and x is *//*  an (n-1)-element real vector.  H is represented in the form *//*        H = I - tau * ( 1 ) * ( 1 v' ) , *//*                      ( v ) *//*  where tau is a real scalar and v is a real (n-1)-element *//*  vector. *//*  If the elements of x are all zero, then tau = 0 and H is taken to be *//*  the unit matrix. *//*  Otherwise  1 <= tau <= 2. *//*  Arguments *//*  ========= *//*  N       (input) INTEGER *//*          The order of the elementary reflector. *//*  ALPHA   (input/output) REAL *//*          On entry, the value alpha. *//*          On exit, it is overwritten with the value beta. *//*  X       (input/output) REAL array, dimension *//*                         (1+(N-2)*abs(INCX)) *//*          On entry, the vector x. *//*          On exit, it is overwritten with the vector v. *//*  INCX    (input) INTEGER *//*          The increment between elements of X. INCX > 0. *//*  TAU     (output) REAL *//*          The value tau. *//*  ===================================================================== */    /* Parameter adjustments */    --x;    /* Function Body */    if (*n <= 0) {	*tau = 0.f;	return 0;    }    i__1 = *n - 1;    xnorm = snrm2_(&i__1, &x[1], incx);    if (xnorm == 0.f) {/*        H  =  [+/-1, 0; I], sign chosen so ALPHA >= 0. */	if (*alpha >= 0.f) {/*           When TAU.eq.ZERO, the vector is special-cased to be *//*           all zeros in the application routines.  We do not need *//*           to clear it. */	    *tau = 0.f;	} else {/*           However, the application routines rely on explicit *//*           zero checks when TAU.ne.ZERO, and we must clear X. */	    *tau = 2.f;	    i__1 = *n - 1;	    for (j = 1; j <= i__1; ++j) {		x[(j - 1) * *incx + 1] = 0.f;	    }	    *alpha = -(*alpha);	}    } else {/*        general case */	r__1 = slapy2_(alpha, &xnorm);	beta = r_sign(&r__1, alpha);	safmin = slamch_("S") / slamch_("E");//.........这里部分代码省略.........

//.........这里部分代码省略.........    Internal Parameters       ===================       MAXITS  INTEGER, default = 5               The maximum number of iterations performed.       EXTRA   INTEGER, default = 2               The number of iterations performed after norm growth               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]

int sfgmr(int n,     void (*smatvec) (float, float[], float, float[]),     void (*spsolve) (int, float[], float[]),     float *rhs, float *sol, double tol, int im, int *itmax, FILE * fits){/*----------------------------------------------------------------------|                 *** Preconditioned FGMRES ***+-----------------------------------------------------------------------| This is a simple version of the ARMS preconditioned FGMRES algorithm.+-----------------------------------------------------------------------| Y. S. Dec. 2000. -- Apr. 2008+-----------------------------------------------------------------------| on entry:|----------|| rhs     = real vector of length n containing the right hand side.| sol     = real vector of length n containing an initial guess to the|           solution on input.| tol     = tolerance for stopping iteration| im      = Krylov subspace dimension| (itmax) = max number of iterations allowed.| fits    = NULL: no output|        != NULL: file handle to output " resid vs time and its"|| on return:|----------| fgmr      int =  0 --> successful return.|           int =  1 --> convergence not achieved in itmax iterations.| sol     = contains an approximate solution (upon successful return).| itmax   = has changed. It now contains the number of steps required|           to converge --+-----------------------------------------------------------------------| internal work arrays:|----------| vv      = work array of length [im+1][n] (used to store the Arnoldi|           basis)| hh      = work array of length [im][im+1] (Householder matrix)| z       = work array of length [im][n] to store preconditioned vectors+-----------------------------------------------------------------------| subroutines called :| matvec - matrix-vector multiplication operation| psolve - (right) preconditionning operation|	   psolve can be a NULL pointer (GMRES without preconditioner)+---------------------------------------------------------------------*/    int maxits = *itmax;    int i, i1, ii, j, k, k1, its, retval, i_1 = 1, i_2 = 2;    float beta, eps1 = 0.0, t, t0, gam;    float **hh, *c, *s, *rs;    float **vv, **z, tt;    float zero = 0.0;    float one = 1.0;    its = 0;    vv = (float **)SUPERLU_MALLOC((im + 1) * sizeof(float *));    for (i = 0; i <= im; i++) vv[i] = floatMalloc(n);    z = (float **)SUPERLU_MALLOC(im * sizeof(float *));    hh = (float **)SUPERLU_MALLOC(im * sizeof(float *));    for (i = 0; i < im; i++)    {	hh[i] = floatMalloc(i + 2);	z[i] = floatMalloc(n);    }    c = floatMalloc(im);    s = floatMalloc(im);    rs = floatMalloc(im + 1);    /*---- outer loop starts here ----*/    do    {	/*---- compute initial residual vector ----*/	smatvec(one, sol, zero, vv[0]);	for (j = 0; j < n; j++)	    vv[0][j] = rhs[j] - vv[0][j];	/* vv[0]= initial residual */	beta = snrm2_(&n, vv[0], &i_1);	/*---- print info if fits != null ----*/	if (fits != NULL && its == 0)	    fprintf(fits, "%8d   %10.2e/n", its, beta);	/*if ( beta <= tol * dnrm2_(&n, rhs, &i_1) )*/	if ( !(beta > tol * snrm2_(&n, rhs, &i_1)) )	    break;	t = 1.0 / beta;	/*---- normalize: vv[0] = vv[0] / beta ----*/	for (j = 0; j < n; j++)	    vv[0][j] = vv[0][j] * t;	if (its == 0)	    eps1 = tol * beta;	/*---- initialize 1-st term of rhs of hessenberg system ----*/	rs[0] = beta;	for (i = 0; i < im; i++)	{	    its++;	    i1 = i + 1;	    /*------------------------------------------------------------	    |  (Right) Preconditioning Operation   z_{j} = M^{-1} v_{j}	    +-----------------------------------------------------------*///.........这里部分代码省略.........

float snrm2( int n, float *x, int incx){    return snrm2_(&n, x, &incx);}

/* Subroutine */ int check1_(real *sfac){    /* Initialized data */    static real sa[10] = { .3f,-1.f,0.f,1.f,.3f,.3f,.3f,.3f,.3f,.3f };    static real dv[80]	/* was [8][5][2] */ = { .1f,2.f,2.f,2.f,2.f,2.f,2.f,	    2.f,.3f,3.f,3.f,3.f,3.f,3.f,3.f,3.f,.3f,-.4f,4.f,4.f,4.f,4.f,4.f,	    4.f,.2f,-.6f,.3f,5.f,5.f,5.f,5.f,5.f,.1f,-.3f,.5f,-.1f,6.f,6.f,	    6.f,6.f,.1f,8.f,8.f,8.f,8.f,8.f,8.f,8.f,.3f,9.f,9.f,9.f,9.f,9.f,	    9.f,9.f,.3f,2.f,-.4f,2.f,2.f,2.f,2.f,2.f,.2f,3.f,-.6f,5.f,.3f,2.f,	    2.f,2.f,.1f,4.f,-.3f,6.f,-.5f,7.f,-.1f,3.f };    static real dtrue1[5] = { 0.f,.3f,.5f,.7f,.6f };    static real dtrue3[5] = { 0.f,.3f,.7f,1.1f,1.f };    static real dtrue5[80]	/* was [8][5][2] */ = { .1f,2.f,2.f,2.f,2.f,	    2.f,2.f,2.f,-.3f,3.f,3.f,3.f,3.f,3.f,3.f,3.f,0.f,0.f,4.f,4.f,4.f,	    4.f,4.f,4.f,.2f,-.6f,.3f,5.f,5.f,5.f,5.f,5.f,.03f,-.09f,.15f,	    -.03f,6.f,6.f,6.f,6.f,.1f,8.f,8.f,8.f,8.f,8.f,8.f,8.f,.09f,9.f,	    9.f,9.f,9.f,9.f,9.f,9.f,.09f,2.f,-.12f,2.f,2.f,2.f,2.f,2.f,.06f,	    3.f,-.18f,5.f,.09f,2.f,2.f,2.f,.03f,4.f,-.09f,6.f,-.15f,7.f,-.03f,	    3.f };    static integer itrue2[5] = { 0,1,2,2,3 };    /* System generated locals */    integer i__1;    real r__1;    /* Builtin functions */    integer s_wsle(cilist *), do_lio(integer *, integer *, char *, ftnlen), 	    e_wsle(void);    /* Subroutine */ int s_stop(char *, ftnlen);    /* Local variables */    integer i__;    real sx[8];    integer np1, len;    extern doublereal snrm2_(integer *, real *, integer *);    extern /* Subroutine */ int sscal_(integer *, real *, real *, integer *);    real stemp[1];    extern doublereal sasum_(integer *, real *, integer *);    real strue[8];    extern /* Subroutine */ int stest_(integer *, real *, real *, real *, 	    real *), itest1_(integer *, integer *), stest1_(real *, real *, 	    real *, real *);    extern integer isamax_(integer *, real *, integer *);    /* Fortran I/O blocks */    static cilist io___32 = { 0, 6, 0, 0, 0 };/*     .. Parameters .. *//*     .. Scalar Arguments .. *//*     .. Scalars in Common .. *//*     .. Local Scalars .. *//*     .. Local Arrays .. *//*     .. External Functions .. *//*     .. External Subroutines .. *//*     .. Intrinsic Functions .. *//*     .. Common blocks .. *//*     .. Data statements .. *//*     .. Executable Statements .. */    for (combla_1.incx = 1; combla_1.incx <= 2; ++combla_1.incx) {	for (np1 = 1; np1 <= 5; ++np1) {	    combla_1.n = np1 - 1;	    len = max(combla_1.n,1) << 1;/*           .. Set vector arguments .. */	    i__1 = len;	    for (i__ = 1; i__ <= i__1; ++i__) {		sx[i__ - 1] = dv[i__ + (np1 + combla_1.incx * 5 << 3) - 49];/* L20: */	    }	    if (combla_1.icase == 7) {/*              .. SNRM2 .. */		stemp[0] = dtrue1[np1 - 1];		r__1 = snrm2_(&combla_1.n, sx, &combla_1.incx);		stest1_(&r__1, stemp, stemp, sfac);	    } else if (combla_1.icase == 8) {/*              .. SASUM .. */		stemp[0] = dtrue3[np1 - 1];		r__1 = sasum_(&combla_1.n, sx, &combla_1.incx);		stest1_(&r__1, stemp, stemp, sfac);	    } else if (combla_1.icase == 9) {/*              .. SSCAL .. */		sscal_(&combla_1.n, &sa[(combla_1.incx - 1) * 5 + np1 - 1], 			sx, &combla_1.incx);		i__1 = len;		for (i__ = 1; i__ <= i__1; ++i__) {		    strue[i__ - 1] = dtrue5[i__ + (np1 + combla_1.incx * 5 << 			    3) - 49];/* L40: */		}		stest_(&len, sx, strue, strue, sfac);	    } else if (combla_1.icase == 10) {/*              .. ISAMAX .. */		i__1 = isamax_(&combla_1.n, sx, &combla_1.incx);		itest1_(&i__1, &itrue2[np1 - 1]);	    } else {		s_wsle(&io___32);		do_lio(&c__9, &c__1, " Shouldn't be here in CHECK1", (ftnlen)			28);//.........这里部分代码省略.........

 int slaqps_(int *m, int *n, int *offset, int 	*nb, int *kb, float *a, int *lda, int *jpvt, float *tau, 	float *vn1, float *vn2, float *auxv, float *f, int *ldf){    /* System generated locals */    int a_dim1, a_offset, f_dim1, f_offset, i__1, i__2;    float r__1, r__2;    /* Builtin functions */    double sqrt(double);    int i_nint(float *);    /* Local variables */    int j, k, rk;    float akk;    int pvt;    float temp, temp2;    extern double snrm2_(int *, float *, int *);    float tol3z;    extern  int sgemm_(char *, char *, int *, int *, 	    int *, float *, float *, int *, float *, int *, float *, 	    float *, int *);    int itemp;    extern  int sgemv_(char *, int *, int *, float *, 	    float *, int *, float *, int *, float *, float *, int *), sswap_(int *, float *, int *, float *, int *);    extern double slamch_(char *);    int lsticc;    extern int isamax_(int *, float *, int *);    extern  int slarfp_(int *, float *, float *, int *, 	    float *);    int lastrk;/*  -- LAPACK auxiliary routine (version 3.2) -- *//*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. *//*     November 2006 *//*     .. Scalar Arguments .. *//*     .. *//*     .. Array Arguments .. *//*     .. *//*  Purpose *//*  ======= *//*  SLAQPS computes a step of QR factorization with column pivoting *//*  of a float M-by-N matrix A by using Blas-3.  It tries to factorize *//*  NB columns from A starting from the row OFFSET+1, and updates all *//*  of the matrix with Blas-3 xGEMM. *//*  In some cases, due to catastrophic cancellations, it cannot *//*  factorize NB columns.  Hence, the actual number of factorized *//*  columns is returned in KB. *//*  Block A(1:OFFSET,1:N) is accordingly pivoted, but not factorized. *//*  Arguments *//*  ========= *//*  M       (input) INTEGER *//*          The number of rows of the matrix A. M >= 0. *//*  N       (input) INTEGER *//*          The number of columns of the matrix A. N >= 0 *//*  OFFSET  (input) INTEGER *//*          The number of rows of A that have been factorized in *//*          previous steps. *//*  NB      (input) INTEGER *//*          The number of columns to factorize. *//*  KB      (output) INTEGER *//*          The number of columns actually factorized. *//*  A       (input/output) REAL array, dimension (LDA,N) *//*          On entry, the M-by-N matrix A. *//*          On exit, block A(OFFSET+1:M,1:KB) is the triangular *//*          factor obtained and block A(1:OFFSET,1:N) has been *//*          accordingly pivoted, but no factorized. *//*          The rest of the matrix, block A(OFFSET+1:M,KB+1:N) has *//*          been updated. *//*  LDA     (input) INTEGER *//*          The leading dimension of the array A. LDA >= MAX(1,M). *//*  JPVT    (input/output) INTEGER array, dimension (N) *//*          JPVT(I) = K <==> Column K of the full matrix A has been *//*          permuted into position I in AP. *//*  TAU     (output) REAL array, dimension (KB) *//*          The scalar factors of the elementary reflectors. *//*  VN1     (input/output) REAL array, dimension (N) *//*          The vector with the partial column norms. *//*  VN2     (input/output) REAL array, dimension (N) *//*          The vector with the exact column norms. *//*  AUXV    (input/output) REAL array, dimension (NB) *///.........这里部分代码省略.........

/* Subroutine */ int slaror_(char *side, char *init, integer *m, integer *n, 	real *a, integer *lda, integer *iseed, real *x, integer *info){    /* System generated locals */    integer a_dim1, a_offset, i__1, i__2;    real r__1;    /* Builtin functions */    double r_sign(real *, real *);    /* Local variables */    static integer kbeg, jcol;    extern /* Subroutine */ int sger_(integer *, integer *, real *, real *, 	    integer *, real *, integer *, real *, integer *);    static integer irow;    extern real snrm2_(integer *, real *, integer *);    static integer j;    extern logical lsame_(char *, char *);    extern /* Subroutine */ int sscal_(integer *, real *, real *, integer *), 	    sgemv_(char *, integer *, integer *, real *, real *, integer *, 	    real *, integer *, real *, real *, integer *);    static integer ixfrm, itype, nxfrm;    static real xnorm;    extern /* Subroutine */ int xerbla_(char *, integer *);    static real factor;    extern doublereal slarnd_(integer *, integer *);    extern /* Subroutine */ int slaset_(char *, integer *, integer *, real *, 	    real *, real *, integer *);    static real xnorms;/*  -- LAPACK auxiliary test routine (version 2.0) --          Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,          Courant Institute, Argonne National Lab, and Rice University          September 30, 1994       Purpose       =======       SLAROR pre- or post-multiplies an M by N matrix A by a random       orthogonal matrix U, overwriting A.  A may optionally be initialized       to the identity matrix before multiplying by U.  U is generated using       the method of G.W. Stewart (SIAM J. Numer. Anal. 17, 1980, 403-409).       Arguments       =========       SIDE    (input) CHARACTER*1               Specifies whether A is multiplied on the left or right by U.               = 'L':         Multiply A on the left (premultiply) by U               = 'R':         Multiply A on the right (postmultiply) by U'               = 'C' or 'T':  Multiply A on the left by U and the right                               by U' (Here, U' means U-transpose.)       INIT    (input) CHARACTER*1               Specifies whether or not A should be initialized to the               identity matrix.               = 'I':  Initialize A to (a section of) the identity matrix                        before applying U.               = 'N':  No initialization.  Apply U to the input matrix A.               INIT = 'I' may be used to generate square or rectangular               orthogonal matrices:               For M = N and SIDE = 'L' or 'R', the rows will be orthogonal               to each other, as will the columns.               If M < N, SIDE = 'R' produces a dense matrix whose rows are               orthogonal and whose columns are not, while SIDE = 'L'               produces a matrix whose rows are orthogonal, and whose first               M columns are orthogonal, and whose remaining columns are               zero.               If M > N, SIDE = 'L' produces a dense matrix whose columns               are orthogonal and whose rows are not, while SIDE = 'R'               produces a matrix whose columns are orthogonal, and whose               first M rows are orthogonal, and whose remaining rows are               zero.       M       (input) INTEGER               The number of rows of A.       N       (input) INTEGER               The number of columns of A.       A       (input/output) REAL array, dimension (LDA, N)               On entry, the array A.               On exit, overwritten by U A ( if SIDE = 'L' ),                or by A U ( if SIDE = 'R' ),                or by U A U' ( if SIDE = 'C' or 'T').       LDA     (input) INTEGER               The leading dimension of the array A.  LDA >= max(1,M).   //.........这里部分代码省略.........

//.........这里部分代码省略.........            = 0:  successful exit.               < 0:  if INFO = -i, the i-th argument had an illegal value.       Further Details       ===============       Based on contributions by          Ming Gu and Ren-Cang Li, Computer Science Division, University of            California at Berkeley, USA          Osni Marques, LBNL/NERSC, USA       =====================================================================          Test the input parameters.          Parameter adjustments */    /* Table of constant values */    static real c_b5 = -1.f;    static integer c__1 = 1;    static real c_b11 = 1.f;    static real c_b13 = 0.f;    static integer c__0 = 0;        /* System generated locals */    integer givcol_dim1, givcol_offset, b_dim1, b_offset, bx_dim1, bx_offset, 	    difr_dim1, difr_offset, givnum_dim1, givnum_offset, poles_dim1, 	    poles_offset, i__1, i__2;    real r__1;    /* Local variables */    static real temp;    extern /* Subroutine */ int srot_(integer *, real *, integer *, real *, 	    integer *, real *, real *);    extern doublereal snrm2_(integer *, real *, integer *);    static integer i__, j, m, n;    static real diflj, difrj, dsigj;    extern /* Subroutine */ int sscal_(integer *, real *, real *, integer *), 	    sgemv_(char *, integer *, integer *, real *, real *, integer *, 	    real *, integer *, real *, real *, integer *), scopy_(	    integer *, real *, integer *, real *, integer *);    extern doublereal slamc3_(real *, real *);    static real dj;    extern /* Subroutine */ int xerbla_(char *, integer *);    static real dsigjp;    extern /* Subroutine */ int slascl_(char *, integer *, integer *, real *, 	    real *, integer *, integer *, real *, integer *, integer *), slacpy_(char *, integer *, integer *, real *, integer *, 	    real *, integer *);    static integer nlp1;#define difr_ref(a_1,a_2) difr[(a_2)*difr_dim1 + a_1]#define b_ref(a_1,a_2) b[(a_2)*b_dim1 + a_1]#define poles_ref(a_1,a_2) poles[(a_2)*poles_dim1 + a_1]#define bx_ref(a_1,a_2) bx[(a_2)*bx_dim1 + a_1]#define givcol_ref(a_1,a_2) givcol[(a_2)*givcol_dim1 + a_1]#define givnum_ref(a_1,a_2) givnum[(a_2)*givnum_dim1 + a_1]    b_dim1 = *ldb;    b_offset = 1 + b_dim1 * 1;    b -= b_offset;    bx_dim1 = *ldbx;    bx_offset = 1 + bx_dim1 * 1;    bx -= bx_offset;    --perm;    givcol_dim1 = *ldgcol;    givcol_offset = 1 + givcol_dim1 * 1;    givcol -= givcol_offset;

 int slaein_(int *rightv, int *noinit, int *n, 	float *h__, int *ldh, float *wr, float *wi, float *vr, float *vi, float 	*b, int *ldb, float *work, float *eps3, float *smlnum, float *bignum, 	int *info){    /* System generated locals */    int b_dim1, b_offset, h_dim1, h_offset, i__1, i__2, i__3, i__4;    float r__1, r__2, r__3, r__4;    /* Builtin functions */    double sqrt(double);    /* Local variables */    int i__, j;    float w, x, y;    int i1, i2, i3;    float w1, ei, ej, xi, xr, rec;    int its, ierr;    float temp, norm, vmax;    extern double snrm2_(int *, float *, int *);    float scale;    extern  int sscal_(int *, float *, float *, int *);    char trans[1];    float vcrit;    extern double sasum_(int *, float *, int *);    float rootn, vnorm;    extern double slapy2_(float *, float *);    float absbii, absbjj;    extern int isamax_(int *, float *, int *);    extern  int sladiv_(float *, float *, float *, float *, float *, float *);    char normin[1];    float nrmsml;    extern  int slatrs_(char *, char *, char *, char *, 	    int *, float *, int *, float *, float *, float *, int *);    float growto;/*  -- LAPACK auxiliary routine (version 3.2) -- *//*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. *//*     November 2006 *//*     .. Scalar Arguments .. *//*     .. *//*     .. Array Arguments .. *//*     .. *//*  Purpose *//*  ======= *//*  SLAEIN uses inverse iteration to find a right or left eigenvector *//*  corresponding to the eigenvalue (WR,WI) of a float upper Hessenberg *//*  matrix H. *//*  Arguments *//*  ========= *//*  RIGHTV   (input) LOGICAL *//*          = .TRUE. : compute right eigenvector; *//*          = .FALSE.: compute left eigenvector. *//*  NOINIT   (input) LOGICAL *//*          = .TRUE. : no initial vector supplied in (VR,VI). *//*          = .FALSE.: initial vector supplied in (VR,VI). *//*  N       (input) INTEGER *//*          The order of the matrix H.  N >= 0. *//*  H       (input) REAL array, dimension (LDH,N) *//*          The upper Hessenberg matrix H. *//*  LDH     (input) INTEGER *//*          The leading dimension of the array H.  LDH >= MAX(1,N). *//*  WR      (input) REAL *//*  WI      (input) REAL *//*          The float and imaginary parts of the eigenvalue of H whose *//*          corresponding right or left eigenvector is to be computed. *//*  VR      (input/output) REAL array, dimension (N) *//*  VI      (input/output) REAL array, dimension (N) *//*          On entry, if NOINIT = .FALSE. and WI = 0.0, VR must contain *//*          a float starting vector for inverse iteration using the float *//*          eigenvalue WR; if NOINIT = .FALSE. and WI.ne.0.0, VR and VI *//*          must contain the float and imaginary parts of a complex *//*          starting vector for inverse iteration using the complex *//*          eigenvalue (WR,WI); otherwise VR and VI need not be set. *//*          On exit, if WI = 0.0 (float eigenvalue), VR contains the *//*          computed float eigenvector; if WI.ne.0.0 (complex eigenvalue), *//*          VR and VI contain the float and imaginary parts of the *//*          computed complex eigenvector. The eigenvector is normalized *//*          so that the component of largest magnitude has magnitude 1; *//*          here the magnitude of a complex number (x,y) is taken to be *//*          |x| + |y|. *//*          VI is not referenced if WI = 0.0. *//*  B       (workspace) REAL array, dimension (LDB,N) *//*  LDB     (input) INTEGER *//*          The leading dimension of the array B.  LDB >= N+1. *///.........这里部分代码省略.........

/*! /brief * <pre> * Purpose * ======= *    ilu_sdrop_row() - Drop some small rows from the previous *    supernode (L-part only). * </pre> */int ilu_sdrop_row(        superlu_options_t *options, /* options */        int    first,       /* index of the first column in the supernode */        int    last,        /* index of the last column in the supernode */        double drop_tol,    /* dropping parameter */        int    quota,       /* maximum nonzero entries allowed */        int    *nnzLj,      /* in/out number of nonzeros in L(:, 1:last) */        double *fill_tol,   /* in/out - on exit, fill_tol=-num_zero_pivots,                             * does not change if options->ILU_MILU != SMILU1 */        GlobalLU_t *Glu,    /* modified */        float swork[],   /* working space                             * the length of swork[] should be no less than                             * the number of rows in the supernode */        float swork2[], /* working space with the same size as swork[],                             * used only by the second dropping rule */        int    lastc        /* if lastc == 0, there is nothing after the                             * working supernode [first:last];                             * if lastc == 1, there is one more column after                             * the working supernode. */ ){    register int i, j, k, m1;    register int nzlc; /* number of nonzeros in column last+1 */    register int xlusup_first, xlsub_first;    int m, n; /* m x n is the size of the supernode */    int r = 0; /* number of dropped rows */    register float *temp;    register float *lusup = Glu->lusup;    register int *lsub = Glu->lsub;    register int *xlsub = Glu->xlsub;    register int *xlusup = Glu->xlusup;    register float d_max = 0.0, d_min = 1.0;    int    drop_rule = options->ILU_DropRule;    milu_t milu = options->ILU_MILU;    norm_t nrm = options->ILU_Norm;    float zero = 0.0;    float one = 1.0;    float none = -1.0;    int i_1 = 1;    int inc_diag; /* inc_diag = m + 1 */    int nzp = 0;  /* number of zero pivots */    float alpha = pow((double)(Glu->n), -1.0 / options->ILU_MILU_Dim);    xlusup_first = xlusup[first];    xlsub_first = xlsub[first];    m = xlusup[first + 1] - xlusup_first;    n = last - first + 1;    m1 = m - 1;    inc_diag = m + 1;    nzlc = lastc ? (xlusup[last + 2] - xlusup[last + 1]) : 0;    temp = swork - n;    /* Quick return if nothing to do. */    if (m == 0 || m == n || drop_rule == NODROP)    {        *nnzLj += m * n;        return 0;    }    /* basic dropping: ILU(tau) */    for (i = n; i <= m1; )    {        /* the average abs value of ith row */        switch (nrm)        {            case ONE_NORM:                temp[i] = sasum_(&n, &lusup[xlusup_first + i], &m) / (double)n;                break;            case TWO_NORM:                temp[i] = snrm2_(&n, &lusup[xlusup_first + i], &m)                    / sqrt((double)n);                break;            case INF_NORM:            default:                k = isamax_(&n, &lusup[xlusup_first + i], &m) - 1;                temp[i] = fabs(lusup[xlusup_first + i + m * k]);                break;        }        /* drop small entries due to drop_tol */        if (drop_rule & DROP_BASIC && temp[i] < drop_tol)        {            r++;            /* drop the current row and move the last undropped row here */            if (r > 1) /* add to last row */            {                /* accumulate the sum (for MILU) */                switch (milu)                {                    case SMILU_1:                    case SMILU_2:                        saxpy_(&n, &one, &lusup[xlusup_first + i], &m,                                &lusup[xlusup_first + m - 1], &m);//.........这里部分代码省略.........

/* 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 | *///.........这里部分代码省略.........

void test08 ( void )/******************************************************************************//*  Purpose:    TEST08 demonstrates SNRM2.  Modified:    29 March 2007  Author:    John Burkardt*/{/*  These parameters illustrate the fact that matrices are typically  dimensioned with more space than the user requires.*/  float *a;  int i;  int inc;  int j;  int lda = 10;  int n = 5;  int ncopy;  float sum1;  float *x;  a = malloc ( lda * lda * sizeof ( float ) );  x = malloc ( n * sizeof ( float ) );  printf ( "/n" );  printf ( "TEST08/n" );  printf ( "  SNRM2 computes the Euclidean norm of a vector./n" );  printf ( "/n" );/*  Compute the euclidean norm of a vector:*/  for ( i = 0; i < n; i++ )  {    x[i] = ( float ) ( i + 1 );  }  printf ( "/n" );  printf ( "  X =/n" );  printf ( "/n" );  for ( i = 0; i < n; i++ )  {    printf ( "  %6d  %14d/n", i + 1, x[i] );  }  printf ( "/n" );  ncopy = n;  inc = 1;  printf ( "  The 2-norm of X is %f/n", snrm2_ ( &ncopy, x, &inc ) );/*  Compute the euclidean norm of a row or column of a matrix:*/  for ( i = 0; i < n; i++ )  {    for ( j = 0; j < n; j++ )    {      a[i+j*lda] = ( float ) ( i + 1 + j + 1 );    }  }  printf ( "/n" );  ncopy = n;  inc = lda;  printf ( "  The 2-norm of row 2 of A is %f/n",     snrm2_ ( &ncopy, a+1+0*lda, &inc ) );  printf ( "/n" );  ncopy = n;  inc = 1;  printf ( "  The 2-norm of column 2 of A is %f/n" ,       snrm2_ ( &ncopy, a+0+1*lda, &inc ) );  free ( a );  free ( x );   return;}

/* Subroutine */int sgebal_(char *job, integer *n, real *a, integer *lda, integer *ilo, integer *ihi, real *scale, integer *info){    /* System generated locals */    integer a_dim1, a_offset, i__1, i__2;    real r__1, r__2;    /* Local variables */    real c__, f, g;    integer i__, j, k, l, m;    real r__, s, ca, ra;    integer ica, ira, iexc;    extern real snrm2_(integer *, real *, integer *);    extern logical lsame_(char *, char *);    extern /* Subroutine */    int sscal_(integer *, real *, real *, integer *), sswap_(integer *, real *, integer *, real *, integer *);    real sfmin1, sfmin2, sfmax1, sfmax2;    extern real slamch_(char *);    extern /* Subroutine */    int xerbla_(char *, integer *);    extern integer isamax_(integer *, real *, integer *);    extern logical sisnan_(real *);    logical noconv;    /* -- LAPACK computational routine (version 3.5.0) -- */    /* -- LAPACK is a software package provided by Univ. of Tennessee, -- */    /* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- */    /* November 2013 */    /* .. 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;    --scale;    /* Function Body */    *info = 0;    if (! lsame_(job, "N") && ! lsame_(job, "P") && ! lsame_(job, "S") && ! lsame_(job, "B"))    {        *info = -1;    }    else if (*n < 0)    {        *info = -2;    }    else if (*lda < max(1,*n))    {        *info = -4;    }    if (*info != 0)    {        i__1 = -(*info);        xerbla_("SGEBAL", &i__1);        return 0;    }    k = 1;    l = *n;    if (*n == 0)    {        goto L210;    }    if (lsame_(job, "N"))    {        i__1 = *n;        for (i__ = 1;                i__ <= i__1;                ++i__)        {            scale[i__] = 1.f;            /* L10: */        }        goto L210;    }    if (lsame_(job, "S"))    {        goto L120;    }    /* Permutation to isolate eigenvalues if possible */    goto L50;    /* Row and column exchange. */L20:    scale[m] = (real) j;    if (j == m)    {        goto L30;    }    sswap_(&l, &a[j * a_dim1 + 1], &c__1, &a[m * a_dim1 + 1], &c__1);//.........这里部分代码省略.........

/* Subroutine */ int sqrt15_(integer *scale, integer *rksel, integer *m,                             integer *n, integer *nrhs, real *a, integer *lda, real *b, integer *                             ldb, real *s, integer *rank, real *norma, real *normb, integer *iseed,                             real *work, integer *lwork){    /* System generated locals */    integer a_dim1, a_offset, b_dim1, b_offset, i__1, i__2;    real r__1;    /* Local variables */    static integer info;    static real temp;    extern doublereal snrm2_(integer *, real *, integer *);    static integer j;    extern /* Subroutine */ int sscal_(integer *, real *, real *, integer *),           slarf_(char *, integer *, integer *, real *, integer *, real *,                  real *, integer *, real *), sgemm_(char *, char *,                          integer *, integer *, integer *, real *, real *, integer *, real *                          , integer *, real *, real *, integer *);    extern doublereal sasum_(integer *, real *, integer *);    static real dummy[1];    static integer mn;    extern doublereal slamch_(char *), slange_(char *, integer *,            integer *, real *, integer *, real *);    extern /* Subroutine */ int xerbla_(char *, integer *);    static real bignum;    extern /* Subroutine */ int slascl_(char *, integer *, integer *, real *,                                        real *, integer *, integer *, real *, integer *, integer *);    extern doublereal slarnd_(integer *, integer *);    extern /* Subroutine */ int slaord_(char *, integer *, real *, integer *), slaset_(char *, integer *, integer *, real *, real *,            real *, integer *), slaror_(char *, char *, integer *,                                        integer *, real *, integer *, integer *, real *, integer *), slarnv_(integer *, integer *, integer *, real *);    static real smlnum, eps;#define a_ref(a_1,a_2) a[(a_2)*a_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           September 30, 1994        Purpose        =======        SQRT15 generates a matrix with full or deficient rank and of various        norms.        Arguments        =========        SCALE   (input) INTEGER                SCALE = 1: normally scaled matrix                SCALE = 2: matrix scaled up                SCALE = 3: matrix scaled down        RKSEL   (input) INTEGER                RKSEL = 1: full rank matrix                RKSEL = 2: rank-deficient matrix        M       (input) INTEGER                The number of rows of the matrix A.        N       (input) INTEGER                The number of columns of A.        NRHS    (input) INTEGER                The number of columns of B.        A       (output) REAL array, dimension (LDA,N)                The M-by-N matrix A.        LDA     (input) INTEGER                The leading dimension of the array A.        B       (output) REAL array, dimension (LDB, NRHS)                A matrix that is in the range space of matrix A.        LDB     (input) INTEGER                The leading dimension of the array B.        S       (output) REAL array, dimension MIN(M,N)                Singular values of A.        RANK    (output) INTEGER                number of nonzero singular values of A.        NORMA   (output) REAL                one-norm of A.        NORMB   (output) REAL                one-norm of B.        ISEED   (input/output) integer array, dimension (4)                seed for random number generator.        WORK    (workspace) REAL array, dimension (LWORK)//.........这里部分代码省略.........

/* Subroutine */ int slasd3_(integer *nl, integer *nr, integer *sqre, integer 	*k, real *d__, real *q, integer *ldq, real *dsigma, real *u, integer *	ldu, real *u2, integer *ldu2, real *vt, integer *ldvt, real *vt2, 	integer *ldvt2, integer *idxc, integer *ctot, real *z__, integer *	info){    /* System generated locals */    integer q_dim1, q_offset, u_dim1, u_offset, u2_dim1, u2_offset, vt_dim1, 	    vt_offset, vt2_dim1, vt2_offset, i__1, i__2;    real r__1, r__2;    /* Builtin functions */    double sqrt(doublereal), r_sign(real *, real *);    /* Local variables */    integer i__, j, m, n, jc;    real rho;    integer nlp1, nlp2, nrp1;    real temp;    extern doublereal snrm2_(integer *, real *, integer *);    integer ctemp;    extern /* Subroutine */ int sgemm_(char *, char *, integer *, integer *, 	    integer *, real *, real *, integer *, real *, integer *, real *, 	    real *, integer *);    integer ktemp;    extern /* Subroutine */ int scopy_(integer *, real *, integer *, real *, 	    integer *);    extern doublereal slamc3_(real *, real *);    extern /* Subroutine */ int slasd4_(integer *, integer *, real *, real *, 	    real *, real *, real *, real *, integer *), xerbla_(char *, 	    integer *), slascl_(char *, integer *, integer *, real *, 	    real *, integer *, integer *, real *, integer *, integer *), slacpy_(char *, integer *, integer *, real *, integer *, 	    real *, integer *);/*  -- LAPACK auxiliary routine (version 3.2) -- *//*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. *//*     November 2006 *//*     .. Scalar Arguments .. *//*     .. *//*     .. Array Arguments .. *//*     .. *//*  Purpose *//*  ======= *//*  SLASD3 finds all the square roots of the roots of the secular *//*  equation, as defined by the values in D and Z.  It makes the *//*  appropriate calls to SLASD4 and then updates the singular *//*  vectors by matrix multiplication. *//*  This code makes very mild assumptions about floating point *//*  arithmetic. It will work on machines with a guard digit in *//*  add/subtract, or on those binary machines without guard digits *//*  which subtract like the Cray XMP, Cray YMP, Cray C 90, or Cray 2. *//*  It could conceivably fail on hexadecimal or decimal machines *//*  without guard digits, but we know of none. *//*  SLASD3 is called from SLASD1. *//*  Arguments *//*  ========= *//*  NL     (input) INTEGER *//*         The row dimension of the upper block.  NL >= 1. *//*  NR     (input) INTEGER *//*         The row dimension of the lower block.  NR >= 1. *//*  SQRE   (input) INTEGER *//*         = 0: the lower block is an NR-by-NR square matrix. *//*         = 1: the lower block is an NR-by-(NR+1) rectangular matrix. *//*         The bidiagonal matrix has N = NL + NR + 1 rows and *//*         M = N + SQRE >= N columns. *//*  K      (input) INTEGER *//*         The size of the secular equation, 1 =< K = < N. *//*  D      (output) REAL array, dimension(K) *//*         On exit the square roots of the roots of the secular equation, *//*         in ascending order. *//*  Q      (workspace) REAL array, *//*                     dimension at least (LDQ,K). *//*  LDQ    (input) INTEGER *//*         The leading dimension of the array Q.  LDQ >= K. *//*  DSIGMA (input/output) REAL array, dimension(K) *//*         The first K elements of this array contain the old roots *//*         of the deflated updating problem.  These are the poles *//*         of the secular equation. *//*  U      (output) REAL array, dimension (LDU, N) *//*         The last N - K columns of this matrix contain the deflated *//*         left singular vectors. *//*  LDU    (input) INTEGER *///.........这里部分代码省略.........

/* ----------------------------------------------------------------------- *//* Subroutine */ int sneupd_(logical *rvec, char *howmny, logical *select, 	real *dr, real *di, real *z__, integer *ldz, real *sigmar, real *	sigmai, real *workev, char *bmat, integer *n, char *which, integer *	nev, real *tol, real *resid, integer *ncv, real *v, integer *ldv, 	integer *iparam, integer *ipntr, real *workd, real *workl, integer *	lworkl, integer *info, ftnlen howmny_len, ftnlen bmat_len, ftnlen 	which_len){    /* System generated locals */    integer v_dim1, v_offset, z_dim1, z_offset, i__1;    real r__1, r__2;    doublereal d__1;    /* Local variables */    static integer j, k, ih, jj, np;    static real vl[1]	/* was [1][1] */;    static integer ibd, ldh, ldq, iri;    static real sep;    static integer irr, wri, wrr, mode;    static real eps23;    extern /* Subroutine */ int sger_(integer *, integer *, real *, real *, 	    integer *, real *, integer *, real *, integer *);    static integer ierr;    static real temp;    static integer iwev;    static char type__[6];    static real temp1;    extern doublereal snrm2_(integer *, real *, integer *);    static integer ihbds, iconj;    extern /* Subroutine */ int sscal_(integer *, real *, real *, integer *);    static real conds;    static logical reord;    extern /* Subroutine */ int sgemv_(char *, integer *, integer *, real *, 	    real *, integer *, real *, integer *, real *, real *, integer *, 	    ftnlen);    static integer nconv, iwork[1];    static real rnorm;    extern /* Subroutine */ int scopy_(integer *, real *, integer *, real *, 	    integer *);    static integer ritzi;    extern /* Subroutine */ int strmm_(char *, char *, char *, char *, 	    integer *, integer *, real *, real *, integer *, real *, integer *	    , ftnlen, ftnlen, ftnlen, ftnlen), ivout_(integer *, integer *, 	    integer *, integer *, char *, ftnlen), smout_(integer *, integer *	    , integer *, real *, integer *, integer *, char *, ftnlen);    static integer ritzr;    extern /* Subroutine */ int svout_(integer *, integer *, real *, integer *	    , char *, ftnlen), sgeqr2_(integer *, integer *, real *, integer *	    , real *, real *, integer *);    static integer nconv2;    extern doublereal slapy2_(real *, real *);    extern /* Subroutine */ int sorm2r_(char *, char *, integer *, integer *, 	    integer *, real *, integer *, real *, real *, integer *, real *, 	    integer *, ftnlen, ftnlen);    static integer iheigi, iheigr, bounds, invsub, iuptri, msglvl, outncv, 	    ishift, numcnv;    extern /* Subroutine */ int slacpy_(char *, integer *, integer *, real *, 	    integer *, real *, integer *, ftnlen), slahqr_(logical *, logical 	    *, integer *, integer *, integer *, real *, integer *, real *, 	    real *, integer *, integer *, real *, integer *, integer *), 	    slaset_(char *, integer *, integer *, real *, real *, real *, 	    integer *, ftnlen), strevc_(char *, char *, logical *, integer *, 	    real *, integer *, real *, integer *, real *, integer *, integer *	    , integer *, real *, integer *, ftnlen, ftnlen), strsen_(char *, 	    char *, logical *, integer *, real *, integer *, real *, integer *	    , real *, real *, integer *, real *, real *, real *, integer *, 	    integer *, integer *, integer *, ftnlen, ftnlen);    extern doublereal slamch_(char *, ftnlen);    extern /* Subroutine */ int sngets_(integer *, char *, integer *, integer 	    *, real *, real *, real *, real *, 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 | *///.........这里部分代码省略.........

/* Subroutine */int slaein_(logical *rightv, logical *noinit, integer *n, real *h__, integer *ldh, real *wr, real *wi, real *vr, real *vi, real *b, integer *ldb, real *work, real *eps3, real *smlnum, real *bignum, integer *info){    /* System generated locals */    integer b_dim1, b_offset, h_dim1, h_offset, i__1, i__2, i__3, i__4;    real r__1, r__2, r__3, r__4;    /* Builtin functions */    double sqrt(doublereal);    /* Local variables */    integer i__, j;    real w, x, y;    integer i1, i2, i3;    real w1, ei, ej, xi, xr, rec;    integer its, ierr;    real temp, norm, vmax;    extern real snrm2_(integer *, real *, integer *);    real scale;    extern /* Subroutine */    int sscal_(integer *, real *, real *, integer *);    char trans[1];    real vcrit;    extern real sasum_(integer *, real *, integer *);    real rootn, vnorm;    extern real slapy2_(real *, real *);    real absbii, absbjj;    extern integer isamax_(integer *, real *, integer *);    extern /* Subroutine */    int sladiv_(real *, real *, real *, real *, real * , real *);    char normin[1];    real nrmsml;    extern /* Subroutine */    int slatrs_(char *, char *, char *, char *, integer *, real *, integer *, real *, real *, real *, integer *);    real growto;    /* -- 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 */    h_dim1 = *ldh;    h_offset = 1 + h_dim1;    h__ -= h_offset;    --vr;    --vi;    b_dim1 = *ldb;    b_offset = 1 + b_dim1;    b -= b_offset;    --work;    /* Function Body */    *info = 0;    /* GROWTO is the threshold used in the acceptance test for an */    /* eigenvector. */    rootn = sqrt((real) (*n));    growto = .1f / rootn;    /* Computing MAX */    r__1 = 1.f;    r__2 = *eps3 * rootn; // , expr subst    nrmsml = max(r__1,r__2) * *smlnum;    /* Form B = H - (WR,WI)*I (except that the subdiagonal elements and */    /* the imaginary parts of the diagonal elements are not stored). */    i__1 = *n;    for (j = 1;            j <= i__1;            ++j)    {        i__2 = j - 1;        for (i__ = 1;                i__ <= i__2;                ++i__)        {            b[i__ + j * b_dim1] = h__[i__ + j * h_dim1];            /* L10: */        }        b[j + j * b_dim1] = h__[j + j * h_dim1] - *wr;        /* L20: */    }    if (*wi == 0.f)    {        /* Real eigenvalue. */        if (*noinit)        {            /* Set initial vector. */            i__1 = *n;            for (i__ = 1;                    i__ <= i__1;//.........这里部分代码省略.........

/* Subroutine */ int slaqp2_(integer *m, integer *n, integer *offset, real *a,	 integer *lda, integer *jpvt, real *tau, real *vn1, real *vn2, real *	work){    /* System generated locals */    integer a_dim1, a_offset, i__1, i__2, i__3;    real r__1, r__2;    /* Builtin functions */    double sqrt(doublereal);    /* Local variables */    static integer i__, j, mn;    static real aii;    static integer pvt;    static real temp, temp2;    extern doublereal snrm2_(integer *, real *, integer *);    static integer offpi;    extern /* Subroutine */ int slarf_(char *, integer *, integer *, real *, 	    integer *, real *, real *, integer *, real *, ftnlen);    static integer itemp;    extern /* Subroutine */ int sswap_(integer *, real *, integer *, real *, 	    integer *), slarfg_(integer *, real *, real *, integer *, real *);    extern integer isamax_(integer *, real *, integer *);/*  -- LAPACK auxiliary routine (version 3.0) -- *//*     Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., *//*     Courant Institute, Argonne National Lab, and Rice University *//*     June 30, 1999 *//*     .. Scalar Arguments .. *//*     .. *//*     .. Array Arguments .. *//*     .. *//*  Purpose *//*  ======= *//*  SLAQP2 computes a QR factorization with column pivoting of *//*  the block A(OFFSET+1:M,1:N). *//*  The block A(1:OFFSET,1:N) is accordingly pivoted, but not factorized. *//*  Arguments *//*  ========= *//*  M       (input) INTEGER *//*          The number of rows of the matrix A. M >= 0. *//*  N       (input) INTEGER *//*          The number of columns of the matrix A. N >= 0. *//*  OFFSET  (input) INTEGER *//*          The number of rows of the matrix A that must be pivoted *//*          but no factorized. OFFSET >= 0. *//*  A       (input/output) REAL array, dimension (LDA,N) *//*          On entry, the M-by-N matrix A. *//*          On exit, the upper triangle of block A(OFFSET+1:M,1:N) is *//*          the triangular factor obtained; the elements in block *//*          A(OFFSET+1:M,1:N) below the diagonal, together with the *//*          array TAU, represent the orthogonal matrix Q as a product of *//*          elementary reflectors. Block A(1:OFFSET,1:N) has been *//*          accordingly pivoted, but no factorized. *//*  LDA     (input) INTEGER *//*          The leading dimension of the array A. LDA >= max(1,M). *//*  JPVT    (input/output) INTEGER array, dimension (N) *//*          On entry, if JPVT(i) .ne. 0, the i-th column of A is permuted *//*          to the front of A*P (a leading column); if JPVT(i) = 0, *//*          the i-th column of A is a free column. *//*          On exit, if JPVT(i) = k, then the i-th column of A*P *//*          was the k-th column of A. *//*  TAU     (output) REAL array, dimension (min(M,N)) *//*          The scalar factors of the elementary reflectors. *//*  VN1     (input/output) REAL array, dimension (N) *//*          The vector with the partial column norms. *//*  VN2     (input/output) REAL array, dimension (N) *//*          The vector with the exact column norms. *//*  WORK    (workspace) REAL array, dimension (N) *//*  Further Details *//*  =============== *//*  Based on contributions by *//*    G. Quintana-Orti, Depto. de Informatica, Universidad Jaime I, Spain *//*    X. Sun, Computer Science Dept., Duke University, USA *//*  ===================================================================== *//*     .. Parameters .. *//*     .. *//*     .. Local Scalars .. *//*     .. *//*     .. External Subroutines .. *///.........这里部分代码省略.........

/*<       SUBROUTINE SLARFG( N, ALPHA, X, INCX, TAU ) >*//* Subroutine */ int slarfg_(integer *n, real *alpha, real *x, integer *incx,        real *tau){    /* System generated locals */    integer i__1;    real r__1;    /* Builtin functions */    double r_sign(real *, real *);    /* Local variables */    integer j, knt;    real beta;    extern doublereal snrm2_(integer *, real *, integer *);    extern /* Subroutine */ int sscal_(integer *, real *, real *, integer *);    real xnorm;    extern doublereal slapy2_(real *, real *), slamch_(char *, ftnlen);    real safmin, rsafmn;/*  -- LAPACK auxiliary routine (version 3.0) -- *//*     Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., *//*     Courant Institute, Argonne National Lab, and Rice University *//*     September 30, 1994 *//*     .. Scalar Arguments .. *//*<       INTEGER            INCX, N >*//*<       REAL               ALPHA, TAU >*//*     .. *//*     .. Array Arguments .. *//*<       REAL               X( * ) >*//*     .. *//*  Purpose *//*  ======= *//*  SLARFG generates a real elementary reflector H of order n, such *//*  that *//*        H * ( alpha ) = ( beta ),   H' * H = I. *//*            (   x   )   (   0  ) *//*  where alpha and beta are scalars, and x is an (n-1)-element real *//*  vector. H is represented in the form *//*        H = I - tau * ( 1 ) * ( 1 v' ) , *//*                      ( v ) *//*  where tau is a real scalar and v is a real (n-1)-element *//*  vector. *//*  If the elements of x are all zero, then tau = 0 and H is taken to be *//*  the unit matrix. *//*  Otherwise  1 <= tau <= 2. *//*  Arguments *//*  ========= *//*  N       (input) INTEGER *//*          The order of the elementary reflector. *//*  ALPHA   (input/output) REAL *//*          On entry, the value alpha. *//*          On exit, it is overwritten with the value beta. *//*  X       (input/output) REAL array, dimension *//*                         (1+(N-2)*abs(INCX)) *//*          On entry, the vector x. *//*          On exit, it is overwritten with the vector v. *//*  INCX    (input) INTEGER *//*          The increment between elements of X. INCX > 0. *//*  TAU     (output) REAL *//*          The value tau. *//*  ===================================================================== *//*     .. Parameters .. *//*<       REAL               ONE, ZERO >*//*<       PARAMETER          ( ONE = 1.0E+0, ZERO = 0.0E+0 ) >*//*     .. *//*     .. Local Scalars .. *//*<       INTEGER            J, KNT >*//*<       REAL               BETA, RSAFMN, SAFMIN, XNORM >*//*     .. *//*     .. External Functions .. *//*<       REAL               SLAMCH, SLAPY2, SNRM2 >*//*<       EXTERNAL           SLAMCH, SLAPY2, SNRM2 >*//*     .. *//*     .. Intrinsic Functions .. *//*<       INTRINSIC          ABS, SIGN >*//*     .. *//*     .. External Subroutines .. *//*<       EXTERNAL           SSCAL >*//*     .. *//*     .. Executable Statements .. *///.........这里部分代码省略.........

/* Subroutine */ int slaed3_(integer *k, integer *n, integer *n1, real *d__, 	real *q, integer *ldq, real *rho, real *dlamda, real *q2, integer *	indx, integer *ctot, real *w, real *s, integer *info){    /* System generated locals */    integer q_dim1, q_offset, i__1, i__2;    real r__1;    /* Builtin functions */    double sqrt(doublereal), r_sign(real *, real *);    /* Local variables */    static real temp;    extern doublereal snrm2_(integer *, real *, integer *);    static integer i__, j;    extern /* Subroutine */ int sgemm_(char *, char *, integer *, integer *, 	    integer *, real *, real *, integer *, real *, integer *, real *, 	    real *, integer *), scopy_(integer *, real *, 	    integer *, real *, integer *);    static integer n2;    extern /* Subroutine */ int slaed4_(integer *, integer *, real *, real *, 	    real *, real *, real *, integer *);    extern doublereal slamc3_(real *, real *);    static integer n12, ii, n23;    extern /* Subroutine */ int xerbla_(char *, integer *), slacpy_(	    char *, integer *, integer *, real *, integer *, real *, integer *	    ), slaset_(char *, integer *, integer *, real *, real *, 	    real *, integer *);    static integer iq2;#define q_ref(a_1,a_2) q[(a_2)*q_dim1 + a_1]/*  -- LAPACK routine (instrumented to count operations, version 3.0) --          Univ. of Tennessee, Oak Ridge National Lab, Argonne National Lab,          Courant Institute, NAG Ltd., and Rice University          June 30, 1999          Common block to return operation count and iteration count          ITCNT is unchanged, OPS is only incremented       Purpose       =======       SLAED3 finds the roots of the secular equation, as defined by the       values in D, W, and RHO, between 1 and K.  It makes the       appropriate calls to SLAED4 and then updates the eigenvectors by       multiplying the matrix of eigenvectors of the pair of eigensystems       being combined by the matrix of eigenvectors of the K-by-K system       which is solved here.       This code makes very mild assumptions about floating point       arithmetic. It will work on machines with a guard digit in       add/subtract, or on those binary machines without guard digits       which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or Cray-2.       It could conceivably fail on hexadecimal or decimal machines       without guard digits, but we know of none.       Arguments       =========       K       (input) INTEGER               The number of terms in the rational function to be solved by               SLAED4.  K >= 0.       N       (input) INTEGER               The number of rows and columns in the Q matrix.               N >= K (deflation may result in N>K).       N1      (input) INTEGER               The location of the last eigenvalue in the leading submatrix.               min(1,N) <= N1 <= N/2.       D       (output) REAL array, dimension (N)               D(I) contains the updated eigenvalues for               1 <= I <= K.       Q       (output) REAL array, dimension (LDQ,N)               Initially the first K columns are used as workspace.               On output the columns 1 to K contain               the updated eigenvectors.       LDQ     (input) INTEGER               The leading dimension of the array Q.  LDQ >= max(1,N).       RHO     (input) REAL               The value of the parameter in the rank one update equation.               RHO >= 0 required.       DLAMDA  (input/output) REAL array, dimension (K)               The first K elements of this array contain the old roots               of the deflated updating problem.  These are the poles               of the secular equation. May be changed on output by               having lowest order bit set to zero on Cray X-MP, Cray Y-MP,               Cray-2, or Cray C-90, as described above.       Q2      (input) REAL array, dimension (LDQ2, N)               The first K columns of this matrix contain the non-deflated               eigenvectors for the split problem.   //.........这里部分代码省略.........

//.........这里部分代码省略.........		    dsigjp = -poles[j + 1 + (poles_dim1 << 1)];		}		if (z__[j] == 0.f || poles[j + (poles_dim1 << 1)] == 0.f) {		    work[j] = 0.f;		} else {		    work[j] = -poles[j + (poles_dim1 << 1)] * z__[j] / diflj /			     (poles[j + (poles_dim1 << 1)] + dj);		}		i__2 = j - 1;		for (i__ = 1; i__ <= i__2; ++i__) {		    if (z__[i__] == 0.f || poles[i__ + (poles_dim1 << 1)] == 			    0.f) {			work[i__] = 0.f;		    } else {			work[i__] = poles[i__ + (poles_dim1 << 1)] * z__[i__] 				/ (slamc3_(&poles[i__ + (poles_dim1 << 1)], &				dsigj) - diflj) / (poles[i__ + (poles_dim1 << 				1)] + dj);		    }		}		i__2 = *k;		for (i__ = j + 1; i__ <= i__2; ++i__) {		    if (z__[i__] == 0.f || poles[i__ + (poles_dim1 << 1)] == 			    0.f) {			work[i__] = 0.f;		    } else {			work[i__] = poles[i__ + (poles_dim1 << 1)] * z__[i__] 				/ (slamc3_(&poles[i__ + (poles_dim1 << 1)], &				dsigjp) + difrj) / (poles[i__ + (poles_dim1 <<				 1)] + dj);		    }		}		work[1] = -1.f;		temp = snrm2_(k, &work[1], &c__1);		sgemv_("T", k, nrhs, &c_b11, &bx[bx_offset], ldbx, &work[1], &			c__1, &c_b13, &b[j + b_dim1], ldb);		slascl_("G", &c__0, &c__0, &temp, &c_b11, &c__1, nrhs, &b[j + 			b_dim1], ldb, info);	    }	}/*        Move the deflated rows of BX to B also. */	if (*k < max(m,n)) {	    i__1 = n - *k;	    slacpy_("A", &i__1, nrhs, &bx[*k + 1 + bx_dim1], ldbx, &b[*k + 1 		    + b_dim1], ldb);	}    } else {/*        Apply back the right orthogonal transformations. *//*        Step (1R): apply back the new right singular vector matrix *//*        to B. */	if (*k == 1) {	    scopy_(nrhs, &b[b_offset], ldb, &bx[bx_offset], ldbx);	} else {	    i__1 = *k;	    for (j = 1; j <= i__1; ++j) {		dsigj = poles[j + (poles_dim1 << 1)];		if (z__[j] == 0.f) {		    work[j] = 0.f;		} else {		    work[j] = -z__[j] / difl[j] / (dsigj + poles[j + 			    poles_dim1]) / difr[j + (difr_dim1 << 1)];

GURLS_EXPORT float nrm2(const int N, const float* X, const int incX){    return snrm2_(const_cast<int*>(&N), const_cast<float*>(X), const_cast<int*>(&incX));}

/* Subroutine */ int slarfg_(integer *n, real *alpha, real *x, integer *incx, 	real *tau){    /* System generated locals */    integer i__1;    real r__1;    /* Builtin functions */    double r_sign(real *, real *);    /* Local variables */    static integer j, knt;    static real beta;    extern doublereal snrm2_(integer *, real *, integer *);    extern /* Subroutine */ int sscal_(integer *, real *, real *, integer *);    static real xnorm;    extern doublereal slapy2_(real *, real *), slamch_(char *, ftnlen);    static real safmin, rsafmn;/*  -- LAPACK auxiliary routine (version 3.0) -- *//*     Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., *//*     Courant Institute, Argonne National Lab, and Rice University *//*     September 30, 1994 *//*     .. Scalar Arguments .. *//*     .. *//*     .. Array Arguments .. *//*     .. *//*  Purpose *//*  ======= *//*  SLARFG generates a real elementary reflector H of order n, such *//*  that *//*        H * ( alpha ) = ( beta ),   H' * H = I. *//*            (   x   )   (   0  ) *//*  where alpha and beta are scalars, and x is an (n-1)-element real *//*  vector. H is represented in the form *//*        H = I - tau * ( 1 ) * ( 1 v' ) , *//*                      ( v ) *//*  where tau is a real scalar and v is a real (n-1)-element *//*  vector. *//*  If the elements of x are all zero, then tau = 0 and H is taken to be *//*  the unit matrix. *//*  Otherwise  1 <= tau <= 2. *//*  Arguments *//*  ========= *//*  N       (input) INTEGER *//*          The order of the elementary reflector. *//*  ALPHA   (input/output) REAL *//*          On entry, the value alpha. *//*          On exit, it is overwritten with the value beta. *//*  X       (input/output) REAL array, dimension *//*                         (1+(N-2)*abs(INCX)) *//*          On entry, the vector x. *//*          On exit, it is overwritten with the vector v. *//*  INCX    (input) INTEGER *//*          The increment between elements of X. INCX > 0. *//*  TAU     (output) REAL *//*          The value tau. *//*  ===================================================================== *//*     .. Parameters .. *//*     .. *//*     .. Local Scalars .. *//*     .. *//*     .. External Functions .. *//*     .. *//*     .. Intrinsic Functions .. *//*     .. *//*     .. External Subroutines .. *//*     .. *//*     .. Executable Statements .. */    /* Parameter adjustments */    --x;    /* Function Body */    if (*n <= 1) {	*tau = 0.f;	return 0;    }    i__1 = *n - 1;    xnorm = snrm2_(&i__1, &x[1], incx);//.........这里部分代码省略.........

/* Subroutine */ int cstein_(integer *n, real *d__, real *e, integer *m, real 	*w, integer *iblock, integer *isplit, complex *z__, integer *ldz, 	real *work, integer *iwork, integer *ifail, integer *info){    /* System generated locals */    integer z_dim1, z_offset, i__1, i__2, i__3, i__4, i__5;    real r__1, r__2, r__3, r__4, r__5;    complex q__1;    /* Builtin functions */    double sqrt(doublereal);    /* Local variables */    integer i__, j, b1, j1, bn, jr;    real xj, scl, eps, ctr, sep, nrm, tol;    integer its;    real xjm, eps1;    integer jblk, nblk, jmax;    extern doublereal snrm2_(integer *, real *, integer *);    integer iseed[4], gpind, iinfo;    extern /* Subroutine */ int sscal_(integer *, real *, real *, integer *);    extern doublereal sasum_(integer *, real *, integer *);    extern /* Subroutine */ int scopy_(integer *, real *, integer *, real *, 	    integer *);    real ortol;    integer indrv1, indrv2, indrv3, indrv4, indrv5;    extern doublereal slamch_(char *);    extern /* Subroutine */ int xerbla_(char *, integer *), slagtf_(	    integer *, real *, real *, real *, real *, real *, real *, 	    integer *, integer *);    integer nrmchk;    extern integer isamax_(integer *, real *, integer *);    extern /* Subroutine */ int slagts_(integer *, integer *, real *, real *, 	    real *, real *, integer *, real *, real *, integer *);    integer blksiz;    real onenrm, pertol;    extern /* Subroutine */ int slarnv_(integer *, integer *, integer *, real 	    *);    real stpcrt;/*  -- LAPACK routine (version 3.1) -- *//*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. *//*     November 2006 *//*     .. Scalar Arguments .. *//*     .. *//*     .. Array Arguments .. *//*     .. *//*  Purpose *//*  ======= *//*  CSTEIN computes the eigenvectors of a real symmetric tridiagonal *//*  matrix T corresponding to specified eigenvalues, using inverse *//*  iteration. *//*  The maximum number of iterations allowed for each eigenvector is *//*  specified by an internal parameter MAXITS (currently set to 5). *//*  Although the eigenvectors are real, they are stored in a complex *//*  array, which may be passed to CUNMTR or CUPMTR for back *//*  transformation to the eigenvectors of a complex Hermitian matrix *//*  which was reduced to tridiagonal form. *//*  Arguments *//*  ========= *//*  N       (input) INTEGER *//*          The order of the matrix.  N >= 0. *//*  D       (input) REAL array, dimension (N) *//*          The n diagonal elements of the tridiagonal matrix T. *//*  E       (input) REAL array, dimension (N-1) *//*          The (n-1) subdiagonal elements of the tridiagonal matrix *//*          T, stored in elements 1 to N-1. *//*  M       (input) INTEGER *//*          The number of eigenvectors to be found.  0 <= M <= N. *//*  W       (input) REAL array, dimension (N) *//*          The first M elements of W contain the eigenvalues for *//*          which eigenvectors are to be computed.  The eigenvalues *//*          should be grouped by split-off block and ordered from *//*          smallest to largest within the block.  ( The output array *//*          W from SSTEBZ with ORDER = 'B' is expected here. ) *//*  IBLOCK  (input) INTEGER array, dimension (N) *//*          The submatrix indices associated with the corresponding *//*          eigenvalues in W; IBLOCK(i)=1 if eigenvalue W(i) belongs to *//*          the first submatrix from the top, =2 if W(i) belongs to *//*          the second submatrix, etc.  ( The output array IBLOCK *//*          from SSTEBZ is expected here. ) *//*  ISPLIT  (input) INTEGER array, dimension (N) *//*          The splitting points, at which T breaks up into submatrices. *//*          The first submatrix consists of rows/columns 1 to *//*          ISPLIT( 1 ), the second of rows/columns ISPLIT( 1 )+1 *///.........这里部分代码省略.........

//.........这里部分代码省略.........    printf("L//U MB %.3f/ttotal MB needed %.3f/n",	   mem_usage.for_lu/1e6, mem_usage.total_needed/1e6);    fflush(stdout);    /* Set the global variables. */    GLOBAL_A = &A;    GLOBAL_L = &L;    GLOBAL_U = &U;    GLOBAL_STAT = &stat;    GLOBAL_PERM_C = perm_c;    GLOBAL_PERM_R = perm_r;    GLOBAL_OPTIONS = &options;    GLOBAL_R = R;    GLOBAL_C = C;    GLOBAL_MEM_USAGE = &mem_usage;    /* Set the options to do solve-only. */    options.Fact = FACTORED;    options.PivotGrowth = NO;    options.ConditionNumber = NO;    /* Set the variables used by GMRES. */    restrt = SUPERLU_MIN(n / 3 + 1, 50);    maxit = 1000;    iter = maxit;    resid = 1e-8;    if (!(x = floatMalloc(n))) ABORT("Malloc fails for x[].");    if (info <= n + 1)    {	int i_1 = 1;	double maxferr = 0.0, nrmA, nrmB, res, t;        float temp;	extern float snrm2_(int *, float [], int *);	extern void saxpy_(int *, float *, float [], int *, float [], int *);	/* Initial guess */	for (i = 0; i < n; i++) x[i] = zero;	t = SuperLU_timer_();	/* Call GMRES */	sfgmr(n, smatvec_mult, spsolve, b, x, resid, restrt, &iter, stdout);	t = SuperLU_timer_() - t;	/* Output the result. */	nrmA = snrm2_(&(Astore->nnz), (float *)((DNformat *)A.Store)->nzval,		&i_1);	nrmB = snrm2_(&m, b, &i_1);	sp_sgemv("N", -1.0, &A, x, 1, 1.0, b, 1);	res = snrm2_(&m, b, &i_1);	resid = res / nrmB;	printf("||A||_F = %.1e, ||B||_2 = %.1e, ||B-A*X||_2 = %.1e, "		"relres = %.1e/n", nrmA, nrmB, res, resid);	if (iter >= maxit)	{	    if (resid >= 1.0) iter = -180;	    else if (resid > 1e-8) iter = -111;	}	printf("iteration: %d/nresidual: %.1e/nGMRES time: %.2f seconds./n",		iter, resid, t);	/* Scale the solution back if equilibration was performed. */	if (*equed == 'C' || *equed == 'B') 

/* Subroutine */ int stgsna_(char *job, char *howmny, logical *select, 	integer *n, real *a, integer *lda, real *b, integer *ldb, real *vl, 	integer *ldvl, real *vr, integer *ldvr, real *s, real *dif, integer *	mm, integer *m, real *work, integer *lwork, integer *iwork, integer *	info){    /* System generated locals */    integer a_dim1, a_offset, b_dim1, b_offset, vl_dim1, vl_offset, vr_dim1, 	    vr_offset, i__1, i__2;    real r__1, r__2;    /* Builtin functions */    double sqrt(doublereal);    /* Local variables */    integer i__, k;    real c1, c2;    integer n1, n2, ks, iz;    real eps, beta, cond;    logical pair;    integer ierr;    real uhav, uhbv;    integer ifst;    real lnrm;    extern doublereal sdot_(integer *, real *, integer *, real *, integer *);    integer ilst;    real rnrm;    extern /* Subroutine */ int slag2_(real *, integer *, real *, integer *, 	    real *, real *, real *, real *, real *, real *);    extern doublereal snrm2_(integer *, real *, integer *);    real root1, root2, scale;    extern logical lsame_(char *, char *);    real uhavi, uhbvi;    extern /* Subroutine */ int sgemv_(char *, integer *, integer *, real *, 	    real *, integer *, real *, integer *, real *, real *, integer *);    real tmpii;    integer lwmin;    logical wants;    real tmpir, tmpri, dummy[1], tmprr;    extern doublereal slapy2_(real *, real *);    real dummy1[1], alphai, alphar;    extern doublereal slamch_(char *);    extern /* Subroutine */ int xerbla_(char *, integer *);    logical wantbh, wantdf;    extern /* Subroutine */ int slacpy_(char *, integer *, integer *, real *, 	    integer *, real *, integer *), stgexc_(logical *, logical 	    *, integer *, real *, integer *, real *, integer *, real *, 	    integer *, real *, integer *, integer *, integer *, real *, 	    integer *, integer *);    logical somcon;    real alprqt, smlnum;    logical lquery;    extern /* Subroutine */ int stgsyl_(char *, integer *, integer *, integer 	    *, real *, integer *, real *, integer *, real *, integer *, real *, integer *, real *, integer *, real *, integer *, real *, real *, 	     real *, integer *, integer *, integer *);/*  -- LAPACK routine (version 3.1) -- *//*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. *//*     November 2006 *//*     .. Scalar Arguments .. *//*     .. *//*     .. Array Arguments .. *//*     .. *//*  Purpose *//*  ======= *//*  STGSNA estimates reciprocal condition numbers for specified *//*  eigenvalues and/or eigenvectors of a matrix pair (A, B) in *//*  generalized real Schur canonical form (or of any matrix pair *//*  (Q*A*Z', Q*B*Z') with orthogonal matrices Q and Z, where *//*  Z' denotes the transpose of Z. *//*  (A, B) must be in generalized real Schur form (as returned by SGGES), *//*  i.e. A is block upper triangular with 1-by-1 and 2-by-2 diagonal *//*  blocks. B is upper triangular. *//*  Arguments *//*  ========= *//*  JOB     (input) CHARACTER*1 *//*          Specifies whether condition numbers are required for *//*          eigenvalues (S) or eigenvectors (DIF): *//*          = 'E': for eigenvalues only (S); *//*          = 'V': for eigenvectors only (DIF); *//*          = 'B': for both eigenvalues and eigenvectors (S and DIF). *//*  HOWMNY  (input) CHARACTER*1 *//*          = 'A': compute condition numbers for all eigenpairs; *//*          = 'S': compute condition numbers for selected eigenpairs *//*                 specified by the array SELECT. *//*  SELECT  (input) LOGICAL array, dimension (N) *//*          If HOWMNY = 'S', SELECT specifies the eigenpairs for which *//*          condition numbers are required. To select condition numbers *//*          for the eigenpair corresponding to a real eigenvalue w(j), *///.........这里部分代码省略.........