C++ rhs函数代码示例

//.........这里部分代码省略.........
    boost::dynamic_bitset<Block> b(std::string("1010"));
    Tests::operator_shift_left(b, pos);
  }
  { // case pos == size()/2
    std::size_t pos = long_string.size() / 2;
    boost::dynamic_bitset<Block> b(long_string);
    Tests::operator_shift_left(b, pos);
  }
  { // case pos >= n
    std::size_t pos = long_string.size();
    boost::dynamic_bitset<Block> b(long_string);
    Tests::operator_shift_left(b, pos);
  }
  //=====================================================================
  // Test b >> pos
  { // case pos == 0
    std::size_t pos = 0;
    boost::dynamic_bitset<Block> b(std::string("1010"));
    Tests::operator_shift_right(b, pos);
  }
  { // case pos == size()/2
    std::size_t pos = long_string.size() / 2;
    boost::dynamic_bitset<Block> b(long_string);
    Tests::operator_shift_right(b, pos);
  }
  { // case pos >= n
    std::size_t pos = long_string.size();
    boost::dynamic_bitset<Block> b(long_string);
    Tests::operator_shift_right(b, pos);
  }
  //=====================================================================
  // Test a & b
  {
    boost::dynamic_bitset<Block> lhs, rhs;
    Tests::operator_and(lhs, rhs);
  }
  {
    boost::dynamic_bitset<Block> lhs(std::string("1")), rhs(std::string("0"));
    Tests::operator_and(lhs, rhs);
  }
  {
    boost::dynamic_bitset<Block> lhs(long_string.size(), 0), rhs(long_string);
    Tests::operator_and(lhs, rhs);
  }
  {
    boost::dynamic_bitset<Block> lhs(long_string.size(), 1), rhs(long_string);
    Tests::operator_and(lhs, rhs);
  }
  //=====================================================================
  // Test a | b
  {
    boost::dynamic_bitset<Block> lhs, rhs;
    Tests::operator_or(lhs, rhs);
  }
  {
    boost::dynamic_bitset<Block> lhs(std::string("1")), rhs(std::string("0"));
    Tests::operator_or(lhs, rhs);
  }
  {
    boost::dynamic_bitset<Block> lhs(long_string.size(), 0), rhs(long_string);
    Tests::operator_or(lhs, rhs);
  }
  {
    boost::dynamic_bitset<Block> lhs(long_string.size(), 1), rhs(long_string);
    Tests::operator_or(lhs, rhs);
  }

//.........这里部分代码省略.........
			// the inverse of mmat is just the inverse of the diagonal components
			// here, we are extracting the inverse diagonal components only
			minv_vec[ee*np+jj] = 1./mmat[ee*np*np+jj*np+jj];
		}
		for(jj=0;jj<np;jj++){
			for(ii=0;ii<np;ii++){
				dv[ee*np*np+jj*np+ii] = minv_vec[ee*np+ii]*smat[jj*np+ii];
			}
		}
		for(jj=0;jj<2;jj++){
			for(ii=0;ii<np;ii++){
				df[ee*np*2+jj*np+ii]  = minv_vec[ee*np+ii]*mf[jj*np+ii];
			}
		}

	}
	
	// initialize qq field
	initialize(qq, xx, xmax, xmin, init_type);
	cfl = 1./(np*np);
	dt = cfl * min_dx / fabs(speed);
	rtime = 0.;
	nstep = 0;

	printf("Start Time Integration\n");

	// Runge-Kutta 4th order Time integration loop
	
	t_sta = clock();

	while(rtime < tend){
		dt = fmin(dt, tend-rtime);

		rhs(qq,	   k1, dv, df, ib, speed);

		for(ii=0;ii<ne*np;ii++)
			qtemp[ii] = qq[ii]+0.5*dt*k1[ii];
		rhs(qtemp, k2, dv, df, ib, speed);

		for(ii=0;ii<ne*np;ii++)
			qtemp[ii] = qq[ii]+0.5*dt*k2[ii];
		rhs(qtemp, k3, dv, df, ib, speed);
		
		for(ii=0;ii<ne*np;ii++)
			qtemp[ii] = qq[ii]+dt*k3[ii];
		rhs(qtemp, k4, dv, df, ib, speed);

		for(ii=0;ii<ne*np;ii++)
			qq[ii] += 1./6.*dt*(k1[ii]+2*k2[ii]+2*k3[ii]+k4[ii]);

		rtime += dt;
		nstep += 1;
		if(nstep%10000 == 0) 
			printf("nstep = %10ld, %5.1f%% complete\n", nstep, rtime/tend*100);
	}

	// timeloop ends here;

	printf("Integration complete\n");

	if(ne > 200){
		eres = 2;
	}
	else if (ne > 60){
		eres = 3;
	}

    void PolynomialFit4<Real>::DoLeastSquaresFit ( int numSamples,
            Real* trgSamples[4] )
    {
        // The matrix and vector for a linear system that determines the
        // coefficients of the fitted polynomial.
        GMatrix<Real> mat( mNumPowers, mNumPowers ); // initially zero
        GVector<Real> rhs( mNumPowers ); // initially zero
        mCoefficients = new1<Real>( mNumPowers );

        int row, col;
        for ( int i = 0; i < numSamples; ++i )
        {
            // Compute relevant powers of x and y.
            Real x = trgSamples[0][i];
            Real y = trgSamples[1][i];
            Real z = trgSamples[2][i];
            Real w = trgSamples[3][i];
            int j;
            for ( j = 1; j <= 2 * mMaxXPower; ++j )
            {
                mXPowers[j] = mXPowers[j - 1] * x;
            }
            for ( j = 1; j <= 2 * mMaxYPower; ++j )
            {
                mYPowers[j] = mYPowers[j - 1] * y;
            }
            for ( j = 1; j <= 2 * mMaxZPower; ++j )
            {
                mZPowers[j] = mZPowers[j - 1] * z;
            }

            for ( row = 0; row < mNumPowers; ++row )
            {
                // Update the upper-triangular portion of the symmetric matrix.
                Real xp, yp, zp;
                for ( col = row; col < mNumPowers; ++col )
                {
                    xp = mXPowers[mPowers[row][0] + mPowers[col][0]];
                    yp = mYPowers[mPowers[row][1] + mPowers[col][1]];
                    zp = mYPowers[mPowers[row][2] + mPowers[col][2]];
                    mat[row][col] += xp * yp * zp;
                }

                // Update the right-hand side of the system.
                xp = mXPowers[mPowers[row][0]];
                yp = mYPowers[mPowers[row][1]];
                zp = mYPowers[mPowers[row][2]];
                rhs[row] += xp * yp * zp * w;
            }
        }

        // Copy the upper-triangular portion of the symmetric matrix to the
        // lower-triangular portion.
        for ( row = 0; row < mNumPowers; ++row )
        {
            for ( col = 0; col < row; ++col )
            {
                mat[row][col] = mat[col][row];
            }
        }

        // Precondition by normalizing the sums.
        Real invNumSamples = ( ( Real )1 ) / ( Real )numSamples;
        for ( row = 0; row < mNumPowers; ++row )
        {
            for ( col = 0; col < mNumPowers; ++col )
            {
                mat[row][col] *= invNumSamples;
            }
            rhs[row] *= invNumSamples;
        }

        if ( LinearSystem<Real>().Solve( mat, rhs, mCoefficients ) )
        {
            mSolved = true;
        }
        else
        {
            memset( mCoefficients, 0, mNumPowers * sizeof( Real ) );
            mSolved = false;
        }
    }

    void NaturalSpline1<Real>::CreatePeriodicSpline ()
    {
        mB = new1<Real>( mNumSegments );
        mC = new1<Real>( mNumSegments );
        mD = new1<Real>( mNumSegments );

#if 1
        // Solving the system using a standard linear solver appears to be
        // numerically stable.
        const int size = 4 * mNumSegments;
        GMatrix<Real> mat( size, size );
        GVector<Real> rhs( size );
        int i, j, k;
        Real delta, delta2, delta3;
        for ( i = 0, j = 0; i < mNumSegments - 1; ++i, j += 4 )
        {
            delta = mTimes[i + 1] - mTimes[i];
            delta2 = delta * delta;
            delta3 = delta * delta2;

            mat[j + 0][j + 0] = ( Real )1;
            mat[j + 0][j + 1] = ( Real )0;
            mat[j + 0][j + 2] = ( Real )0;
            mat[j + 0][j + 3] = ( Real )0;
            mat[j + 1][j + 0] = ( Real )1;
            mat[j + 1][j + 1] = delta;
            mat[j + 1][j + 2] = delta2;
            mat[j + 1][j + 3] = delta3;
            mat[j + 2][j + 0] = ( Real )0;
            mat[j + 2][j + 1] = ( Real )1;
            mat[j + 2][j + 2] = ( ( Real )2 ) * delta;
            mat[j + 2][j + 3] = ( ( Real )3 ) * delta2;
            mat[j + 3][j + 0] = ( Real )0;
            mat[j + 3][j + 1] = ( Real )0;
            mat[j + 3][j + 2] = ( Real )1;
            mat[j + 3][j + 3] = ( ( Real )3 ) * delta;

            k = j + 4;
            mat[j + 0][k + 0] = ( Real )0;
            mat[j + 0][k + 1] = ( Real )0;
            mat[j + 0][k + 2] = ( Real )0;
            mat[j + 0][k + 3] = ( Real )0;
            mat[j + 1][k + 0] = ( Real ) - 1;
            mat[j + 1][k + 1] = ( Real )0;
            mat[j + 1][k + 2] = ( Real )0;
            mat[j + 1][k + 3] = ( Real )0;
            mat[j + 2][k + 0] = ( Real )0;
            mat[j + 2][k + 1] = ( Real ) - 1;
            mat[j + 2][k + 2] = ( Real )0;
            mat[j + 2][k + 3] = ( Real )0;
            mat[j + 3][k + 0] = ( Real )0;
            mat[j + 3][k + 1] = ( Real )0;
            mat[j + 3][k + 2] = ( Real ) - 1;
            mat[j + 3][k + 3] = ( Real )0;
        }

        delta = mTimes[i + 1] - mTimes[i];
        delta2 = delta * delta;
        delta3 = delta * delta2;

        mat[j + 0][j + 0] = ( Real )1;
        mat[j + 0][j + 1] = ( Real )0;
        mat[j + 0][j + 2] = ( Real )0;
        mat[j + 0][j + 3] = ( Real )0;
        mat[j + 1][j + 0] = ( Real )1;
        mat[j + 1][j + 1] = delta;
        mat[j + 1][j + 2] = delta2;
        mat[j + 1][j + 3] = delta3;
        mat[j + 2][j + 0] = ( Real )0;
        mat[j + 2][j + 1] = ( Real )1;
        mat[j + 2][j + 2] = ( ( Real )2 ) * delta;
        mat[j + 2][j + 3] = ( ( Real )3 ) * delta2;
        mat[j + 3][j + 0] = ( Real )0;
        mat[j + 3][j + 1] = ( Real )0;
        mat[j + 3][j + 2] = ( Real )1;
        mat[j + 3][j + 3] = ( ( Real )3 ) * delta;

        k = 0;
        mat[j + 0][k + 0] = ( Real )0;
        mat[j + 0][k + 1] = ( Real )0;
        mat[j + 0][k + 2] = ( Real )0;
        mat[j + 0][k + 3] = ( Real )0;
        mat[j + 1][k + 0] = ( Real ) - 1;
        mat[j + 1][k + 1] = ( Real )0;
        mat[j + 1][k + 2] = ( Real )0;
        mat[j + 1][k + 3] = ( Real )0;
        mat[j + 2][k + 0] = ( Real )0;
        mat[j + 2][k + 1] = ( Real ) - 1;
        mat[j + 2][k + 2] = ( Real )0;
        mat[j + 2][k + 3] = ( Real )0;
        mat[j + 3][k + 0] = ( Real )0;
        mat[j + 3][k + 1] = ( Real )0;
        mat[j + 3][k + 2] = ( Real ) - 1;
        mat[j + 3][k + 3] = ( Real )0;

        for ( i = 0, j = 0; i < mNumSegments; ++i, j += 4 )
        {
            rhs[j + 0] = mA[i];
            rhs[j + 1] = ( Real )0;
            rhs[j + 2] = ( Real )0;
//.........这里部分代码省略.........

vector<Grid> Planner::dStarLite(Map &map)
{
    // this is the implementation of the D star lite algorithm
    
    vector<Grid> wayPoint;

    
    // Initialization
    vector<vector<int> > g(map._rows,vector<int>(map._cols,map._rows*map._cols+1));
    vector<vector<int> > rhs(map._rows,vector<int>(map._cols,map._rows*map._cols+1));
    km =0;
    rhs[map.goal.y][map.goal.x] =0;
    Key n = calculateKey(map.goal,g,rhs,map.start,0);
    U.push_back({map.goal,n});
    
    
    // Computer current shortest path
    // Update the states of the map if the first element in the priority list can be updated or the goal is not reached
    while(U.front().key<calculateKey(map.goal,g,rhs,map.start,0)||(rhs[map.start.y][map.start.x]!=g[map.start.y][map.start.x]))
          {
              // take the key value and postion of the first element in the priority queue
              Key kold = U.front().key;
              Grid u = U.front().point;
              
              U.pop_front();
              Key knew = calculateKey(u,g,rhs,map.start,km); // calculate the new key value
              
              // update map if old value is different from the new value
              if (kold<knew)
              {
                  // if the new key is larger, the cost of the edge of the grid might be change
                  // the current grid should be updated and re-expanded
                  // insert it in the priority queue
                  auto pt =find_if(U.begin(),U.end(),[knew](Uelem &u){return u.key<knew;});
                  U.insert(pt,{u,knew});
              }
              else if (g[u.y][u.x]>rhs[u.y][u.x])
              {
                  // if the grid is overconstraint, there are new shorter paths detected
                  g[u.y][u.x] = rhs[u.y][u.x];
                  
                  // update all its neighbour value
                  vector<Grid> neightbour = findNeighbour(map, u, 8);
                  for(auto &n:neightbour)
                  {
                      updateVertex(U,n,g,rhs,map,km);
                  }
              }
              else
              {
                  // if the grid is underconstraint, the grid it self and its neightbour should all be updated
                  g[u.y][u.x] = map._cols*map._rows;
                  vector<Grid> neightbour = findNeighbour(map, u, 8);
                  for(auto &n:neightbour)
                  {
                      updateVertex(U,n,g,rhs,map,km);
                  }
                  updateVertex(U,u,g,rhs,map,km);
              }
              
                  
              
              
              
          }
    
    return wayPoint;
}

void check_params( struct user_parameters* params, int matrix_size, 
                   int block_size, double dx, double dy, double *f_,
                   int niter, double *u_, double *unew_) 
{
    double x, y;
    int i, j;
    double *udiff_ =(double*)malloc(matrix_size * matrix_size * sizeof(double));
    double (*udiff)[matrix_size][matrix_size] = (double (*)[matrix_size][matrix_size])udiff_;
    double (*unew)[matrix_size][matrix_size] = (double (*)[matrix_size][matrix_size])unew_;
    double (*u)[matrix_size][matrix_size] = (double (*)[matrix_size][matrix_size])u_;
    double (*f)[matrix_size][matrix_size] = (double (*)[matrix_size][matrix_size])f_;

    // Check for convergence.
    for (j = 0; j < matrix_size; j++) {
        y = (double) (j) / (double) (matrix_size - 1);
        for (i = 0; i < matrix_size; i++) {
            x = (double) (i) / (double) (matrix_size - 1);
            (*udiff)[i][j] = (*unew)[i][j] - u_exact(x, y);
            if( (*udiff)[i][j] > 1.0E-6 ) { 
                printf("error: %d, %d: %f\n", i, j, (*udiff)[i][j]);
            }
        }
    }
    double error = r8mat_rms(matrix_size, matrix_size, udiff_);

    double error1;
    // Set the right hand side array F.
    rhs(matrix_size, matrix_size, f_, block_size);

    for (j = 0; j < matrix_size; j++) {
        for (i = 0; i < matrix_size; i++) {
            if (i == 0 || i == matrix_size - 1 || j == 0 || j == matrix_size - 1) {
                (*unew)[i][j] = (*f)[i][j];
                (*u)[i][j] = (*f)[i][j];
            } else {
                (*unew)[i][j] = 0.0;
                (*u)[i][j] = 0.0;
            }
        }
    }

    sweep_seq(matrix_size, matrix_size, dx, dy, f_, 0, niter, u_, unew_);

    // Check for convergence.
    for (j = 0; j < matrix_size; j++) {
        y = (double) (j) / (double) (matrix_size - 1);
        for (i = 0; i < matrix_size; i++) {
            x = (double) (i) / (double) (matrix_size - 1);
            (*udiff)[i][j] = (*unew)[i][j] - u_exact(x, y);
            if( (*udiff)[i][j] > 1.0E-6 ) { 
                printf("error: %d, %d: %f\n", i, j, (*udiff)[i][j]);
            }
        }
    }
    error1 = r8mat_rms(matrix_size, matrix_size, udiff_);
    params->succeed = fabs(error - error1) < 1.0E-6;
    if(!params->succeed) {
        printf("error = %f, error1 = %f\n", error, error1);
    }
    free(udiff_);
}

double aug (cholmod_sparse *A)
{
    double r, maxerr = 0, bnorm, anorm ;
    cholmod_sparse *S, *Im, *In, *At, *A1, *A2, *Sup ;
    cholmod_dense *Alpha, *B, *Baug, *X, *W1, *W2, *R, *X2, X2mat ;
    cholmod_factor *L ;
    double *b, *baug, *rx, *w, *x ;
    Int nrow, ncol, nrhs, i, j, d, d2, save, save2, save3 ;

    if (A == NULL)
    {
	ERROR (CHOLMOD_INVALID, "cm: no A for aug") ;
	return (1) ;
    }

    if (A->xtype != CHOLMOD_REAL)
    {
	return (0) ;
    }

    /* ---------------------------------------------------------------------- */
    /* A is m-by-n, B must be m-by-nrhs */
    /* ---------------------------------------------------------------------- */

    nrow = A->nrow ;
    ncol = A->ncol ;
    B = rhs (A, 5, A->nrow + 7) ;

    /* ---------------------------------------------------------------------- */
    /* create scalars */
    /* ---------------------------------------------------------------------- */

    bnorm = CHOLMOD(norm_dense) (B, 0, cm) ;
    anorm = CHOLMOD(norm_sparse) (A, 1, cm) ;

    Alpha = CHOLMOD(eye) (1, 1, CHOLMOD_REAL, cm) ;
    if (Alpha != NULL)
    {
	((double *) (Alpha->x)) [0] = anorm ;
    }

    CHOLMOD(print_dense) (M1, "MinusOne", cm) ;
    CHOLMOD(print_dense) (Alpha, "Alpha = norm(A)", cm) ;

    /* ---------------------------------------------------------------------- */
    /* create augmented system, S = [-I A' ; A anorm*I] */
    /* ---------------------------------------------------------------------- */

    Im = CHOLMOD(speye) (nrow, nrow, CHOLMOD_REAL, cm) ;
    In = CHOLMOD(speye) (ncol, ncol, CHOLMOD_REAL, cm) ;
    CHOLMOD(scale) (Alpha, CHOLMOD_SCALAR, Im, cm) ;
    CHOLMOD(scale) (M1, CHOLMOD_SCALAR, In, cm) ;
    At = CHOLMOD(transpose) (A, 2, cm) ;

    /* use one of two equivalent methods */
    if (nrow % 2)
    {
	/* S = [[-In A'] ; [A alpha*Im]] */
	A1 = CHOLMOD(horzcat) (In, At, TRUE, cm) ;
	A2 = CHOLMOD(horzcat) (A,  Im, TRUE, cm) ;
	S = CHOLMOD(vertcat) (A1, A2, TRUE, cm) ;
    }
    else
    {
	/* S = [[-In ; A] [A' ; alpha*Im]] */
	A1 = CHOLMOD(vertcat) (In, A, TRUE, cm) ;
	A2 = CHOLMOD(vertcat) (At, Im, TRUE, cm) ;
	S = CHOLMOD(horzcat) (A1, A2, TRUE, cm) ;
    }

    CHOLMOD(free_sparse) (&Im, cm) ;
    CHOLMOD(free_sparse) (&In, cm) ;

    CHOLMOD(print_sparse) (S, "S, augmented system", cm) ;

    /* make a symmetric (upper) copy of S */
    Sup = CHOLMOD(copy) (S, 1, 1, cm) ;

    CHOLMOD(print_sparse) (S, "S, augmented system (upper)", cm) ;
    CHOLMOD(print_sparse) (Sup, "Sup", cm) ;

    /* ---------------------------------------------------------------------- */
    /* create augmented right-hand-side, Baug = [ zeros(ncol,nrhs) ; B ] */
    /* ---------------------------------------------------------------------- */

    b = NULL ;
    d = 0 ;
    nrhs = 0 ;
    d2 = 0 ;
    if (B != NULL)
    {
	nrhs = B->ncol ;
	d = B->d ;
	b = B->x ;
	Baug = CHOLMOD(zeros) (nrow+ncol, nrhs, CHOLMOD_REAL, cm) ;
	if (Baug != NULL)
	{
	    d2 = Baug->d ;
	    baug = Baug->x ;
	    for (j = 0 ; j < nrhs ; j++)
//.........这里部分代码省略.........

// This routine simply glues together many of the routines that are already
// written in the Poisson solver library
//
// phi( 1:SubNumPhysNodes       ) is a scalar quantity.  
//
// E1 ( 1:NumElems, 1:kmax2d ) is a vector quantity.
// E2 ( 1:NumElems, 1:kmax2d ) is a vector quantity.
//
// See also: ConvertEfieldOntoDGbasis
void ComputeElectricField( const double t, const mesh& Mesh, const dTensorBC5& q,
    dTensor2& E1, dTensor2& E2)
{

    //
    const int       mx = q.getsize(1);   assert_eq(mx,dogParamsCart2.get_mx());
    const int       my = q.getsize(2);   assert_eq(my,dogParamsCart2.get_my());
    const int NumElems = q.getsize(3);
    const int     meqn = q.getsize(4);
    const int     kmax = q.getsize(5);

    const int space_order = dogParams.get_space_order();

    // unstructured parameters:
    const int kmax2d    = E2.getsize(2);
    const int NumBndNodes  = Mesh.get_NumBndNodes();
    const int NumPhysNodes = Mesh.get_NumPhysNodes();

    // Quick error check
    if( !Mesh.get_is_submesh() )
    {
        printf("ERROR: mesh needs to have subfactor set to %d\n", space_order);
        printf("Go to Unstructured mesh and remesh the problem\n");
        exit(-1);
    }
    const int SubFactor    = Mesh.get_SubFactor();

    assert_eq( NumElems, Mesh.get_NumElems() );

    // -- Step 1: Compute rho -- //
    dTensor3 rho(NumElems, 1, kmax2d );
    void ComputeDensity( const mesh& Mesh, const dTensorBC5& q, dTensor3& rho );
    ComputeDensity( Mesh, q, rho );

    // -- Step 2: Figure out how large phi needs to be
    int SubNumPhysNodes = 0;
    int SubNumBndNodes  = 0;
    switch( dogParams.get_space_order() )
    {
        case 1:
            SubNumPhysNodes = NumPhysNodes;
            SubNumBndNodes  = NumBndNodes;
            break;

        case 2:
            SubNumPhysNodes = Mesh.get_SubNumPhysNodes();
            SubNumBndNodes  = Mesh.get_SubNumBndNodes();
            if(SubFactor!=2)
            {
                printf("\n");
                printf(" Error: for space_order = %i, need SubFactor = %i\n",space_order,2);
                printf("      SubFactor = %i\n",SubFactor);
                printf("\n");
                exit(1);
            }
            break;

        case 3:
            SubNumPhysNodes = Mesh.get_SubNumPhysNodes();
            SubNumBndNodes  = Mesh.get_SubNumBndNodes();
            if(SubFactor!=3)
            {
                printf("\n");
                printf(" Error: for space_order = %i, need SubFactor = %i\n",space_order,3);
                printf("      SubFactor = %i\n",SubFactor);
                printf("\n");
                exit(1);
            }
            break;

        default:
            printf("\n");
            printf(" ERROR in RunDogpack_unst.cpp: space_order value not supported.\n");
            printf("       space_order = %i\n",space_order);
            printf("\n");
            exit(1);
    }

    // local storage:
    dTensor1 rhs(SubNumPhysNodes);
    dTensor1 phi(SubNumPhysNodes);

    // Get Cholesky factorization matrix R
    //
    // TODO - this should be saved earlier in the code rather than reading
    // from file every time we with to run a Poisson solve!
    //
    SparseCholesky R(SubNumPhysNodes);
    string outputdir = dogParams.get_outputdir();
    R.init(outputdir);
    R.read(outputdir);
//.........这里部分代码省略.........

/* Cholesky update/downdate */
cs_long_t demo3 (problem *Prob)
{
    cs_cl *A, *C, *W = NULL, *WW, *WT, *E = NULL, *W2 ;
    cs_long_t n, k, *Li, *Lp, *Wi, *Wp, p1, p2, *p = NULL, ok ;
    cs_complex_t *b, *x, *resid, *y = NULL, *Lx, *Wx, s ;
    double t, t1 ;
    cs_cls *S = NULL ;
    cs_cln *N = NULL ;
    if (!Prob || !Prob->sym || Prob->A->n == 0) return (0) ;
    A = Prob->A ; C = Prob->C ; b = Prob->b ; x = Prob->x ; resid = Prob->resid;
    n = A->n ;
    if (!Prob->sym || n == 0) return (1) ;
    rhs (x, b, n) ;                             /* compute right-hand side */
    printf ("\nchol then update/downdate ") ;
    print_order (1) ;
    y = cs_cl_malloc (n, sizeof (cs_complex_t)) ;
    t = tic () ;
    S = cs_cl_schol (1, C) ;                       /* symbolic Chol, amd(A+A') */
    printf ("\nsymbolic chol time %8.2f\n", toc (t)) ;
    t = tic () ;
    N = cs_cl_chol (C, S) ;                        /* numeric Cholesky */
    printf ("numeric  chol time %8.2f\n", toc (t)) ;
    if (!S || !N || !y) return (done3 (0, S, N, y, W, E, p)) ;
    t = tic () ;
    cs_cl_ipvec (S->pinv, b, y, n) ;               /* y = P*b */
    cs_cl_lsolve (N->L, y) ;                       /* y = L\y */
    cs_cl_ltsolve (N->L, y) ;                      /* y = L'\y */
    cs_cl_pvec (S->pinv, y, x, n) ;                /* x = P'*y */
    printf ("solve    chol time %8.2f\n", toc (t)) ;
    printf ("original: ") ;
    print_resid (1, C, x, b, resid) ;           /* print residual */
    k = n/2 ;                                   /* construct W  */
    W = cs_cl_spalloc (n, 1, n, 1, 0) ;
    if (!W) return (done3 (0, S, N, y, W, E, p)) ;
    Lp = N->L->p ; Li = N->L->i ; Lx = N->L->x ;
    Wp = W->p ; Wi = W->i ; Wx = W->x ;
    Wp [0] = 0 ;
    p1 = Lp [k] ;
    Wp [1] = Lp [k+1] - p1 ;
    s = Lx [p1] ;
    srand (1) ;
    for ( ; p1 < Lp [k+1] ; p1++)
    {
        p2 = p1 - Lp [k] ;
        Wi [p2] = Li [p1] ;
        Wx [p2] = s * rand () / ((double) RAND_MAX) ;
    }
    t = tic () ;
    ok = cs_cl_updown (N->L, +1, W, S->parent) ;   /* update: L*L'+W*W' */
    t1 = toc (t) ;
    printf ("update:   time: %8.2f\n", t1) ;
    if (!ok) return (done3 (0, S, N, y, W, E, p)) ;
    t = tic () ;
    cs_cl_ipvec (S->pinv, b, y, n) ;               /* y = P*b */
    cs_cl_lsolve (N->L, y) ;                       /* y = L\y */
    cs_cl_ltsolve (N->L, y) ;                      /* y = L'\y */
    cs_cl_pvec (S->pinv, y, x, n) ;                /* x = P'*y */
    t = toc (t) ;
    p = cs_cl_pinv (S->pinv, n) ;
    W2 = cs_cl_permute (W, p, NULL, 1) ;           /* E = C + (P'W)*(P'W)' */
    WT = cs_cl_transpose (W2,1) ;
    WW = cs_cl_multiply (W2, WT) ;
    cs_cl_spfree (WT) ;
    cs_cl_spfree (W2) ;
    E = cs_cl_add (C, WW, 1, 1) ;
    cs_cl_spfree (WW) ;
    if (!E || !p) return (done3 (0, S, N, y, W, E, p)) ;
    printf ("update:   time: %8.2f (incl solve) ", t1+t) ;
    print_resid (1, E, x, b, resid) ;           /* print residual */
    cs_cl_nfree (N) ;                              /* clear N */
    t = tic () ;
    N = cs_cl_chol (E, S) ;                        /* numeric Cholesky */
    if (!N) return (done3 (0, S, N, y, W, E, p)) ;
    cs_cl_ipvec (S->pinv, b, y, n) ;               /* y = P*b */
    cs_cl_lsolve (N->L, y) ;                       /* y = L\y */
    cs_cl_ltsolve (N->L, y) ;                      /* y = L'\y */
    cs_cl_pvec (S->pinv, y, x, n) ;                /* x = P'*y */
    t = toc (t) ;
    printf ("rechol:   time: %8.2f (incl solve) ", t) ;
    print_resid (1, E, x, b, resid) ;           /* print residual */
    t = tic () ;
    ok = cs_cl_updown (N->L, -1, W, S->parent) ;   /* downdate: L*L'-W*W' */
    t1 = toc (t) ;
    if (!ok) return (done3 (0, S, N, y, W, E, p)) ;
    printf ("downdate: time: %8.2f\n", t1) ;
    t = tic () ;
    cs_cl_ipvec (S->pinv, b, y, n) ;               /* y = P*b */
    cs_cl_lsolve (N->L, y) ;                       /* y = L\y */
    cs_cl_ltsolve (N->L, y) ;                      /* y = L'\y */
    cs_cl_pvec (S->pinv, y, x, n) ;                /* x = P'*y */
    t = toc (t) ;
    printf ("downdate: time: %8.2f (incl solve) ", t1+t) ;
    print_resid (1, C, x, b, resid) ;           /* print residual */
    return (done3 (1, S, N, y, W, E, p)) ;
} 

TEST(uri_comparison_test, equality_test_capitalized_scheme_with_case_normalization) {
  network::uri lhs("http://www.example.com/");
  network::uri rhs("HTTP://www.example.com/");
  ASSERT_EQ(lhs.compare(rhs, network::uri_comparison_level::syntax_based), 0);
}

void DisparityProc::interpolate_natural_cubic_spline (
	// input
	std::vector<double> const & x,
	std::vector<double> const & y,
	int l_clamped,
	int r_clamped,
	std::vector<double> const & x_vis,
	// output
	std::vector<double> & y_vis )
{
  // Assemble the system; note that it is a tridiagonal one, so we store it as 3 vectors 

  int N = x.size();

  std::vector<double> lower(N-1, 0);
  std::vector<double> diag(N, 0);
  std::vector<double> upper(N-1, 0);
  std::vector<double> rhs(N);

  // Compute difference between data points
  std::vector<double> delta(N-1, 0);
  for (unsigned int i = 0; i < delta.size(); ++i)
    delta[i] = x[i+1] - x[i];

  // Initialize the system
  if (l_clamped) {
    upper[0] = (delta[0])/6;
    diag[0] = (delta[0])/3;
    rhs[0] = (y[1] - y[0])/(delta[0]);
  }
  else
    diag[0] = 1;

  for (int i = 1; i < N - 1; ++i) {
    upper[i] = (delta[i])/6;
    lower[i-1] = (delta[i-1])/6;
    diag[i] = (x[i+1] - x[i-1])/3;
    rhs[i] = ((y[i+1] - y[i])/delta[i]) - ((y[i] - y[i-1])/delta[i-1]);
  }

  if (r_clamped) {
    diag[N-1] = -(delta[N-2])/3;
    lower[N-2] = -(delta[N-2])/6;
    rhs[N-1] = (y[N-1] - y[N-2])/(delta[N-2]);
  }
  else
    diag[N-1] = 1;

  // Solve the system and store the result in d
  std::vector<double> d(N,0);
  solve_thomas(lower,diag,upper,rhs,d);

  // Evaluate the interploated curve at x_vis and store data in y_vis
  // We assume a sorted data list -> maybe implement later
  // The counter keeps track of which cubic polynomial the current x_vis values are
  int counter = 0;
  for (unsigned int i = 0; i < x_vis.size(); ++i) {
    while (x_vis[i] >= x[counter]) {
      counter++;
      if (counter == N) {
	counter--;
	break;
      }
    }
    y_vis[i] =	d[counter-1] * ((pow(x[counter] - x_vis[i], 3))/(6*delta[counter-1])) + 
		d[counter]*((pow(x_vis[i] - x[counter-1], 3))/(6*delta[counter-1])) + 
		((y[counter] - y[counter-1])/delta[counter-1] - (d[counter] -
		d[counter-1])*(delta[counter-1])/6)*(x_vis[i] - x[counter-1]) + (y[counter-1] - 
		d[counter-1]*(delta[counter-1]*delta[counter-1]/6));
  }
}

RegressionTreeNode* RegressionTree::buildTree(const RegressionData &trainingData,RegressionTreeNode *parent,Vector< UINT > features,UINT nodeID){
    
    const UINT M = trainingData.getNumSamples();
    const UINT N = trainingData.getNumInputDimensions();
    const UINT T = trainingData.getNumTargetDimensions();
    VectorFloat regressionData(T);
    
    //Update the nodeID
    
    //Get the depth
    UINT depth = 0;
    
    if( parent != NULL )
        depth = parent->getDepth() + 1;
    
    //If there are no training data then return NULL
    if( trainingData.getNumSamples() == 0 )
        return NULL;
    
    //Create the new node
    RegressionTreeNode *node = new RegressionTreeNode;
    
    if( node == NULL )
        return NULL;
    
    //Set the parent
    node->initNode( parent, depth, nodeID );
    
    //If there are no features left then create a leaf node and return
    if( features.size() == 0 || M < minNumSamplesPerNode || depth >= maxDepth ){
        
        //Flag that this is a leaf node
        node->setIsLeafNode( true );
        
        //Compute the regression data that will be stored at this node
        computeNodeRegressionData( trainingData, regressionData );
        
        //Set the node
        node->set( trainingData.getNumSamples(), 0, 0, regressionData );
        
        Regressifier::trainingLog << "Reached leaf node. Depth: " << depth << " NumSamples: " << trainingData.getNumSamples() << std::endl;
        
        return node;
    }
    
    //Compute the best spilt point
    UINT featureIndex = 0;
    Float threshold = 0;
    Float minError = 0;
    if( !computeBestSpilt( trainingData, features, featureIndex, threshold, minError ) ){
        delete node;
        return NULL;
    }
    
    Regressifier::trainingLog << "Depth: " << depth << " FeatureIndex: " << featureIndex << " Threshold: " << threshold << " MinError: " << minError << std::endl;
    
    //If the minError is below the minRMSError then create a leaf node and return
    if( minError <= minRMSErrorPerNode ){
        //Compute the regression data that will be stored at this node
        computeNodeRegressionData( trainingData, regressionData );
        
        //Set the node
        node->set( trainingData.getNumSamples(), featureIndex, threshold, regressionData );
        
        Regressifier::trainingLog << "Reached leaf node. Depth: " << depth << " NumSamples: " << M << std::endl;
        
        return node;
    }
    
    //Set the node
    node->set( trainingData.getNumSamples(), featureIndex, threshold, regressionData );
    
    //Remove the selected feature so we will not use it again
    if( removeFeaturesAtEachSpilt ){
        for(UINT i=0; i<features.getSize(); i++){
            if( features[i] == featureIndex ){
                features.erase( features.begin()+i );
                break;
            }
        }
    }
    
    //Split the data
    RegressionData lhs(N,T);
    RegressionData rhs(N,T);
    
    for(UINT i=0; i<M; i++){
        if( node->predict( trainingData[i].getInputVector() ) ){
            rhs.addSample(trainingData[i].getInputVector(), trainingData[i].getTargetVector());
        }else lhs.addSample(trainingData[i].getInputVector(), trainingData[i].getTargetVector());
    }
    
    //Run the recursive tree building on the children
    node->setLeftChild( buildTree( lhs, node, features, nodeID ) );
    node->setRightChild( buildTree( rhs, node, features, nodeID ) );
    
    return node;
}

double run(struct user_parameters* params)
{
    int matrix_size = params->matrix_size;
    if (matrix_size <= 0) {
        matrix_size = 512;
        params->matrix_size = matrix_size;
    }
    int block_size = params->blocksize;
    if (block_size <= 0) {
        block_size = 128;
        params->blocksize = block_size;
    }
    if ( (matrix_size % block_size) || (matrix_size % block_size) ) {
        params->succeed = 0;
        params->string2display = "*****ERROR: blocsize must divide NX and NY";
        return 0;
    }
    int niter = params->titer;
    if (niter <= 0) {
        niter = 4;
        params->titer = niter;
    }
    int ii,i,jj,j;
    double *f_ = (double*)malloc(matrix_size * matrix_size * sizeof(double));
    double (*f)[matrix_size][matrix_size] = (double (*)[matrix_size][matrix_size])f_;
    double *u_ = (double*)malloc(matrix_size * matrix_size * sizeof(double));
    double *unew_ = (double*)malloc(matrix_size * matrix_size * sizeof(double));
    double (*unew)[matrix_size][matrix_size] = (double (*)[matrix_size][matrix_size])unew_;
    double (*u)[matrix_size][matrix_size] = (double (*)[matrix_size][matrix_size])u_;

    double dx = 1.0 / (double) (matrix_size - 1);
    double dy = 1.0 / (double) (matrix_size - 1);

    rhs(matrix_size, matrix_size, f_, block_size);

    //Set the initial solution estimate UNEW.
    //We are "allowed" to pick up the boundary conditions exactly.
#pragma omp parallel
#pragma omp master
    for (j = 0; j < matrix_size; j+= block_size)
        for (i = 0; i < matrix_size; i+= block_size)
#pragma omp task firstprivate(i,j) private(ii,jj)
            for (jj=j; jj<j+block_size; ++jj)
                for (ii=i; ii<i+block_size; ++ii) {
                    if (ii == 0 || ii == matrix_size - 1 || jj == 0 || jj == matrix_size - 1) {
                        (*unew)[ii][jj] = (*f)[ii][jj];
                        (*u)[ii][jj] = (*f)[ii][jj];
                    } else {
                        (*unew)[ii][jj] = 0.0;
                        (*u)[ii][jj] = 0.0;
                    }
                }
    /// KERNEL INTENSIVE COMPUTATION
    START_TIMER;
#ifndef _OPENMP
    sweep_seq(matrix_size, matrix_size, dx, dy, f_, 0, niter, u_, unew_);
#else
    sweep(matrix_size, matrix_size, dx, dy, f_, 0, niter, u_, unew_, block_size);
#endif
    END_TIMER;

#ifdef _OPENMP
    if(params->check) {
        check_params(params, matrix_size, block_size, dx, dy, f_, niter, u_, unew_) ;
    }
#else
    params->succeed = 1;
#endif
    free(f_);
    free(u_);
    free(unew_);
    return TIMER;
}

//.........这里部分代码省略.........
    boost::dynamic_bitset<Block> b(std::string("1010"));
    Tests::operator_shift_left(b, pos);
  }
  { // case pos == size()/2
    std::size_t pos = long_string.size() / 2;
    boost::dynamic_bitset<Block> b(long_string);
    Tests::operator_shift_left(b, pos);
  }
  { // case pos >= n
    std::size_t pos = long_string.size();
    boost::dynamic_bitset<Block> b(long_string);
    Tests::operator_shift_left(b, pos);
  }
  //=====================================================================
  // Test b >> pos  
  { // case pos == 0
    std::size_t pos = 0;
    boost::dynamic_bitset<Block> b(std::string("1010"));
    Tests::operator_shift_right(b, pos);
  }
  { // case pos == size()/2
    std::size_t pos = long_string.size() / 2;
    boost::dynamic_bitset<Block> b(long_string);
    Tests::operator_shift_right(b, pos);
  }
  { // case pos >= n
    std::size_t pos = long_string.size();
    boost::dynamic_bitset<Block> b(long_string);
    Tests::operator_shift_right(b, pos);
  }
  //=====================================================================
  // Test a & b
  {
    boost::dynamic_bitset<Block> lhs, rhs;
    Tests::operator_and(lhs, rhs);
  }
  {
    boost::dynamic_bitset<Block> lhs(std::string("1")), rhs(std::string("0"));
    Tests::operator_and(lhs, rhs);
  }
  {
    boost::dynamic_bitset<Block> lhs(long_string.size(), 0), rhs(long_string);
    Tests::operator_and(lhs, rhs);
  }
  {
    boost::dynamic_bitset<Block> lhs(long_string.size(), 1), rhs(long_string);
    Tests::operator_and(lhs, rhs);
  }
  //=====================================================================
  // Test a | b
  {
    boost::dynamic_bitset<Block> lhs, rhs;
    Tests::operator_or(lhs, rhs);
  }
  {
    boost::dynamic_bitset<Block> lhs(std::string("1")), rhs(std::string("0"));
    Tests::operator_or(lhs, rhs);
  }
  {
    boost::dynamic_bitset<Block> lhs(long_string.size(), 0), rhs(long_string);
    Tests::operator_or(lhs, rhs);
  }
  {
    boost::dynamic_bitset<Block> lhs(long_string.size(), 1), rhs(long_string);
    Tests::operator_or(lhs, rhs);
  }

    void ContinuousArithmeticAsianVecerEngine::calculate() const {
        Real expectedAverage;

        QL_REQUIRE(arguments_.averageType == Average::Arithmetic,
                   "not an Arithmetic average option");
        QL_REQUIRE(arguments_.exercise->type() == Exercise::European,
                   "not an European Option");

        DayCounter rfdc  = process_->riskFreeRate()->dayCounter();
        DayCounter divdc = process_->dividendYield()->dayCounter();
        DayCounter voldc = process_->blackVolatility()->dayCounter();
        Real S_0 = process_->stateVariable()->value();

        // payoff
        ext::shared_ptr<StrikedTypePayoff> payoff =
            ext::dynamic_pointer_cast<StrikedTypePayoff>(arguments_.payoff);
        QL_REQUIRE(payoff, "non-plain payoff given");

        // original time to maturity
        Date maturity = arguments_.exercise->lastDate();

        Real X = payoff->strike();
        QL_REQUIRE(z_min_<=0 && z_max_>=0,
                   "strike (0 for vecer fixed strike asian)  not on Grid");

        Volatility sigma =
            process_->blackVolatility()->blackVol(maturity, X);

        Rate r = process_->riskFreeRate()->
            zeroRate(maturity, rfdc, Continuous, NoFrequency);
        Rate q = process_->dividendYield()->
            zeroRate(maturity, divdc, Continuous, NoFrequency);

        Date today(Settings::instance().evaluationDate());

        QL_REQUIRE(startDate_>=today,
                   "Seasoned Asian not yet implemented");

        // Expiry in Years
        Time T = rfdc.yearFraction(today,
                                   arguments_.exercise->lastDate());
        Time T1 = rfdc.yearFraction(today,
                                    startDate_ );            // Average Begin
        Time T2 = T;            // Average End (In this version only Maturity...)

        if ((T2 - T1) < 0.001) {
            // its a vanilla option. Use vanilla engine
            VanillaOption europeanOption(payoff, arguments_.exercise);
            europeanOption.setPricingEngine(
                        ext::make_shared<AnalyticEuropeanEngine>(process_));
            results_.value = europeanOption.NPV();

        } else {
            Real Theta = 0.5;        // Mixed Scheme: 0.5 = Crank Nicolson
            Real Z_0 = cont_strategy(0,T1,T2,q,r) - std::exp(-r*T) * X /S_0;

            QL_REQUIRE(Z_0>=z_min_ && Z_0<=z_max_,
                       "spot not on grid");

            Real h = (z_max_ - z_min_) / assetSteps_; // Space step size
            Real k = T / timeSteps_;         // Time Step size

            Real sigma2 = sigma * sigma, vecerTerm;

            Array SVec(assetSteps_+1),u_initial(assetSteps_+1),
                  u(assetSteps_+1),rhs(assetSteps_+1);

            for (Natural i= 0; i<= SVec.size()-1;i++) {
                SVec[i] = z_min_ + i * h;     // Value of Underlying on the grid
            }

            // Begin gamma construction
            TridiagonalOperator gammaOp = DPlusDMinus(assetSteps_+1,h);

            Array upperD = gammaOp.upperDiagonal();
            Array lowerD = gammaOp.lowerDiagonal();
            Array Dia    = gammaOp.diagonal();

            // Construct Vecer operator
            TridiagonalOperator explicit_part(gammaOp.size());
            TridiagonalOperator implicit_part(gammaOp.size());

            for (Natural i= 0; i<= SVec.size()-1;i++) {
                u_initial[i] = std::max<Real>(SVec[i] , 0.0); // Call Payoff
            }

            u = u_initial;

            // Start Time Loop

            for (Natural j = 1; j<=timeSteps_;j++) {
                if (Theta != 1.0) { // Explicit Part
                    for (Natural i = 1; i<= SVec.size()-2;i++) {
                        vecerTerm = SVec[i] - std::exp(-q * (T-(j-1)*k))
                                  * cont_strategy(T-(j-1)*k,T1,T2,q,r);
                        gammaOp.setMidRow(i,
                            0.5 * sigma2 * vecerTerm * vecerTerm  * lowerD[i-1],
                            0.5 * sigma2 * vecerTerm * vecerTerm  * Dia[i],
                            0.5 * sigma2 *  vecerTerm * vecerTerm * upperD[i]);
                    }
//.........这里部分代码省略.........