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;
}