double ProtoMol::Lapack::dnrm2(int *n, double *x, int *incx) {
  FAHCheckIn();
#if defined(HAVE_LAPACK)
  return dnrm2_(n, x, incx);
#elif defined(HAVE_SIMTK_LAPACK)
  return dnrm2_(*n, x, *incx);
#elif defined(HAVE_MKL_LAPACK)
  return DNRM2(n, x, incx);
#else
  THROW(std::string(__func__) + " not supported");
#endif
}

void dngeres(int *n,        int *nev,   int *aptr, int *aind, 
             double *avals, int *bptr,  int *bind, double *bvals,
             double *dr,    double *di, double *z, int *ldz,
             double *res)
{
   int     i, first=1, ione=1, j;
   double  md, rnrm, *ax, *bx;
 
   ax = (double*)malloc((*n)*sizeof(double)); 
   bx = (double*)malloc((*n)*sizeof(double)); 

   for (i = 0; i<(*nev); i++) {
      if (di[i] == 0.0) {
         dmvm_(n, avals, aind, aptr, &z(1,i+1), ax, &ione);
         dmvm_(n, bvals, bind, bptr, &z(1,i+1), bx, &ione);
         md = -dr[i];
         daxpy_(n, &md, bx, &ione, ax, &ione);
         res[i] = dnrm2_(n, ax, &ione);
         res[i] = res[i]/abs(dr[i]); 
      }
      else if (first) {
         dmvm_(n, avals, aind, aptr, &z(1,i+1), ax, &ione);
         dmvm_(n, bvals, bind, bptr, &z(1,i+1), bx, &ione);
         md = -dr[i];
         daxpy_(n, &md, bx, &ione, ax, &ione);
         dmvm_(n, bvals, bind, bptr, &z(1,i+2), bx, &ione);
         daxpy_(n, &di[i], bx, &ione, ax, &ione);
         rnrm = dnrm2_(n, ax, &ione);
         res[i] = rnrm*rnrm;

         dmvm_(n, avals, aind, aptr, &z(1,i+2), ax, &ione);
         dmvm_(n, bvals, bind, bptr, &z(1,i+2), bx, &ione);
         md = -dr[i];
         daxpy_(n, &md, bx, &ione, ax, &ione);
         dmvm_(n, bvals, bind, bptr, &z(1,i+1), bx, &ione);
         md = -di[i];
         daxpy_(n, &md, bx, &ione, ax, &ione);
         rnrm = dnrm2_(n, ax, &ione);
         res[i] = dlapy2_(&res[i], &rnrm);
         res[i] = res[i]/dlapy2_(&dr[i],&di[i]);  
         res[i+1] = res[i];
         first = 0; 
      }
      else {
         first = 1;
      }
   }
   free(ax);
   free(bx);
}

/*---------------------------------------------------------
  fem_vec_norm - to compute a norm of global vector (in parallel)
---------------------------------------------------------*/
double fem_vec_norm( /* returns: L2 norm of global Vector */
  int Solver_id,        /* in: solver data structure to be used */
  int Level_id,         /* in: level number */
  int Nrdof,            /* in: number of vector components */
  double* Vector        /* in: local part of global Vector */
  )
{

  int IONE=1; int problem_id; double vec_norm;

/*++++++++++++++++ executable statements ++++++++++++++++*/
  
  if(siv_solver[Solver_id].parallel){

    problem_id = siv_solver[Solver_id].problem_id;
    vec_norm = pdr_vec_norm(problem_id, Solver_id, Level_id, Nrdof, Vector);

  } else {

    vec_norm = dnrm2_(&Nrdof, Vector, &IONE);

  }

  return(vec_norm);
}

  /**
   * @param[in] space The list of points in space that need to be iterated.
   */
  void operator() (const pfunc::space_1D& space) const {

    for (size_t i=space.begin(); i<space.end(); ++i) {
      EdgeListType edge_list = (*set_map)(i);
      RandomAccessContainer A (M*edge_list.size());
      RandomAccessContainer X (edge_list.size());

      /* Response is common across all the vertices */
      int response;

      /* Set up A, X, and R */
      EdgeListIteratorType list_iter = edge_list.begin ();
      int index=0;
      while (list_iter != edge_list.end()) {
        edge_type current_entry = (*list_iter);
        const snp_type snp = (*snp_map)(current_entry.source);
        response = current_entry.target;
        const float_type weight = current_entry.weight;

        /* populate A */
        (*materializer)(snp, A.begin()+(index*M));

        /* populate X */
        X[index] = weight;

        /* populate R */
        int int_M = static_cast<int>(M);
        double* Y_ptr = const_cast<double*>(&((*Y)[0])+(response*M));
        double* R_ptr = const_cast<double*>(&((*R)[0])+(response*M));
        dcopy_ (&int_M, Y_ptr, &ONE_STEP, R_ptr, &ONE_STEP);

        ++list_iter;
        ++index;
      }

      /* Now, compute R = AX-R */
      double* A_ptr = &(A[0]);
      double* X_ptr = &(X[0]);
      double* R_ptr = &((*R)[0])+(response*M);
      int LDA = static_cast<int>(M);
      int LDX = static_cast<int>(edge_list.size());
      int LDR = static_cast<int>(M);

      dgemm_ (&NO_TRANS,
              &NO_TRANS,
              &LDR,
              &ONE_STEP,
              &LDX,
              &PLUS_ONE,
              A_ptr,
              &LDA,
              X_ptr,
              &LDX,
              &MINUS_ONE,
              R_ptr,
              &LDR);
      (*l2_norm)[response] = dnrm2_(&LDR, R_ptr, &ONE_STEP);
    }
  }

  /**
   * @param[in] space The list of points in space that need to be iterated.
   */
  void operator() (const pfunc::space_1D& space) const {

    double* R_ptr = &((*R)[0]);
    int LDR = M;

    /* Loop over each element in the range that is given to us */
    for (size_t i=space.begin(); i<space.end(); ++i) 
      (*l2)[i] = dnrm2_ (&LDR, R_ptr+(i*LDR), &ONE_STEP);
  }

int TRON::trcg(double *g, double *s, double *r, int num_sample, long int *num_data)
{
	int i, inc = 1;
	int n = fun_obj->get_nr_variable();
	double one = 1;
	double *d = new double[n];
	double *Hd = new double[n];
	double rTr, rnewTrnew, alpha, beta, cgtol;

	for (i=0; i<n; i++)
	{
		Hd[i] = 0.0;
		s[i] = 0.0;
		r[i] = -g[i];
		d[i] = r[i];
	}
	
	cgtol = 0.1*dnrm2_(&n, g, &inc);

	int cg_iter = 0;
	rTr = ddot_(&n, r, &inc, r, &inc);
	while (1)
	{
		if (dnrm2_(&n, r, &inc) <= cgtol || cg_iter >= max_cg_iter)
			break;
		cg_iter++;
		fun_obj->sample_Hv(d, Hd, num_data);

		alpha = rTr/ddot_(&n, d, &inc, Hd, &inc);
		daxpy_(&n, &alpha, d, &inc, s, &inc);
		alpha = -alpha;
		daxpy_(&n, &alpha, Hd, &inc, r, &inc);
		rnewTrnew = ddot_(&n, r, &inc, r, &inc);
		beta = rnewTrnew/rTr;
		dscal_(&n, &beta, d, &inc);
		daxpy_(&n, &one, r, &inc, d, &inc);
		rTr = rnewTrnew;
	}
	delete[] d;
	delete[] Hd;

	return(cg_iter);
}

	double SPDTensorDL::f(Variable *x) const
	{
		const double *xptr = x->ObtainReadData();
		/*each slice of Xaj is B alpha_j in [(6), CS15]*/
		SharedSpace *Xalpha = new SharedSpace(3, dim, dim, N);
		double *Xalphaptr = Xalpha->ObtainWriteEntireData();

		integer dd = dim * dim, nnum = num, NN = N;
		/*Xalpha <-- \mathbb{B} alpha*/
		dgemm_(GLOBAL::N, GLOBAL::N, &dd, &NN, &nnum, &GLOBAL::DONE, const_cast<double *> (xptr), &dd,
			alpha, &nnum, &GLOBAL::DZERO, Xalphaptr, &dd);
		x->AddToTempData("Xalpha", Xalpha);
		/*compute cholesky decomposition for all slices in Xalpha*/
		SPDTensor *Mani = dynamic_cast<SPDTensor *> (Domain);
		Mani->CholeskyRepresentation(x);

		const SharedSpace *SharedL = x->ObtainReadTempData("XaL");
		const double *L = SharedL->ObtainReadData();

		SharedSpace *SharedlogLXL = new SharedSpace(3, dim, dim, N);
		double *logLXL = SharedlogLXL->ObtainWriteEntireData();
		double *Ltmp = new double[dim * dim];
		integer length = dim * dim, ddim = dim, info;
		for (integer i = 0; i < N; i++)
		{
			dcopy_(&length, const_cast<double *> (L) + i * length, &GLOBAL::IONE, Ltmp, &GLOBAL::IONE);
			/*Solve the linear system Ls X = Li, i.e., X = Ls^{-1} Li. The solution X is stored in Li.
			Note that Li is a lower triangular matrix.
			Details: http://www.netlib.org/lapack/explore-html/d6/d6f/dtrtrs_8f.html */
			dtrtrs_(GLOBAL::L, GLOBAL::N, GLOBAL::N, &ddim, &ddim, Ls + dim * dim * i, &ddim, Ltmp, &ddim, &info);
			if (info != 0)
			{
				std::cout << "Warning: Solving linear system in SPDTensorDL::f failed with info:" << info << "!" << std::endl;
			}
			dgemm_(GLOBAL::N, GLOBAL::T, &ddim, &ddim, &ddim, &GLOBAL::DONE, Ltmp, &ddim, Ltmp, &ddim, &GLOBAL::DZERO, logLXL + ddim * ddim * i, &ddim);
			Matrix MMt(logLXL + ddim * ddim * i, ddim, ddim);
			Matrix::LogSymmetricM(GLOBAL::L, MMt, MMt);
		}
		delete[] Ltmp;

		length = dim * dim * N;
		double result = dnrm2_(&length, logLXL, &GLOBAL::IONE);
		x->AddToTempData("logLXL", SharedlogLXL);
		result *= result;
		result /= 2.0;
		/*add \Omega(X) = \sum \tr(X_i)*/
		for (integer i = 0; i < num; i++)
		{
			for (integer j = 0; j < dim; j++)
			{
				result += lambdaX * xptr[i * dim * dim + j * dim + j];
			}
		}
		return result;
	};

/*! return its Euclidean norm */
inline double nrm2(const _drovector& vec)
{
#ifdef  CPPL_VERBOSE
    std::cerr << "# [MARK] nrm2(const _drovector&)"
              << std::endl;
#endif//CPPL_VERBOSE

    double val( dnrm2_(vec.L, vec.Array, 1) );
    vec.destroy();
    return val;
}

int toScalarR(int code, KDVEC(x), DVEC(r)) {
    REQUIRES(rn==1,BAD_SIZE);
    DEBUGMSG("toScalarR");
    double res;
    integer one = 1;
    integer n = xn;
    switch(code) {
        case 0: { res = dnrm2_(&n,xp,&one); break; }
        case 1: { res = dasum_(&n,xp,&one);  break; }
        case 2: { res = vector_max_index(V(x));  break; }
        case 3: { res = vector_max(V(x));  break; }
        case 4: { res = vector_min_index(V(x)); break; }
        case 5: { res = vector_min(V(x)); break; }
        default: ERROR(BAD_CODE);
    }
    rp[0] = res;
    OK
}

//.........这里部分代码省略.........
            The leading dimension of S.  LDS >= max(1,K).   

    INFO    (output) INTEGER   
            = 0:  successful exit.   
            < 0:  if INFO = -i, the i-th argument had an illegal value.   
            > 0:  if INFO = 1, an eigenvalue did not converge   

    Further Details   
    ===============   

    Based on contributions by   
       Jeff Rutter, Computer Science Division, University of California   
       at Berkeley, USA   
    Modified by Francoise Tisseur, University of Tennessee.   

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


       Test the input parameters.   

       Parameter adjustments */
    /* Table of constant values */
    static integer c__1 = 1;
    static doublereal c_b22 = 1.;
    static doublereal c_b23 = 0.;
    
    /* System generated locals */
    integer q_dim1, q_offset, i__1, i__2;
    doublereal d__1;
    /* Builtin functions */
    double sqrt(doublereal), d_sign(doublereal *, doublereal *);
    /* Local variables */
    static doublereal temp;
    extern doublereal dnrm2_(integer *, doublereal *, integer *);
    static integer i__, j;
    extern /* Subroutine */ int dgemm_(char *, char *, integer *, integer *, 
	    integer *, doublereal *, doublereal *, integer *, doublereal *, 
	    integer *, doublereal *, doublereal *, integer *),
	     dcopy_(integer *, doublereal *, integer *, doublereal *, integer 
	    *), dlaed4_(integer *, integer *, doublereal *, doublereal *, 
	    doublereal *, doublereal *, doublereal *, integer *);
    static integer n2;
    extern doublereal dlamc3_(doublereal *, doublereal *);
    static integer n12, ii, n23;
    extern /* Subroutine */ int dlacpy_(char *, integer *, integer *, 
	    doublereal *, integer *, doublereal *, integer *), 
	    dlaset_(char *, integer *, integer *, doublereal *, doublereal *, 
	    doublereal *, integer *), xerbla_(char *, integer *);
    static integer iq2;
#define q_ref(a_1,a_2) q[(a_2)*q_dim1 + a_1]


    --d__;
    q_dim1 = *ldq;
    q_offset = 1 + q_dim1 * 1;
    q -= q_offset;
    --dlamda;
    --q2;
    --indx;
    --ctot;
    --w;
    --s;

    /* Function Body */
    *info = 0;

double dnrm2( int N, double *a, int inca ){
	return dnrm2_( &N, a, &inca );
};

/* Subroutine */ int dstein_(integer * n, doublereal * d__, doublereal * e,
                             integer * m, doublereal * w, integer * iblock,
                             integer * isplit, doublereal * z__, integer * ldz,
                             doublereal * work, integer * iwork,
                             integer * ifail, integer * info)
{
        /* System generated locals */
        integer z_dim1, z_offset, i__1, i__2, i__3;
        doublereal d__1, d__2, d__3, d__4, d__5;

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

        /* Local variables */
        static integer jblk, nblk;
        extern doublereal ddot_(integer *, doublereal *, integer *,
                                doublereal *, integer *);
        static integer jmax;
        extern doublereal dnrm2_(integer *, doublereal *, integer *);
        static integer i__, j;
        extern /* Subroutine */ int dscal_(integer *, doublereal *,
                                           doublereal *,
                                           integer *);
        static integer iseed[4], gpind, iinfo;
        extern doublereal dasum_(integer *, doublereal *, integer *);
        extern /* Subroutine */ int dcopy_(integer *, doublereal *, integer *,
                                           doublereal *, integer *);
        static integer b1;
        extern /* Subroutine */ int daxpy_(integer *, doublereal *,
                                           doublereal *,
                                           integer *, doublereal *, integer *);
        static integer j1;
        static doublereal ortol;
        static integer indrv1, indrv2, indrv3, indrv4, indrv5, bn;
        extern doublereal dlamch_(char *);
        extern /* Subroutine */ int dlagtf_(integer *, doublereal *,
                                            doublereal *,
                                            doublereal *, doublereal *,
                                            doublereal *, doublereal *,
                                            integer *, integer *);
        static doublereal xj;
        extern integer idamax_(integer *, doublereal *, integer *);
        extern /* Subroutine */ int xerbla_(char *, integer *), dlagts_(
                                                                               integer
                                                                               *,
                                                                               integer
                                                                               *,
                                                                               doublereal
                                                                               *,
                                                                               doublereal
                                                                               *,
                                                                               doublereal
                                                                               *,
                                                                               doublereal
                                                                               *,
                                                                               integer
                                                                               *,
                                                                               doublereal
                                                                               *,
                                                                               doublereal
                                                                               *,
                                                                               integer
                                                                               *);
        static integer nrmchk;
        extern /* Subroutine */ int dlarnv_(integer *, integer *, integer *,
                                            doublereal *);
        static integer blksiz;
        static doublereal onenrm, dtpcrt, pertol, scl, eps, sep, nrm, tol;
        static integer its;
        static doublereal xjm, ztr, eps1;

#define z___ref(a_1,a_2) z__[(a_2)*z_dim1 + a_1]

/*  -- LAPACK routine (instrumented to count operations, version 3.0) --   
       Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,   
       Courant Institute, Argonne National Lab, and Rice University   
       September 30, 1994   

       Common block to return operation count and iteration count   
       ITCNT is initialized to 0, OPS is only incremented   

    Purpose   
    =======   

    DSTEIN 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).   

    Arguments   
    =========   

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

    D       (input) DOUBLE PRECISION array, dimension (N)   
            The n diagonal elements of the tridiagonal matrix T.   

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

void TRON::tron(double *w)
{
	// Parameters for updating the iterates.
	double eta0 = 1e-4, eta1 = 0.25, eta2 = 0.75;

	// Parameters for updating the trust region size delta.
	double sigma1 = 0.25, sigma2 = 0.5, sigma3 = 4;

	Int n = fun_obj->get_nr_variable();
	Int i, cg_iter;
	double delta, snorm, one=1.0;
	double alpha, f, fnew, prered, actred, gs;
	Int search = 1, iter = 1, inc = 1;
	double *s = new double[n];
	double *r = new double[n];
	double *g = new double[n];

	// calculate gradient norm at w=0 for stopping condition.
	double *w0 = new double[n];
	for (i=0; i<n; i++)
		w0[i] = 0;
	fun_obj->fun(w0);
	fun_obj->grad(w0, g);
	double gnorm0 = dnrm2_(&n, g, &inc);
	delete [] w0;

	f = fun_obj->fun(w);
	fun_obj->grad(w, g);
	delta = dnrm2_(&n, g, &inc);
	double gnorm = delta;

	if (gnorm <= eps*gnorm0)
		search = 0;

	iter = 1;

	double *w_new = new double[n];
	while (iter <= max_iter && search)
	{
		cg_iter = trcg(delta, g, s, r);

		memcpy(w_new, w, sizeof(double)*n);
		daxpy_(&n, &one, s, &inc, w_new, &inc);

		gs = ddot_(&n, g, &inc, s, &inc);
		prered = -0.5*(gs-ddot_(&n, s, &inc, r, &inc));
		fnew = fun_obj->fun(w_new);

		// Compute the actual reduction.
		actred = f - fnew;

		// On the first iteration, adjust the initial step bound.
		snorm = dnrm2_(&n, s, &inc);
		if (iter == 1)
			delta = min(delta, snorm);

		// Compute prediction alpha*snorm of the step.
		if (fnew - f - gs <= 0)
			alpha = sigma3;
		else
			alpha = max(sigma1, -0.5*(gs/(fnew - f - gs)));

		// Update the trust region bound according to the ratio of actual to predicted reduction.
		if (actred < eta0*prered)
			delta = min(max(alpha, sigma1)*snorm, sigma2*delta);
		else if (actred < eta1*prered)
			delta = max(sigma1*delta, min(alpha*snorm, sigma2*delta));
		else if (actred < eta2*prered)
			delta = max(sigma1*delta, min(alpha*snorm, sigma3*delta));
		else
			delta = max(delta, min(alpha*snorm, sigma3*delta));

		info("iter %2d act %5.3e pre %5.3e delta %5.3e f %5.3e |g| %5.3e CG %3d\n", iter, actred, prered, delta, f, gnorm, cg_iter);

		if (actred > eta0*prered)
		{
			iter++;
			memcpy(w, w_new, sizeof(double)*n);
			f = fnew;
			fun_obj->grad(w, g);

			gnorm = dnrm2_(&n, g, &inc);
			if (gnorm <= eps*gnorm0)
				break;
		}
		if (f < -1.0e+32)
		{
			info("WARNING: f < -1.0e+32\n");
			break;
		}
		if (fabs(actred) <= 0 && prered <= 0)
		{
			info("WARNING: actred and prered <= 0\n");
			break;
		}
		if (fabs(actred) <= 1.0e-12*fabs(f) &&
		    fabs(prered) <= 1.0e-12*fabs(f))
		{
			info("WARNING: actred and prered too small\n");
			break;
//.........这里部分代码省略.........

/* Subroutine */ int dqrt15_(integer *scale, integer *rksel, integer *m, 
	integer *n, integer *nrhs, doublereal *a, integer *lda, doublereal *b,
	 integer *ldb, doublereal *s, integer *rank, doublereal *norma, 
	doublereal *normb, integer *iseed, doublereal *work, integer *lwork)
{
    /* System generated locals */
    integer a_dim1, a_offset, b_dim1, b_offset, i__1, i__2;
    doublereal d__1;

    /* Local variables */
    static integer info;
    static doublereal temp;
    extern doublereal dnrm2_(integer *, doublereal *, integer *);
    static integer j;
    extern /* Subroutine */ int dscal_(integer *, doublereal *, doublereal *, 
	    integer *), dlarf_(char *, integer *, integer *, doublereal *, 
	    integer *, doublereal *, doublereal *, integer *, doublereal *), dgemm_(char *, char *, integer *, integer *, integer *, 
	    doublereal *, doublereal *, integer *, doublereal *, integer *, 
	    doublereal *, doublereal *, integer *);
    extern doublereal dasum_(integer *, doublereal *, integer *);
    static doublereal dummy[1];
    extern doublereal dlamch_(char *), dlange_(char *, integer *, 
	    integer *, doublereal *, integer *, doublereal *);
    static integer mn;
    extern /* Subroutine */ int dlascl_(char *, integer *, integer *, 
	    doublereal *, doublereal *, integer *, integer *, doublereal *, 
	    integer *, integer *);
    extern doublereal dlarnd_(integer *, integer *);
    extern /* Subroutine */ int dlaord_(char *, integer *, doublereal *, 
	    integer *), dlaset_(char *, integer *, integer *, 
	    doublereal *, doublereal *, doublereal *, integer *), 
	    xerbla_(char *, integer *);
    static doublereal bignum;
    extern /* Subroutine */ int dlaror_(char *, char *, integer *, integer *, 
	    doublereal *, integer *, integer *, doublereal *, integer *), dlarnv_(integer *, integer *, integer *, 
	    doublereal *);
    static doublereal 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   
    =======   

    DQRT15 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) DOUBLE PRECISION array, dimension (LDA,N)   
            The M-by-N matrix A.   

    LDA     (input) INTEGER   
            The leading dimension of the array A.   

    B       (output) DOUBLE PRECISION 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) DOUBLE PRECISION array, dimension MIN(M,N)   
            Singular values of A.   

    RANK    (output) INTEGER   
            number of nonzero singular values of A.   

    NORMA   (output) DOUBLE PRECISION   
            one-norm of A.   

    NORMB   (output) DOUBLE PRECISION   
            one-norm of B.   

    ISEED   (input/output) integer array, dimension (4)   
//.........这里部分代码省略.........

doublereal dqrt12_(integer *m, integer *n, doublereal *a, integer *lda, 
	doublereal *s, doublereal *work, integer *lwork)
{
    /* System generated locals */
    integer a_dim1, a_offset, i__1, i__2;
    doublereal ret_val;

    /* Local variables */
    integer i__, j, mn, iscl, info;
    doublereal anrm;
    extern doublereal dnrm2_(integer *, doublereal *, integer *), dasum_(
	    integer *, doublereal *, integer *);
    extern /* Subroutine */ int daxpy_(integer *, doublereal *, doublereal *, 
	    integer *, doublereal *, integer *), dgebd2_(integer *, integer *, 
	     doublereal *, integer *, doublereal *, doublereal *, doublereal *
, doublereal *, doublereal *, integer *);
    doublereal dummy[1];
    extern /* Subroutine */ int dlabad_(doublereal *, doublereal *);
    extern doublereal dlamch_(char *), dlange_(char *, integer *, 
	    integer *, doublereal *, integer *, doublereal *);
    extern /* Subroutine */ int dlascl_(char *, integer *, integer *, 
	    doublereal *, doublereal *, integer *, integer *, doublereal *, 
	    integer *, integer *), dlaset_(char *, integer *, integer 
	    *, doublereal *, doublereal *, doublereal *, integer *), 
	    xerbla_(char *, integer *), dbdsqr_(char *, integer *, 
	    integer *, integer *, integer *, doublereal *, doublereal *, 
	    doublereal *, integer *, doublereal *, integer *, doublereal *, 
	    integer *, doublereal *, integer *);
    doublereal bignum, 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 */
/*  ======= */

/*  DQRT12 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) DOUBLE PRECISION 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) DOUBLE PRECISION array, dimension (min(M,N)) */
/*          The singular values of the matrix A. */

/*  WORK    (workspace) DOUBLE PRECISION 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.;

/*     Test that enough workspace is supplied */

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

/* Subroutine */ int dlarge_(integer *n, doublereal *a, integer *lda, integer
        *iseed, doublereal *work, integer *info)
{
    /* System generated locals */
    integer a_dim1, a_offset, i__1;
    doublereal d__1;

    /* Builtin functions */
    double d_sign(doublereal *, doublereal *);

    /* Local variables */
    extern /* Subroutine */ int dger_(integer *, integer *, doublereal *,
            doublereal *, integer *, doublereal *, integer *, doublereal *,
            integer *);
    extern doublereal dnrm2_(integer *, doublereal *, integer *);
    static integer i;
    extern /* Subroutine */ int dscal_(integer *, doublereal *, doublereal *,
            integer *), dgemv_(char *, integer *, integer *, doublereal *,
            doublereal *, integer *, doublereal *, integer *, doublereal *,
            doublereal *, integer *);
    static doublereal wa, wb, wn;
    extern /* Subroutine */ int xerbla_(char *, integer *), dlarnv_(
            integer *, integer *, integer *, doublereal *);
    static doublereal tau;


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


    Purpose
    =======

    DLARGE pre- and post-multiplies a real general n by n matrix A
    with a random orthogonal matrix: A = U*D*U'.

    Arguments
    =========

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

    A       (input/output) DOUBLE PRECISION array, dimension (LDA,N)
            On entry, the original n by n matrix A.
            On exit, A is overwritten by U*A*U' for some random
            orthogonal matrix U.

    LDA     (input) INTEGER
            The leading dimension of the array A.  LDA >= N.

    ISEED   (input/output) INTEGER array, dimension (4)
            On entry, the seed of the random number generator; the array

            elements must be between 0 and 4095, and ISEED(4) must be
            odd.
            On exit, the seed is updated.

    WORK    (workspace) DOUBLE PRECISION array, dimension (2*N)

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

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



       Test the input arguments

       Parameter adjustments */
    a_dim1 = *lda;
    a_offset = a_dim1 + 1;
    a -= a_offset;
    --iseed;
    --work;

    /* Function Body */
    *info = 0;
    if (*n < 0) {
        *info = -1;
    } else if (*lda < max(1,*n)) {
        *info = -3;
    }
    if (*info < 0) {
        i__1 = -(*info);
        xerbla_("DLARGE", &i__1);
        return 0;
    }

/*     pre- and post-multiply A by random orthogonal matrix */

    for (i = *n; i >= 1; --i) {

/*        generate random reflection */

        i__1 = *n - i + 1;
        dlarnv_(&c__3, &iseed[1], &i__1, &work[1]);
        i__1 = *n - i + 1;
//.........这里部分代码省略.........

/* Subroutine */
int dlaein_(logical *rightv, logical *noinit, integer *n, doublereal *h__, integer *ldh, doublereal *wr, doublereal *wi, doublereal *vr, doublereal *vi, doublereal *b, integer *ldb, doublereal *work, doublereal *eps3, doublereal *smlnum, doublereal * bignum, integer *info)
{
    /* System generated locals */
    integer b_dim1, b_offset, h_dim1, h_offset, i__1, i__2, i__3, i__4;
    doublereal d__1, d__2, d__3, d__4;
    /* Builtin functions */
    double sqrt(doublereal);
    /* Local variables */
    integer i__, j;
    doublereal w, x, y;
    integer i1, i2, i3;
    doublereal w1, ei, ej, xi, xr, rec;
    integer its, ierr;
    doublereal temp, norm, vmax;
    extern doublereal dnrm2_(integer *, doublereal *, integer *);
    extern /* Subroutine */
    int dscal_(integer *, doublereal *, doublereal *, integer *);
    doublereal scale;
    extern doublereal dasum_(integer *, doublereal *, integer *);
    char trans[1];
    doublereal vcrit, rootn, vnorm;
    extern doublereal dlapy2_(doublereal *, doublereal *);
    doublereal absbii, absbjj;
    extern integer idamax_(integer *, doublereal *, integer *);
    extern /* Subroutine */
    int dladiv_(doublereal *, doublereal *, doublereal *, doublereal *, doublereal *, doublereal *), dlatrs_( char *, char *, char *, char *, integer *, doublereal *, integer * , doublereal *, doublereal *, doublereal *, integer *);
    char normin[1];
    doublereal nrmsml, 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((doublereal) (*n));
    growto = .1 / rootn;
    /* Computing MAX */
    d__1 = 1.;
    d__2 = *eps3 * rootn; // , expr subst
    nrmsml = max(d__1,d__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.)
    {
        /* Real eigenvalue. */
        if (*noinit)
        {
            /* Set initial vector. */
            i__1 = *n;
            for (i__ = 1;
                    i__ <= i__1;
                    ++i__)
            {
                vr[i__] = *eps3;
                /* L30: */
//.........这里部分代码省略.........

int fgmr(int n,
     void (*matvec) (double, double[], double, double[]),
     void (*psolve) (int, double[], double[]),
     double *rhs, double *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;
    double **hh, *c, *s, *rs, t, t0;
    double beta, eps1 = 0.0, gam, **vv, **z;

    its = 0;
    vv = (double **)SUPERLU_MALLOC((im + 1) * sizeof(double *));
    for (i = 0; i <= im; i++)
	vv[i] = doubleMalloc(n);
    z = (double **)SUPERLU_MALLOC(im * sizeof(double *));
    hh = (double **)SUPERLU_MALLOC(im * sizeof(double *));
    for (i = 0; i < im; i++)
    {
	hh[i] = doubleMalloc(i + 2);
	z[i] = doubleMalloc(n);
    }
    c = doubleMalloc(im);
    s = doubleMalloc(im);
    rs = doubleMalloc(im + 1);

    /*---- outer loop starts here ----*/
    do
    {
	/*---- compute initial residual vector ----*/
	matvec(1.0, sol, 0.0, vv[0]);
	for (j = 0; j < n; j++)
	    vv[0][j] = rhs[j] - vv[0][j];	/* vv[0]= initial residual */
	beta = dnrm2_(&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 * dnrm2_(&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}
	    +-----------------------------------------------------------*/
	    if (psolve)
		psolve(n, z[i], vv[i]);
//.........这里部分代码省略.........

/*! return its Euclidean norm */
inline double nrm2(const _dcovector& vec)
{   VERBOSE_REPORT;
    double val( dnrm2_(vec.l, vec.array, 1) );
    vec.destroy();
    return val;
}