C++ MDArray类代码示例

/// sequentially read nxn array in row major order
void sequentialRead()
{
   MDArray *array = new MDArray(fileName);
   PagedStorageContainer::resetPerfCounts();
   gettimeofday(&tim, NULL);
   t1 = tim.tv_sec + tim.tv_usec/1000000.0;
   MDIterator *it = array->createIterator(Dense, row);
   while (it->moveNext())
   {
      it->get(coord, datum);
   }
   delete it;
   delete array;

   gettimeofday(&tim, NULL);
   t2 = tim.tv_sec + tim.tv_usec/1000000.0;
   readCount += PagedStorageContainer::readCount;
   writeCount += PagedStorageContainer::writeCount;
   accessTime += PagedStorageContainer::accessTime;

   execTime += (t2 - t1);

   fprintf(pFile, "%-14s %8i %8i %11.3f %11.3f\n", "seq read", PagedStorageContainer::readCount, PagedStorageContainer::writeCount, PagedStorageContainer::accessTime, t2-t1);

}

    void StringFunction::operator()(MDArray& input_phy_points, MDArray& output_values,
                                    const stk::mesh::Bucket& bucket, const MDArray& parametric_coordinates, double time_value_optional) 
    {
      PRINT("tmp srk StringFunction::operator(bucket) getName()= " << getName() << " input_phy_points= " << input_phy_points << " output_values= " << output_values);

      VERIFY_OP(input_phy_points.rank(), ==, 3, "StringFunction::operator() must pass in input_phy_points(numCells, numPointsPerCell, spaceDim)");
      int nPoints = input_phy_points.dimension(1);
      int nSpaceDim = input_phy_points.dimension(2);
      int nOutDim = output_values.dimension(2);
      MDArray input_phy_points_one(1, nPoints, nSpaceDim);
      MDArray output_values_one   (1, nPoints, nOutDim);
      const unsigned num_elements_in_bucket = bucket.size();
      for (unsigned iElement = 0; iElement < num_elements_in_bucket; iElement++)
        {
          stk::mesh::Entity& element = bucket[iElement];
          for (int iPoint = 0; iPoint<nPoints; iPoint++)
            {
              for (int iSpaceDim=0; iSpaceDim < nSpaceDim; iSpaceDim++)
                {
                  input_phy_points_one(0, iPoint, iSpaceDim) = input_phy_points(iElement, iPoint, iSpaceDim);
                }
            }
          (*this)(input_phy_points_one, output_values_one, element, parametric_coordinates, time_value_optional);
          for (int iPoint = 0; iPoint<nPoints; iPoint++)
            {
              for (int iDOF = 0; iDOF < nOutDim; iDOF++)
                {
                  output_values(iElement, iPoint, iDOF) = output_values_one(0, iPoint, iDOF);
                }
            }
        }
    }

/// strided insertion, using the iterator
void stridedIteratorInsert()
{
    MDArray *array;
    Btree::MSplitter<Datum_t> leafSp;
    Btree::MSplitter<PID_t> intSp;
    if (type == DMA)
        array = new MDArray(dim, type, row, fileName);
    else if (type == BTREE)
        array = new MDArray(dim, row, &leafSp, &intSp, fileName);
   PagedStorageContainer::resetPerfCounts();
   
   gettimeofday(&tim, NULL);
   t1 = tim.tv_sec + tim.tv_usec/1000000.0;
   MDIterator *it = array->createIterator(Dense, col);
   while (it->moveNext())
   {
      it->put(12345);
   }
   delete it;
   delete array;

   gettimeofday(&tim, NULL);
   t2 = tim.tv_sec + tim.tv_usec/1000000.0;
   readCount += PagedStorageContainer::readCount;
   writeCount += PagedStorageContainer::writeCount;
   accessTime += PagedStorageContainer::accessTime;
   
   execTime += (t2 - t1);

   fprintf(pFile, "%-14s %8i %8i %11.3f %11.3f\n", "batch stride", PagedStorageContainer::readCount, PagedStorageContainer::writeCount, PagedStorageContainer::accessTime, t2-t1);

}

/// strided insertion into nxn array (insert going down columns first into row
//major array
void stridedInsert()
{
    MDArray *array;
    Btree::MSplitter<Datum_t> leafSp;
    Btree::MSplitter<PID_t> intSp;
    if (type == DMA)
        array = new MDArray(dim, type, row, fileName);
    else if (type == BTREE)
        array = new MDArray(dim, row, &leafSp, &intSp, fileName);
   PagedStorageContainer::resetPerfCounts();

   gettimeofday(&tim, NULL);
   t1 = tim.tv_sec + tim.tv_usec/1000000.0;
   for (i64 j = 0; j < n; j++)
   {
      for (i64 i = 0; i < n; i++)
      {
         coord = MDCoord(2, i, j);
         array->put(coord, 12345);
      }
   }
   delete array;

   gettimeofday(&tim, NULL);
   t2 = tim.tv_sec + tim.tv_usec/1000000.0;
   readCount += PagedStorageContainer::readCount;
   writeCount += PagedStorageContainer::writeCount;
   accessTime += PagedStorageContainer::accessTime;
   execTime += (t2 - t1);

   fprintf(pFile, "%-14s %8i %8i %11.3f %11.3f\n", "stride insert", PagedStorageContainer::readCount, PagedStorageContainer::writeCount, PagedStorageContainer::accessTime, t2-t1);
   
}

/// random read from nxn array. all indices are read once
void randomRead()
{
   i64 k[n*n];
   for (int i = 0; i < n*n; i++)
   {
      k[i] = i;
   }
   MDArray *array = new MDArray(fileName);
   PagedStorageContainer::resetPerfCounts();
   permute(k, n*n);

   gettimeofday(&tim, NULL);
   t1 = tim.tv_sec + tim.tv_usec/1000000.0;
   for (int i = 0; i < n*n; i++)
   {
      coord = MDCoord(2, k[i]/n, k[i]%n);
      array->get(coord, datum);
   }
   delete array;

   gettimeofday(&tim, NULL);
   t2 = tim.tv_sec + tim.tv_usec/1000000.0;
   readCount += PagedStorageContainer::readCount;
   writeCount += PagedStorageContainer::writeCount;
   accessTime += PagedStorageContainer::accessTime;
   execTime += (t2 - t1);

   fprintf(pFile, "%-14s %8i %8i %11.3f %11.3f\n", "rand read", PagedStorageContainer::readCount, PagedStorageContainer::writeCount, PagedStorageContainer::accessTime, t2-t1);

}

void MultinomialDensities::estimate(vector<MDArray<double> > & ess)
{
	assert(!ess.empty());

    MDArray<double> e = ess[0];
    e.add_inplace(pseudo_count);
	set_parameters(e);
}

TEST(MDArray, Read)
{
    i64 rows = 20L;
    i64 cols = 20L;
    MDCoord dim(2, rows, cols);
    StorageType type = DMA;
    i64 blockDims[] = {10L, 10L};
    u8 blockOrders[] = {0, 1};
    u8 microOrders[] = {0, 1};
    //permute(blockOrders, dim.nDim);
    //permute(microOrders, dim.nDim);
    Linearization *block = new BlockBased(dim.nDim, dim.coords, blockDims,
                                          blockOrders, microOrders);
    const char *fileName = "test1.bin";
    MDCoord key;
    Datum_t datum = 1.0;

    MDArray *array = new MDArray(dim, type, block, fileName);
    
    MDIterator *it = array->createIterator(Dense, block);

    srand(time(NULL));
    const int num=100;
    std::map<int, MDCoord> keys;
    std::map<int, Datum_t> data;
    int total = rows*cols;
    for(int i=0; i<num; i++) {
        int k = rand() % total;
        if (data.find(k) == data.end()) {
            data[k] = rand();
            it->setIndexRange(k, k+1);
            it->moveNext();
            it->put(data[k]);
            MDCoord key;
            it->get(key, datum);
            keys[k] = key;
        }
        else
            i--;
    }
    
    delete it;

    delete array;

    array = new MDArray(fileName);
    Datum_t datum1;
    for (std::map<int,MDCoord>::iterator s = keys.begin();
         s != keys.end();
         ++s) {
        //cout<<"checking "<<s->second.toString()<<endl;
        array->get(s->second, datum1);
        ASSERT_DOUBLE_EQ(data[s->first], datum1);
    }
    delete array;
    
    delete block;
}

int main(int argc, char **argv)
{
	const int required = 2;
	char fileName[100] = "/riot/mb";
	unsigned int tm;
	if (argc >= required+1)
		tm = atoi(argv[required]);
	else
		tm = time(NULL);
	srand(tm);
	cerr<<"seed = "<<tm<<endl;
	if (argc < required) {
		cerr<<"Usage: "<<argv[0]<<"<splitter type>"<<endl
			<<"splitter type: M,A,D"<<endl;
		return 0;
	}

	char splitterType = argv[1][0];

    i64 arraydims[] = {20000,20000};
    i64 blockdims[] = {31,31};
    u8 orders[] = {0, 1};
    Linearization<2> *lin = new BlockBased<2>(arraydims, blockdims, orders,
                                              orders);
    MDArray<2> *array;
    StorageParam sp;
    sp.fileName = fileName;
    sp.intSp = 'M';
    sp.leafSp = splitterType;
    if (splitterType == 'D') {
        // dma
        sp.type = DMA;
    }
    else if (splitterType == 'M') {
        sp.type = BTREE;
        sp.useDenseLeaf = false;
    }
    else if (splitterType == 'A') {
        sp.type = BTREE;
        sp.useDenseLeaf = true;
    }
    else {
        cerr<<"wrong splitter type"<<endl;
        exit(1);
    }
    
    array = new MDArray<2>(&sp, MDCoord<2>(arraydims), lin);

    for (i64 i=0; i<20000; ++i) {
        for (i64 j=0; j<20000; ++j) {
            MDCoord<2> c(i,j);
            array->put(c, 1.0);
        }
    }
	delete lin;
    delete array;
}

void ParentMap::set(uint h, MDArray<double> & v) {
	vector<uint> & v1 = indices[h];
	uint i = v1.back();
	assert(v.get_shape().size() == 1);
	uint dim = v.get_shape().front();
	uint k = 0;
	for (uint j = i; j < i + dim; j++, k++) {
		*(seq[j]) = v[k];
	}
}

MDArray<double> Shrinkage::combine(MDArray<double> target,MDArray<double> MLE,double lambda)
{
	assert(target.get_shape() == MLE.get_shape());

	//        assert(lambda>=0 && lambda<=1);
	lambda = std::max(lambda, 0.0);
	lambda = std::min(lambda, 1.0);

	return add_matrices(mult_num_to_matrix(target,lambda),mult_num_to_matrix(MLE,1-lambda));
}

Shrinkage::Shrinkage(MDArray<double> data, MDArray<double> MLE)
{	
        N=data.get_shape()[0];
        P=data.get_shape()[1];
	assert(MLE.get_shape()[0]==P && MLE.get_shape()[1]==P);
        this->data=data;
        wk=wk_StDM=vec(N,P,P);
        S_MLE=vec(P,P);
        StDM=vec(N,P);
	set_S_MLE(MLE);
}

    void FieldFunction::operator()(MDArray& input_phy_points, MDArray& output_field_values,
                                   const stk_classic::mesh::Bucket& bucket, const MDArray& parametric_coordinates, double time_value_optional)
    {
      EXCEPTWATCH;
#ifndef NDEBUG
      int num_elements_in_bucket   = bucket.size();
      VERIFY_OP(input_phy_points.dimension(0), ==, num_elements_in_bucket, "FieldFunction::operator() mismatch in input_phy_points and num_elements_in_bucket");
      VERIFY_OP(output_field_values.dimension(0), ==, num_elements_in_bucket, "FieldFunction::operator() mismatch in input_phy_points and num_elements_in_bucket");
#endif
      helper(input_phy_points, output_field_values, bucket, parametric_coordinates, time_value_optional);
    }

void printMDArray(const MDArray<nDim> &array)
{
	MDCoord<nDim> dim = array.getDims();
	Datum_t datum;
    for (int i=0; i<dim[0]; i++) {
        for (int j=0; j<dim[1]; j++) {
            MDCoord<2> c(i,j);
            array.get(c, datum);
			cout<<datum<<"\t";
           }
        cout<<endl;
    }
}

void IntrepidSideCell<MDArray>::mapToCellPhysicalFrame( 
    const MDArray& parametric_coords, MDArray& physical_coords )
{
    DTK_REQUIRE( 2 == parametric_coords.rank() );
    DTK_REQUIRE( 3 == physical_coords.rank() );
    DTK_REQUIRE( parametric_coords.dimension(1) ==
		   Teuchos::as<int>(this->d_topology.getDimension()) );
    DTK_REQUIRE( physical_coords.dimension(0) ==
		   this->d_cell_node_coords.dimension(0) );
    DTK_REQUIRE( physical_coords.dimension(1) ==
		   parametric_coords.dimension(0) );
    DTK_REQUIRE( physical_coords.dimension(2) ==
		   Teuchos::as<int>(d_parent_topology.getDimension()) );

    MDArray mapped_coords( parametric_coords.dimension(0),
			   d_parent_topology.getDimension() );

    Intrepid::CellTools<Scalar>::mapToReferenceSubcell(
	mapped_coords, parametric_coords, this->d_topology.getDimension(),
	d_side_id, d_parent_topology );

    Intrepid::CellTools<Scalar>::mapToPhysicalFrame( 
	physical_coords, mapped_coords, 
	this->d_cell_node_coords, d_parent_topology );
}

void checkRect(const MDArray<2> &array, const MDCoord<2> &begin, const MDCoord<2> &end, const Datum_t *data)
{
	int k = 0;
	Datum_t datum;
	for (i64 j=begin[1]; j<=end[1]; ++j) 
		for (i64 i=begin[0]; i<=end[0]; ++i) {
			array.get(MDCoord<2>(i,j), datum);
			ASSERT_EQ(data[k++], datum);
		}
}