Commit 79a5e447 authored by Alvarez, Gonzalo's avatar Alvarez, Gonzalo
Browse files

ChemicalH: working on creation API

parent 3d2bfc1f
/*
Copyright (c) 2009-2012-2019, UT-Battelle, LLC
Copyright (c) 2009-2012-2019-2020, UT-Battelle, LLC
All rights reserved
[DMRG++, Version 5.]
......@@ -172,14 +172,15 @@ public:
// set this basis to the outer product of
// basis2 and basis3 or basis3 and basis2 depending on dir
template<typename SomeSuperOperatorHelperType>
void setToProduct(const ThisType& basis2,
const ThisType& basis3,
ProgramGlobals::DirectionEnum dir)
const SomeSuperOperatorHelperType& someSuperOpHelper)
{
if (dir == ProgramGlobals::DirectionEnum::EXPAND_SYSTEM)
setToProduct(basis2, basis3);
if (someSuperOpHelper.dir() == ProgramGlobals::DirectionEnum::EXPAND_SYSTEM)
setToProductInternal(basis2, basis3, someSuperOpHelper);
else
setToProduct(basis3, basis2);
setToProductInternal(basis3, basis2, someSuperOpHelper);
}
//! transform this basis by transform
......@@ -306,8 +307,10 @@ private:
//! set this basis to the outer product of basis2 and basis3
//!PTEX_LABEL{setToProductOps}
void setToProduct(const ThisType& basis2,
const ThisType& basis3)
template<typename SomeSuperOperatorHelperType>
void setToProductInternal(const ThisType& basis2,
const ThisType& basis3,
const SomeSuperOperatorHelperType& someSuperOpHelper)
{
// reorder the basis
BasisType::setToProduct(basis2, basis3);
......@@ -317,7 +320,8 @@ private:
basis2.operators_,
basis3,
basis3.operators_,
BaseType::permutationInverse());
BaseType::permutationInverse(),
someSuperOpHelper);
//! Calc. hamiltonian
......
/*
Copyright (c) 2009-2014, UT-Battelle, LLC
Copyright (c) 2009-2014-2020, UT-Battelle, LLC
All rights reserved
[DMRG++, Version 5.]
......@@ -82,6 +82,7 @@ DISCLOSED WOULD NOT INFRINGE PRIVATELY OWNED RIGHTS.
#include "ProgressIndicator.h"
#include "KroneckerDumper.h"
#include "Io/IoNg.h"
#include "SuperOperatorHelper.h"
namespace Dmrg {
......@@ -364,7 +365,8 @@ private:
BasisWithOperatorsType Xbasis("Xbasis");
typedef LeftRightSuper<BasisWithOperatorsType, BasisType> LeftRightSuper2Type;
Xbasis.setOneSite(X, model, time);
leftOrRight.setToProduct(pS, Xbasis, dir);
SuperOperatorHelper<SomeModelType> superOperatorHelper(model, pS.block(), X, dir);
leftOrRight.setToProduct(pS, Xbasis, superOperatorHelper);
SparseMatrixType matrix = leftOrRight.hamiltonian().getCRS();
......
......@@ -355,6 +355,21 @@ public:
permutationFull);
}
void outerProduct(const Operator& A,
const Operator& B,
const VectorRealType& signs,
bool order,
const VectorSizeType& permutationFull)
{
externalProduct2(data_,
A.getStorage(),
B.getStorage(),
signs,
order,
permutationFull);
}
SizeType metaDiff(const Operator& op2) const
{
const Operator& op1 = *this;
......
......@@ -249,6 +249,24 @@ public:
throw PsimagLite::RuntimeError("OperatorStorage: externalProduct\n");
}
friend void externalProduct2(OperatorStorage& C,
const OperatorStorage& A,
const OperatorStorage& B,
const VectorRealType& signs,
bool order,
const VectorSizeType& permutationFull)
{
if (A.justCRS() && B.justCRS() && C.justCRS())
return externalProduct(C.crs_,
A.getCRS(),
B.getCRS(),
signs,
order,
permutationFull);
throw PsimagLite::RuntimeError("OperatorStorage: externalProduct\n");
}
friend void fullMatrixToCrsMatrix(OperatorStorage& dest,
const PsimagLite::Matrix<ComplexOrRealType>& src)
{
......
......@@ -302,14 +302,16 @@ public:
ChangeOfBasisType::changeBasis(hamiltonian_, ftransform);
}
template<typename SomeSuperOperatorHelperType>
void setToProduct(const BasisType& basis1,
const ThisType& ops1,
const BasisType& basis2,
const ThisType& ops2,
const VectorSizeType& permutationInverse)
const VectorSizeType& permutationInverse,
const SomeSuperOperatorHelperType& someSuperOpHelper)
{
setToProductLocal(basis1, ops1, basis2, ops2, permutationInverse);
// ChemicalH: Set super ops here
setToProductSuper(basis1, ops1, basis2, ops2, permutationInverse, someSuperOpHelper);
}
void outerProductHamiltonian(const StorageType& h2,
......@@ -378,6 +380,12 @@ public:
hamiltonian_.clear();
}
const OperatorType& getSuperByIndex(SizeType ind) const
{
assert(ind < superOps_.size());
return superOps_[ind];
}
SizeType superIndices(const VectorSizeType&, SizeType) const
{
PsimagLite::String msg(__FILE__);
......@@ -439,6 +447,42 @@ private:
}
}
template<typename SomeSuperOperatorHelperType>
void setToProductSuper(const BasisType& basis2,
const ThisType& ops2,
const BasisType& basis3,
const ThisType& ops3,
const VectorSizeType& permutationInverse,
const SomeSuperOperatorHelperType& someSuperOpHelper)
{
typename PsimagLite::Vector<RealType>::Type fermionicSigns;
SizeType nSuperOps = someSuperOpHelper.size();
superOps_.resize(nSuperOps);
ProgramGlobals::FermionOrBosonEnum savedSign = ProgramGlobals::FermionOrBosonEnum::BOSON;
typedef typename SomeSuperOperatorHelperType::PairBoolSizeType PairBoolSizeType;
const bool option = (basis3.block().size() == 1);
for (SizeType i = 0; i < nSuperOps; ++i) {
const PairBoolSizeType op2Index = someSuperOpHelper.leftOperatorIndex(i);
const PairBoolSizeType op3Index = someSuperOpHelper.rightOperatorIndex(i);
const OperatorType& op1 = (!op2Index.first) ? ops2.getLocalByIndex(op2Index.second)
: ops2.getSuperByIndex(op2Index.
second);
const OperatorType& op3 = (!op3Index.first) ? ops3.getLocalByIndex(op3Index.second)
: ops3.getSuperByIndex(op3Index.
second);
bool isFermion = (op3.fermionOrBoson() == ProgramGlobals::FermionOrBosonEnum::FERMION);
if (savedSign != op3.fermionOrBoson() || fermionicSigns.size() == 0) {
utils::fillFermionicSigns(fermionicSigns, basis2.signs(), (isFermion) ? -1 : 1);
savedSign = op3.fermionOrBoson();
}
superOps_[i].outerProduct(op1, op3, fermionicSigns, option, permutationInverse);
}
}
/* PSIDOC OperatorsExternalProduct
I will know explain how the full outer product between two operators
is implemented. If local operator $A$ lives in Hilbert space
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment