C++ Array::size方法代码示例

void IceModel_Decode::run_timestep(double time_s,
	blitz::Array<int,1> const &indices,
	std::map<IceField, blitz::Array<double,1>> const &vals2)
{
printf("BEGIN IceModel_Decode::run_timestep(%f) size=%ld\n", time_s, indices.size());
	std::map<IceField, blitz::Array<double,1>> vals2d;	/// Decoded fields

	// Loop through the fields we require
	std::set<IceField> fields;
	get_required_fields(fields);
	for (auto field = fields.begin(); field != fields.end(); ++field) {
printf("Looking for required field %s\n", field->str());

		// Look up the field we require
		auto ii = vals2.find(*field);
		if (ii == vals2.end()) {
			fprintf(stderr, "Cannot find required ice field = %s\n", field->str());
			throw std::exception();
		}
		blitz::Array<double,1> vals(ii->second);

		// Decode the field!
		blitz::Array<double,1> valsd(ndata());
		valsd = nan;
		int n = indices.size();
		for (int i=0; i < n; ++i) {
			int ix = indices(i);
			// Do our own bounds checking!
			if (ix < 0 || ix >= ndata()) {
				fprintf(stderr, "IceModel: index %d out of range [0, %d)\n", ix, ndata());
				throw std::exception();
			}

#if 0
			// Sanity check for NaN coming through
			if (std::isnan(vals(i))) {
				fprintf(stderr, "IceModel::decode: vals[%d] (index=%d) is NaN!\n", i, ix);
				throw std::exception();
			}
#endif

			// Add this value to existing field
			double &oval = valsd(ix);
			if (std::isnan(oval)) oval = vals(i);
			else oval += vals(i);
		}

		// Store decoded field in our output
		vals2d.insert(std::make_pair(*field, valsd));
printf("Done decoding required field, %s\n", field->str());
	}

	// Pass decoded fields on to subclass
	run_decoded(time_s, vals2d);
printf("END IceModel_Decode::run_timestep(%ld)\n", time_s);
}

double VectorInnerProduct(const blitz::Array<double, Rank> &u, const blitz::Array<double, Rank> &v)
{
	if (u.size() != v.size())
	{
		cout << "Vector u and v is of different size: " << u.size() << " != " << v.size() << endl;
		throw std::runtime_error("invalid vector sizes for inner product");
	}
	if (!u.isStorageContiguous())
	{
		throw std::runtime_error("Vector u is not contiguous");
	}
	if (!v.isStorageContiguous())
	{
		throw std::runtime_error("Vector v is not contiguous");
	}
	int N = u.size();
	int uStride = 1; //u.stride(0);
	int vStride = 1; //v.stride(0);
	return BLAS_NAME(ddot)(N, (double*)u.data(), uStride, (double*)v.data(), vStride);
}

/*
 * Sets the radial Coulomb wave F_l(k*r, eta), with eta = Z/k into data for all radial
 * grid points specified by r
 */
void SetRadialCoulombWave(int Z, int l, double k, blitz::Array<double, 1> r, blitz::Array<double, 1> data)
{
	double eta = Z / k;

	for (int i=0; i<r.size(); i++)
	{	
		double x = k * r(i);
		gsl_sf_result F, Fp, G, Gp;
		double exp_F, exp_G;
		int error = gsl_sf_coulomb_wave_FG_e(eta, x, (double)l, 0., &F, &Fp, &G, &Gp, &exp_F, &exp_G);
		if (error == GSL_EOVRFLW)
		{
			cout << "WARNING: Overflow in SetCoulombWave(" << Z << ", " << l << ", " << k << ", r=" << r(i) << ");" << endl;
			cout << "         exp_F = " << exp_F << ", exp_G = " << exp_G << endl;
		}

		data(i) = F.val;
	}
}

void CopyTensorPotentialToEpetraMatrix(Epetra_FECrsMatrix_Ptr epetraMatrix, blitz::Array<cplx, Rank> potentialData, list pyLocalBasisPairs, blitz::TinyVector<int, Rank> globalStrides, double cutoff)
{
	blitz::Array<int, 2> indexArray;
	indexArray.resize(4 * potentialData.size(), 2);
	indexArray = -100;

	double sqrCutoff = sqr(cutoff);

	//Setup structures for calculating matrix row/col indices from the 
	//basis pairs in the tensor potential
	blitz::TinyVector< blitz::Array<int, 2>, Rank > localBasisPairs;
	for (int rank=0; rank<Rank; rank++)
	{
		localBasisPairs(rank).reference( boost::python::extract< blitz::Array<int, 2> >(pyLocalBasisPairs[rank]) );
	}

	//Iterate over all items in potentialData
	typename blitz::Array<cplx, Rank>::iterator it = potentialData.begin();
	for (int linearCount=0; linearCount<potentialData.size(); linearCount++)
	{
		int globalRow = 0;
		int globalCol = 0;
		for (int rank=0; rank<Rank; rank++)
		{
			int rankPos = it.position()(rank);
			globalRow += globalStrides(rank) * localBasisPairs(rank)(rankPos, 0);
			globalCol += globalStrides(rank) * localBasisPairs(rank)(rankPos, 1);
		}

		double realVal = real(*it);
		double imagVal = imag(*it);

		//Skip padded elements (they have negative row/col index)
		if ((globalRow < 0) || (globalCol < 0))
		{
			it++;
			continue;
		}
		
		/*
		 * Because epetra does not support complex natively,
		 * each matrix element is a 2x2 block
		 *
		 * (A_r  -A_i )  (c_r)  =  (A_r + i A_i) * (c_r + i c_i)  =  A * c
		 * (A_i   A_r )  (c_i)
		 *
		 * Detect if A_i or A_r is zero to avoid redundant elements
		 */

		//Insert values into matrix
		if (sqr(realVal) > sqrCutoff)
		{
			int r,c;
			r = 2*globalRow; c = 2*globalCol;
			indexArray(4*linearCount, 0) = r;
			indexArray(4*linearCount, 1) = c;
			epetraMatrix->InsertGlobalValues(r, 1, &realVal, &c);
			r++; c++;
			epetraMatrix->InsertGlobalValues(r, 1, &realVal, &c);
			indexArray(4*linearCount+1, 0) = r;
			indexArray(4*linearCount+1, 1) = c;
		}
		if (sqr(imagVal) > sqrCutoff)
		{
			int r,c;
			//Upper row, A_i with minus sign
			r = 2*globalRow; c = 2*globalCol+1;
			indexArray(4*linearCount+2, 0) = r;
			indexArray(4*linearCount+2, 1) = c;
			imagVal = -imagVal;
			epetraMatrix->InsertGlobalValues(r, 1, &imagVal, &c);

			//Lower row, A_i without minus sign
			imagVal = -imagVal;
			epetraMatrix->InsertGlobalValues(c, 1, &imagVal, &r);
			indexArray(4*linearCount+3, 0) = c;
			indexArray(4*linearCount+3, 1) = r;
		}
		++it;	
	}
}