PyObject* multi_dotu(PyObject *self, PyObject *args)
{
  PyArrayObject* a;
  PyArrayObject* b;
  PyArrayObject* c;
  if (!PyArg_ParseTuple(args, "OOO", &a, &b, &c)) 
    return NULL;
  int n0 = PyArray_DIMS(a)[0];
  int n = PyArray_DIMS(a)[1];
  for (int i = 2; i < PyArray_NDIM(a); i++)
    n *= PyArray_DIMS(a)[i];
  int incx = 1;
  int incy = 1;
  if (PyArray_DESCR(a)->type_num == NPY_DOUBLE)
    {
      double *ap = DOUBLEP(a);
      double *bp = DOUBLEP(b);
      double *cp = DOUBLEP(c);

      for (int i = 0; i < n0; i++)
        {
          cp[i] = ddot_(&n, (void*)ap, 
             &incx, (void*)bp, &incy);
          ap += n;
          bp += n;
        }
    }
  else
    {
      double_complex* ap = COMPLEXP(a);
      double_complex* bp = COMPLEXP(b);
      double_complex* cp = COMPLEXP(c);
      for (int i = 0; i < n0; i++)
        {
          cp[i] = 0.0;
          for (int j = 0; j < n; j++)
              cp[i] += ap[j] * bp[j];
          ap += n;
          bp += n;
        }
    }
  Py_RETURN_NONE;
}

void QToCurve(const double *Q, integer d, integer n, double *C, bool isclosed)
{
	double *q2n = new double[n + n * d];
	double *q2nTimesQ = q2n + n;
	integer inc = n;

	for (integer i = 0; i < n; i++)
	{
		q2n[i] = sqrt(ddot_(&d, const_cast<double *> (Q + i), &inc, const_cast<double *> (Q + i), &inc));
	}

	ElasticCurvesRO::PointwiseQProdl(Q, q2n, d, n, q2nTimesQ);
	for (integer i = 0; i < d; i++)
	{
		ElasticCurvesRO::CumTrapz(q2nTimesQ + i * n, n, 1.0 / (n - 1), C);
	}

	delete[] q2n;
};

double CheMPS2::Excitation::third_middle( const int ikappa, const SyBookkeeper * book_up, const SyBookkeeper * book_down, const double gamma, Sobject * S_up, Sobject * S_down, TensorO * Lovlp, TensorO * Rovlp, double * workmem1, double * workmem2 ){

   const int index = S_up->gIndex();
   const int TwoSL = S_up->gTwoSL( ikappa );
   const int TwoSR = S_up->gTwoSR( ikappa );
   const int TwoJ = S_up->gTwoJ( ikappa );
   const int NL = S_up->gNL( ikappa );
   const int NR = S_up->gNR( ikappa );
   const int IL = S_up->gIL( ikappa );
   const int IR = S_up->gIR( ikappa );
   const int N1 = S_up->gN1( ikappa );
   const int N2 = S_up->gN2( ikappa );

   int dimLup   = book_up  ->gCurrentDim( index,     NL, TwoSL, IL );
   int dimRup   = book_up  ->gCurrentDim( index + 2, NR, TwoSR, IR );
   int dimLdown = book_down->gCurrentDim( index,     NL, TwoSL, IL );
   int dimRdown = book_down->gCurrentDim( index + 2, NR, TwoSR, IR );

   double inproduct = 0.0;
   if (( dimLdown > 0 ) && ( dimRdown > 0 )){
      double * block_down  = S_down->gStorage( NL, TwoSL, IL, N1, N2, TwoJ, NR, TwoSR, IR );
      double * block_left  = Lovlp ->gStorage( NL, TwoSL, IL, NL, TwoSL, IL );
      double * block_right = Rovlp ->gStorage( NR, TwoSR, IR, NR, TwoSR, IR );
      char trans = 'T';
      char notrans = 'N';
      double one = 1.0;
      double set = 0.0;
      dgemm_( &notrans, &notrans, &dimLup, &dimRdown, &dimLdown, &one, block_left, &dimLup, block_down,  &dimLdown, &set, workmem1, &dimLup );
      dgemm_( &notrans, &trans,   &dimLup, &dimRup,   &dimRdown, &one, workmem1,   &dimLup, block_right, &dimRup,   &set, workmem2, &dimLup );

      double * block_up = S_up->gStorage() + S_up->gKappa2index( ikappa );
      int size = dimLup * dimRup;
      int inc1 = 1;
      if ( fabs( gamma ) > 0.0 ){
         double factor = gamma;
         daxpy_( &size, &factor, workmem2, &inc1, block_up, &inc1 );
      }
      inproduct = ddot_( &size, workmem2, &inc1, block_up, &inc1 );
   }
   return inproduct;

}

void CheMPS2::Heff::addDiagramExcitations(const int ikappa, double * memS, double * memHeff, const Sobject * denS, int nLower, double ** VeffTilde) const{

   int dimTotal = denS->gKappa2index(denS->gNKappa());
   int ptr = denS->gKappa2index(ikappa);
   int dimBlock = denS->gKappa2index(ikappa+1) - ptr;
   int inc = 1;
   #ifdef CHEMPS2_MPI_COMPILATION
   const int MPIRANK = MPIchemps2::mpi_rank();
   #endif
   
   for (int state=0; state<nLower; state++){
      #ifdef CHEMPS2_MPI_COMPILATION
      if ( MPIchemps2::owner_specific_excitation( Prob->gL(), state ) == MPIRANK )
      #endif
      {
         double alpha = ddot_(&dimTotal, memS, &inc, VeffTilde[state], &inc);
         daxpy_(&dimBlock,&alpha,VeffTilde[state]+ptr,&inc,memHeff+ptr,&inc);
      }
   }

}

int main(int argc, char** argv) {
    /* You can define arrays on the heap */
    double *a = (double*) malloc( 3 * sizeof(double) );
    a[0] = 1.0;
    a[1] = 2.0;
    a[2] = 3.0;
    
    /* or on the stack */
    double b[3] = { 4.0, 5.0, 6.0 };
    
    printf("A=\tB=\n");
    for (int i = 0; i < 3; i++) {
        printf("%2.1f\t%2.1f\n", a[i], b[i]);
    }
    
    int N = 3, one = 1;
    double dot_product = ddot_(&N, a, &one, b, &one);
    printf("\nA dot B = %2.1f\n", dot_product);
    
    return 0;
};

/*! drovector*dcovector operator */
inline double operator*(const drovector& rovec, const dcovector& covec)
{
#ifdef  CPPL_VERBOSE
  std::cerr << "# [MARK] operator*(const drovector&, const dcovector&)"
            << std::endl;
#endif//CPPL_VERBOSE
  
#ifdef  CPPL_DEBUG
  if(rovec.L!=covec.L){
    std::cerr << "[ERROR] operator*(const drovector&, const dcovector&)"
              << std::endl
              << "These two vectors can not make a product." << std::endl
              << "Your input was (" << rovec.L << ") * (" << covec.L << ")."
              << std::endl;
    exit(1);
  }
#endif//CPPL_DEBUG
  
  double val( ddot_( rovec.L, rovec.Array, 1, covec.Array, 1 ) );
  
  return val;
}

//.........这里部分代码省略.........
      res2 = cblas_ddot (n, r, 1, r, 1);
      if (it->debug == 2)
	{
	  fprintf (it->out, "libiter-cgs %d %e\n", i, res2 / b2);
	}
      if (res2 <= eps2)
	{
	  ret = 0; // success
	  break;
	}

      rho1 = cblas_ddot (n, r0, 1, r, 1); // rho = (r0*, r)
      beta = rho1 / rho;
      rho = rho1;

      // here t is not used so that this is used for working area.
      double *qu = t;
      cblas_dcopy (n, q, 1, qu, 1); // qu = q
      cblas_daxpy (n, -1.0, u, 1, qu, 1); // qu = q - u
      cblas_dcopy (n, r, 1, u, 1); // u = r
      cblas_daxpy (n, beta, qu, 1, u, 1); // u = r + beta (q - u)

      cblas_daxpy (n, beta, p, 1, qu, 1); // qu = q - u + beta * p
      cblas_dcopy (n, u, 1, p, 1); // p = u
      cblas_daxpy (n, beta, qu, 1, p, 1); // p = u + beta (q - u + b * p)
    }

#else // !HAVE_CBLAS_H
# ifdef HAVE_BLAS_H
  /**
   * BLAS version
   */

  double b2 = ddot_ (&n, b, &i_1, b, &i_1); // (b,b)
  eps2 *= b2;

  // initial residue
  atimes (n, x, r, atimes_param); // r = A.x
  daxpy_ (&n, &d_m1, b, &i_1, r, &i_1); // r = A.x - b

  dcopy_ (&n, r, &i_1, r0, &i_1); // r0* = r
  dcopy_ (&n, r, &i_1, p, &i_1); // p = r
  dcopy_ (&n, r, &i_1, u, &i_1); // u = r

  rho = ddot_ (&n, r0, &i_1, r, &i_1); // rho = (r0*, r)

  int i;
  for (i = 0; i < itmax; i ++)
    {
      atimes (n, p, ap, atimes_param); // ap = A.p
      r0ap = ddot_ (&n, r0, &i_1, ap, &i_1); // r0ap = (r0*, A.p)
      delta = - rho / r0ap;

      dcopy_ (&n, u, &i_1, q, &i_1); // q = u
      dscal_ (&n, &d_2, q, &i_1); // q = 2 u
      daxpy_ (&n, &delta, ap, &i_1, q, &i_1); // q = 2 u + delta A.p

      atimes (n, q, t, atimes_param); // t = A.q

      daxpy_ (&n, &delta, t, &i_1, r, &i_1); // r = r + delta t
      daxpy_ (&n, &delta, q, &i_1, x, &i_1); // x = x + delta q

      res2 = ddot_ (&n, r, &i_1, r, &i_1);
      if (it->debug == 2)
	{
	  fprintf (it->out, "libiter-cgs %d %e\n", i, res2 / b2);

double CheMPS2::Excitation::neighbours( const int ikappa, const SyBookkeeper * book_up, const SyBookkeeper * book_down, const double alpha, const double beta, const double gamma, Sobject * S_up, Sobject * S_down ){

   const int index = S_up->gIndex();
   const int TwoSL = S_up->gTwoSL( ikappa );
   const int TwoSR = S_up->gTwoSR( ikappa );
   const int TwoJ = S_up->gTwoJ( ikappa );
   const int NL = S_up->gNL( ikappa );
   const int NR = S_up->gNR( ikappa );
   const int IL = S_up->gIL( ikappa );
   const int IR = S_up->gIR( ikappa );
   const int N1 = S_up->gN1( ikappa );
   const int N2 = S_up->gN2( ikappa );

   const int dimLup   = book_up  ->gCurrentDim( index,     NL, TwoSL, IL );
   const int dimRup   = book_up  ->gCurrentDim( index + 2, NR, TwoSR, IR );
   const int dimLdown = book_down->gCurrentDim( index,     NL, TwoSL, IL );
   const int dimRdown = book_down->gCurrentDim( index + 2, NR, TwoSR, IR );
   assert( dimLup == dimLdown );
   assert( dimRup == dimRdown );
   assert( book_up->gIrrep( index ) == book_up->gIrrep( index + 1 ) );

   double * block_up = S_up->gStorage() + S_up->gKappa2index( ikappa );
   int size = dimLup * dimRup;
   int inc1 = 1;

   // Add a^+ a
   if ( fabs( alpha ) > 0.0 ){
      if (( TwoJ == 0 ) && ( N1 == 1 ) && ( N2 == 1 )){
         double * block_down = S_down->gStorage( NL, TwoSL, IL, 0, 2, TwoJ, NR, TwoSR, IR );
         assert( block_down != NULL );
         double factor = sqrt( 2.0 ) * alpha;
         daxpy_( &size, &factor, block_down, &inc1, block_up, &inc1 );
      }
      if (( TwoJ == 1 ) && ( N1 == 2 ) && ( N2 == 1 )){
         double * block_down = S_down->gStorage( NL, TwoSL, IL, 1, 2, TwoJ, NR, TwoSR, IR );
         assert( block_down != NULL );
         double factor = -alpha;
         daxpy_( &size, &factor, block_down, &inc1, block_up, &inc1 );
      }
      if (( TwoJ == 1 ) && ( N1 == 1 ) && ( N2 == 0 )){
         double * block_down = S_down->gStorage( NL, TwoSL, IL, 0, 1, TwoJ, NR, TwoSR, IR );
         assert( block_down != NULL );
         double factor = alpha;
         daxpy_( &size, &factor, block_down, &inc1, block_up, &inc1 );
      }
      if (( TwoJ == 0 ) && ( N1 == 2 ) && ( N2 == 0 )){
         double * block_down = S_down->gStorage( NL, TwoSL, IL, 1, 1, TwoJ, NR, TwoSR, IR );
         assert( block_down != NULL );
         double factor = sqrt( 2.0 ) * alpha;
         daxpy_( &size, &factor, block_down, &inc1, block_up, &inc1 );
      }
   }

   // Add a a^+
   if ( fabs( beta ) > 0.0 ){
      if (( TwoJ == 0 ) && ( N1 == 1 ) && ( N2 == 1 )){
         double * block_down = S_down->gStorage( NL, TwoSL, IL, 2, 0, TwoJ, NR, TwoSR, IR );
         assert( block_down != NULL );
         double factor = sqrt( 2.0 ) * beta;
         daxpy_( &size, &factor, block_down, &inc1, block_up, &inc1 );
      }
      if (( TwoJ == 1 ) && ( N1 == 0 ) && ( N2 == 1 )){
         double * block_down = S_down->gStorage( NL, TwoSL, IL, 1, 0, TwoJ, NR, TwoSR, IR );
         assert( block_down != NULL );
         double factor = beta;
         daxpy_( &size, &factor, block_down, &inc1, block_up, &inc1 );
      }
      if (( TwoJ == 1 ) && ( N1 == 1 ) && ( N2 == 2 )){
         double * block_down = S_down->gStorage( NL, TwoSL, IL, 2, 1, TwoJ, NR, TwoSR, IR );
         assert( block_down != NULL );
         double factor = -beta;
         daxpy_( &size, &factor, block_down, &inc1, block_up, &inc1 );
      }
      if (( TwoJ == 0 ) && ( N1 == 0 ) && ( N2 == 2 )){
         double * block_down = S_down->gStorage( NL, TwoSL, IL, 1, 1, TwoJ, NR, TwoSR, IR );
         assert( block_down != NULL );
         double factor = sqrt( 2.0 ) * beta;
         daxpy_( &size, &factor, block_down, &inc1, block_up, &inc1 );
      }
   }

   // Add the constant part
   double * block_down = S_down->gStorage( NL, TwoSL, IL, N1, N2, TwoJ, NR, TwoSR, IR );
   assert( block_down != NULL );
   if ( fabs( gamma ) > 0.0 ){
      double factor = gamma;
      daxpy_( &size, &factor, block_down, &inc1, block_up, &inc1 );
   }
   const double inproduct = ddot_( &size, block_down, &inc1, block_up, &inc1 );
   return inproduct;

}

/* Subroutine */
int dpptri_(char *uplo, integer *n, doublereal *ap, integer * info)
{
    /* System generated locals */
    integer i__1, i__2;
    /* Local variables */
    integer j, jc, jj;
    doublereal ajj;
    integer jjn;
    extern doublereal ddot_(integer *, doublereal *, integer *, doublereal *, integer *);
    extern /* Subroutine */
    int dspr_(char *, integer *, doublereal *, doublereal *, integer *, doublereal *), dscal_(integer *, doublereal *, doublereal *, integer *);
    extern logical lsame_(char *, char *);
    extern /* Subroutine */
    int dtpmv_(char *, char *, char *, integer *, doublereal *, doublereal *, integer *);
    logical upper;
    extern /* Subroutine */
    int xerbla_(char *, integer *), dtptri_( char *, char *, integer *, doublereal *, integer *);
    /* -- LAPACK computational routine (version 3.4.0) -- */
    /* -- LAPACK is a software package provided by Univ. of Tennessee, -- */
    /* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- */
    /* November 2011 */
    /* .. Scalar Arguments .. */
    /* .. */
    /* .. Array Arguments .. */
    /* .. */
    /* ===================================================================== */
    /* .. Parameters .. */
    /* .. */
    /* .. Local Scalars .. */
    /* .. */
    /* .. External Functions .. */
    /* .. */
    /* .. External Subroutines .. */
    /* .. */
    /* .. Executable Statements .. */
    /* Test the input parameters. */
    /* Parameter adjustments */
    --ap;
    /* Function Body */
    *info = 0;
    upper = lsame_(uplo, "U");
    if (! upper && ! lsame_(uplo, "L"))
    {
        *info = -1;
    }
    else if (*n < 0)
    {
        *info = -2;
    }
    if (*info != 0)
    {
        i__1 = -(*info);
        xerbla_("DPPTRI", &i__1);
        return 0;
    }
    /* Quick return if possible */
    if (*n == 0)
    {
        return 0;
    }
    /* Invert the triangular Cholesky factor U or L. */
    dtptri_(uplo, "Non-unit", n, &ap[1], info);
    if (*info > 0)
    {
        return 0;
    }
    if (upper)
    {
        /* Compute the product inv(U) * inv(U)**T. */
        jj = 0;
        i__1 = *n;
        for (j = 1;
                j <= i__1;
                ++j)
        {
            jc = jj + 1;
            jj += j;
            if (j > 1)
            {
                i__2 = j - 1;
                dspr_("Upper", &i__2, &c_b8, &ap[jc], &c__1, &ap[1]);
            }
            ajj = ap[jj];
            dscal_(&j, &ajj, &ap[jc], &c__1);
            /* L10: */
        }
    }
    else
    {
        /* Compute the product inv(L)**T * inv(L). */
        jj = 1;
        i__1 = *n;
        for (j = 1;
                j <= i__1;
                ++j)
        {
            jjn = jj + *n - j + 1;
            i__2 = *n - j + 1;
            ap[jj] = ddot_(&i__2, &ap[jj], &c__1, &ap[jj], &c__1);
//.........这里部分代码省略.........

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

/*  UPLO    (input) CHARACTER*1 */
/*          = 'U':  Upper triangular factor is stored in AP; */
/*          = 'L':  Lower triangular factor is stored in AP. */

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

/*  AP      (input/output) DOUBLE PRECISION array, dimension (N*(N+1)/2) */
/*          On entry, the triangular factor U or L from the Cholesky */
/*          factorization A = U**T*U or A = L*L**T, packed columnwise as */
/*          a linear array.  The j-th column of U or L is stored in the */
/*          array AP as follows: */
/*          if UPLO = 'U', AP(i + (j-1)*j/2) = U(i,j) for 1<=i<=j; */
/*          if UPLO = 'L', AP(i + (j-1)*(2n-j)/2) = L(i,j) for j<=i<=n. */

/*          On exit, the upper or lower triangle of the (symmetric) */
/*          inverse of A, overwriting the input factor U or L. */

/*  INFO    (output) INTEGER */
/*          = 0:  successful exit */
/*          < 0:  if INFO = -i, the i-th argument had an illegal value */
/*          > 0:  if INFO = i, the (i,i) element of the factor U or L is */
/*                zero, and the inverse could not be computed. */

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

/*     Test the input parameters. */

    /* Parameter adjustments */
    --ap;

    /* Function Body */
    *info = 0;
    upper = lsame_(uplo, "U");
    if (! upper && ! lsame_(uplo, "L")) {
	*info = -1;
    } else if (*n < 0) {
	*info = -2;
    }
    if (*info != 0) {
	i__1 = -(*info);
	xerbla_("DPPTRI", &i__1);
	return 0;
    }

/*     Quick return if possible */

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

/*     Invert the triangular Cholesky factor U or L. */

    dtptri_(uplo, "Non-unit", n, &ap[1], info);
    if (*info > 0) {
	return 0;
    }

    if (upper) {

/*        Compute the product inv(U) * inv(U)'. */

	jj = 0;
	i__1 = *n;
	for (j = 1; j <= i__1; ++j) {
	    jc = jj + 1;
	    jj += j;
	    if (j > 1) {
		i__2 = j - 1;
		dspr_("Upper", &i__2, &c_b8, &ap[jc], &c__1, &ap[1]);
	    }
	    ajj = ap[jj];
	    dscal_(&j, &ajj, &ap[jc], &c__1);
	}

    } else {

/*        Compute the product inv(L)' * inv(L). */

	jj = 1;
	i__1 = *n;
	for (j = 1; j <= i__1; ++j) {
	    jjn = jj + *n - j + 1;
	    i__2 = *n - j + 1;
	    ap[jj] = ddot_(&i__2, &ap[jj], &c__1, &ap[jj], &c__1);
	    if (j < *n) {
		i__2 = *n - j;
		dtpmv_("Lower", "Transpose", "Non-unit", &i__2, &ap[jjn], &ap[
			jj + 1], &c__1);
	    }
	    jj = jjn;
	}
    }

    return 0;

/*     End of DPPTRI */

} /* dpptri_ */

// BLAS Level 1
double ddot( int N, double *a, int inca, double *b, int incb ){
	return ddot_( &N, a, &inca, b, &incb );
};

/* Subroutine */ int dlaqtr_(logical *ltran, logical *lreal, integer *n, 
	doublereal *t, integer *ldt, doublereal *b, doublereal *w, doublereal 
	*scale, doublereal *x, doublereal *work, integer *info)
{
    /* System generated locals */
    integer t_dim1, t_offset, i__1, i__2;
    doublereal d__1, d__2, d__3, d__4, d__5, d__6;

    /* Local variables */
    doublereal d__[4]	/* was [2][2] */;
    integer i__, j, k;
    doublereal v[4]	/* was [2][2] */, z__;
    integer j1, j2, n1, n2;
    doublereal si, xj, sr, rec, eps, tjj, tmp;
    extern doublereal ddot_(integer *, doublereal *, integer *, doublereal *, 
	    integer *);
    integer ierr;
    doublereal smin, xmax;
    extern /* Subroutine */ int dscal_(integer *, doublereal *, doublereal *, 
	    integer *);
    extern doublereal dasum_(integer *, doublereal *, integer *);
    extern /* Subroutine */ int daxpy_(integer *, doublereal *, doublereal *, 
	    integer *, doublereal *, integer *);
    integer jnext;
    doublereal sminw, xnorm;
    extern /* Subroutine */ int dlaln2_(logical *, integer *, integer *, 
	    doublereal *, doublereal *, doublereal *, integer *, doublereal *, 
	     doublereal *, doublereal *, integer *, doublereal *, doublereal *
, doublereal *, integer *, doublereal *, doublereal *, integer *);
    extern doublereal dlamch_(char *), dlange_(char *, integer *, 
	    integer *, doublereal *, integer *, doublereal *);
    extern integer idamax_(integer *, doublereal *, integer *);
    doublereal scaloc;
    extern /* Subroutine */ int dladiv_(doublereal *, doublereal *, 
	    doublereal *, doublereal *, doublereal *, doublereal *);
    doublereal bignum;
    logical notran;
    doublereal smlnum;


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

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

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

/*  DLAQTR solves the real quasi-triangular system */

/*               op(T)*p = scale*c,               if LREAL = .TRUE. */

/*  or the complex quasi-triangular systems */

/*             op(T + iB)*(p+iq) = scale*(c+id),  if LREAL = .FALSE. */

/*  in real arithmetic, where T is upper quasi-triangular. */
/*  If LREAL = .FALSE., then the first diagonal block of T must be */
/*  1 by 1, B is the specially structured matrix */

/*                 B = [ b(1) b(2) ... b(n) ] */
/*                     [       w            ] */
/*                     [           w        ] */
/*                     [              .     ] */
/*                     [                 w  ] */

/*  op(A) = A or A', A' denotes the conjugate transpose of */
/*  matrix A. */

/*  On input, X = [ c ].  On output, X = [ p ]. */
/*                [ d ]                  [ q ] */

/*  This subroutine is designed for the condition number estimation */
/*  in routine DTRSNA. */

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

/*  LTRAN   (input) LOGICAL */
/*          On entry, LTRAN specifies the option of conjugate transpose: */
/*             = .FALSE.,    op(T+i*B) = T+i*B, */
/*             = .TRUE.,     op(T+i*B) = (T+i*B)'. */

/*  LREAL   (input) LOGICAL */
/*          On entry, LREAL specifies the input matrix structure: */
/*             = .FALSE.,    the input is complex */
/*             = .TRUE.,     the input is real */

/*  N       (input) INTEGER */
/*          On entry, N specifies the order of T+i*B. N >= 0. */

/*  T       (input) DOUBLE PRECISION array, dimension (LDT,N) */
/*          On entry, T contains a matrix in Schur canonical form. */
/*          If LREAL = .FALSE., then the first diagonal block of T mu */
/*          be 1 by 1. */

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

//.........这里部分代码省略.........
    part of the matrix, using a symmetric rank-2k update of the form:   
    A := A - V*W' - W*V'.   

    The contents of A on exit are illustrated by the following examples   
    with n = 5 and nb = 2:   

    if UPLO = 'U':                       if UPLO = 'L':   

      (  a   a   a   v4  v5 )              (  d                  )   
      (      a   a   v4  v5 )              (  1   d              )   
      (          a   1   v5 )              (  v1  1   a          )   
      (              d   1  )              (  v1  v2  a   a      )   
      (                  d  )              (  v1  v2  a   a   a  )   

    where d denotes a diagonal element of the reduced matrix, a denotes   
    an element of the original matrix that is unchanged, and vi denotes   
    an element of the vector defining H(i).   

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


       Quick return if possible   

       Parameter adjustments */
    /* Table of constant values */
    static doublereal c_b5 = -1.;
    static doublereal c_b6 = 1.;
    static integer c__1 = 1;
    static doublereal c_b16 = 0.;
    
    /* System generated locals */
    integer a_dim1, a_offset, w_dim1, w_offset, i__1, i__2, i__3;
    /* Local variables */
    extern doublereal ddot_(integer *, doublereal *, integer *, doublereal *, 
	    integer *);
    static integer i__;
    static doublereal alpha;
    extern /* Subroutine */ HYPRE_Int dscal_(integer *, doublereal *, doublereal *, 
	    integer *);
    extern logical lsame_(const char *,const char *);
    extern /* Subroutine */ HYPRE_Int dgemv_(const char *, integer *, integer *, 
	    doublereal *, doublereal *, integer *, doublereal *, integer *, 
	    doublereal *, doublereal *, integer *), daxpy_(integer *, 
	    doublereal *, doublereal *, integer *, doublereal *, integer *), 
	    dsymv_(const char *, integer *, doublereal *, doublereal *, integer *, 
	    doublereal *, integer *, doublereal *, doublereal *, integer *), dlarfg_(integer *, doublereal *, doublereal *, integer *,
	     doublereal *);
    static integer iw;
#define a_ref(a_1,a_2) a[(a_2)*a_dim1 + a_1]
#define w_ref(a_1,a_2) w[(a_2)*w_dim1 + a_1]


    a_dim1 = *lda;
    a_offset = 1 + a_dim1 * 1;
    a -= a_offset;
    --e;
    --tau;
    w_dim1 = *ldw;
    w_offset = 1 + w_dim1 * 1;
    w -= w_offset;

    /* Function Body */
    if (*n <= 0) {
	return 0;
    }

Int TRON::trcg(double delta, double *g, double *s, double *r)
{
	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++)
	{
		s[i] = 0;
		r[i] = -g[i];
		d[i] = r[i];
	}
	cgtol = eps_cg*dnrm2_(&n, g, &inc);

	Int cg_iter = 0;
	rTr = ddot_(&n, r, &inc, r, &inc);
	while (1)
	{
		if (dnrm2_(&n, r, &inc) <= cgtol)
			break;
		cg_iter++;
		fun_obj->Hv(d, Hd);

		alpha = rTr/ddot_(&n, d, &inc, Hd, &inc);
		daxpy_(&n, &alpha, d, &inc, s, &inc);
		if (dnrm2_(&n, s, &inc) > delta)
		{
			info("cg reaches trust region boundary\n");
			alpha = -alpha;
			daxpy_(&n, &alpha, d, &inc, s, &inc);

			double std = ddot_(&n, s, &inc, d, &inc);
			double sts = ddot_(&n, s, &inc, s, &inc);
			double dtd = ddot_(&n, d, &inc, d, &inc);
			double dsq = delta*delta;
			double rad = sqrt(std*std + dtd*(dsq-sts));
			if (std >= 0)
				alpha = (dsq - sts)/(std + rad);
			else
				alpha = (rad - std)/dtd;
			daxpy_(&n, &alpha, d, &inc, s, &inc);
			alpha = -alpha;
			daxpy_(&n, &alpha, Hd, &inc, r, &inc);
			break;
		}
		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);
}

/* DECK DCGS */
/* Subroutine */ int dcgs_(integer *n, doublereal *b, doublereal *x, integer *
	nelt, integer *ia, integer *ja, doublereal *a, integer *isym, S_fp 
	matvec, S_fp msolve, integer *itol, doublereal *tol, integer *itmax, 
	integer *iter, doublereal *err, integer *ierr, integer *iunit, 
	doublereal *r__, doublereal *r0, doublereal *p, doublereal *q, 
	doublereal *u, doublereal *v1, doublereal *v2, doublereal *rwork, 
	integer *iwork)
{
    /* System generated locals */
    integer i__1, i__2;
    doublereal d__1;

    /* Local variables */
    static integer i__, k;
    static doublereal ak, bk, akm;
    extern doublereal ddot_(integer *, doublereal *, integer *, doublereal *, 
	    integer *);
    static doublereal bnrm, rhon, fuzz, sigma;
    extern /* Subroutine */ int daxpy_(integer *, doublereal *, doublereal *, 
	    integer *, doublereal *, integer *);
    extern doublereal d1mach_(integer *);
    static doublereal rhonm1;
    extern integer isdcgs_(integer *, doublereal *, doublereal *, integer *, 
	    integer *, integer *, doublereal *, integer *, S_fp, S_fp, 
	    integer *, doublereal *, integer *, integer *, doublereal *, 
	    integer *, integer *, doublereal *, doublereal *, doublereal *, 
	    doublereal *, doublereal *, doublereal *, doublereal *, 
	    doublereal *, integer *, doublereal *, doublereal *, doublereal *,
	     doublereal *);
    static doublereal tolmin, solnrm;

/* ***BEGIN PROLOGUE  DCGS */
/* ***PURPOSE  Preconditioned BiConjugate Gradient Squared Ax=b Solver. */
/*            Routine to solve a Non-Symmetric linear system  Ax = b */
/*            using the Preconditioned BiConjugate Gradient Squared */
/*            method. */
/* ***LIBRARY   SLATEC (SLAP) */
/* ***CATEGORY  D2A4, D2B4 */
/* ***TYPE      DOUBLE PRECISION (SCGS-S, DCGS-D) */
/* ***KEYWORDS  BICONJUGATE GRADIENT, ITERATIVE PRECONDITION, */
/*             NON-SYMMETRIC LINEAR SYSTEM, SLAP, SPARSE */
/* ***AUTHOR  Greenbaum, Anne, (Courant Institute) */
/*           Seager, Mark K., (LLNL) */
/*             Lawrence Livermore National Laboratory */
/*             PO BOX 808, L-60 */
/*             Livermore, CA 94550 (510) 423-3141 */
/*             [email protected] */
/* ***DESCRIPTION */

/* *Usage: */
/*      INTEGER N, NELT, IA(NELT), JA(NELT), ISYM, ITOL, ITMAX */
/*      INTEGER ITER, IERR, IUNIT, IWORK(USER DEFINED) */
/*      DOUBLE PRECISION B(N), X(N), A(NELT), TOL, ERR, R(N), R0(N), P(N) */
/*      DOUBLE PRECISION Q(N), U(N), V1(N), V2(N), RWORK(USER DEFINED) */
/*      EXTERNAL MATVEC, MSOLVE */

/*      CALL DCGS(N, B, X, NELT, IA, JA, A, ISYM, MATVEC, */
/*     $     MSOLVE, ITOL, TOL, ITMAX, ITER, ERR, IERR, IUNIT, */
/*     $     R, R0, P, Q, U, V1, V2, RWORK, IWORK) */

/* *Arguments: */
/* N      :IN       Integer */
/*         Order of the Matrix. */
/* B      :IN       Double Precision B(N). */
/*         Right-hand side vector. */
/* X      :INOUT    Double Precision X(N). */
/*         On input X is your initial guess for solution vector. */
/*         On output X is the final approximate solution. */
/* NELT   :IN       Integer. */
/*         Number of Non-Zeros stored in A. */
/* IA     :IN       Integer IA(NELT). */
/* JA     :IN       Integer JA(NELT). */
/* A      :IN       Double Precision A(NELT). */
/*         These arrays contain the matrix data structure for A. */
/*         It could take any form.  See "Description", below, */
/*         for more details. */
/* ISYM   :IN       Integer. */
/*         Flag to indicate symmetric storage format. */
/*         If ISYM=0, all non-zero entries of the matrix are stored. */
/*         If ISYM=1, the matrix is symmetric, and only the upper */
/*         or lower triangle of the matrix is stored. */
/* MATVEC :EXT      External. */
/*         Name of a routine which  performs the matrix vector multiply */
/*         operation  Y = A*X  given A and X.  The  name of  the MATVEC */
/*         routine must  be declared external  in the  calling program. */
/*         The calling sequence of MATVEC is: */
/*             CALL MATVEC( N, X, Y, NELT, IA, JA, A, ISYM ) */
/*         Where N is the number of unknowns, Y is the product A*X upon */
/*         return,  X is an input  vector.  NELT, IA,  JA,  A and  ISYM */
/*         define the SLAP matrix data structure: see Description,below. */
/* MSOLVE :EXT      External. */
/*         Name of a routine which solves a linear system MZ = R  for Z */
/*         given R with the preconditioning matrix M (M is supplied via */
/*         RWORK  and IWORK arrays).   The name  of  the MSOLVE routine */
/*         must be declared  external  in the  calling   program.   The */
/*         calling sequence of MSOLVE is: */
/*             CALL MSOLVE(N, R, Z, NELT, IA, JA, A, ISYM, RWORK, IWORK) */
/*         Where N is the number of unknowns, R is  the right-hand side */
/*         vector, and Z is the solution upon return.  NELT,  IA, JA, A */
//.........这里部分代码省略.........

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

    where tau is a real scalar, and v is a real vector with   
    v(1:i) = 0 and v(i+1) = 1; v(i+2:n) is stored on exit in A(i+2:n,i),   
    and tau in TAU(i).   

    The contents of A on exit are illustrated by the following examples   
    with n = 5:   

    if UPLO = 'U':                       if UPLO = 'L':   

      (  d   e   v2  v3  v4 )              (  d                  )   
      (      d   e   v3  v4 )              (  e   d              )   
      (          d   e   v4 )              (  v1  e   d          )   
      (              d   e  )              (  v1  v2  e   d      )   
      (                  d  )              (  v1  v2  v3  e   d  )   

    where d and e denote diagonal and off-diagonal elements of T, and vi   
    denotes an element of the vector defining H(i).   

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


       Test the input parameters   

       Parameter adjustments */
    /* Table of constant values */
    static integer c__1 = 1;
    static doublereal c_b8 = 0.;
    static doublereal c_b14 = -1.;
    
    /* System generated locals */
    integer a_dim1, a_offset, i__1, i__2, i__3;
    /* Local variables */
    extern doublereal ddot_(integer *, doublereal *, integer *, doublereal *, 
	    integer *);
    static doublereal taui;
    extern /* Subroutine */ int dsyr2_(char *, integer *, doublereal *, 
	    doublereal *, integer *, doublereal *, integer *, doublereal *, 
	    integer *);
    static integer i__;
    static doublereal alpha;
    extern logical lsame_(char *, char *);
    extern /* Subroutine */ int daxpy_(integer *, doublereal *, doublereal *, 
	    integer *, doublereal *, integer *);
    static logical upper;
    extern /* Subroutine */ int dsymv_(char *, integer *, doublereal *, 
	    doublereal *, integer *, doublereal *, integer *, doublereal *, 
	    doublereal *, integer *), dlarfg_(integer *, doublereal *,
	     doublereal *, integer *, doublereal *), xerbla_(char *, integer *
	    );
#define a_ref(a_1,a_2) a[(a_2)*a_dim1 + a_1]


    a_dim1 = *lda;
    a_offset = 1 + a_dim1 * 1;
    a -= a_offset;
    --d__;
    --e;
    --tau;

    /* Function Body */
    *info = 0;
    upper = lsame_(uplo, "U");
    if (! upper && ! lsame_(uplo, "L")) {
	*info = -1;
    } else if (*n < 0) {

/* Subroutine */ int dpotf2_(char *uplo, integer *n, doublereal *a, integer *
                             lda, integer *info)
{
    /* System generated locals */
    integer a_dim1, a_offset, i__1, i__2, i__3;
    doublereal d__1;

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

    /* Local variables */
    integer j;
    doublereal ajj;
    extern doublereal ddot_(integer *, doublereal *, integer *, doublereal *,
                            integer *);
    extern /* Subroutine */ int dscal_(integer *, doublereal *, doublereal *,
                                       integer *);
    extern logical lsame_(char *, char *);
    extern /* Subroutine */ int dgemv_(char *, integer *, integer *,
                                       doublereal *, doublereal *, integer *, doublereal *, integer *,
                                       doublereal *, doublereal *, integer *);
    logical upper;
    extern logical disnan_(doublereal *);
    extern /* Subroutine */ int xerbla_(char *, integer *);


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

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

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

    /*  DPOTF2 computes the Cholesky factorization of a real symmetric */
    /*  positive definite matrix A. */

    /*  The factorization has the form */
    /*     A = U' * U ,  if UPLO = 'U', or */
    /*     A = L  * L',  if UPLO = 'L', */
    /*  where U is an upper triangular matrix and L is lower triangular. */

    /*  This is the unblocked version of the algorithm, calling Level 2 BLAS. */

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

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

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

    /*  A       (input/output) DOUBLE PRECISION array, dimension (LDA,N) */
    /*          On entry, the symmetric matrix A.  If UPLO = 'U', the leading */
    /*          n by n upper triangular part of A contains the upper */
    /*          triangular part of the matrix A, and the strictly lower */
    /*          triangular part of A is not referenced.  If UPLO = 'L', the */
    /*          leading n by n lower triangular part of A contains the lower */
    /*          triangular part of the matrix A, and the strictly upper */
    /*          triangular part of A is not referenced. */

    /*          On exit, if INFO = 0, the factor U or L from the Cholesky */
    /*          factorization A = U'*U  or A = L*L'. */

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

    /*  INFO    (output) INTEGER */
    /*          = 0: successful exit */
    /*          < 0: if INFO = -k, the k-th argument had an illegal value */
    /*          > 0: if INFO = k, the leading minor of order k is not */
    /*               positive definite, and the factorization could not be */
    /*               completed. */

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

    /*     .. Parameters .. */
    /*     .. */
    /*     .. Local Scalars .. */
    /*     .. */
    /*     .. External Functions .. */
    /*     .. */
    /*     .. External Subroutines .. */
    /*     .. */
    /*     .. Intrinsic Functions .. */
    /*     .. */
    /*     .. Executable Statements .. */

    /*     Test the input parameters. */

    /* Parameter adjustments */
    a_dim1 = *lda;
//.........这里部分代码省略.........

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, mi