C++ DoubleMatrix::nr方法代码示例

static DoubleMatrix funExactLambda_b(const DoubleMatrix &dvecAlp, 
                                     const DoubleMatrix &dvecB, 
                                     const DoubleMatrix &dvecY,
                                     const bool withD){
	
	// Updated 1-8-01 Alyssa
	// fixed for const correctness
    DoubleMatrix term1(1, dvecB.nr()), term2(1, dvecB.nr());
    const double *pdAlp   = dvecAlp.data();
    const double *pdB     = dvecB.data();
    const double *pdY     = dvecY.data();
    double *pdTerm1 = term1.data();
    double *pdTerm2 = term2.data();

    term1.fill(0.0);
    term2.fill(0.0);

    pdTerm1[0] = 1.0/pdB[0]  
                -0.5 * (pow(pdB[0],-2.0)*pow((pdY[0]-pdAlp[1]-pdB[1]),2.0)
                        +pow(pdB[0],-2.0)*pow((pdY[1]-pdAlp[1]-pdB[1]),2.0));
    pdTerm1[1] = -(pdY[0]+pdY[1]-2.0*pdAlp[1]-2.0*pdB[1])/pdB[0] ;

    if( withD ){
        //pdAns[0] = 1.0/pdB[0] - pow((1.0 - pdAlp[1] - pdB[1]*pdB[1]), 2.0) + pdB[0]/pdAlp[0];
        //pdAns[1] = -2.0 * (1.0 - pdAlp[1] - pdB[1]) / pdB[0] + pdB[1]/pdAlp[0];
        pdTerm2[0] = pdB[0]/pdAlp[0];
        pdTerm2[1] = pdB[1]/pdAlp[0];
        return add( term1, term2 );
    }
    else{
        return term1;
    }
}

void symmetrize(const DoubleMatrix& L, DoubleMatrix& Sym)
{
    int m = L.nr();
    int n = L.nc();
    assert(m==n);
    assert(Sym.nr() == m);
    assert(Sym.nc() == n);

    const double* pL = L.data();
    double*       pS = Sym.data();

    if( &L != &Sym )
    {
        std::copy(pL, pL+m*n, pS);
    }
    int j,i;
    for( j=0; j<n; j++ )
    {
        // Copy the lower triangle data into the resulting matrix's upper triangle
        for( i=0; i<j; i++ )
        {
            pS[i+j*m] = pL[j+i*m];
        }        
    }
}

void AkronItimesCTest::testAkronItimesCVal()
{
    const int m = 2;
    const int n = 3;
    const int k = 1;
    int i;
    double seq[] = {0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25,26,27,28,29,30};
    DoubleMatrix A(m,n);
    DoubleMatrix I = identity(n);
    DoubleMatrix C(n*n,k);

    std::copy(seq, seq+m*n, A.data());
    std::copy(seq, seq+n*n*k, C.data());

    DoubleMatrix AIC = AkronItimesC(A,I,C);
    CPPUNIT_ASSERT_EQUAL(m*n, AIC.nr());
    CPPUNIT_ASSERT_EQUAL(k, AIC.nc());

    DoubleMatrix ABC = AkronBtimesC(A,I,C);
    CPPUNIT_ASSERT_EQUAL(m*n, ABC.nr());
    CPPUNIT_ASSERT_EQUAL(k, ABC.nc());
    //std::cout << "ABC=" << ABC << std::endl;
    //std::cout << "AIC=" << AIC << std::endl;
    for( i = 0; i < m*n*k; i++ )
    {
        CPPUNIT_ASSERT_EQUAL(ABC.data()[i], AIC.data()[i]);
    }

}

void multiply(const DoubleMatrix &X, const DoubleMatrix &Y, DoubleMatrix& Z)
{
  assert(&X!=&Z);
  assert(&Y!=&Z);
  if( X.isEmpty() || Y.isEmpty() )
    return;

  //
  // Row Major Matrices
  //
  //    C = alpha * A * B + beta * C  --- (1)
  //
  // A : m by k         lda (stride) = k
  // B : k by n         ldb (stride) = n
  // C : m by n         ldc (stride) = n
  //
  // Column Major Matrices
  //
  //    Z = alpha * X * Y + beta * C  --- (2)
  //    Z = C^t 
  //      = alpha * B^t * A^t + beta * C^t  --- (3)
  //
  // X = B^t : n by k   ldx (stride) = n
  // Y = A^t : k by m   ldy (stride) = k
  // Z = C^t : n by m   ldz (stride) = n
  //
  int m = Y.nc();
  int k = X.nc();
  int n = X.nr();

  Z.resize( n, m );
  
  assert( X.nr() == n );
  assert( X.nc() == k );
  assert( Y.nr() == k );
  assert( Y.nc() == m );
  assert( Z.nr() == n );
  assert( Z.nc() == m );
  
  const double * pX = X.data();
  const double * pY = Y.data();
  double       * pZ = Z.data(); 
  
  int lda = n;
  int ldb = k;
  int ldc = n;

  cblas_dgemm( CblasColMajor, CblasNoTrans, CblasNoTrans, n, m, k, ALPHA, pX, lda, pY, ldb , BETA, pZ, ldc );

}

static void replaceIth( DoubleMatrix &target, int ith, const DoubleMatrix &data )
{
    assert(target.nc() == data.nc());
    assert( ith >= 0 && ith < target.nr());

    double *pTarget = target.data();
    const double *pData   = data.data();
    int j;

    for( j=0; j<target.nc(); j++ )
    {
        pTarget[j*target.nr()+ith] = pData[j];
    }
    return;
}

void lambdaTest::test(
    SpkModel<double> &model, 
    const DoubleMatrix &dvecY, 
    const DoubleMatrix &dvecAlp, 
    const DoubleMatrix &dvecB, 
    double lambdaOut, 
    DoubleMatrix &lambda_alpOut, 
    DoubleMatrix &lambda_bOut,
    const bool withD)
{
    DoubleMatrix exactTemp;

    double       dblexactTemp;
    int i;

    lambdaOut     = lambda(model, dvecY, dvecAlp, dvecB, withD);
    lambda_alpOut = lambda_alp(model, dvecY, dvecAlp, dvecB, withD);
    lambda_bOut   = lambda_b(model, dvecY, dvecAlp, dvecB, withD);

    dblexactTemp = funExactLambda(dvecAlp, dvecB, dvecY, withD);
    CPPUNIT_ASSERT_DOUBLES_EQUAL( dblexactTemp, lambdaOut, 0.0001);

    exactTemp = funExactLambda_alp(dvecAlp, dvecB, dvecY, withD);
    for( i=0; i<exactTemp.nr()*exactTemp.nc(); i++ )
    {
        CPPUNIT_ASSERT_DOUBLES_EQUAL(
            exactTemp.data()[i], 
            lambda_alpOut.data()[i], 0.0001
            );
    }
    
    exactTemp = funExactLambda_b(dvecAlp, dvecB, dvecY, withD);
    for( i=0; i<exactTemp.nr()*exactTemp.nc(); i++ )
    {
        CPPUNIT_ASSERT_DOUBLES_EQUAL( 
            exactTemp.data()[i],
            lambda_bOut.data()[i], 0.0001
            );
    }
}

static double compareAgainst(
            const DoubleMatrix &dvecZ,     // column vector
            const DoubleMatrix &dvecH,     // column vector
            const DoubleMatrix &dmatQ,     // symmetric matrix Q(x)
            const DoubleMatrix &dmatInvQ   // inverse of Q
            )
{

    assert( dvecZ.nc() == 1 );
    assert( dvecH.nc() == 1 );
    assert( dmatQ.nr() == dmatQ.nc() );
    assert( hasPosDet( dmatQ ) );

    DoubleMatrix dmatQ2PI( dmatQ.nr(), dmatQ.nc() );
    
	const double *pQ    = dmatQ.data();
    double *pQ2PI = dmatQ2PI.data();

    for( int i=0; i<dmatQ.nr()*dmatQ.nc(); i++ )
        pQ2PI[i] = pQ[i] * 2.0 * PI;

    // Compute det(Q2PI) = b * 2^c .
    double b;
    long c;
    det( dmatQ2PI , &b, &c );

    // Compute log(det(Q2PI)).
    double dTerm1 = log( b ) + c * log( 2.0 );

    DoubleMatrix  dmatR = subtract( dvecZ, dvecH );

    DoubleMatrix  term2 = multiply(multiply(transpose(dmatR),dmatInvQ), dmatR);
    double        dTerm2 = oneByOneToScalar(term2);

    return (dTerm1 + dTerm2) / 2.0;
}

static DoubleMatrix elsq_x_x( 
    const DoubleMatrix &dvecX,          // parameter
    const DoubleMatrix &dvecZ,          // data
    const DoubleMatrix &dvecH,          // h(x)
    const DoubleMatrix &dmatH_x,    
    const DoubleMatrix &dmatH_x_x,
    const DoubleMatrix &dmatQ,          // Q(x)
    const DoubleMatrix &dmatQ_x,
    const DoubleMatrix &dmatQinv,       // Q(x)inv
    const DoubleMatrix &dmatQinv_x,
    const DoubleMatrix &dmatQinv_x_x
    )
{
    using namespace std;

    DoubleMatrix    term1, term2, term1sub2div2, term3, term4, term5, term4add5, term6, term6div2, 
                    term7, term8, term9, term7add8div2, term7add8div2Trans;
    DoubleMatrix ZsubH, tranZsubH, identX;
    DoubleMatrix dmatElsq_x_x;

    ZsubH     = subtract(dvecZ, dvecH);
    tranZsubH = transpose(ZsubH);
    identX    = identity(dvecX.nr());

    term1 = AkronBtimesC(transpose(rvec(dmatQ)), identX, dmatQinv_x_x);
    term2 = multiply(transpose(dmatQinv_x), dmatQ_x);
    term1sub2div2 = divByScalar( add(term1, term2), -2.0 );

    term3 = mulByScalar( AkronBtimesC(multiply(tranZsubH,dmatQinv), identX, dmatH_x_x), -1.0);

    term4 = AkronBtimesC(tranZsubH, transpose(dmatH_x), dmatQinv_x);
    term5 = multiply(multiply(transpose(dmatH_x),dmatQinv),dmatH_x);
    term4add5 = subtract( term5, term4 );

    term6 = AkronBtimesC( transpose(rvec(multiply(ZsubH,tranZsubH))), identX, dmatQinv_x_x );
    term6div2 = divByScalar(term6, 2.0);


    term7 = AkronBtimesC( transpose(dmatH_x), transpose(ZsubH), dmatQinv_x );
    term8 = AkronBtimesC( transpose(ZsubH), transpose(dmatH_x), dmatQinv_x );
    term7add8div2 = divByScalar( add(term7, term8), -2.0 );
    term7add8div2Trans = transpose(term7add8div2);
    
    dmatElsq_x_x = add( term1sub2div2, add( term3, add(term4add5, add(term6div2, term7add8div2Trans))));
    
    return dmatElsq_x_x;
}

static bool check(DoubleMatrix dmatB,
                  int rows,
                  int cols,
                  int d1, int d2, int d3, 
                  int d4, int d5, int d6,
                  int d7, int d8, int d9)
{
    using namespace std;
    bool isOkay = true;
    int m = dmatB.nr();
    int n = dmatB.nc();
    double *pdB = dmatB.data();

    if( m != rows || n != n ){
        isOkay = false;
        return isOkay;
    }

    int d[9];
    d[0] = d1;
    d[1] = d2;
    d[2] = d3;
    d[3] = d4;
    d[4] = d5;
    d[5] = d6;
    d[6] = d7;
    d[7] = d8;
    d[8] = d9;
    
    for( int i=0; i<rows*cols && isOkay; i++ ){
        if( d[i] != pdB[i] ){
            cout << "d[" << i << "] received was " << d[i] << endl;
            cout << "mat[" << i << "] was " << pdB[i] << endl;
            dmatB.print();
            isOkay = false;
        }
    }
    return isOkay;
}

const DoubleMatrix backDiv(const DoubleMatrix &dmatA, const DoubleMatrix &dmatB)
{
  // A is assumed to be square.
  int m = dmatA.nr();
  int n = dmatA.nc();
  assert( m == n );

  // B is m by l matrix, where l is the number of right hand sides.
  int l = dmatB.nc();
  assert( dmatB.nr() == m );

  if( dmatA.isEmpty() || dmatB.isEmpty() )
    return DoubleMatrix( 0, 0 );

  //==============================================================
  // First decompose A into LU such that A = P * L * U, 
  // where P is the permutation matrix,
  // L is the lower triangle and the U the upper triangle.
  //
  // We use CLAPACK's DGETRF() which does LU decomposition
  // with partial (ie. row interchanges only) pivoting.
  //==============================================================
  
  // enum CBLAS_ORDER order =: (CblasColMajor | CblasRowMajor)
  //
  // If order = CblasColMajor, the array, a, is assumed to 
  // hold each matrix A's column in the contiguous manner
  // in memory (ie. A is said to be in the column major order).
  // If order = CblasRowMajor, the array, a, is assumed to
  // hold each matrix A's row in the contiguous manner
  // in memory (ie. A is said to be in the row major order).
  enum CBLAS_ORDER order = CblasColMajor;

  // double *a
  //
  // (on entry) a points to the elements of matrix A(m,n) 
  // in the column major order if "order" = CblasColMajor, 
  // or in the row major order if "order" = CblasRowMajor.
  //
  // (on exit) The lower triangle (j<=i) is replaced by L
  // and the upper triangle (j>i) is replaced by U.
  double a[m*n];
  copy( dmatA.data(), dmatA.data()+m*n, a );

  // int lda
  // 
  // The leading dimension of A.  
  // If A is in the column major order, lda = m.
  // If A is in the row major order, lda = n.
  int lda = m;

  // int ipiv(m)
  //
  // (on exit) The i-th row in A was interchanged with the row
  // indicated by the value in ipiv[i].
  int ipiv[m];

  int info = clapack_dgetrf( order, m, n, a, lda, ipiv );
  if( info < 0 )
    {
      char mess[ SpkError::maxMessageLen() ];
      snprintf( mess, SpkError::maxMessageLen(), "Solution of a system of linear equations using the LU decomposition failed: \n the %s argument to the function that performs the LU decomposition  had an illegal value.", 
               intToOrdinalString( -info, ONE_IS_FIRST_INT ).c_str() );
      throw SpkException( SpkError::SPK_UNKNOWN_ERR, mess, __LINE__, __FILE__ );
    }
  else if( info > 0 )
    {
      char mess[ SpkError::maxMessageLen() ];
      snprintf( mess, SpkError::maxMessageLen(), "Solution of a system of linear equations using the LU decomposition failed: \nthe %s diagonal element of U is exactly zero.", 
               intToOrdinalString( info, ONE_IS_FIRST_INT ).c_str() );
      throw SpkException( SpkError::SPK_NOT_POS_DEF_ERR, mess, __LINE__, __FILE__ );
    }

  //==============================================================
  // Solve A x = B for x using the LU computed in the previous
  // step.
  // Note that A is now assumed to be square: m = n.
  //==============================================================
  
  // int rhs
  //
  // The number of right hand sides (ie. the number of columns of B).
  int nrhs = l;

  // int ldb
  // The leading dimension of B.
  // If B is in the column major order, ldb = m.
  // If B is in the row major order, ldb = l.
  int ldb = m;

  // double *x
  //
  // (on entry) x points to the elements of B in the column major
  // order if "order" = CblasColMajor or in the row major otherwise.
  // (on exit) x points to the solution matrix, x (m=n by l).
  DoubleMatrix X( dmatB );
  double * x = X.data();// This points to b on entry and contains the solution x upon exit

  info = clapack_dgetrs( order, CblasNoTrans, n, nrhs, a, lda, ipiv, x, ldb );
  if( info < 0 )
//.........这里部分代码省略.........

void ppdOpt( SpdModel&             model,
             Optimizer&            optInfo,
             const DoubleMatrix&   dvecXLow,
             const DoubleMatrix&   dvecXUp,
             const DoubleMatrix&   dvecXIn,
             DoubleMatrix*         pdvecXOut,
             const DoubleMatrix&   dvecXStep,
             const DoubleMatrix&   dvecAlp,
             const DoubleMatrix&   dvecAlpStep,
             double*               pdPhiTildeOut,
             DoubleMatrix*         pdrowPhiTilde_xOut,
             DoubleMatrix*         pdmatPhiTilde_x_xOut )
{
  //------------------------------------------------------------
  // Preliminaries.
  //------------------------------------------------------------

  using namespace std;

  // If no evaluation is requested, return immediately.
  if ( !pdvecXOut && 
       !pdPhiTildeOut &&
       !pdrowPhiTilde_xOut &&
       !pdmatPhiTilde_x_xOut )
  {
    return;
  }

  // Get the number of design parameters.
  int nX = dvecXIn.nr();


  //------------------------------------------------------------
  // Validate the inputs (Debug mode).
  //------------------------------------------------------------

  // Check the length of the design parameter vector.
  assert( nX == model.getNDesPar() );


  //------------------------------------------------------------
  // Validate the inputs (All modes).
  //------------------------------------------------------------

  // [Revisit - Warm Start Disabled - Mitch]
  // In order to simplify the testing of this function,
  // warm starts have been disabled for now.
  if ( optInfo.getIsWarmStart() )
  {
    throw SpkException(
      SpkError::SPK_USER_INPUT_ERR,
      "The input Optimizer object requested a warm start, which have been disabled for now.",
      __LINE__,
      __FILE__ );
  }


  //------------------------------------------------------------
  // Prepare the values required to calculate the output values.
  //------------------------------------------------------------

  // If this value is required to calculate the output values, 
  // initialize it so that it wil be calculated.  Otherwise, 
  // set its pointer to zero so that it won't be calculated.
  DoubleMatrix dvecXOutTemp;
  DoubleMatrix* pdvecXOutTemp = &dvecXOutTemp;
  if ( pdvecXOut || pdmatPhiTilde_x_xOut )
  {
    dvecXOutTemp.resize( nX, 1 );
  }
  else
  {
    pdvecXOutTemp = 0;
  }
  
  // If this value is required to calculate the output values, 
  // initialize it so that it wil be calculated.  Otherwise, 
  // set its pointer to zero so that it won't be calculated.
  double dNegLogOut;
  double* pdNegLogOut = &dNegLogOut;
  if ( pdPhiTildeOut || pdrowPhiTilde_xOut || pdmatPhiTilde_x_xOut )
  {
    dNegLogOut = 0.0;
  }
  else
  {
    pdNegLogOut = 0;
  }

  // If this value is required to calculate the output values, 
  // initialize it so that it wil be calculated.  Otherwise, 
  // set its pointer to zero so that it won't be calculated.
  DoubleMatrix drowNegLog_xOut;
  DoubleMatrix* pdrowNegLog_xOut = &drowNegLog_xOut;
  if ( pdrowPhiTilde_xOut || pdmatPhiTilde_x_xOut )
  {
    drowNegLog_xOut.resize( 1, nX );
  }
  else
  {
//.........这里部分代码省略.........

extern void spk_non_par(
	size_t                             level       ,
	SpkModel< CppAD::AD<double> >     &admodel     ,
	SpkModel<double>                  &model       ,
	const DoubleMatrix                &N           ,
	const DoubleMatrix                &y           ,
	const DoubleMatrix                &max_itr     ,
	const DoubleMatrix                &epsilon     ,
	const DoubleMatrix                &blow        ,
	const DoubleMatrix                &bup         ,
	const DoubleMatrix                &Bin         ,
	DoubleMatrix                      &Bout        ,
	DoubleMatrix                      &lamout      ,
	DoubleMatrix                      &Pout  )
{
	// temporary indices
	size_t i, j, k;

	// temporary double pointer
	double       *ptr;
	const double *ptr_c;

	// number of discrete measure points
	size_t J = Bin.nc();

	// number of random effects
	size_t n = blow.nr();
	
	// ------------ Arguments to non_par::opt_measure --------------------
	assert( max_itr.nr()  == 2 );
	mat2cpp::matrix<size_t> maxitr(2, 1);
	maxitr(0, 0) = size_t( *(max_itr.data() + 0) );
	maxitr(1, 0) = size_t( *(max_itr.data() + 1) );

	assert( epsilon.nr()  == 5 );
	mat2cpp::matrix<double> eps(5, 1);
	eps(0, 0)    = *(epsilon.data() + 0);
	eps(1, 0)    = *(epsilon.data() + 1);
	eps(2, 0)    = *(epsilon.data() + 2);
	eps(3, 0)    = *(epsilon.data() + 3);
	eps(4, 0)    = *(epsilon.data() + 4);

	// input likelihood function
	Like like(admodel, model, N, y, n);

	// input number of individuals in the population
	size_t M = N.nr();

	// input lower limit for the random effects
	mat2cpp::matrix<double> xLow(1, n);
	ptr_c = blow.data();
	for(k = 0; k < n; k++)
		xLow(0, k) = ptr_c[k];

	// input upper limit for the random effects
	mat2cpp::matrix<double> xUp(1, n);
	ptr_c = bup.data();
	for(k = 0; k < n; k++)
		xUp(0, k) = ptr_c[k];

	// input and return discrete measure points
	mat2cpp::matrix<double> X(J, n);
	ptr_c = Bin.data();
	for(j = 0; j < J; j++)
	{	for(k = 0; k < n; k++)
			X(j, k) = ptr_c[k + j * n];
	}

	// return weight corresponding to each measure oint
	mat2cpp::matrix<double> lambda(J, 1);

	// return convergence information
	mat2cpp::matrix<double> info; 

	// return status message
	const char *msg;

	// -----------------------------------------------------------------
	msg = non_par::opt_measure( level, maxitr, eps, 
		&like, M, xLow, xUp, X, lambda, info
	);
	// -----------------------------------------------------------------

	if( strcmp(msg, "ok") != 0 )
	{	throw SpkException(
      			SpkError::SPK_NON_PAR_ERR,
			msg,
      			__LINE__,
      			__FILE__ 
		);
	}
	// determine number of discrete measure points
	assert( n == X.size2() );
	J = X.size1();

	// dimension the return matrices
	Bout.resize(n, J);
	lamout.resize(J, 1);
	Pout.resize(M, J);

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

static double oneByOneToScalar( const DoubleMatrix &dmatA ){
    assert( dmatA.nr() == 1 );
    assert( dmatA.nc() == 1 );

    return (dmatA.data())[0];
}

const DoubleMatrix elsq_x(const DoubleMatrix &dvecR,      // m size vector, z - h
                    const DoubleMatrix &dmatQ,      // m by m symmetric, positive definite
                    const DoubleMatrix &dmatQinv,   // m by m symmetric, positive definite
                    const DoubleMatrix &dmatH_x,    // m by n matrix
                    const DoubleMatrix &dmatQ_x     // m*m by n matrix
                    )
{
    using namespace std;
    // x:   m by 1
    // z:   m by 1
    // h:   m by 1
    // Q:   m by m
    // invQ:m by m
    // h_x: m by n
    // Q_x: m*m by n
    const int m = dvecR.nr();
    const int n = dmatQ_x.nc();
    double  val;

    assert( dmatQ.nr() == m );
    assert( dmatQ.nc() == m );
    assert( dmatQinv.nr() == m );
    assert( dmatQinv.nc() == m );
    assert( dmatH_x.nr() == m );
    assert( dmatH_x.nc() == n );
    assert( dmatQ_x.nr() == m*m );
    assert( dmatQ_x.nc() == n );

    transpose(dvecR, dvecRTrans);
    multiply(dvecRTrans, dmatQinv, drowW);
    transpose( rvec( dmatQinv ), rvecQinvTrans );

    multiply(rvecQinvTrans, dmatQ_x, drowTerm1);
    assert(drowTerm1.nr()== 1);
    assert(drowTerm1.nc()== n);

    multiply(drowW, dmatH_x, drowTerm2);
    assert(drowTerm2.nr()== 1);
    assert(drowTerm2.nc()== n);

    //
    // The loop below untangles the following statement
    // for eliminating matrix object constructions and taking advantage of
    // Q being symmetric.
    //
    //    DoubleMatrix dmatTerm3  = AkronBtimesC( drowW, rvecQinvTrans, dmatQ_x );
    // 
    // Let w be (z-h).
    // 
    // The orignal formula
    //    0.5 [w^T kron w^T] partial_x(Qinv) --- (eq. 1)
    //
    // is equivalent to:
    //    0.5 partial_x(k) [ w^T Qinv w ] --- (eq. 2)
    //   =0.5 SUM [ w(i) w(j) partial_x(k)(Q(i,j)) ], over 1<=i<=m and 1<=j<=m.
    //
    // For Q being symmetric, the following is true:
    //    w(i) w(j) partial_x(k)(Q(i,j)) == w(j) w(i) partial_x(k)(Q(j,i))
    //
    //
    const double* pW   = drowW.data();
    const double* pQ_x = dmatQ_x.data();
    drowTerm3.resize(1,n);
    drowTerm3.fill(0.0);
    double* pTerm3 = drowTerm3.data();

    for( int k=0; k<n; k++ )
    {
        for( int j=0; j<m; j++ )
        {
            for( int i=0; i<m; i++ )
            {
                if( i<=j )
                {
                    val = pW[i] * pW[j] * pQ_x[ j*m+i+k*m*m ] / 2.0;
                    drowTerm3.data()[k] += val;
                    if( i<j )
                        pTerm3[k] += val;
                }
            }
        }
    }
    
    subtract( subtract(mulByScalar(drowTerm1, 0.5), drowTerm2 )
                                        ,drowTerm3, drowAns );


    // elsq_x returns a 1 by n matrix
    assert(drowAns.nr()==1);
    assert(drowAns.nc()==n);
    return drowAns;

}

void ppedOpt( SpdModel&             model,
              Optimizer&            xOptInfo,
              const DoubleMatrix&   dvecXLow,
              const DoubleMatrix&   dvecXUp,
              const DoubleMatrix&   dvecXIn,
              DoubleMatrix*         pdvecXOut,
              const DoubleMatrix&   dvecXStep,
              Optimizer&            alpOptInfo,
              const DoubleMatrix&   dvecAlpLow,
              const DoubleMatrix&   dvecAlpUp,
              const DoubleMatrix&   dvecAlpIn,
              DoubleMatrix*         pdvecAlpOut,
              const DoubleMatrix&   dvecAlpStep,
              double*               pdPhiTildeOut,
              DoubleMatrix*         pdrowPhiTilde_xOut,
              DoubleMatrix*         pdmatPhiTilde_x_xOut )
{
  //------------------------------------------------------------
  // Preliminaries.
  //------------------------------------------------------------

  using namespace std;

  // If no evaluation is requested, return immediately.
  if ( !pdvecXOut && 
       !pdPhiTildeOut &&
       !pdrowPhiTilde_xOut &&
       !pdmatPhiTilde_x_xOut )
  {
    return;
  }

  // Get the number of design and fixed population parameters.
  int nX   = dvecXIn.nr();
  int nAlp = dvecAlpIn.nr();


  //------------------------------------------------------------
  // Validate the inputs (Debug mode).
  //------------------------------------------------------------

  // Check the length of the design parameter vector.
  assert( nX == model.getNDesPar() );


  //------------------------------------------------------------
  // Validate the inputs (All modes).
  //------------------------------------------------------------

  // [Revisit - Warm Start Disabled - Mitch]
  // In order to simplify the testing of this function,
  // warm starts have been disabled for now.
  if ( xOptInfo.getIsWarmStart() || alpOptInfo.getIsWarmStart() )
  {
    throw SpkException(
      SpkError::SPK_USER_INPUT_ERR,
      "One of the input Optimizer objects requested a warm start, which have been disabled for now.",
      __LINE__,
      __FILE__ );
  }


  //------------------------------------------------------------
  // Prepare the values required to calculate the output values.
  //------------------------------------------------------------

  // If this value is required to calculate the output values, 
  // initialize it so that it wil be calculated.  Otherwise, 
  // set its pointer to zero so that it won't be calculated.
  DoubleMatrix dvecXOutTemp;
  DoubleMatrix* pdvecXOutTemp = &dvecXOutTemp;
  if ( pdvecXOut )
  {
    dvecXOutTemp.resize( nX, 1 );
  }
  else
  {
    pdvecXOutTemp = 0;
  }

  // If this value is required to calculate the output values, 
  // initialize it so that it wil be calculated.  Otherwise, 
  // set its pointer to zero so that it won't be calculated.
  DoubleMatrix dvecAlpOutTemp;
  DoubleMatrix* pdvecAlpOutTemp = &dvecAlpOutTemp;
  if ( pdvecAlpOut )
  {
    dvecAlpOutTemp.resize( nAlp, 1 );
  }
  else
  {
    pdvecAlpOutTemp = 0;
  }

  // If this value is required to calculate the output values, 
  // initialize it so that it wil be calculated.  Otherwise, 
  // set its pointer to zero so that it won't be calculated.
  double dLTildeOut;
  double* pdLTildeOut = &dLTildeOut;
  if ( pdPhiTildeOut || pdrowPhiTilde_xOut || pdmatPhiTilde_x_xOut )
//.........这里部分代码省略.........