21 #ifndef _libint2_src_bin_libint_hrr_h_ 22 #define _libint2_src_bin_libint_hrr_h_ 34 #include <prefactors.h> 35 #include <default_params.h> 56 template<
class IntType,
class BFSet,
int part, FunctionPosition loc_a,
57 unsigned int pos_a, FunctionPosition loc_b,
unsigned int pos_b>
64 typedef IntType TargetType;
65 typedef IntType ChildType;
76 static SafePtr<ThisType> Instance(
const SafePtr<TargetType>&,
unsigned int dir = 0);
81 if (loc_b == InBra && loc_a == InKet)
82 return BraketDirection::BraToKet;
83 else if (loc_b == InKet && loc_a == InBra)
84 return BraketDirection::KetToBra;
86 return BraketDirection::None;
93 if (boost::is_same<BasisFunctionType,CGShell>::value)
101 SafePtr<TargetType>
target()
const {
return target_;};
103 SafePtr<ChildType> child(
unsigned int i)
const;
105 SafePtr<DGVertex>
rr_target()
const {
return static_pointer_cast<DGVertex,TargetType>(target());}
107 SafePtr<DGVertex>
rr_child(
unsigned int i)
const {
return static_pointer_cast<DGVertex,ChildType>(child(i));}
113 std::string spfunction_call(
const SafePtr<CodeContext>& context,
114 const SafePtr<ImplicitDimensions>& dims)
const;
122 HRR(
const SafePtr<TargetType>&,
unsigned int dir);
125 SafePtr<TargetType> target_;
126 static const unsigned int max_nchildren_ = 8;
127 SafePtr<ChildType> children_[max_nchildren_];
128 unsigned int nchildren_;
130 void oper_checks()
const;
133 std::string generate_label()
const;
135 SafePtr<ImplicitDimensions> adapt_dims_(
const SafePtr<ImplicitDimensions>& dims)
const;
137 bool register_with_rrstack()
const;
142 bool expl_high_dim()
const;
143 bool expl_low_dim()
const;
146 template <
class IntType,
class F,
int part,
147 FunctionPosition loc_a,
unsigned int pos_a,
148 FunctionPosition loc_b,
unsigned int pos_b>
149 SafePtr< HRR<IntType,F,part,loc_a,pos_a,loc_b,pos_b> >
152 SafePtr<ThisType> this_ptr(
new ThisType(Tint,dir));
154 if (this_ptr->num_children() != 0) {
155 this_ptr->register_with_rrstack();
158 return SafePtr<ThisType>();
161 template <
class IntType,
class F,
int part,
162 FunctionPosition loc_a,
unsigned int pos_a,
163 FunctionPosition loc_b,
unsigned int pos_b>
165 dir_(dir), target_(Tint), nchildren_(0)
167 using namespace libint2::algebra;
168 using namespace libint2::prefactor;
171 assert(loc_a != loc_b);
174 const typename IntType::AuxQuantaType& aux = Tint->aux();
175 const typename IntType::OperType& oper = Tint->oper();
178 if (loc_a != loc_b && oper.hermitian(part) != +1) {
182 typedef typename IntType::BraType IBraType;
183 typedef typename IntType::KetType IKetType;
184 IBraType* bra =
new IBraType(Tint->bra());
185 IKetType* ket =
new IKetType(Tint->ket());
190 if (loc_a == InKet && loc_b == InBra) {
191 F a(ket->member(part,pos_a));
192 F b(bra->member(part,pos_b));
195 F bm1(b); bm1.dec(dir_);
201 bra->set_member(bm1,part,pos_b);
204 F ap1(a); ap1.inc(dir_);
205 ket->set_member(ap1,part,pos_a);
206 children_[nchildren_++] = IntType::Instance(*bra,*ket,aux,oper);
207 ket->set_member(a,part,pos_a);
210 children_[nchildren_++] = IntType::Instance(*bra,*ket,aux,oper);
216 for(
unsigned int xyz=0; xyz<3; ++xyz) {
218 if (a.deriv().d(xyz) > 0) {
220 a_der_m1.deriv().dec(xyz);
221 ket->set_member(a_der_m1,part,pos_a);
222 children_[nchildren_++] = IntType::Instance(*bra,*ket,aux,oper);
223 ket->set_member(a,part,pos_a);
226 if (bm1.deriv().d(xyz) > 0) {
228 bm1_der_m1.deriv().dec(xyz);
229 bra->set_member(bm1_der_m1,part,pos_b);
230 children_[nchildren_++] = IntType::Instance(*bra,*ket,aux,oper);
231 bra->set_member(bm1,part,pos_b);
239 expr_ = children_[0] + prefactors.
Y_X[part][dir] * children_[1];
242 const bool aderiv = a.deriv().d(dir_) > 0;
245 a_der_m1.deriv().dec(dir_);
246 ket->set_member(a_der_m1,part,pos_a);
247 children_[nchildren_++] = IntType::Instance(*bra,*ket,aux,oper);
248 ket->set_member(a,part,pos_a);
252 const bool bderiv = bm1.deriv().d(dir_) > 0;
255 bm1_der_m1.deriv().dec(dir_);
256 bra->set_member(bm1_der_m1,part,pos_b);
257 children_[nchildren_++] = IntType::Instance(*bra,*ket,aux,oper);
258 bra->set_member(bm1,part,pos_b);
262 expr_ += Vector(a.deriv())[dir_] * children_[2];
264 expr_ -= Vector(b.deriv())[dir_] * children_[aderiv ? 3 : 2];
268 if (loc_a == InBra && loc_b == InKet) {
269 F a(bra->member(part,pos_a));
270 F b(ket->member(part,pos_b));
273 F bm1(b); bm1.dec(dir_);
279 ket->set_member(bm1,part,pos_b);
282 F ap1(a); ap1.inc(dir_);
283 bra->set_member(ap1,part,pos_a);
284 children_[nchildren_++] = IntType::Instance(*bra,*ket,aux,oper);
285 bra->set_member(a,part,pos_a);
288 children_[nchildren_++] = IntType::Instance(*bra,*ket,aux,oper);
294 for(
unsigned int xyz=0; xyz<3; ++xyz) {
296 if (a.deriv().d(xyz) > 0) {
298 a_der_m1.deriv().dec(xyz);
299 bra->set_member(a_der_m1,part,pos_a);
300 children_[nchildren_++] = IntType::Instance(*bra,*ket,aux,oper);
301 bra->set_member(a,part,pos_a);
304 if (bm1.deriv().d(xyz) > 0) {
306 bm1_der_m1.deriv().dec(xyz);
307 ket->set_member(bm1_der_m1,part,pos_b);
308 children_[nchildren_++] = IntType::Instance(*bra,*ket,aux,oper);
309 ket->set_member(bm1,part,pos_b);
316 expr_ = children_[0] + prefactors.
X_Y[part][dir] * children_[1];
319 const bool aderiv = a.deriv().d(dir_) > 0;
322 a_der_m1.deriv().dec(dir_);
323 bra->set_member(a_der_m1,part,pos_a);
324 children_[nchildren_++] = IntType::Instance(*bra,*ket,aux,oper);
325 bra->set_member(a,part,pos_a);
329 const bool bderiv = bm1.deriv().d(dir_) > 0;
332 bm1_der_m1.deriv().dec(dir_);
333 ket->set_member(bm1_der_m1,part,pos_b);
334 children_[nchildren_++] = IntType::Instance(*bra,*ket,aux,oper);
335 ket->set_member(bm1,part,pos_b);
339 expr_ += Vector(a.deriv())[dir_] * children_[2];
341 expr_ -= Vector(b.deriv())[dir_] * children_[aderiv ? 3 : 2];
350 template <
class IntType,
class F,
int part,
351 FunctionPosition loc_a,
unsigned int pos_a,
352 FunctionPosition loc_b,
unsigned int pos_b>
354 HRR<IntType,F,part,loc_a,pos_a,loc_b,pos_b>::register_with_rrstack()
const 361 if (TrivialBFSet<F>::result)
363 typedef typename IntType::BraType IBraType;
364 typedef typename IntType::KetType IKetType;
365 const IBraType& bra = target_->bra();
366 const IKetType& ket = target_->ket();
369 bool nonzero_quanta =
false;
370 unsigned const int npart = IntType::OperatorType::Properties::np;
371 for(
unsigned int p=0; p<npart; p++) {
374 int nfbra = bra.num_members(p);
376 for(
int f=0; f<nfbra; f++)
377 if (!bra.member(p,f).zero() || !bra.member(p,f).deriv().zero())
378 nonzero_quanta =
true;
379 int nfket = ket.num_members(p);
381 for(
int f=0; f<nfket; f++)
382 if (!ket.member(p,f).zero() || !ket.member(p,f).deriv().zero())
383 nonzero_quanta =
true;
386 if (!nonzero_quanta) {
388 SafePtr<ThisType> this_ptr =
389 const_pointer_cast<ThisType,const ThisType>(
390 static_pointer_cast<const ThisType, const ParentType>(
391 EnableSafePtrFromThis<ParentType>::SafePtr_from_this()
394 rrstack->find(this_ptr);
403 IBraType bra_zero(bra);
404 IKetType ket_zero(ket);
405 for(
unsigned int p=0; p<npart; p++) {
408 int nfbra = bra_zero.num_members(p);
409 for(
int f=0; f<nfbra; f++) {
410 typedef typename IBraType::bfs_type bfs_type;
411 typedef typename IBraType::bfs_ref bfs_ref;
412 bfs_ref bfs = bra_zero.member(p,f);
413 if (!bfs.zero() || !bfs.deriv().zero()) {
418 int nfket = ket_zero.num_members(p);
419 for(
int f=0; f<nfket; f++) {
420 typedef typename IKetType::bfs_type bfs_type;
421 typedef typename IKetType::bfs_ref bfs_ref;
422 bfs_ref bfs = ket_zero.member(p,f);
423 if (!bfs.zero() || !bfs.deriv().zero()) {
431 typedef GenOper< GenMultSymmOper_Descr<IntType::OperatorType::Properties::np> > DummyOper;
432 typedef typename IBraType::bfs_type bfs_type;
434 typedef GenIntegralSet<DummyOper, IncableBFSet, IBraType, IKetType, DummyQuanta> DummyIntegral;
435 DummyOper dummy_oper;
436 DummyQuanta dummy_quanta(std::vector<int>(0,0));
437 SafePtr<DummyIntegral> dummy_integral = DummyIntegral::Instance(bra_zero,ket_zero,dummy_quanta,dummy_oper);
440 typedef HRR<DummyIntegral,F,part,loc_a,pos_a,loc_b,pos_b> DummyHRR;
441 SafePtr<DummyHRR> dummy_hrr = DummyHRR::Instance(dummy_integral,dir_);
443 rrstack->find(dummy_hrr);
447 template <
class IntType,
class F,
int part,
448 FunctionPosition loc_a,
unsigned int pos_a,
449 FunctionPosition loc_b,
unsigned int pos_b>
450 HRR<IntType,F,part,loc_a,pos_a,loc_b,pos_b>::~HRR()
455 template <
class IntType,
class F,
int part,
456 FunctionPosition loc_a,
unsigned int pos_a,
457 FunctionPosition loc_b,
unsigned int pos_b>
459 HRR<IntType,F,part,loc_a,pos_a,loc_b,pos_b>::oper_checks()
const 467 typedef typename IntType::OperatorType Oper;
468 if (part < 0 || part >= Oper::Properties::np) {
473 if (loc_a == loc_b && pos_a == pos_b) {
479 template <
class IntType,
class F,
int part,
480 FunctionPosition loc_a,
unsigned int pos_a,
481 FunctionPosition loc_b,
unsigned int pos_b>
482 SafePtr<typename HRR<IntType,F,part,loc_a,pos_a,loc_b,pos_b>::ChildType>
485 assert(i>=0 && i<nchildren_);
488 for(
unsigned int c=0; c<max_nchildren_; c++) {
498 template <
class IntType,
class F,
int part,
499 FunctionPosition loc_a,
unsigned int pos_a,
500 FunctionPosition loc_b,
unsigned int pos_b>
506 os <<
"HRR Part " << part <<
" " 507 << (loc_a == InBra ?
"bra" :
"ket") <<
" " << pos_a <<
" " 508 << (loc_b == InBra ?
"bra" :
"ket") <<
" " << pos_b <<
" ";
510 if (loc_a == InBra) {
511 F sh_a(target_->bra(part,pos_a));
515 os << sh_a.
label() <<
" ";
517 if (loc_b == InBra) {
518 F sh_b(target_->bra(part,pos_b));
523 F sh_b(target_->ket(part,pos_b));
529 F sh_a(target_->ket(part,pos_a));
531 os << sh_a.label() <<
" ";
533 if (loc_b == InBra) {
534 F sh_b(target_->bra(part,pos_b));
539 F sh_b(target_->ket(part,pos_b));
548 template <
class IntType,
class F,
int part,
549 FunctionPosition loc_a,
unsigned int pos_a,
550 FunctionPosition loc_b,
unsigned int pos_b>
553 const SafePtr<CodeContext>& context,
const SafePtr<ImplicitDimensions>& dims)
const 560 << context->value_to_pointer(
rr_target()->symbol());
563 for(
unsigned int c=0; c<nchildren; c++) {
564 os <<
", " << context->value_to_pointer(
rr_child(c)->symbol());
567 unsigned int hsr = 1;
572 for(
int p=0; p<part; p++) {
573 unsigned int nbra = target_->bra().num_members(p);
574 for(
unsigned int i=0; i<nbra; i++) {
575 SubIterator* iter = target_->bra().member_subiter(p,i);
579 unsigned int nket = target_->ket().num_members(p);
580 for(
unsigned int i=0; i<nket; i++) {
581 SubIterator* iter = target_->ket().member_subiter(p,i);
588 taskmgr.
current().params()->max_hrr_hsrank(hsr);
592 if (loc_a == loc_b && pos_a != 0 && pos_b != 0)
593 throw CodeDoesNotExist(
"HRR::spfunction_call -- has not been generalized yet");
596 unsigned int lsr = 1;
597 unsigned int np = IntType::OperType::Properties::np;
598 for(
unsigned int p=part+1; p<np; p++) {
599 unsigned int nbra = target_->bra().num_members(p);
600 for(
unsigned int i=0; i<nbra; i++) {
601 SubIterator* iter = target_->bra().member_subiter(p,i);
605 unsigned int nket = target_->ket().num_members(p);
606 for(
unsigned int i=0; i<nket; i++) {
607 SubIterator* iter = target_->ket().member_subiter(p,i);
613 taskmgr.
current().params()->max_hrr_hsrank(hsr);
619 os <<
")" << context->end_of_stat() << endl;
623 template <
class IntType,
class F,
int part,
624 FunctionPosition loc_a,
unsigned int pos_a,
625 FunctionPosition loc_b,
unsigned int pos_b>
635 template <
class IntType,
class F,
int part,
636 FunctionPosition loc_a,
unsigned int pos_a,
637 FunctionPosition loc_b,
unsigned int pos_b>
639 HRR<IntType,F,part,loc_a,pos_a,loc_b,pos_b>::expl_low_dim()
const 641 unsigned int np = IntType::OperType::Properties::np;
652 template <
class IntType,
class F,
int part,
653 FunctionPosition loc_a,
unsigned int pos_a,
654 FunctionPosition loc_b,
unsigned int pos_b>
655 SafePtr<ImplicitDimensions>
656 HRR<IntType,F,part,loc_a,pos_a,loc_b,pos_b>::adapt_dims_(
const SafePtr<ImplicitDimensions>& dims)
const 658 bool high_rank = expl_high_dim();
659 bool low_rank = expl_low_dim();
661 SafePtr<Entity> high_dim, low_dim;
666 high_dim = dims->high();
672 low_dim = dims->low();
675 SafePtr<ImplicitDimensions> localdims(
new ImplicitDimensions(high_dim,low_dim,dims->vecdim()));
BraketDirection braket_direction() const
overrides RecurrenceRelation::braket_direction()
Definition: hrr.h:80
std::string label_to_funcname(const std::string &label)
Converts a label, e.g. name of the target node, to the name of the function to compute it.
Definition: default_params.cc:216
Manages tasks. This is a Singleton.
Definition: task.h:63
TrivialBFSet<T> defines static member result, which is true if T is a basis function set consisting o...
Definition: bfset.h:892
SafePtr< rdouble > Y_X[np][3]
Cartesian components of Y_X vectors.
Definition: prefactors.h:66
Defaults definitions for various parameters assumed by Libint.
Definition: algebra.cc:24
SafePtr< ChildType > child(unsigned int i) const
child(i) returns pointer to i-th child
Definition: hrr.h:483
std::string spfunction_call(const SafePtr< CodeContext > &context, const SafePtr< ImplicitDimensions > &dims) const
Implementation of RecurrenceRelation::spfunction_call()
Definition: hrr.h:552
void current(const std::string &task_label)
Makes this task current (must have been added already)
Definition: task.cc:76
RecurrenceRelation::ExprType ExprType
A short alias.
Definition: hrr.h:67
bool is_simple() const
Implementation of RecurrenceRelation::is_simple()
Definition: hrr.h:109
AlgebraicOperator is an algebraic operator that acts on objects of type T.
Definition: algebra.h:48
const std::string & label() const
label() returns a unique, short, descriptive label of this RR (e.g.
Definition: rr.cc:292
bool exists(const IncableBFSet &A)
Return true if A is valid.
Definition: bfset.h:92
SafePtr< DGVertex > rr_target() const
Implementation of RecurrenceRelation::target()
Definition: hrr.h:105
Iterator provides a base class for all object iterator classes.
Definition: iter.h:45
virtual unsigned int num_iter() const =0
Returns a number of iterations (number of elements in a set over which to iterate).
unsigned int num_children() const
Implementation of RecurrenceRelation::num_children()
Definition: hrr.h:99
RecurrenceRelation describes all recurrence relations.
Definition: rr.h:101
SafePtr< DGVertex > rr_child(unsigned int i) const
Implementation of RecurrenceRelation::child()
Definition: hrr.h:107
Set of basis functions.
Definition: bfset.h:42
static LibraryTaskManager & Instance()
LibraryTaskManager is a Singleton.
Definition: task.cc:34
static bool directional()
is this recurrence relation parameterized by a direction.
Definition: hrr.h:92
A generic Horizontal Recurrence Relation:
Definition: hrr.h:58
SafePtr< TargetType > target() const
returns pointer to the target
Definition: hrr.h:101
SafePtr< rdouble > X_Y[np][3]
Cartesian components of X-Y vectors.
Definition: prefactors.h:58
static SafePtr< RRStackBase< RR > > & Instance()
Obtain the unique Instance of RRStack.
Definition: rr.h:71
DefaultQuantumNumbers< int, 0 >::Result EmptySet
EmptySet is the type that describes null set of auxiliary indices.
Definition: quanta.h:407
This exception used to indicate that some code hasn't been developed or generalized yet.
Definition: exception.h:87