C++ Wavefunction::GetData方法代码示例

void DistributedOverlapMatrix<Rank>::WavefunctionToMultiVector(Wavefunction<Rank> &psi, Epetra_MultiVector_Ptr vec, int opRank)
{
	//Map wavefunction to 3D array (compress before- and after-ranks)
	blitz::Array<cplx, 3> psiData = MapToRank3(psi.GetData(), opRank, 1);
	int beforeSize = psiData.extent(0);
	int opSize = psiData.extent(1);
	int afterSize = psiData.extent(2);
	int otherSize = beforeSize*afterSize;

	//Copy real and imag part of wavefunction into multivector
	double realVal, imagVal;
	for (int i=0; i<beforeSize; i++)
	{
		for (int j=0; j<opSize; j++)
		{
			for (int k=0; k<afterSize; k++)
			{
				//cout << "cur vec = " << i * afterSize + k << ", cur idx = " << j << endl;
				//cout << "cur vec = " << i * afterSize + k + otherSize << ", cur idx = " << j << endl;
				
				realVal = real( psiData(i,j,k) );
				imagVal = imag( psiData(i,j,k) );
				vec->ReplaceMyValue(j, i*afterSize + k, realVal);
				vec->ReplaceMyValue(j, i*afterSize + k + otherSize, imagVal);
			}
		}
	}
}

void DistributedOverlapMatrix<Rank>::MultiVectorToWavefunction(Wavefunction<Rank> &psi, Epetra_MultiVector_Ptr vec, int opRank)
{
	//Put result into psi
	blitz::Array<cplx, 3> data = MapToRank3(psi.GetData(), opRank, 1);
	int beforeSize = data.extent(0);
	int opSize = data.extent(1);
	int afterSize = data.extent(2);
	int otherSize = beforeSize*afterSize;
	double *destVecView=0;
	vec->ExtractView(&destVecView, &opSize);
	blitz::Array<double, 2> data2(destVecView, shape(2 * otherSize, opSize), shape(opSize, 1), neverDeleteData);
	//blitz::Array<double, 2> data2;
	//data2.resize(2 * otherSize, opSize);
	//vec->ExtractCopy(data2.data(), opSize);
	for (int i=0; i<beforeSize; i++)
	{
		for (int j=0; j<opSize; j++)
		{
			for (int k=0; k<afterSize; k++)
			{
				double realVal = data2(i*afterSize + k, j);
				double imagVal = data2(i*afterSize + k + otherSize, j);
				data(i,j,k) = realVal + I * imagVal;
			}
		}
	}
}

/*
 * Maps the given wavefunction to one where the particles are exchanged
 * psi(1,2) -> psi(2,1)
 */
Wavefunction<3>::Ptr GetWavefunctionParticleExchange(Wavefunction<3>::Ptr psi, list angularSymmetrizationPairs)
{
	typedef Array<cplx, 3> ArrayType;
	ArrayType data = psi->GetData();

	int countr = data.extent(1);
	typedef stl_input_iterator<tuple> Iterator;
	Iterator begin(angularSymmetrizationPairs);
	Iterator end;

	Wavefunction<3>::Ptr exchgPsi = psi->Copy();
	ArrayType exchgData = exchgPsi->GetData();

	for (Iterator i=begin; i!=end; i++)
	{
		int a1 = extract<int>((*i)[0]);
		int a2 = extract<int>((*i)[1]);

		for (int r1=0; r1<countr; r1++)
		{
			for (int r2=0; r2<countr; r2++)
			{
				exchgData(a1, r1, r2) = data(a2, r2, r1);
			}
		}
	}
	
	return exchgPsi;
}

void CombinedRepresentation<Rank>::MultiplyOverlap(Wavefunction<Rank> &srcPsi, Wavefunction<Rank> &dstPsi, int rank)
{
	using namespace blitz;

	if (this->IsOrthogonalBasis(rank))
	{
		dstPsi.GetData() = srcPsi.GetData();
	}
	else
	{
		DistributedOverlap->MultiplyOverlapRank(srcPsi, dstPsi, rank, true);
		//Map the data to a 3D array, where the b-spline part is the middle rank
		//Array<cplx, 3> srcData = MapToRank3(srcPsi.Data, rank, 1);
		//Array<cplx, 3> dstData = MapToRank3(dstPsi.Data, rank, 1);

		//this->GetGlobalOverlapMatrix(rank)->MultiplyOverlapTensor(srcData, dstData);
	}

}

void DistributedOverlapMatrix<Rank>::SetupRank(Wavefunction<Rank> &srcPsi, int opRank)
{
	//Need a copy of srcPsi for future reference (only on first call to this function)
	if (!HasPsi)
	{
		Psi = srcPsi.Copy();
		HasPsi = true;
	}

	//Check that distribution for opRank has not changed since last call.
	//Also check that typeID of representation is the same
/*	int curDistribOpRank = Psi->GetRepresentation()->GetDistributedModel()->GetDistribution()(opRank);
	int srcDistribOpRank = srcPsi.GetRepresentation()->GetDistributedModel()->GetDistribution()(opRank);
	if ( (curDistribOpRank != srcDistribOpRank) )
	{
		Psi = srcPsi.Copy();

		//NB: We reset IsSetup flag for _all_ ranks!
		IsSetupRank = false;
	}
*/
	if (!IsSetupRank(opRank))
	{
		//Sanity check: operation rank should be less than rank of wavefunction (and nonzero, duh)
		assert(opRank < Rank);
		assert(opRank > -1);

		//Assert non-orthogonal rank opRank
		assert (!srcPsi.GetRepresentation()->IsOrthogonalBasis(opRank));

		//Create Epetra map for this rank
		WavefunctionMaps(opRank) = CreateWavefunctionMultiVectorEpetraMap<Rank>(Psi, opRank);

		//Setup overlap matrix
		SetupOverlapMatrixRank(srcPsi, opRank);

		//Setup work multivectors
		blitz::Array<cplx, 3> psiData = MapToRank3(srcPsi.GetData(), opRank, 1);
		int numVectors = 2 * psiData.extent(0) * psiData.extent(2);
		InputVector(opRank) = Epetra_MultiVector_Ptr( new Epetra_MultiVector(*WavefunctionMaps(opRank), numVectors, false) );
		OutputVector(opRank) = Epetra_MultiVector_Ptr( new Epetra_MultiVector(*WavefunctionMaps(opRank), numVectors, false) );

		//Allocate mem for multivectors
		//InputData.resize(srcPsi.GetData().size(), 1);
		//OutputData.resize(srcPsi.GetData().size(), 1);
	
		//Setup Amesos solvers
		SetupOverlapSolvers(srcPsi, opRank);
		
		//Flag this rank as set up
		IsSetupRank(opRank) = true;
	}
}

void CombinedRepresentation<Rank>::MultiplyOverlap(cplx sourceScaling, Wavefunction<Rank> &srcPsi, cplx destScaling, Wavefunction<Rank> &dstPsi, int rank)
{
	//Not implemented correctly atm.
	throw std::runtime_error("Not implemented properly yet!");

	using namespace blitz;

	if (this->IsOrthogonalBasis(rank))
	{
		dstPsi.GetData() = srcPsi.GetData();
	}
	else
	{
		//Map the data to a 3D array, where the b-spline part is the middle rank
		//Array<cplx, 3> srcData = MapToRank3(srcPsi.Data, rank, 1);
		//Array<cplx, 3> dstData = MapToRank3(dstPsi.Data, rank, 1);

		//this->GetGlobalOverlapMatrix(rank)->MultiplyOverlapTensor(srcData, dstData);
		DistributedOverlap->MultiplyOverlapRank(srcPsi, dstPsi, rank, true);
	}

}

void BSplineTransform<Rank>::InverseTransform(Wavefunction<Rank> &psi)
{
	using namespace blitz;

	if(psi.GetActiveBufferName() != BSplineDataName)
	{
		throw std::runtime_error("Active databuffer is not what is should be...");
	}

	//Project input and output into 3D arrays with bsplineRank in the middle. 
	Array<cplx, Rank> inputData(psi.GetData());
	Array<cplx, Rank> outputData(psi.GetData(BSplineGridDataName));
	Array<cplx, 3> input3d = MapToRank3(inputData, BaseRank, 1);
	Array<cplx, 3> output3d = MapToRank3(outputData, BaseRank, 1);

	int psiSliceExtent = input3d.extent(1);
	Array<cplx, 1> psiSlice(psiSliceExtent);

	output3d = 0;
	int preCount = input3d.extent(0);
	int postCount = input3d.extent(2);
	for (int i=0; i<preCount; i++)
	{
		for (int j=0; j<postCount; j++)
		{
			/*
			 * Copy wavefunction slice along bspline rank to temp
			 * array.
			 */
			psiSlice = input3d(i, Range::all(), j).copy();

			// Call on BSpline function to perform expansion
			output3d(i, Range::all() , j) = 
				BSplineObject->ConstructFunctionFromBSplineExpansion(psiSlice);
		}
	}
	psi.SetActiveBuffer(BSplineGridDataName);
}

void BSplineTransform<Rank>::SetupStep(Wavefunction<Rank> &psi, BSpline::Ptr bsplineObject, int baseRank)
{
	using namespace blitz;
	
	SetBaseRank(baseRank);
	BSplineObject = bsplineObject;

	/*
	 * Get shape of wavefunction (we are in the bspline repr). 
	 * Then allocate bspline grid representation buffer of wavefunction
	 * and store buffer names on object.
	 */
	TinyVector<int, Rank> gridShape = psi.GetData().shape();
	gridShape(baseRank) = BSplineObject->GetQuadratureGridGlobal().extent(0);
	BSplineGridDataName = psi.AllocateData(gridShape);
	BSplineDataName = psi.GetActiveBufferName();

	//TempData.resize(BSplineObject->GetQuadratureGridGlobal().extent(0));
	TempData.resize(BSplineObject->NumberOfBSplines);
}

void CombinedRepresentation<Rank>::MultiplyIntegrationWeights(Wavefunction<Rank> &psi, int rank)
{
	blitz::Array<cplx, 3> data = MapToRank3(psi.GetData(), rank, 1);
	if (this->IsOrthogonalBasis(rank))
	{
		blitz::Array<double, 1> weights = this->GetLocalWeights(rank);
		data *= weights(blitz::tensor::j) + 0*blitz::tensor::k;
	}
	else
	{
		if (this->GetDistributedModel()->IsDistributedRank(rank))
		{
			DistributedOverlap->MultiplyOverlapRank(psi, rank, true);
		}
		else
		{
			this->GetGlobalOverlapMatrix(rank)->MultiplyOverlapTensor(data);
		}
	}
}

cplx CombinedRepresentation<Rank>::InnerProduct(const Wavefunction<Rank>& w1, const Wavefunction<Rank>& w2)
{
	blitz::Array<cplx, Rank> d1(w1.GetData());
	blitz::Array<cplx, Rank> d2(w2.GetData());

	/*
	 * Algorithm1 is faster for orthogonal basises
	 * Algorithm2 is faster for nonorthogonal basies
	 * Algorithm3 is the only one working for parallel problems, but require
	 *    more memory
	 *
	 * For a combination of orthogonal and non-orthogonal
	 * basises, 1 and 2 are most likely almost equally fast 
	 *
	 * Conclusion: Algo 3 is default
	 */

	if (Algorithm == 1)
	{
		return InnerProductImpl_Algo1(d1, d2);
	}
	else if (Algorithm == 2)
	{
		return InnerProductImpl_Algo2(d1, d2);
	}
	else if (Algorithm == 3)
	{
		blitz::TinyVector<int, Rank> shape = d1.shape();
		
		blitz::Array<cplx, Rank> temp1;
		blitz::Array<cplx, Rank> temp2;
		
		int tempName[2];
		int tempNamePsi[2];
		Wavefunction<Rank>* psiList[2];
		for (int i=0; i<2; i++)
		{
			tempName[i] = -1;
			tempNamePsi[i] = -1;
		}
		psiList[0] = const_cast<Wavefunction<Rank>*>(&w1);
		psiList[1] = const_cast<Wavefunction<Rank>*>(&w2);

		//Find any available buffers of correct size on any of the wavefunctions
		for (int i=0; i<2; i++)
		{
			//See if there is an available buffer in psi j
			for (int j=0; j<2; j++)
			{
				int name = psiList[j]->GetAvailableDataBufferName(shape);
				if (name != -1)
				{
					tempName[i] = name;
					tempNamePsi[i] = j;
					psiList[j]->LockBuffer(name);
					break;
				}
			}
		}

		//If we didnt find two available buffers, we must allocate
		//We'll allocate on w2
		for (int i=0; i<2; i++)
		{
			if (tempName[i] == -1)
			{
				tempName[i] = psiList[1]->AllocateData(shape);
				tempNamePsi[i] = 1;
				psiList[1]->LockBuffer(tempName[i]);
			}
		}

		//Get the actual data buffers
		temp1.reference(psiList[tempNamePsi[0]]->GetData(tempName[0]));
		temp2.reference(psiList[tempNamePsi[1]]->GetData(tempName[1]));

		//Perform MatrixVector multiplication
		//first step
		//
		
		for (int i=0; i<Rank; i++)
		{
			if (this->GetDistributedModel()->IsDistributedRank(i) && !this->IsOrthogonalBasis(i))
			{
				throw std::runtime_error("This inner product only supports parallelization for orthogonal ranks");
			}

			if (this->IsOrthogonalBasis(i))
			{
				if (i == 0)
				{
					temp1 = d2;
				}
				//TODO: Make this faster by moving it to TensorMultiply
				blitz::Array<double, 1> weights = this->GetLocalWeights(i);
				blitz::Array<cplx, 3> temp3d = MapToRank3(temp1, i, 1);
				temp3d *= weights(blitz::tensor::j) + 0*blitz::tensor::k;
			}
			else
			{
//.........这里部分代码省略.........