20 #ifndef _libint2_src_bin_libint_hrr_h_ 21 #define _libint2_src_bin_libint_hrr_h_ 33 #include <prefactors.h> 34 #include <default_params.h> 55 template<
class IntType,
class BFSet,
int part, FunctionPosition loc_a,
56 unsigned int pos_a, FunctionPosition loc_b,
unsigned int pos_b>
63 typedef IntType TargetType;
64 typedef IntType ChildType;
75 static SafePtr<ThisType> Instance(
const SafePtr<TargetType>&,
unsigned int dir = 0);
82 if (boost::is_same<BasisFunctionType,CGShell>::value)
90 SafePtr<TargetType>
target()
const {
return target_;};
92 SafePtr<ChildType> child(
unsigned int i)
const;
94 SafePtr<DGVertex>
rr_target()
const {
return static_pointer_cast<
DGVertex,TargetType>(target());}
96 SafePtr<DGVertex>
rr_child(
unsigned int i)
const {
return static_pointer_cast<
DGVertex,ChildType>(child(i));}
102 std::string spfunction_call(
const SafePtr<CodeContext>& context,
103 const SafePtr<ImplicitDimensions>& dims)
const;
111 HRR(
const SafePtr<TargetType>&,
unsigned int dir);
114 SafePtr<TargetType> target_;
115 static const unsigned int max_nchildren_ = 8;
116 SafePtr<ChildType> children_[max_nchildren_];
117 unsigned int nchildren_;
119 void oper_checks()
const;
122 std::string generate_label()
const;
124 SafePtr<ImplicitDimensions> adapt_dims_(
const SafePtr<ImplicitDimensions>& dims)
const;
126 bool register_with_rrstack()
const;
131 bool expl_high_dim()
const;
132 bool expl_low_dim()
const;
135 template <
class IntType,
class F,
int part,
136 FunctionPosition loc_a,
unsigned int pos_a,
137 FunctionPosition loc_b,
unsigned int pos_b>
138 SafePtr< HRR<IntType,F,part,loc_a,pos_a,loc_b,pos_b> >
141 SafePtr<ThisType> this_ptr(
new ThisType(Tint,dir));
143 if (this_ptr->num_children() != 0) {
144 this_ptr->register_with_rrstack();
147 return SafePtr<ThisType>();
150 template <
class IntType,
class F,
int part,
151 FunctionPosition loc_a,
unsigned int pos_a,
152 FunctionPosition loc_b,
unsigned int pos_b>
154 dir_(dir), target_(Tint), nchildren_(0)
160 assert(loc_a != loc_b);
163 const typename IntType::AuxQuantaType& aux = Tint->aux();
164 const typename IntType::OperType& oper = Tint->oper();
167 if (loc_a != loc_b && oper.hermitian(part) != +1) {
171 typedef typename IntType::BraType IBraType;
172 typedef typename IntType::KetType IKetType;
173 IBraType* bra =
new IBraType(Tint->bra());
174 IKetType* ket =
new IKetType(Tint->ket());
179 if (loc_a == InKet && loc_b == InBra) {
180 F a(ket->member(part,pos_a));
181 F b(bra->member(part,pos_b));
184 F bm1(b); bm1.dec(dir_);
190 bra->set_member(bm1,part,pos_b);
193 F ap1(a); ap1.inc(dir_);
194 ket->set_member(ap1,part,pos_a);
195 children_[nchildren_++] = IntType::Instance(*bra,*ket,aux,oper);
196 ket->set_member(a,part,pos_a);
199 children_[nchildren_++] = IntType::Instance(*bra,*ket,aux,oper);
205 for(
unsigned int xyz=0; xyz<3; ++xyz) {
207 if (a.deriv().d(xyz) > 0) {
209 a_der_m1.deriv().dec(xyz);
210 ket->set_member(a_der_m1,part,pos_a);
211 children_[nchildren_++] = IntType::Instance(*bra,*ket,aux,oper);
212 ket->set_member(a,part,pos_a);
215 if (bm1.deriv().d(xyz) > 0) {
217 bm1_der_m1.deriv().dec(xyz);
218 bra->set_member(bm1_der_m1,part,pos_b);
219 children_[nchildren_++] = IntType::Instance(*bra,*ket,aux,oper);
220 bra->set_member(bm1,part,pos_b);
228 expr_ = children_[0] + prefactors.
Y_X[part][dir] * children_[1];
231 const bool aderiv = a.deriv().d(dir_) > 0;
234 a_der_m1.deriv().dec(dir_);
235 ket->set_member(a_der_m1,part,pos_a);
236 children_[nchildren_++] = IntType::Instance(*bra,*ket,aux,oper);
237 ket->set_member(a,part,pos_a);
241 const bool bderiv = bm1.deriv().d(dir_) > 0;
244 bm1_der_m1.deriv().dec(dir_);
245 bra->set_member(bm1_der_m1,part,pos_b);
246 children_[nchildren_++] = IntType::Instance(*bra,*ket,aux,oper);
247 bra->set_member(bm1,part,pos_b);
251 expr_ += Vector(a.deriv())[dir_] * children_[2];
253 expr_ -= Vector(b.deriv())[dir_] * children_[aderiv ? 3 : 2];
257 if (loc_a == InBra && loc_b == InKet) {
258 F a(bra->member(part,pos_a));
259 F b(ket->member(part,pos_b));
262 F bm1(b); bm1.dec(dir_);
268 ket->set_member(bm1,part,pos_b);
271 F ap1(a); ap1.inc(dir_);
272 bra->set_member(ap1,part,pos_a);
273 children_[nchildren_++] = IntType::Instance(*bra,*ket,aux,oper);
274 bra->set_member(a,part,pos_a);
277 children_[nchildren_++] = IntType::Instance(*bra,*ket,aux,oper);
283 for(
unsigned int xyz=0; xyz<3; ++xyz) {
285 if (a.deriv().d(xyz) > 0) {
287 a_der_m1.deriv().dec(xyz);
288 bra->set_member(a_der_m1,part,pos_a);
289 children_[nchildren_++] = IntType::Instance(*bra,*ket,aux,oper);
290 bra->set_member(a,part,pos_a);
293 if (bm1.deriv().d(xyz) > 0) {
295 bm1_der_m1.deriv().dec(xyz);
296 ket->set_member(bm1_der_m1,part,pos_b);
297 children_[nchildren_++] = IntType::Instance(*bra,*ket,aux,oper);
298 ket->set_member(bm1,part,pos_b);
305 expr_ = children_[0] + prefactors.
X_Y[part][dir] * children_[1];
308 const bool aderiv = a.deriv().d(dir_) > 0;
311 a_der_m1.deriv().dec(dir_);
312 bra->set_member(a_der_m1,part,pos_a);
313 children_[nchildren_++] = IntType::Instance(*bra,*ket,aux,oper);
314 bra->set_member(a,part,pos_a);
318 const bool bderiv = bm1.deriv().d(dir_) > 0;
321 bm1_der_m1.deriv().dec(dir_);
322 ket->set_member(bm1_der_m1,part,pos_b);
323 children_[nchildren_++] = IntType::Instance(*bra,*ket,aux,oper);
324 ket->set_member(bm1,part,pos_b);
328 expr_ += Vector(a.deriv())[dir_] * children_[2];
330 expr_ -= Vector(b.deriv())[dir_] * children_[aderiv ? 3 : 2];
339 template <
class IntType,
class F,
int part,
340 FunctionPosition loc_a,
unsigned int pos_a,
341 FunctionPosition loc_b,
unsigned int pos_b>
352 typedef typename IntType::BraType IBraType;
353 typedef typename IntType::KetType IKetType;
354 const IBraType& bra = target_->bra();
355 const IKetType& ket = target_->ket();
358 bool nonzero_quanta =
false;
359 unsigned const int npart = IntType::OperatorType::Properties::np;
360 for(
unsigned int p=0; p<npart; p++) {
363 int nfbra = bra.num_members(p);
365 for(
int f=0; f<nfbra; f++)
366 if (!bra.member(p,f).zero() || !bra.member(p,f).deriv().zero())
367 nonzero_quanta =
true;
368 int nfket = ket.num_members(p);
370 for(
int f=0; f<nfket; f++)
371 if (!ket.member(p,f).zero() || !ket.member(p,f).deriv().zero())
372 nonzero_quanta =
true;
375 if (!nonzero_quanta) {
376 SafePtr<RRStack> rrstack = RRStack::Instance();
377 SafePtr<ThisType> this_ptr =
378 const_pointer_cast<
ThisType,
const ThisType>(
379 static_pointer_cast<
const ThisType,
const ParentType>(
380 EnableSafePtrFromThis<ParentType>::SafePtr_from_this()
383 rrstack->find(this_ptr);
392 IBraType bra_zero(bra);
393 IKetType ket_zero(ket);
394 for(
unsigned int p=0; p<npart; p++) {
397 int nfbra = bra_zero.num_members(p);
398 for(
int f=0; f<nfbra; f++) {
399 typedef typename IBraType::bfs_type bfs_type;
400 typedef typename IBraType::bfs_ref bfs_ref;
401 bfs_ref bfs = bra_zero.member(p,f);
402 if (!bfs.zero() || !bfs.deriv().zero()) {
407 int nfket = ket_zero.num_members(p);
408 for(
int f=0; f<nfket; f++) {
409 typedef typename IKetType::bfs_type bfs_type;
410 typedef typename IKetType::bfs_ref bfs_ref;
411 bfs_ref bfs = ket_zero.member(p,f);
412 if (!bfs.zero() || !bfs.deriv().zero()) {
421 typedef typename IBraType::bfs_type bfs_type;
424 DummyOper dummy_oper;
425 DummyQuanta dummy_quanta(std::vector<int>(0,0));
426 SafePtr<DummyIntegral> dummy_integral = DummyIntegral::Instance(bra_zero,ket_zero,dummy_quanta,dummy_oper);
430 SafePtr<DummyHRR> dummy_hrr = DummyHRR::Instance(dummy_integral,dir_);
431 SafePtr<RRStack> rrstack = RRStack::Instance();
432 rrstack->find(dummy_hrr);
436 template <
class IntType,
class F,
int part,
437 FunctionPosition loc_a,
unsigned int pos_a,
438 FunctionPosition loc_b,
unsigned int pos_b>
444 template <
class IntType,
class F,
int part,
445 FunctionPosition loc_a,
unsigned int pos_a,
446 FunctionPosition loc_b,
unsigned int pos_b>
456 typedef typename IntType::OperatorType
Oper;
457 if (part < 0 || part >= Oper::Properties::np) {
462 if (loc_a == loc_b && pos_a == pos_b) {
468 template <
class IntType,
class F,
int part,
469 FunctionPosition loc_a,
unsigned int pos_a,
470 FunctionPosition loc_b,
unsigned int pos_b>
471 SafePtr<typename HRR<IntType,F,part,loc_a,pos_a,loc_b,pos_b>::ChildType>
474 assert(i>=0 && i<nchildren_);
477 for(
unsigned int c=0; c<max_nchildren_; c++) {
487 template <
class IntType,
class F,
int part,
488 FunctionPosition loc_a,
unsigned int pos_a,
489 FunctionPosition loc_b,
unsigned int pos_b>
495 os <<
"HRR Part " << part <<
" " 496 << (loc_a == InBra ?
"bra" :
"ket") <<
" " << pos_a <<
" " 497 << (loc_b == InBra ?
"bra" :
"ket") <<
" " << pos_b <<
" ";
499 if (loc_a == InBra) {
500 F sh_a(target_->bra(part,pos_a));
504 os << sh_a.
label() <<
" ";
506 if (loc_b == InBra) {
507 F sh_b(target_->bra(part,pos_b));
512 F sh_b(target_->ket(part,pos_b));
518 F sh_a(target_->ket(part,pos_a));
520 os << sh_a.label() <<
" ";
522 if (loc_b == InBra) {
523 F sh_b(target_->bra(part,pos_b));
528 F sh_b(target_->ket(part,pos_b));
537 template <
class IntType,
class F,
int part,
538 FunctionPosition loc_a,
unsigned int pos_a,
539 FunctionPosition loc_b,
unsigned int pos_b>
542 const SafePtr<CodeContext>& context,
const SafePtr<ImplicitDimensions>& dims)
const 545 os << context->label_to_name(
label_to_funcname(context->cparams()->api_prefix() + label()))
549 << context->value_to_pointer(rr_target()->symbol());
551 const unsigned int nchildren = num_children();
552 for(
unsigned int c=0; c<nchildren; c++) {
553 os <<
", " << context->value_to_pointer(rr_child(c)->symbol());
556 unsigned int hsr = 1;
561 for(
int p=0; p<part; p++) {
562 unsigned int nbra = target_->bra().num_members(p);
563 for(
unsigned int i=0; i<nbra; i++) {
564 SubIterator* iter = target_->bra().member_subiter(p,i);
568 unsigned int nket = target_->ket().num_members(p);
569 for(
unsigned int i=0; i<nket; i++) {
570 SubIterator* iter = target_->ket().member_subiter(p,i);
577 taskmgr.
current().params()->max_hrr_hsrank(hsr);
581 if (loc_a == loc_b && pos_a != 0 && pos_b != 0)
582 throw CodeDoesNotExist(
"HRR::spfunction_call -- has not been generalized yet");
585 unsigned int lsr = 1;
586 unsigned int np = IntType::OperType::Properties::np;
587 for(
unsigned int p=part+1; p<np; p++) {
588 unsigned int nbra = target_->bra().num_members(p);
589 for(
unsigned int i=0; i<nbra; i++) {
590 SubIterator* iter = target_->bra().member_subiter(p,i);
594 unsigned int nket = target_->ket().num_members(p);
595 for(
unsigned int i=0; i<nket; i++) {
596 SubIterator* iter = target_->ket().member_subiter(p,i);
602 taskmgr.
current().params()->max_hrr_hsrank(hsr);
608 os <<
")" << context->end_of_stat() << endl;
612 template <
class IntType,
class F,
int part,
613 FunctionPosition loc_a,
unsigned int pos_a,
614 FunctionPosition loc_b,
unsigned int pos_b>
624 template <
class IntType,
class F,
int part,
625 FunctionPosition loc_a,
unsigned int pos_a,
626 FunctionPosition loc_b,
unsigned int pos_b>
630 unsigned int np = IntType::OperType::Properties::np;
641 template <
class IntType,
class F,
int part,
642 FunctionPosition loc_a,
unsigned int pos_a,
643 FunctionPosition loc_b,
unsigned int pos_b>
644 SafePtr<ImplicitDimensions>
647 bool high_rank = expl_high_dim();
648 bool low_rank = expl_low_dim();
650 SafePtr<Entity> high_dim, low_dim;
655 high_dim = dims->high();
661 low_dim = dims->low();
664 SafePtr<ImplicitDimensions> localdims(
new ImplicitDimensions(high_dim,low_dim,dims->vecdim()));
SafePtr< DGVertex > rr_child(unsigned int i) const
Implementation of RecurrenceRelation::child()
Definition: hrr.h:96
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:215
Manages tasks. This is a Singleton.
Definition: task.h:62
TrivialBFSet<T> defines static member result, which is true if T is a basis function set consisting o...
Definition: bfset.h:879
SafePtr< rdouble > Y_X[np][3]
Cartesian components of Y_X vectors.
Definition: prefactors.h:65
Defaults definitions for various parameters assumed by Libint.
Definition: algebra.cc:23
GenIntegralSet is a set of integrals over functions derived from BFS.
Definition: integral.h:89
Definition: stdarray.h:18
Definition: prefactors.h:159
GenOper is a single operator described by descriptor Descr.
Definition: oper.h:152
This is a vertex of a Directed Graph (DG)
Definition: dgvertex.h:42
void current(const std::string &task_label)
Makes this task current (must have been added already)
Definition: task.cc:64
RecurrenceRelation::ExprType ExprType
A short alias.
Definition: hrr.h:66
QuantumNumbersA<T,N> is a set of N quantum numbers of type T implemented in terms of a C-style array...
Definition: quanta.h:198
RTimeEntity is an Entity of type T that exists at runtime of the generated code (hence has no value k...
Definition: entity.h:99
bool is_simple() const
Implementation of RecurrenceRelation::is_simple()
Definition: hrr.h:98
AlgebraicOperator is an algebraic operator that acts on objects of type T.
Definition: algebra.h:47
Oper is OperSet characterized by properties Props.
Definition: oper.h:82
SafePtr< DGVertex > rr_target() const
Implementation of RecurrenceRelation::target()
Definition: hrr.h:94
ImplicitDimensions describes basis functions or other "degrees of freedom" not actively engaged in a ...
Definition: dims.h:43
bool exists(const IncableBFSet &A)
Return true if A is valid.
Definition: bfset.h:91
Iterator provides a base class for all object iterator classes.
Definition: iter.h:44
RecurrenceRelation describes all recurrence relations.
Definition: rr.h:100
Set of basis functions.
Definition: bfset.h:41
virtual unsigned int num_iter() const =0
Returns a number of iterations (number of elements in a set over which to iterate).
const std::string & label() const
label() returns a unique, short, descriptive label of this RR (e.g.
Definition: rr.cc:290
unsigned int num_children() const
Implementation of RecurrenceRelation::num_children()
Definition: hrr.h:88
static bool directional()
is this recurrence relation parameterized by a direction.
Definition: hrr.h:81
A generic Horizontal Recurrence Relation:
Definition: hrr.h:57
SafePtr< rdouble > X_Y[np][3]
Cartesian components of X-Y vectors.
Definition: prefactors.h:57
SafePtr< TargetType > target() const
returns pointer to the target
Definition: hrr.h:90
This exception used to indicate that some code hasn't been developed or generalized yet...
Definition: exception.h:86