C++ SpinBlock类代码示例

  template<> void Op_component<CreCreDesComp>::build_iterators(SpinBlock& b)
    {
      if (b.get_sites().size () == 0) return; // blank construction (used in unset_initialised() Block copy construction, for use with STL)
      const double screen_tol = dmrginp.oneindex_screen_tol();
      vector< int > screened_cdd_ix = (dmrginp.hamiltonian() == BCS) ?
        screened_cddcomp_indices(b.get_complementary_sites(), b.get_sites(), v_1, *b.get_twoInt(), v_cc, v_cccc, v_cccd, screen_tol) :
        screened_cddcomp_indices(b.get_complementary_sites(), b.get_sites(), v_1, *b.get_twoInt(), screen_tol);
      m_op.set_indices(screened_cdd_ix, dmrginp.last_site());      
      std::vector<int> orbs(1);
      for (int i = 0; i < m_op.local_nnz(); ++i)
	{
	  orbs[0] = m_op.get_local_indices()[i];
	  m_op.get_local_element(i).resize(1);
	  m_op.get_local_element(i)[0]=boost::shared_ptr<CreCreDesComp>(new CreCreDesComp);
	  SparseMatrix& op = *m_op.get_local_element(i)[0];
	  op.set_orbs() = orbs;
	  op.set_initialised() = true;
	  op.set_fermion() = true;
	  //op.set_deltaQuantum() = SpinQuantum(1, SpinOf(orbs[0]), SymmetryOfSpatialOrb(orbs[0]) );
      if (dmrginp.hamiltonian() == BCS) {
        op.resize_deltaQuantum(4);
        SpinQuantum qorb = getSpinQuantum(orbs[0]);
        op.set_deltaQuantum(0) = qorb;
        op.set_deltaQuantum(1) = SpinQuantum(3, qorb.get_s(), qorb.get_symm());
        op.set_deltaQuantum(2) = SpinQuantum(-1, qorb.get_s(), qorb.get_symm());
        op.set_deltaQuantum(3) = SpinQuantum(-3, qorb.get_s(), qorb.get_symm());
      } else {
	    op.set_deltaQuantum(1, getSpinQuantum(orbs[0]));
      }
	}
    }

void do_4index_2_2_tensor_products( bool forwards, const opTypes& optype, const opTypes& rhsType, const opTypes& lhsType,
                                    SpinBlock& big, SpinBlock* rhsBlock, SpinBlock* lhsBlock, std::ofstream& ofs,
                                    const std::vector<Matrix>& rotateMatrix, const StateInfo *stateinfo )
{
  Op_component_base& rhs_array = rhsBlock->get_op_array(rhsType);
  Op_component_base& lhs_array = lhsBlock->get_op_array(lhsType);
  assert ( (rhs_array.get_size() == 1) || (lhs_array.get_size() == 1) );

  // Loop over all rhs operator indices //FIXME copy or reference?
  std::vector<boost::shared_ptr<SparseMatrix> > rhs_ops;
  for (int iidx = 0; iidx < rhs_array.get_size(); ++iidx) {
    rhs_ops = rhs_array.get_local_element(iidx);

    // Loop over all lhs operator indices
    std::vector<boost::shared_ptr<SparseMatrix> > lhs_ops;
    for (int idx = 0; idx < lhs_array.get_size(); ++idx) {
      lhs_ops = lhs_array.get_local_element(idx);
      int i = rhs_ops[0]->get_orbs()[0]; int j = rhs_ops[0]->get_orbs()[1];
      int k = lhs_ops[0]->get_orbs()[0]; int l = lhs_ops[0]->get_orbs()[1];
      // In parallel calculations not all operators are built on each proc
      if ( ! big.get_op_array(optype).has_local_index(i,j,k,l) ) continue;
//pout << "p" << mpigetrank() << "; i,j,k,l = " << i << "," << j << "," << k << "," << l << endl;
      std::vector<boost::shared_ptr<SparseMatrix> > vec = big.get_op_array(optype).get_element(i,j,k,l);

      // Loop over rhs spin-op components
      for (int jjdx=0; jjdx < rhs_ops.size(); jjdx++) {
        boost::shared_ptr<SparseMatrix>& rhs_op = rhs_ops[jjdx];
        assert( rhs_op->get_built() );
        const SpinQuantum& spin_12 = rhs_op->get_deltaQuantum(0);

        // Loop over lhs spin-op components
        for (int jdx=0; jdx < lhs_ops.size(); jdx++) {
          boost::shared_ptr<SparseMatrix>& lhs_op = lhs_ops[jdx];
          assert( lhs_op->get_built() );
          std::string build_12 = rhs_op->get_build_pattern();
          std::string build_34 = lhs_op->get_build_pattern();
          std::string build_pattern = "(" + build_12 + build_34 + ")";
          const SpinQuantum& spin_34 = lhs_op->get_deltaQuantum(0);

          // Allocate and build new operator
          for (int sx=0; sx < vec.size(); sx++) {
            boost::shared_ptr<SparseMatrix>& op = vec[sx];
            std::vector<SpinQuantum> s1 = { op->get_quantum_ladder().at(build_pattern).at(0), op->get_quantum_ladder().at(build_pattern).at(1) };
            std::vector<SpinQuantum> s2 = { spin_12, spin_34 };
            // Select relevant spin components
            if ( s1 == s2 ) {
              finish_tensor_product( big, rhsBlock, *rhs_op, *lhs_op, *op, forwards, build_pattern );
              // Renormalise operator
              op->renormalise_transform( rotateMatrix, stateinfo );
            }
          }
        }
      }
      // Store spin-batch on disk 
      if ( ! dmrginp.do_npdm_in_core() ) store_ops_on_disk( ofs, vec );
    }
  }

}

void SpinAdapted::InitBlocks::InitNewSystemBlock(SpinBlock &system, SpinBlock &systemDot, SpinBlock &newSystem, int leftState, int rightState, const int& sys_add, const bool &direct, int integralIndex, const Storagetype &storage, bool haveNormops, bool haveCompops, int constraint)
{
  newSystem.set_integralIndex() = integralIndex;
  newSystem.default_op_components(direct, system, systemDot, haveNormops, haveCompops, leftState==rightState);
  newSystem.setstoragetype(storage);
  newSystem.BuildSumBlock (constraint, system, systemDot);

  p2out << "\t\t\t NewSystem block " << endl << newSystem << endl;
  newSystem.printOperatorSummary();
}

void finish_tensor_trace( SpinBlock& b, SpinBlock* sysdot, SparseMatrix& sysdot_op, SparseMatrix& op, std::string& build_pattern )
{
  // Build and store new operator
  assert( ! op.get_built() );
  op.set_built() = true;
  op.set_build_pattern() = build_pattern;
//FIXME magic number 2
  op.set_deltaQuantum(1, op.get_quantum_ladder().at( build_pattern ).at(2) );
  op.allocate(b.get_stateInfo());

  SpinAdapted::operatorfunctions::TensorTrace(sysdot, sysdot_op, &b, &(b.get_stateInfo()), op);
}

void do_4index_tensor_trace( const opTypes& optype, SpinBlock& big, SpinBlock* sysdot, std::ofstream& ofs,
                             const std::vector<Matrix>& rotateMatrix, const StateInfo *stateinfo )
{
  // Get pointer to sparse operator array
  Op_component_base& sysdot_array = sysdot->get_op_array(optype);
  // Open filesystem if necessary
  std::ifstream ifs;
  if ( (! dmrginp.do_npdm_in_core()) && sysdot->size() > 1 ) ifs.open( sysdot_array.get_filename().c_str(), std::ios::binary );

//FIXME need reference?  don't want to copy?
  // Loop over all operator indices
  std::vector<boost::shared_ptr<SparseMatrix> > sysdot_ops;
  for (int idx = 0; idx < sysdot_array.get_size(); ++idx) {
    if ( dmrginp.do_npdm_in_core() || sysdot->size() <= 1) 
      sysdot_ops = sysdot_array.get_local_element(idx);
    else
//FIXME size
      sysdot_ops = get_ops_from_disk( ifs, sysdot_array.get_local_element(0).size() );

    // Loop over spin-op components
    int i = sysdot_ops[0]->get_orbs()[0]; int j = sysdot_ops[0]->get_orbs()[1]; 
    int k = sysdot_ops[0]->get_orbs()[2]; int l = sysdot_ops[0]->get_orbs()[3];
    // In parallel calculations not all operators are built on each proc
    if ( ! big.get_op_array(optype).has_local_index(i,j,k,l) ) continue;
//pout << "p" << mpigetrank() << "; i,j,k,l = " << i << "," << j << "," << k << "," << l << endl;
    std::vector<boost::shared_ptr<SparseMatrix> > new_ops = big.get_op_array(optype).get_element(i,j,k,l);
    for (int jdx=0; jdx < sysdot_ops.size(); jdx++) {
      boost::shared_ptr<SparseMatrix>& sysdot_op = sysdot_ops[jdx];
      assert( sysdot_op->get_built() );
      std::string build_pattern = sysdot_op->get_build_pattern();

      // Allocate and build new operator
      for (int sx=0; sx < new_ops.size(); sx++) {
        boost::shared_ptr<SparseMatrix>& op = new_ops[sx];
        std::vector<SpinQuantum> s1 = sysdot_op->get_quantum_ladder().at(build_pattern);
        std::vector<SpinQuantum> s2 = op->get_quantum_ladder().at(build_pattern);
        // Store spin component in correct location
        if ( s1 == s2 ) {
          finish_tensor_trace( big, sysdot, *sysdot_op, *op, build_pattern );
          // Renormalise operator
          op->renormalise_transform( rotateMatrix, stateinfo );
        }
      }
    }
    // Store spin-batch on disk 
    if ( ! dmrginp.do_npdm_in_core() ) store_ops_on_disk( ofs, new_ops );
  }

  // Close filesystem if necessary
  if ( ifs.is_open() ) ifs.close();

}

void SpinBlock::store (bool forward, const vector<int>& sites, SpinBlock& b, int left, int right, char *name)
{
  Timer disktimer;
  std::string file;
  if(dmrginp.spinAdapted()) {
    if (forward)
      file = str(boost::format("%s%s%d%s%d%s%d%s%d%s%d%s%d%s") % dmrginp.save_prefix() % "/SpinBlock-forward-"% sites[0] % "-" % sites[sites.size()-1] % "." % left % "." % right % "." %b.integralIndex % "." % mpigetrank() % ".tmp" );
    else
      file = str(boost::format("%s%s%d%s%d%s%d%s%d%s%d%s%d%s") % dmrginp.save_prefix() % "/SpinBlock-backward-"% sites[0] % "-" % sites[sites.size()-1] % "." % left % "." % right % "." %b.integralIndex % "." % mpigetrank() % ".tmp" );
  }
  else {
    if (forward)
      file = str(boost::format("%s%s%d%s%d%s%d%s%d%s%d%s%d%s") % dmrginp.save_prefix() % "/SpinBlock-forward-"% (sites[0]/2) % "-" % (sites[sites.size()-1]/2) % "." % left % "." % right % "." %b.integralIndex % "." % mpigetrank() % ".tmp" );
    else
      file = str(boost::format("%s%s%d%s%d%s%d%s%d%s%d%s%d%s") % dmrginp.save_prefix() % "/SpinBlock-backward-"% (sites[0]/2) % "-" % (sites[sites.size()-1]/2) % "." % left % "." % right % "." %b.integralIndex % "." % mpigetrank() % ".tmp" );
  }
  
  p1out << "\t\t\t Saving block file :: " << file << endl;
  
  
  std::ofstream ofs(file.c_str(), std::ios::binary);
  
  int lstate =  left;
  int rstate =  right;
  
  if (mpigetrank()==0) {
    StateInfo::store(forward, sites, b.braStateInfo, lstate);
    StateInfo::store(forward, sites, b.ketStateInfo, rstate);
  }

  b.Save (ofs);
  ofs.close(); 
  //p1out << "\t\t\t block save disk time " << disktimer.elapsedwalltime() << " " << disktimer.elapsedcputime() << endl;
}

double SweepOnepdm::do_one(SweepParams &sweepParams, const bool &warmUp, const bool &forward, const bool &restart, const int &restartSize)
{

  SpinBlock system;
  const int nroots = dmrginp.nroots();
  std::vector<double> finalEnergy(nroots,0.);
  std::vector<double> finalEnergy_spins(nroots,0.);
  double finalError = 0.;

  Matrix onepdm(2*dmrginp.last_site(), 2*dmrginp.last_site());onepdm=0.0;
  for (int i=0; i<nroots; i++)
    for (int j=0; j<=i; j++)      
      save_onepdm_binary(onepdm, i ,j);

  sweepParams.set_sweep_parameters();
  // a new renormalisation sweep routine
  pout << ((forward) ? "\t\t\t Starting renormalisation sweep in forwards direction" : "\t\t\t Starting renormalisation sweep in backwards direction") << endl;
  pout << "\t\t\t ============================================================================ " << endl;
  
  InitBlocks::InitStartingBlock (system,forward, sweepParams.get_forward_starting_size(), sweepParams.get_backward_starting_size(), restartSize, restart, warmUp);

  sweepParams.set_block_iter() = 0;
 
  pout << "\t\t\t Starting block is :: " << endl << system << endl;

  SpinBlock::store (forward, system.get_sites(), system); // if restart, just restoring an existing block --
  sweepParams.savestate(forward, system.get_sites().size());
  bool dot_with_sys = true;

  sweepParams.set_guesstype() = TRANSPOSE;
	  
  SpinBlock newSystem;
  BlockAndDecimate (sweepParams, system, newSystem, warmUp, dot_with_sys);
  pout.precision(12);
  pout << "\t\t\t The lowest sweep energy : "<< sweepParams.get_lowest_energy()[0]+dmrginp.get_coreenergy()<<endl;
  pout << "\t\t\t ============================================================================ " << endl;

  for (int i=0; i<nroots; i++)
    for (int j=0; j<=i; j++) {
      load_onepdm_binary(onepdm, i ,j);
      accumulate_onepdm(onepdm);
      save_onepdm_spatial_text(onepdm, i ,j);
      save_onepdm_text(onepdm, i ,j);
      save_onepdm_spatial_binary(onepdm, i ,j);
    }
  return sweepParams.get_lowest_energy()[0];
}

void build_3index_ops( const opTypes& optype, SpinBlock& big, 
                       const opTypes& lhsType1, const opTypes& lhsType2,
                       const opTypes& rhsType1, const opTypes& rhsType2,
                       const std::vector<Matrix>& rotateMatrix, const StateInfo *stateinfo )
{
  // 3-index output file
//pout << "build_3index_op, ofs =" <<  big.get_op_array(optype).get_filename() << endl;
  std::ofstream ofs;
  if ( ! dmrginp.do_npdm_in_core() ) ofs.open( big.get_op_array(optype).get_filename().c_str(), std::ios::binary );

  SpinBlock* sysBlock = big.get_leftBlock();
  SpinBlock* dotBlock = big.get_rightBlock();

  // All 3 orbitals on sys or dot block
  do_3index_tensor_trace( optype, big, sysBlock, ofs, rotateMatrix, stateinfo );
  do_3index_tensor_trace( optype, big, dotBlock, ofs, rotateMatrix, stateinfo );

  bool forwards = ! ( sysBlock->get_sites().at(0) > dotBlock->get_sites().at(0) );

  // 2,1 partitioning
  if ( forwards ) {
    do_3index_1_2_tensor_products( forwards, optype, lhsType1, rhsType2, big, dotBlock, sysBlock, ofs, rotateMatrix, stateinfo );
    do_3index_2_1_tensor_products( forwards, optype, lhsType2, rhsType1, big, dotBlock, sysBlock, ofs, rotateMatrix, stateinfo );
  } else {
    do_3index_1_2_tensor_products( forwards, optype, lhsType1, rhsType2, big, sysBlock, dotBlock, ofs, rotateMatrix, stateinfo );
    do_3index_2_1_tensor_products( forwards, optype, lhsType2, rhsType1, big, sysBlock, dotBlock, ofs, rotateMatrix, stateinfo );
  }

  if ( ofs.is_open() ) ofs.close();

}

double SparseMatrix::memoryUsed(const SpinBlock& b)
{
  StateInfo stateinfo = b.get_stateInfo();
  double memory = 0.0;

  for (int i=0; i < stateinfo.quanta.size(); i++)
    for (int j=0; j<stateinfo.quanta.size(); j++)
      if (allowedQuantaMatrix(i,j) ) {
	memory += 8.0*operatorMatrix(i,j).Storage();
      }
  return memory;
}

//-------------------------------------------------------------------------------------------------------------------------------------------------------------
//  (Cre,Cre,Cre,Cre)
//-------------------------------------------------------------------------------------------------------------------------------------------------------------
void SpinAdapted::CreCreCreCre::build(const SpinBlock& b) { 
  dmrginp.makeopsT -> start();
  built = true;
  allocate(b.get_braStateInfo(), b.get_ketStateInfo());

  const int i = get_orbs()[0];
  const int j = get_orbs()[1];
  const int k = get_orbs()[2];
  const int l = get_orbs()[3];
  SpinBlock* leftBlock = b.get_leftBlock();
  SpinBlock* rightBlock = b.get_rightBlock();

  if (leftBlock->get_op_array(CRE_CRE_CRE_CRE).has(i,j,k,l))
  {      
    const boost::shared_ptr<SparseMatrix>& op = leftBlock->get_op_rep(CRE_CRE_CRE_CRE, quantum_ladder, i,j,k,l);
    if (rightBlock->get_sites().size() == 0) 
      SpinAdapted::operatorfunctions::TensorTrace(leftBlock, *op, &b, &(b.get_stateInfo()), *this);
    dmrginp.makeopsT -> stop();
    return;
  }
  assert(false && "Only build CRECRECRECRE in the starting block when spin-embeding is used");
}

/*
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();

  Wavefunction u;
  u.resize(leftBraOpSz*leftKetOpSz, rightKetOpSz);

  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))
	  {
	    int lindex = lQ*leftKetOpSz+lQPrime;
	    u.allowed(lindex, rQPrime) = true;
            
	    u(lindex,rQPrime).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',
			    u.operator_element(lindex, rQPrime), factor, 0.);	      

	  }
    }
  }

  pout << "after first step in tensormultiply"<<endl;
      mcheck("before davidson but after all blocks are built");

  {
    for (int lQrQ = 0; lQrQ<leftBraOpSz*rightBraOpSz; ++lQrQ)
    {
      int rQ = lQrQ%rightBraOpSz, lQ=lQrQ/rightBraOpSz;
	if (v.allowed(lQ, rQ))
	  for (int rQPrime = 0; rQPrime < rightKetOpSz; rQPrime++)
	    if (rightOp.allowed(rQ, rQPrime))
	      for (int lQPrime = 0; lQPrime < leftKetOpSz; lQPrime++)
		if (leftOp.allowed(lQ, lQPrime) && u.allowed(lQ*leftKetOpSz+lQPrime, rQPrime))
		{
		  int lindex = lQ*leftKetOpSz+lQPrime;
		  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 (u.operator_element(lindex, rQPrime), 'n',
				  rightOp(rQ, rQPrime), TransposeOf(rightOp.conjugacy()), v.operator_element(lQ, rQ), factor*parity);
		}
    }
  }
	      
}
*/
void SpinAdapted::operatorfunctions::OperatorScaleAdd(double scaleV, const SpinBlock& b, const Baseoperator<Matrix>& op1, Baseoperator<Matrix>& op2)
{
  const StateInfo& s = b.get_stateInfo();
  for (int lQ = 0; lQ< op2.nrows(); lQ++)
    for (int rQ = 0; rQ<op2.ncols(); rQ++)
      if (op2.allowed(lQ, rQ) && op1.allowed(lQ,rQ))
      {
	double factor = op1.get_scaling(s.quanta[lQ], s.quanta[rQ]);
	if (op1.conjugacy() == 't')
	  MatrixScaleAdd(scaleV*factor, op1.operator_element(lQ,rQ).t(), op2.operator_element(lQ,rQ));
	else
	  MatrixScaleAdd(scaleV*factor, op1.operator_element(lQ,rQ), op2.operator_element(lQ,rQ));
      }

}

void SparseMatrix::buildUsingCsf(const SpinBlock& b, vector< vector<Csf> >& ladders, std::vector< Csf >& s) 
{
  StateInfo stateinfo = b.get_stateInfo();
  built = true;
  allocate(stateinfo);

  for (int i=0; i < stateinfo.quanta.size(); i++)
    for (int j=0; j<stateinfo.quanta.size(); j++)
      if (allowedQuantaMatrix(i,j) )
	for (int jq =stateinfo.unBlockedIndex[j]; jq < stateinfo.unBlockedIndex[j]+stateinfo.quantaStates[j]; jq++) 
	{
	  for (int iq =stateinfo.unBlockedIndex[i]; iq < stateinfo.unBlockedIndex[i]+stateinfo.quantaStates[i]; iq++) 
	    operatorMatrix(i,j)(iq-stateinfo.unBlockedIndex[i]+1, jq-stateinfo.unBlockedIndex[j]+1) = redMatrixElement(s[iq], ladders[jq], &b);
	}

}

void finish_tensor_product( SpinBlock& b, SpinBlock* sysdot, 
                            const SparseMatrix& sysdot_op1, const SparseMatrix& sysdot_op2, SparseMatrix& op, 
                            bool include_parity, std::string& build_pattern )
{
  // Build and store new operator
  assert( ! op.get_built() );
  op.set_built() = true;
  op.set_build_pattern() = build_pattern;
  op.set_deltaQuantum(1, op.get_quantum_ladder().at( build_pattern ).at(2) );
  op.allocate(b.get_stateInfo());

  // Do tensor product
  double parity = 1.0;
  if ( include_parity ) parity = getCommuteParity( sysdot_op1.get_deltaQuantum(0), sysdot_op2.get_deltaQuantum(0), op.get_deltaQuantum(0) );
  SpinAdapted::operatorfunctions::TensorProduct(sysdot, sysdot_op1, sysdot_op2, &b, &(b.get_stateInfo()), op, parity);
}

void SpinBlock::store (bool forward, const vector<int>& sites, SpinBlock& b)
{
  Timer disktimer;
  std::string file;
  if (forward)
    file = str(boost::format("%s%s%d%s%d%s%d%s") % dmrginp.save_prefix() % "/SpinBlock-forward-"% sites[0] % "-" % sites[sites.size()-1] % "." % mpigetrank() % ".tmp" );
  else
    file = str(boost::format("%s%s%d%s%d%s%d%s") % dmrginp.save_prefix() % "/SpinBlock-backward-"% sites[0] % "-" % sites[sites.size()-1] % "." % mpigetrank() % ".tmp" );
  if (dmrginp.outputlevel() > 0) 
    pout << "\t\t\t Saving block file :: " << file << endl;


  std::ofstream ofs(file.c_str(), std::ios::binary);
  b.Save (ofs);
  ofs.close();
  //pout << "\t\t\t block save disk time " << disktimer.elapsedwalltime() << " " << disktimer.elapsedcputime() << endl;
}

void SparseMatrix::allocate(const SpinBlock& b)
{
  allocate(b.get_stateInfo());
}