C++ SpinBlock::addAdditionalCompOps方法代码示例

//.........这里部分代码省略.........

      if (nactiveSites > environmentSites.size()) {
        nactiveSites = environmentSites.size();
      }
      ncoreSites = environmentSites.size() - nactiveSites;

      // figure out what sites are in the active and core sites
      int environmentActiveEnd = forward ? environmentStart + nactiveSites - 1 : environmentStart - nactiveSites + 1;
      int environmentCoreStart = forward ? environmentActiveEnd + 1 : environmentActiveEnd - 1;
      
      std::vector<int> activeSites(nactiveSites), coreSites(ncoreSites);
      for (int i = 0; i < nactiveSites; ++i) {
        activeSites[i] = min(environmentStart,environmentActiveEnd) + i;
      }
      for (int i = 0; i < ncoreSites; ++i) {
        coreSites[i] = min(environmentCoreStart,environmentEnd) + i;
      }

      SpinBlock environmentActive, environmentCore;
      environmentActive.nonactive_orb() = system.nonactive_orb();
      environmentCore.nonactive_orb() = system.nonactive_orb();
      if (coreSites.size() > 0) {
	environmentActive.set_integralIndex() = integralIndex;
	environmentCore.set_integralIndex() = integralIndex;
        environmentActive.default_op_components(!forward, leftState==rightState);
        environmentActive.setstoragetype(DISTRIBUTED_STORAGE);
        environmentCore.default_op_components(!forward, leftState==rightState);      
        environmentCore.setstoragetype(DISTRIBUTED_STORAGE);

        environmentActive.BuildTensorProductBlock(activeSites);
        environmentCore.BuildSingleSlaterBlock(coreSites);

        dmrginp.datatransfer -> start();
        environmentCore.addAdditionalCompOps();
        environmentActive.addAdditionalCompOps();
        dmrginp.datatransfer -> stop();

        if ((!dot_with_sys && onedot) || !onedot) {
	  environment.set_integralIndex() = integralIndex;
          environment.default_op_components(!forward, leftState == rightState);
          environment.setstoragetype(DISTRIBUTED_STORAGE);
          environment.BuildSumBlock(constraint, environmentCore, environmentActive,braquanta,ketquanta);
        } else {
	  newEnvironment.set_integralIndex() = integralIndex;
          newEnvironment.default_op_components(direct, environmentCore, environmentActive, haveNormops, haveCompops, leftState == rightState);
          newEnvironment.setstoragetype(DISTRIBUTED_STORAGE);
          newEnvironment.BuildSumBlock(constraint, environmentCore, environmentActive,braquanta,ketquanta);
          if (dmrginp.outputlevel() > 0) {
	    pout << "\t\t\t NewEnvironment block " << endl << newEnvironment << endl;
	    newEnvironment.printOperatorSummary();
          }
        }
      } else { // no core
        if ((!dot_with_sys && onedot) || !onedot) {
	  environment.set_integralIndex() = integralIndex;
          environment.default_op_components(!forward, leftState==rightState);
          environment.setstoragetype(DISTRIBUTED_STORAGE);
          environment.BuildTensorProductBlock(environmentSites); // exact block
        } else {
	  newEnvironment.set_integralIndex() = integralIndex;
          newEnvironment.default_op_components(!forward, leftState==rightState);
          newEnvironment.setstoragetype(DISTRIBUTED_STORAGE);
          newEnvironment.BuildTensorProductBlock(environmentSites);
        }
      }
    } else { //used for warmup guess environemnt

void SweepOnepdm::BlockAndDecimate (SweepParams &sweepParams, SpinBlock& system, SpinBlock& newSystem, const bool &useSlater, const bool& dot_with_sys)
{
  //mcheck("at the start of block and decimate");
  // figure out if we are going forward or backwards
  dmrginp.guessgenT -> start();
  bool forward = (system.get_sites() [0] == 0);
  SpinBlock systemDot;
  SpinBlock envDot;
  int systemDotStart, systemDotEnd;
  int systemDotSize = sweepParams.get_sys_add() - 1;
  if (forward)
  {
    systemDotStart = *system.get_sites().rbegin () + 1;
    systemDotEnd = systemDotStart + systemDotSize;
  }
  else
  {
    systemDotStart = system.get_sites() [0] - 1;
    systemDotEnd = systemDotStart - systemDotSize;
  }
  vector<int> spindotsites(2); 
  spindotsites[0] = systemDotStart;
  spindotsites[1] = systemDotEnd;
  systemDot = SpinBlock(systemDotStart, systemDotEnd);

  SpinBlock environment, environmentDot, newEnvironment;

  int environmentDotStart, environmentDotEnd, environmentStart, environmentEnd;

  const int nexact = forward ? sweepParams.get_forward_starting_size() : sweepParams.get_backward_starting_size();

  system.addAdditionalCompOps();
  InitBlocks::InitNewSystemBlock(system, systemDot, newSystem, sweepParams.get_sys_add(), dmrginp.direct(), DISTRIBUTED_STORAGE, true, true);
  
  InitBlocks::InitNewEnvironmentBlock(environment, systemDot, newEnvironment, system, systemDot,
				      sweepParams.get_sys_add(), sweepParams.get_env_add(), forward, dmrginp.direct(),
				      sweepParams.get_onedot(), nexact, useSlater, true, true, true);
  SpinBlock big;
  newSystem.set_loopblock(true);
  system.set_loopblock(false);
  newEnvironment.set_loopblock(false);
  InitBlocks::InitBigBlock(newSystem, newEnvironment, big); 

  const int nroots = dmrginp.nroots();
  std::vector<Wavefunction> solutions(nroots);

  
  for(int i=0;i<nroots;++i)
    {
      StateInfo newInfo;
      solutions[i].LoadWavefunctionInfo (newInfo, newSystem.get_sites(), i);
    }
  
#ifndef SERIAL
  mpi::communicator world;
  mpi::broadcast(world,solutions,0);
#endif

#ifdef SERIAL
  const int numprocs = 1;
#endif
#ifndef SERIAL
  const int numprocs = world.size();
#endif


  compute_onepdm(solutions, system, systemDot, newSystem, newEnvironment, big, numprocs);

}

void SpinAdapted::mps_nevpt::type1::Startup(const SweepParams &sweepParams, const bool &forward, perturber& pb, int baseState) {

#ifndef SERIAL
  mpi::communicator world;
#endif
  assert(forward);
  SpinBlock system;
  system.nonactive_orb() =pb.orb();
  bool restart=false, warmUp = false;
  int forward_starting_size=1, backward_starting_size=0, restartSize =0;
  InitBlocks::InitStartingBlock(system, forward, pb.wavenumber(), baseState, forward_starting_size, backward_starting_size, restartSize, restart, warmUp, 0,pb.braquanta, pb.ketquanta); 

  SpinBlock::store (forward, system.get_sites(), system, pb.wavenumber(), baseState); // if restart, just restoring an existing block --

  for (int i=0; i<mps_nevpt::sweepIters; i++) {
    SpinBlock newSystem;
    SpinBlock dotSystem(i+1,i+1,pb.orb(),false);

    system.addAdditionalCompOps();
    //newSystem.default_op_components(true, system, dotSystem, true, true, false);
    newSystem.perturb_op_components(false, system, dotSystem, pb);
    newSystem.setstoragetype(DISTRIBUTED_STORAGE);
    newSystem.BuildSumBlock(LessThanQ, system, dotSystem, pb.braquanta, pb.ketquanta);
    newSystem.printOperatorSummary();
    //SpinBlock Environment, big;
    //SpinBlock::restore (!forward, newSystem.get_complementary_sites() , Environment, baseState, baseState);
    //TODO
    //SpinBlock::restore (!forward, newSystem.get_complementary_sites() , Environment,sweepParams.current_root(),sweepParams.current_root());

    //big.BuildSumBlock(PARTICLE_SPIN_NUMBER_CONSTRAINT, newSystem, Environment, pb.braquanta, pb.ketquanta);

    //StateInfo envStateInfo;
    StateInfo ketStateInfo;
    StateInfo braStateInfo;
    StateInfo halfbraStateInfo;// It has the same left and right StateInfo as braStateInfo. However, its total quanta is pb.ketquanta.
    // It is used to project solution into to braStateInfo.

    std::vector<Wavefunction> solution; solution.resize(1);
    std::vector<Wavefunction> outputState; outputState.resize(1);
    std::vector<Wavefunction> solutionprojector; solutionprojector.resize(1);
    solution[0].LoadWavefunctionInfo(ketStateInfo, newSystem.get_sites(), baseState);
    #ifndef SERIAL
      broadcast(world, ketStateInfo, 0);
      broadcast(world, solution, 0);
    #endif
    outputState[0].AllowQuantaFor(newSystem.get_braStateInfo(), *(ketStateInfo.rightStateInfo), pb.braquanta);
    outputState[0].set_onedot(solution[0].get_onedot());
    outputState[0].Clear();
    solutionprojector[0].AllowQuantaFor(newSystem.get_braStateInfo(), *(ketStateInfo.rightStateInfo), pb.ketquanta);
    solutionprojector[0].set_onedot(solution[0].get_onedot());
    solutionprojector[0].Clear();
    //TensorProduct (newSystem.get_braStateInfo(), *(ketStateInfo.rightStateInfo), pb.braquanta[0], EqualQ, braStateInfo);
    //TODO
    //TensorProduct do not support const StateInfo&
    TensorProduct (newSystem.set_braStateInfo(), *(ketStateInfo.rightStateInfo), pb.braquanta[0], EqualQ, braStateInfo);
    TensorProduct (newSystem.set_braStateInfo(), *(ketStateInfo.rightStateInfo), pb.ketquanta[0], EqualQ, halfbraStateInfo);

    //StateInfo::restore(forward, environmentsites, envStateInfo, baseState);

    //DiagonalMatrix e;
    //if(i == 0)
    //  GuessWave::guess_wavefunctions(solution, e, big, TRANSPOSE, true, true, 0.0, baseState); 
    //else
    //  GuessWave::guess_wavefunctions(solution, e, big, TRANSFORM, true, true, 0.0, baseState); 


    //SpinAdapted::operatorfunctions::Product(&newSystem, ccd, solution[0], &ketStateInfo, stateb.getw(), temp, SpinQuantum(0, SpinSpace(0), IrrepSpace(0)), true, 1.0);

    

    boost::shared_ptr<SparseMatrix> O;
    if (pb.type() == TwoPerturbType::Va)
      O = newSystem.get_op_array(CDD_SUM).get_local_element(0)[0]->getworkingrepresentation(&newSystem);
    if (pb.type() == TwoPerturbType::Vi)
      O = newSystem.get_op_array(CCD_SUM).get_local_element(0)[0]->getworkingrepresentation(&newSystem);
    boost::shared_ptr<SparseMatrix> overlap = newSystem.get_op_array(OVERLAP).get_local_element(0)[0]->getworkingrepresentation(&newSystem);
    SpinAdapted::operatorfunctions::TensorMultiply(*O, &braStateInfo, &ketStateInfo , solution[0], outputState[0], pb.delta, true, 1.0);
    SpinAdapted::operatorfunctions::TensorMultiply(*overlap, &halfbraStateInfo, &ketStateInfo , solution[0], solutionprojector[0], overlap->get_deltaQuantum(0), true, 1.0);
    DensityMatrix bratracedMatrix(newSystem.get_braStateInfo());
    bratracedMatrix.allocate(newSystem.get_braStateInfo());
    double norm = DotProduct(outputState[0], outputState[0]);
    if(norm > NUMERICAL_ZERO)
      SpinAdapted::operatorfunctions::MultiplyProduct(outputState[0], Transpose(const_cast<Wavefunction&> (outputState[0])), bratracedMatrix, 0.5/norm);
    SpinAdapted::operatorfunctions::MultiplyProduct(solutionprojector[0], Transpose(const_cast<Wavefunction&> (solutionprojector[0])), bratracedMatrix, 0.5);
    std::vector<Matrix> brarotateMatrix, ketrotateMatrix;
    LoadRotationMatrix (newSystem.get_sites(), ketrotateMatrix, baseState);
    double error;
    if (!mpigetrank())
      error = makeRotateMatrix(bratracedMatrix, brarotateMatrix, sweepParams.get_keep_states(), sweepParams.get_keep_qstates());
    #ifndef SERIAL
      broadcast(world, ketrotateMatrix, 0);
      broadcast(world, brarotateMatrix, 0);
    #endif

    SaveRotationMatrix (newSystem.get_sites(), brarotateMatrix, pb.wavenumber());
    newSystem.transform_operators(brarotateMatrix,ketrotateMatrix);
    SpinBlock::store (forward, newSystem.get_sites(), newSystem, pb.wavenumber(), baseState); // if restart, just restoring an existing block --
    system=newSystem;
  }
  //TODO
//.........这里部分代码省略.........

void SweepGenblock::BlockAndDecimate (SweepParams &sweepParams, SpinBlock& system, SpinBlock& newSystem, const bool &useSlater, const bool& dot_with_sys, int stateA, int stateB)
{
  if (dmrginp.outputlevel() > 0) 
    mcheck("at the start of block and decimate");
  p1out << "\t\t\t Performing Blocking"<<endl;
  dmrginp.guessgenT -> start();
  // figure out if we are going forward or backwards  
  bool forward = (system.get_sites() [0] == 0);
  SpinBlock systemDot;
  int systemDotStart, systemDotEnd;
  int systemDotSize = sweepParams.get_sys_add() - 1;
  if (forward)
  {
    systemDotStart = dmrginp.spinAdapted() ? *system.get_sites().rbegin () + 1 : (*system.get_sites().rbegin ())/2 + 1 ;
    systemDotEnd = systemDotStart + systemDotSize;
  }
  else
  {
    systemDotStart = dmrginp.spinAdapted() ? system.get_sites()[0] - 1 : (system.get_sites()[0])/2 - 1 ;
    systemDotEnd = systemDotStart - systemDotSize;
  }
  vector<int> spindotsites(2); 
  spindotsites[0] = systemDotStart;
  spindotsites[1] = systemDotEnd;
  dmrginp.sysdotmake->start();
  systemDot = SpinBlock(systemDotStart, systemDotEnd, system.get_integralIndex(), stateA==stateB);
  dmrginp.sysdotmake->stop();

  const int nexact = forward ? sweepParams.get_forward_starting_size() : sweepParams.get_backward_starting_size();

  dmrginp.guessgenT -> stop();
  dmrginp.datatransfer -> start();
  system.addAdditionalCompOps();
  dmrginp.datatransfer -> stop();
  dmrginp.initnewsystem->start();
  InitBlocks::InitNewSystemBlock(system, systemDot, newSystem, stateA, stateB, sweepParams.get_sys_add(), dmrginp.direct(), system.get_integralIndex(), DISTRIBUTED_STORAGE, dot_with_sys, true);
  dmrginp.initnewsystem->stop();

  pout << "\t\t\t System  Block"<<newSystem;
  newSystem.printOperatorSummary();

  std::vector<Matrix> leftrotateMatrix, rightrotateMatrix;

  LoadRotationMatrix (newSystem.get_sites(), leftrotateMatrix, stateA);
  LoadRotationMatrix (newSystem.get_sites(), rightrotateMatrix, stateB);

#ifndef SERIAL
  mpi::communicator world;
  broadcast(world, leftrotateMatrix, 0);
  broadcast(world, rightrotateMatrix, 0);
#endif

  p1out <<"\t\t\t Performing Renormalization "<<endl<<endl;

  dmrginp.operrotT->start();
  if (stateB == stateA)
    newSystem.transform_operators(leftrotateMatrix);
  else
    newSystem.transform_operators(leftrotateMatrix, rightrotateMatrix);
  dmrginp.operrotT->stop();


  if (dmrginp.outputlevel() > 0) 
    //mcheck("after rotation and transformation of block");
  p2out <<newSystem<<endl;
  newSystem.printOperatorSummary();
  //mcheck("After renorm transform");

  p2out << *dmrginp.guessgenT<<" "<<*dmrginp.multiplierT<<" "<<*dmrginp.operrotT<< "  "<<globaltimer.totalwalltime()<<" timer "<<endl;
  p2out << *dmrginp.makeopsT<<"  "<<*dmrginp.initnewsystem<<"  "<<*dmrginp.sysdotmake<<"  "<<*dmrginp.buildcsfops<<" makeops "<<endl;
  p2out << *dmrginp.datatransfer<<" datatransfer "<<endl;
  p2out <<"oneindexopmult   twoindexopmult   Hc  couplingcoeff"<<endl;  
  p2out << *dmrginp.oneelecT<<" "<<*dmrginp.twoelecT<<" "<<*dmrginp.hmultiply<<" "<<*dmrginp.couplingcoeff<<" hmult"<<endl;
  p2out << *dmrginp.buildsumblock<<" "<<*dmrginp.buildblockops<<" build block"<<endl;
  p2out << *dmrginp.blockintegrals<<"  "<<*dmrginp.blocksites<<"  "<<*dmrginp.statetensorproduct<<"  "<<*dmrginp.statecollectquanta<<"  "<<*dmrginp.buildsumblock<<" "<<*dmrginp.builditeratorsT<<"  "<<*dmrginp.diskio<<" build sum block"<<endl;
  p2out << "addnoise  S_0_opxop  S_1_opxop   S_2_opxop"<<endl;
  p3out << *dmrginp.addnoise<<" "<<*dmrginp.s0time<<" "<<*dmrginp.s1time<<" "<<*dmrginp.s2time<<endl;

}

void SweepGenblock::BlockAndDecimate (SweepParams &sweepParams, SpinBlock& system, SpinBlock& newSystem, const bool &useSlater, const bool& dot_with_sys, int stateA, int stateB)
{
  if (dmrginp.outputlevel() > 0) 
    mcheck("at the start of block and decimate");
  pout << "\t\t\t Performing Blocking"<<endl;
  dmrginp.guessgenT -> start();
  // figure out if we are going forward or backwards  
  bool forward = (system.get_sites() [0] == 0);
  SpinBlock systemDot;
  int systemDotStart, systemDotEnd;
  int systemDotSize = sweepParams.get_sys_add() - 1;
  if (forward)
  {
    systemDotStart = dmrginp.spinAdapted() ? *system.get_sites().rbegin () + 1 : (*system.get_sites().rbegin ())/2 + 1 ;
    systemDotEnd = systemDotStart + systemDotSize;
  }
  else
  {
    systemDotStart = dmrginp.spinAdapted() ? system.get_sites()[0] - 1 : (system.get_sites()[0])/2 - 1 ;
    systemDotEnd = systemDotStart - systemDotSize;
  }
  vector<int> spindotsites(2); 
  spindotsites[0] = systemDotStart;
  spindotsites[1] = systemDotEnd;
  systemDot = SpinBlock(systemDotStart, systemDotEnd, stateA==stateB);

  const int nexact = forward ? sweepParams.get_forward_starting_size() : sweepParams.get_backward_starting_size();

  system.addAdditionalCompOps();
  InitBlocks::InitNewSystemBlock(system, systemDot, newSystem, stateA, stateB, sweepParams.get_sys_add(), dmrginp.direct(), DISTRIBUTED_STORAGE, dot_with_sys, true);

  pout << "\t\t\t System  Block"<<newSystem;
  if (dmrginp.outputlevel() > 0)
    newSystem.printOperatorSummary();

  std::vector<Matrix> leftrotateMatrix, rightrotateMatrix;

  LoadRotationMatrix (newSystem.get_sites(), leftrotateMatrix, stateA);
  LoadRotationMatrix (newSystem.get_sites(), rightrotateMatrix, stateB);

#ifndef SERIAL
  mpi::communicator world;
  broadcast(world, leftrotateMatrix, 0);
  broadcast(world, rightrotateMatrix, 0);
#endif

  pout <<"\t\t\t Performing Renormalization "<<endl<<endl;

  if (stateB == stateA)
    newSystem.transform_operators(leftrotateMatrix);
  else
    newSystem.transform_operators(leftrotateMatrix, rightrotateMatrix);


  if (dmrginp.outputlevel() > 0) 
    //mcheck("after rotation and transformation of block");
  if (dmrginp.outputlevel() > 0) 
    pout <<newSystem<<endl;
  if (dmrginp.outputlevel() > 0)
    newSystem.printOperatorSummary();
  //mcheck("After renorm transform");
}