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

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

}

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

/*------------------------------------------------------------------------
 * Function definition
 *------------------------------------------------------------------------*/
 static DoubleMatrix getElements( const DoubleMatrix &original, int offset, int howMany )
 {
     DoubleMatrix dvecPart(howMany, 1);
     const double *pdOriginal = original.data();
     double *pdPart     = dvecPart.data();
     for( int i=0; i<howMany; i++ )
         pdPart[i] = pdOriginal[offset+i];
     return dvecPart;
 }

    void doDataMean( valarray<double>& ret ) const
    {
      //
      // h(x) = [ x(3) ]
      //        [ x(3) ]
      //
      ret.resize(_nY);
      const double *pX = _b.data();

      ret[0] = 1.0 * pX[2];
      ret[1] = 1.0 * pX[2];
    }

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) Exact Lambda, Labmda_alp, and Lambda_b 
 **************************************************************/
static double funExactLambda(const DoubleMatrix &dvecAlp, 
                             const DoubleMatrix &dvecB, 
                             const DoubleMatrix &dvecY,
                             const bool withD){
	// Updated 1-8-01 Alyssa
	// fixed for const correctness
    const double *pdAlp = dvecAlp.data();
    const double *pdB   = dvecB.data();
    const double *pdY   = dvecY.data();
    double term1, term2;

    term1 = 0.5 * log( pow( (2.0*PI*pdB[0]), 2.0 ) )
          + 0.5 * (pow(pdY[0]-pdAlp[1]-pdB[1], 2.0) + pow(pdY[1]-pdAlp[1]-pdB[1], 2.0))/pdB[0];

    if( withD ){
        term2 = 0.5 * log( pow(2.0*PI*pdAlp[0],2.0) )
              + 0.5 * (pow(pdB[0], 2.0) + pow(pdB[1], 2.0))/pdAlp[0];
        return term1 + term2;
    }
    else{
        return term1;
    }
}

    void doDataVarianceInv( valarray<double>& ret ) const
    {
      //
      // Q(x)^-1 = [ 1.0 / (2.0*x(1)),   0.0        ]
      //           [ 0.0                 1.0 / x(2) ]
      //
      ret.resize(_nY * _nY);
      const double *x = _b.data();

      ret[0] = 1.0 / (2.0 * x[0]);
      ret[1] = 0.0;
      ret[2] = 0.0;
      ret[3] = 1.0 / x[1];
    }

    void doDataVariance( valarray<double>& ret ) const
    {
      //
      // Q(x) = [ 2 * x(1) ,   0   ]
      //        [     0    ,  x(2) ]
      //
      
      ret.resize(_nY * _nY);
      const double *x = _b.data();

      ret[0] = 2.0 * x[0];
      ret[1] = 0.0;
      ret[2] = 0.0;
      ret[3] = x[1];
    }

    bool doDataVarianceInv_indPar( valarray<double>& ret ) const
    {
      //
      // Q^-1_x(x) = [ -1.0 / (2.0 * x(1))^2  0.0             0.0    ]
      //             [ 0.0                    0.0             0.0    ]
      //             [ 0.0                    0.0             0.0    ]
      //             [ 0.0                    -1.0 / x(2)^2   0.0    ]
      //
      ret.resize(_nY*_nY*_nB);
      const double * x = _b.data();
      for( int i=0; i<_nY*_nY*_nB; i++ )
        ret[i] = 0.0;

      ret[0] = -1.0 / ( (2.0 * x[0])*(2.0 * x[0]) );
      ret[7] = -1.0 / (x[1] * x[1]);
      return true;
    }

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

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

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 )
//.........这里部分代码省略.........