C++ dnrm2_函数代码示例

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

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