void SwappedSystem::reinitialize_phibetas (const Wavefunction<amplitude_t>::Amplitude &phialpha1, const Wavefunction<amplitude_t>::Amplitude &phialpha2)
{
    assert(subsystem_particle_counts_match());

#if defined(DEBUG_VMC_SWAPPED_SYSTEM) || defined(DEBUG_VMC_ALL)
    for (unsigned int species = 0; species < copy1_subsystem_indices.size(); ++species)
        std::cerr << "swapping " << copy1_subsystem_indices[species].size() << " particles of species " << species << std::endl;
    std::cerr << std::endl;
#endif

    PositionArguments swapped_r1(phialpha1.get_positions()), swapped_r2(phialpha2.get_positions());
    swap_positions(swapped_r1, swapped_r2);

    phibeta1 = phialpha1.clone();
    phibeta1->reset(swapped_r1);
    phibeta1_dirty = false;

    phibeta2 = phialpha2.clone();
    phibeta2->reset(swapped_r2);
    phibeta2_dirty = false;

#ifdef VMC_CAREFUL
    verify_phibetas(phialpha1, phialpha2);
#endif
}

//------------------------------------------------------------------------------
void BasisHarmonicOscillator::createInitalDiscretization()
{
    C = zeros<cx_mat>(nGrid, nSpatialOrbitals);
    mat Ctmp(nGrid, nSpatialOrbitals);
    Wavefunction *wf;

    for(int j=0; j<nSpatialOrbitals; j++){
        wf = new HarmonicOscillator(cfg, states[2*j]);
        for(int i=0; i<nGrid; i++){
            Ctmp(i,j) = wf->evaluate(grid->at(i));
        }
        delete wf;
    }

    C.set_real(Ctmp);
    // Forcing orthogonality by performing a SVD decomposition
    SVD(C);
#ifdef DEBUG
    cout << "BasisHarmonicOscillator::discretization1d()" << endl;
    for(int i=0; i<nSpatialOrbitals; i++){
        for(int j=0; j<nSpatialOrbitals; j++){
            cout <<  "|C(" << i << ", " << j << ")| = " << sqrt(cdot(C.col(i),C.col(j))) << endl ;
        }
    }
#endif
}

/*
 * 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 SwappedSystem::verify_phibetas (const Wavefunction<amplitude_t>::Amplitude &phialpha1, const Wavefunction<amplitude_t>::Amplitude &phialpha2) const
{
#ifdef NDEBUG
    (void) phialpha1;
    (void) phialpha2;
#else
    const PositionArguments &r1 = phialpha1.get_positions();
    const PositionArguments &r2 = phialpha2.get_positions();

    assert(r1.get_N_species() == r2.get_N_species());
    assert(r1.get_N_sites() == r2.get_N_sites());

    assert(copy1_subsystem_indices.size() == r1.get_N_species());
    assert(copy2_subsystem_indices.size() == r1.get_N_species());

    const Lattice &lattice = phialpha1.get_lattice();

    for (unsigned int species = 0; species < r1.get_N_species(); ++species) {
        const unsigned int N = r1.get_N_filled(species);
        assert(N == r2.get_N_filled(species));

        // verify that the subsystem index arrays have everything they need (and no duplicates!)
        unsigned int c1 = 0, c2 = 0;
        for (unsigned int i = 0; i < N; ++i) {
            const Particle particle(i, species);
            const bool b1 = vector_find(copy1_subsystem_indices[species], i) != -1;
            const bool b2 = vector_find(copy2_subsystem_indices[species], i) != -1;
            if (b1)
                ++c1;
            if (b2)
                ++c2;
            assert(b1 == subsystem->position_is_within(r1[particle], lattice));
            assert(b2 == subsystem->position_is_within(r2[particle], lattice));
        }
        assert(c1 == c2);
        assert(c1 == copy1_subsystem_indices[species].size());
        assert(c2 == copy2_subsystem_indices[species].size());
    }

    assert(phibeta1 != 0);
    assert(phibeta2 != 0);

    // verify that the positions in the phibeta's are correct
    PositionArguments swapped_r1(phialpha1.get_positions()), swapped_r2(phialpha2.get_positions());
    swap_positions(swapped_r1, swapped_r2);

    for (unsigned int species = 0; species < r1.get_N_species(); ++species) {
        for (unsigned int i = 0; i < r1.get_N_filled(species); ++i) {
            const Particle particle(i, species);
            assert(swapped_r1[particle] == phibeta1->get_positions()[particle]);
            assert(swapped_r2[particle] == phibeta2->get_positions()[particle]);
        }
    }
#endif
}

//solves the equation (H-E)|psi_1> = -QV|\Psi_0>  where lowerState[0] contains |\Psi_0> to enforce orthogonality (Q), and lowerState[1] contains V|\Psi_0> so we can calculate its projection in the krylov space
//the algorithm is taken from wikipedia page
void SpinAdapted::Linear::ConjugateGradient(Wavefunction& xi, DiagonalMatrix& h_diag, double E0, double normtol, Davidson_functor& h_multiply, std::vector<Wavefunction> &lowerStates)
{

    pout.precision (12);
#ifndef SERIAL
    mpi::communicator world;
#endif
    int iter = 0;
    double levelshift = 0.0;

    Wavefunction b = lowerStates[1];
    //make b[0] orthogonal to lowerStates[0]
    double overlap = DotProduct(b, lowerStates[0]);
    double overlap2 = DotProduct(lowerStates[0], lowerStates[0]);
    ScaleAdd(-overlap/overlap2, lowerStates[0], b);

    //normalise all the guess roots
    if(mpigetrank() == 0) {
        Normalise(xi);
    }

    Wavefunction pi, ri;
    ri=xi;
    ri.Clear();
    h_multiply(xi, ri);
    ScaleAdd(-E0, xi, ri); // (H-E0)|psi>

    ScaleAdd(-1.0, lowerStates[1], ri);
    Scale(-1.0, ri);
    pi = ri;

    double oldError = DotProduct(ri, ri);
    while(true) {
        Wavefunction Hp = pi;
        Hp.Clear();
        h_multiply(pi, Hp);
        ScaleAdd(-E0, pi, Hp); // (H-E0)|psi>

        double alpha = oldError/DotProduct(pi, Hp);

        ScaleAdd(alpha, pi, xi);
        ScaleAdd(-alpha, Hp, ri);

        double Error = DotProduct(ri, ri);
        if (Error > normtol)
            return;
        else {
            double beta = Error/oldError;
            oldError = Error;
            ScaleAdd(1.0/beta, ri, pi);
            Scale(beta, pi);
        }
    }
}

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 SpinAdapted::operatorfunctions::TensorMultiply(const SpinBlock *ablock, const Baseoperator<Matrix>& a, const Baseoperator<Matrix>& b, const SpinBlock *cblock, Wavefunction& c, Wavefunction& v, const SpinQuantum opQ, double scale)
{
  // can be used for situation with different bra and ket
  const int leftBraOpSz = cblock->get_leftBlock()->get_braStateInfo().quanta.size ();
  const int leftKetOpSz = cblock->get_leftBlock()->get_ketStateInfo().quanta.size ();
  const int rightBraOpSz = cblock->get_rightBlock()->get_braStateInfo().quanta.size ();
  const int rightKetOpSz = cblock->get_rightBlock()->get_ketStateInfo().quanta.size ();

  const StateInfo* lbraS = cblock->get_braStateInfo().leftStateInfo, *rbraS = cblock->get_braStateInfo().rightStateInfo;
  const StateInfo* lketS = cblock->get_ketStateInfo().leftStateInfo, *rketS = cblock->get_ketStateInfo().rightStateInfo;

  const char conjC = (cblock->get_leftBlock() == ablock) ? 'n' : 't';

  const Baseoperator<Matrix>& leftOp = (conjC == 'n') ? a : b; // an ugly hack to support the release memory optimisation
  const Baseoperator<Matrix>& rightOp = (conjC == 'n') ? b : a;
  const char leftConj = (conjC == 'n') ? a.conjugacy() : b.conjugacy();
  const char rightConj = (conjC == 'n') ? b.conjugacy() : a.conjugacy();


  int totalmem =0;

  for (int lQrQPrime = 0; lQrQPrime<leftBraOpSz*rightKetOpSz; ++lQrQPrime)
  {
    int rQPrime = lQrQPrime%rightKetOpSz, lQ = lQrQPrime/rightKetOpSz;
    for (int lQPrime = 0; lQPrime < leftKetOpSz; lQPrime++)
      if (leftOp.allowed(lQ, lQPrime) && c.allowed(lQPrime, rQPrime))
      {	    
	Matrix m; m.ReSize(lbraS->getquantastates(lQ), rketS->getquantastates(rQPrime));
        
	double factor = leftOp.get_scaling(lbraS->quanta[lQ], lketS->quanta[lQPrime]);
	MatrixMultiply (leftOp.operator_element(lQ, lQPrime), leftConj, c.operator_element(lQPrime, rQPrime), 'n',
			m, factor, 0.);	      
	
	for (int rQ = 0; rQ<rightBraOpSz; rQ++) {
	  if (v.allowed(lQ, rQ) && rightOp.allowed(rQ, rQPrime)) {
	    double factor = scale;
	    
	    factor *= dmrginp.get_ninej()(lketS->quanta[lQPrime].get_s().getirrep(), rketS->quanta[rQPrime].get_s().getirrep() , c.get_deltaQuantum(0).get_s().getirrep(), 
					  leftOp.get_spin().getirrep(), rightOp.get_spin().getirrep(), opQ.get_s().getirrep(),
					  lbraS->quanta[lQ].get_s().getirrep(), rbraS->quanta[rQ].get_s().getirrep() , v.get_deltaQuantum(0).get_s().getirrep());
	    factor *= Symmetry::spatial_ninej(lketS->quanta[lQPrime].get_symm().getirrep() , rketS->quanta[rQPrime].get_symm().getirrep(), c.get_symm().getirrep(), 
					      leftOp.get_symm().getirrep(), rightOp.get_symm().getirrep(), opQ.get_symm().getirrep(),
					      lbraS->quanta[lQ].get_symm().getirrep() , rbraS->quanta[rQ].get_symm().getirrep(), v.get_symm().getirrep());
	    int parity = rightOp.get_fermion() && IsFermion(lketS->quanta[lQPrime]) ? -1 : 1;
	    factor *=  rightOp.get_scaling(rbraS->quanta[rQ], rketS->quanta[rQPrime]);
	    MatrixMultiply (m, 'n', rightOp(rQ, rQPrime), TransposeOf(rightOp.conjugacy()), v.operator_element(lQ, rQ), factor*parity);
	  }
	}
	
      }
  }

}

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

shared_ptr<Epetra_MultiVector> DistributedOverlapMatrix<Rank>::SetupMultivector(Wavefunction<Rank> &srcPsi, int opRank)
{
	//Setup Epetra stuff for opRank
	SetupRank(srcPsi, opRank);

	//Map wavefunction to 3D array (compress before- and after-ranks)
	blitz::Array<cplx, 3> psiData = MapToRank3(srcPsi.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 2D array
	blitz::Array<double, 2> data;
	data.resize(2 * otherSize, opSize);
	for (int i=0; i<beforeSize; i++)
	{
		for (int j=0; j<opSize; j++)
		{
			for (int k=0; k<afterSize; k++)
			{
				data(i*afterSize + k, j) = real( psiData(i,j,k) );
				data(i*afterSize + k + otherSize, j) = imag( psiData(i,j,k) );
			}
		}
	}

	//Number of vectors (in multivector)
	int numVectors = 2 * otherSize; 

	//Create Epetra multivector (view of data)
	Epetra_MultiVector_Ptr srcVec = Epetra_MultiVector_Ptr(new Epetra_MultiVector(Copy, *WavefunctionMaps(opRank), data.data(), opSize, numVectors));

	return srcVec;
}

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

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

bool count_subsystem_particle_counts_for_match (const Wavefunction<amplitude_t>::Amplitude &wf1, const Wavefunction<amplitude_t>::Amplitude &wf2,
                                                const Subsystem &subsystem)
{
    assert(subsystem.lattice_makes_sense(wf1.get_lattice()));
    assert(subsystem.lattice_makes_sense(wf2.get_lattice()));
    // (we are also assuming that the lattices are in fact equivalent)

    const PositionArguments &r1 = wf1.get_positions();
    const PositionArguments &r2 = wf2.get_positions();

    assert(r1.get_N_species() == r2.get_N_species());
    assert(r1.get_N_sites() == r2.get_N_sites());

    for (unsigned int species = 0; species < r1.get_N_species(); ++species) {
        assert(r1.get_N_filled(species) == r2.get_N_filled(species));

        unsigned int count1 = 0, count2 = 0;

        for (unsigned int i = 0; i < r1.get_N_filled(species); ++i) {
            const Particle particle(i, species);
            if (subsystem.position_is_within(r1[particle], wf1.get_lattice()))
                ++count1;
            if (subsystem.position_is_within(r2[particle], wf2.get_lattice()))
                ++count2;
        }

        if (count1 != count2)
            return false;
    }

    return 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 SwappedSystem::initialize (const Wavefunction<amplitude_t>::Amplitude &phialpha1, const Wavefunction<amplitude_t>::Amplitude &phialpha2)
{
    assert(current_state == UNINITIALIZED);

    const PositionArguments &r1 = phialpha1.get_positions();
    const PositionArguments &r2 = phialpha2.get_positions();

    // FIXME: we need a way to assert that phialpha1 and phialpha2 represent
    // the same wave function, just with different amplitudes.  Then again,
    // we're only calling this function from two places in the code where this
    // can be easily verified ...

#ifndef NDEBUG
    assert(r1.get_N_species() == r2.get_N_species());
    for (unsigned int i = 0; i < r1.get_N_species(); ++i)
        assert(r1.get_N_filled(i) == r2.get_N_filled(i));
#endif
    assert(r1.get_N_sites() == r2.get_N_sites());
    assert(subsystem->lattice_makes_sense(phialpha1.get_lattice()));
    assert(subsystem->lattice_makes_sense(phialpha2.get_lattice()));
    assert(&phialpha1.get_lattice() == &phialpha2.get_lattice());

    const unsigned int N_species = r1.get_N_species();
    copy1_subsystem_indices.resize(N_species);
    copy2_subsystem_indices.resize(N_species);
    for (unsigned int species = 0; species < N_species; ++species) {
        const unsigned int N = r1.get_N_filled(species);
        for (unsigned int i = 0; i < N; ++i) {
            const Particle particle(i, species);
            if (subsystem->position_is_within(r1[particle], phialpha1.get_lattice()))
                copy1_subsystem_indices[species].push_back(i);
            if (subsystem->position_is_within(r2[particle], phialpha2.get_lattice()))
                copy2_subsystem_indices[species].push_back(i);
        }
    }

    assert(subsystem_particle_counts_match());
    reinitialize_phibetas(phialpha1, phialpha2);

    current_state = READY;
}

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

void Propagator<Rank>::Setup(const Parameter &param, const cplx &dt, const Wavefunction<Rank> &psi, int rank)
{
	firstIndex i;
	secondIndex j;
	thirdIndex k;

	//Set class parameters
	N = psi.GetRepresentation()->GetFullShape()(rank);
	PropagateRank = rank;
	Param = param;

	//create some temporary arrays
	Array<double, 2> evInvReal; //inverse eigenvector matrix

	//Call setup routines to get the eigenvector decomposition
	setup(N, param, Eigenvalues, Eigenvectors, evInvReal);

	//We need the eigenvalues and vectors in complex format
	Array<cplx, 2> evExp(Eigenvectors.shape());     //eigenvectors scaled by exp(eigenvalues)
	Array<cplx, 2> evDiff(Eigenvectors.shape());    //eigenvectors scaled by eigenvalues
	Array<cplx, 2> evInv(Eigenvectors.shape());  //inverse eigenvector matrix

	evInv = evInvReal(tensor::i, tensor::j);

	//scale eigenvectors by complex rotation
	//The missing minus sign in the exponent is included in the matrix.
	evExp = Eigenvectors(i,j) * exp( I * dt * Eigenvalues(j) / (2.0 * Mass));
	evDiff = - Eigenvectors(i,j) * Eigenvalues(j) / (2.0 * Mass);

	//Create full matrix to propagate wavefunction
	PropagationMatrix.resize(Eigenvectors.shape());
	MatrixMatrixMultiply(evExp, evInv, PropagationMatrix);

	//Create full differentiation matrix
	DiffMatrix.resize(Eigenvectors.shape());
	MatrixMatrixMultiply(evDiff, evInv, DiffMatrix);

	//Allocate temp data
	TempData.resize(Eigenvalues.extent(0));
}

void DistributedOverlapMatrix<Rank>::SetupOverlapSolvers(Wavefunction<Rank> &psi, int opRank)
{
	//Set up Epetra LinearProblem with overlap for this rank and input/output multivectors
	EpetraProblems(opRank) = Epetra_LinearProblem_Ptr( new Epetra_LinearProblem(OverlapMatrices(opRank).get(), OutputVector(opRank).get(), InputVector(opRank).get()) );

	//Determine solver type. Use SuperLU_dist if opRank is distributed, else KLU
	Amesos Factory;
	std::string SolverType;
	if (psi.GetRepresentation()->GetDistributedModel()->IsDistributedRank(opRank))
	{
		SolverType = "Amesos_Superludist";
	}
	else
	{
		SolverType = "Amesos_Klu";
	}

	//Create Amesos solver for this rank
	Solvers(opRank) = Amesos_BaseSolver_Ptr( Factory.Create(SolverType, *EpetraProblems(opRank)) );

	//Check that requested solver exists in Amesos build
	if (Solvers(opRank) == 0)
	{
		throw std::runtime_error("Specified Amesos solver not available");
	}

	//Setup the parameter list for the solver (TODO: get these from Python config section)
	Teuchos::ParameterList List;
	List.set("MatrixType", "symmetric");
	List.set("PrintTiming", false);
	List.set("PrintStatus", false);
	List.set("ComputeTrueResidual", false);
	Solvers(opRank)->SetParameters(List);

	//Perform symbolic factorization
	Solvers(opRank)->SymbolicFactorization();
	Solvers(opRank)->NumericFactorization();
}

void SpinAdapted::Linear::precondition(Wavefunction& op, double e, DiagonalMatrix& diagonal, double levelshift)
{
    if (!mpigetrank())
    {
        int index = 1;
        for (int lQ = 0; lQ < op.nrows (); ++lQ)
            for (int rQ = 0; rQ < op.ncols (); ++rQ)
                if (op.allowed(lQ, rQ))
                    for (int lQState = 0; lQState < op.operator_element(lQ, rQ).Nrows (); ++lQState)
                        for (int rQState = 0; rQState < op.operator_element(lQ, rQ).Ncols (); ++rQState)
                        {
                            if (fabs(e - diagonal(index)) > 1.e-12)
                                op.operator_element(lQ, rQ).element(lQState, rQState) /= (e - diagonal(index)+levelshift);
                            ++index;
                        }
    }
#ifndef SERIAL
    //mpi::communicator world;
    //mpi::broadcast(world, op, 0);
#endif
}