LIBINT  2.6.0
algebra.h
1 /*
2  * Copyright (C) 2004-2019 Edward F. Valeev
3  *
4  * This file is part of Libint.
5  *
6  * Libint is free software: you can redistribute it and/or modify
7  * it under the terms of the GNU General Public License as published by
8  * the Free Software Foundation, either version 3 of the License, or
9  * (at your option) any later version.
10  *
11  * Libint is distributed in the hope that it will be useful,
12  * but WITHOUT ANY WARRANTY; without even the implied warranty of
13  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14  * GNU General Public License for more details.
15  *
16  * You should have received a copy of the GNU General Public License
17  * along with Libint. If not, see <http://www.gnu.org/licenses/>.
18  *
19  */
20 
21 #ifndef _libint2_src_bin_libint_algebra_h_
22 #define _libint2_src_bin_libint_algebra_h_
23 
24 #include <smart_ptr.h>
25 #include <rr.h>
26 #include <exception.h>
27 #include <global_macros.h>
28 #include <dgvertex.h>
29 #include <class_registry.h>
30 
31 namespace libint2 {
32 
33  namespace algebra {
34  struct OperatorTypes {
35  typedef enum {Plus, Minus, Times, Divide} OperatorType;
36  };
37  static const char OperatorSymbol[][2] = { "+", "-", "*", "/" };
38  };
39 
40 
47  template <class T>
49  public DGVertex
50  {
51 
52  public:
54  typedef algebra::OperatorTypes::OperatorType OperatorType;
55 
56  AlgebraicOperator(OperatorType OT,
57  const SafePtr<T>& left,
58  const SafePtr<T>& right) :
59  DGVertex(ClassInfo<AlgebraicOperator>::Instance().id()), OT_(OT), left_(left), right_(right),
60  label_(algebra::OperatorSymbol[OT_])
61  {
62  }
63  virtual ~AlgebraicOperator() {}
64 
66  AlgebraicOperator(const SafePtr<AlgebraicOperator>& A,
67  const SafePtr<T>& left,
68  const SafePtr<T>& right) :
69  DGVertex(static_cast<DGVertex&>(*A)), OT_(A->OT_),
70  left_(left), right_(right), label_(A->label_)
71  {
72 #if DEBUG
73  if (num_exit_arcs() != 2)
74  cout << "AlgebraicOperator<DGVertex> copy constructor: number of children != 2" << endl;
75  else {
76 
77  typedef DGVertex::ArcSetType::const_iterator aciter;
78  aciter a = this->first_exit_arc();
79  auto left_arg = (*a)->dest(); ++a;
80  auto right_arg = (*a)->dest();
81 
82  if (left_ != left_arg && left_ != right_arg)
83  cout << "AlgebraicOperator<DGVertex> copy constructor: invalid left operand given" << endl;
84  if (right_ != left_arg && right_ != right_arg)
85  cout << "AlgebraicOperator<DGVertex> copy constructor: invalid right operand given" << endl;
86  }
87 #endif
88  }
89 
91  OperatorType type() const { return OT_; }
93  const SafePtr<T>& left() const { return left_; }
95  const SafePtr<T>& right() const { return right_; }
96 
98  void add_exit_arc(const SafePtr<DGArc>& a);
100  unsigned int size() const { return 1; }
102  bool equiv(const SafePtr<DGVertex>& a) const
103  {
104  if (typeid_ == a->typeid_) {
105 #if ALGEBRAICOPERATOR_USE_KEY_TO_COMPARE
106 #if USE_INT_KEY_TO_COMPARE
107  if (key() == a->key())
108  return *this == static_pointer_cast<AlgebraicOperator,DGVertex>(a);
109  else
110  return false;
111 #else
112  return description() == a->description();
113 #endif
114 #else
115  return *this == static_pointer_cast<AlgebraicOperator,DGVertex>(a);
116 #endif
117  }
118  else
119  return false;
120  }
121 
123  bool operator==(const SafePtr<AlgebraicOperator>& a) const {
124 #if ALGEBRAICOPERATOR_USE_SAFEPTR
125  // Find out why sometimes equivalent left_ and a->left_ have non-equivalent pointers
126  if (left_->equiv(a->left()) && left_ != a->left_) {
127  cout << "Left arguments are equivalent but pointers differ!" << endl;
128  cout << left_->description() << endl;
129  cout << a->left_->description() << endl;
130  }
131  // Find out why sometimes equivalent right_ and a->right_ have non-equivalent pointers
132  if (right_->equiv(a->right()) && right_ != a->right_) {
133  cout << "Left arguments are equivalent but pointers differ!" << endl;
134  cout << right_->description() << endl;
135  cout << a->right_->description() << endl;
136  }
137 #endif
138  if (OT_ == a->OT_) {
139 #if ALGEBRAICOPERATOR_USE_SAFEPTR
140  if (left_ == a->left_ && right_ == a->right_)
141 #else
142  if (left_->equiv(a->left()) && right_->equiv(a->right()))
143 #endif
144  return true;
145  else
146  return false;
147  }
148  else
149  return false;
150  }
151 
153  const std::string& label() const
154  {
155  return label_;
156  }
158  const std::string& id() const
159  {
160  return label();
161  }
163  std::string description() const
164  {
165  ostringstream os;
166  os << "( ( " << left_->description() << " ) "
167  << algebra::OperatorSymbol[OT_] << " ( "
168  << right_->description() << " ) )";
169  const std::string descr = os.str();
170  return descr;
171  }
172 
175  {
176  throw CannotAddArc("AlgebraicOperator::del_exit_arcs() -- cannot safely use del_exit_arcs on operator vertices.");
177  }
178 
180  typename DGVertex::KeyReturnType key() const { return 0; }
181 
182  void print(std::ostream& os) const {
183  DGVertex::print(os);
184  using std::endl;
185  std::string prefix("AlgebraicOperator::print: ");
186  os << prefix << "this = " << this << endl;
187  os << prefix << "left_ = " << left_ << endl;
188  os << prefix << "right_ = " << right_ << endl;
189  }
190 
191  private:
192  OperatorType OT_;
193  SafePtr<T> left_;
194  SafePtr<T> right_;
195 
197  bool this_precomputed() const
198  {
199  return false;
200  }
201 
202  std::string label_;
203  };
204 
205  /*
206  template <>
207  void
208  AlgebraicOperator<DGVertex>::add_exit_arc(const SafePtr<DGArc>& a)
209  {
210  DGVertex::add_exit_arc(a);
211  if (left_->equiv(a->dest()))
212  left_ = a->dest();
213  else if (right_->equiv(a->dest()))
214  right_ = a->dest();
215  else
216  throw std::runtime_error("AlgebraicOperator<DGVertex>::add_exit_arc -- trying to add an arc to a vertex not equivalent to either argument.");
217  }
218  */
219 
222  template <class C, class T>
224  public:
225  typedef std::pair<C,T> term_t;
226  typedef std::vector<term_t> data_t;
227 
228  LinearCombination() {}
229  // shallow copy constructor -- used by operator^
230  LinearCombination(data_t* data) { swap(*data,data_); delete data; }
231 
232  LinearCombination& operator+=(const term_t& t) {
233  data_.push_back(t);
234  return *this;
235  }
236  size_t size() const { return data_.size(); }
237  const term_t& operator[](unsigned int i) const { return data_[i]; }
238 
239  private:
240  data_t data_;
241  };
242 
243  namespace algebra {
245  template <class L, class R>
246  struct Wedge {
247  Wedge(const L& l, const R& r) : left(l), right(r) {}
248  L left;
249  R right;
250  };
251  template <class L, class R> Wedge<L,R> make_wedge(const L& l, const R& r) { return Wedge<L,R>(l,r); }
252 
254  template <class C, class Tl, class Tr>
255  typename LinearCombination<C, Wedge<Tl,Tr> >::term_t
256  wedge(const std::pair<C,Tl>& L,
257  const std::pair<C,Tr>& R) {
258  return make_pair(L.first*R.first,
259  L.second ^ R.second
260  );
261  }
262  }
263 
264  template <class C, class Tl, class Tr>
265  typename LinearCombination<C, algebra::Wedge<Tl,Tr> >::data_t*
266  operator^(const LinearCombination<C,Tl>& L,
267  const LinearCombination<C,Tr>& R) {
268  typedef typename LinearCombination<C, algebra::Wedge<Tl,Tr> >::data_t data_t;
269  data_t* result = new data_t;
270  const size_t nL = L.size();
271  const size_t nR = R.size();
272  for(unsigned int l=0; l<nL; ++l)
273  for(unsigned int r=0; r<nR; ++r) {
274  result->push_back(algebra::wedge(L[l],R[r]));
275  }
276  return result;
277  }
278 
279  template <class C, class Tl, class Tr>
280  typename LinearCombination<C, algebra::Wedge<Tl,Tr> >::data_t*
281  operator^(const Tl& L,
282  const LinearCombination<C,Tr>& R) {
283  typedef typename LinearCombination<C, algebra::Wedge<Tl,Tr> >::data_t data_t;
284  data_t* result = new data_t;
285  const size_t nR = R.size();
286  for(unsigned int r=0; r<nR; ++r) {
287  result->push_back(L ^ R[r]);
288  }
289  return result;
290  }
291 
292 };
293 
294 #endif
295 
virtual void print(std::ostream &os) const
print(os) prints vertex info to os
Definition: dgvertex.cc:457
unsigned int size() const
Implements DGVertex::size()
Definition: algebra.h:100
const std::string & id() const
Implements DGVertex::id()
Definition: algebra.h:158
void print(std::ostream &os) const
print(os) prints vertex info to os
Definition: algebra.h:182
bool equiv(const SafePtr< DGVertex > &a) const
Implements DGVertex::equiv()
Definition: algebra.h:102
Defaults definitions for various parameters assumed by Libint.
Definition: algebra.cc:24
const SafePtr< T > & right() const
Returns the right argument.
Definition: algebra.h:95
Wedge is a typeholder for the result of a wedge product.
Definition: algebra.h:246
This is a vertex of a Directed Graph (DG)
Definition: dgvertex.h:43
DGVertex(ClassID tid)
Sets typeid to tid.
Definition: dgvertex.cc:31
std::string description() const
Implements DGVertex::description()
Definition: algebra.h:163
void del_exit_arcs()
Overloads DGVertex::del_exit_arcs()
Definition: algebra.h:174
Definition: algebra.h:34
represents linear combination of objects of type T with coefficients of type C
Definition: algebra.h:223
OperatorType type() const
Returns the OperatorType.
Definition: algebra.h:91
AlgebraicOperator is an algebraic operator that acts on objects of type T.
Definition: algebra.h:48
DGVertex::KeyReturnType key() const
Implements Hashable::key()
Definition: algebra.h:180
const std::string & label() const
Implements DGVertex::label()
Definition: algebra.h:153
bool operator==(const SafePtr< AlgebraicOperator > &a) const
laboriously compare 2 operators element by element
Definition: algebra.h:123
Objects of this type provide limited information about the class at runtime.
Definition: class_registry.h:44
Definition: exception.h:40
void add_exit_arc(const SafePtr< DGArc > &a)
Overloads DGVertex::add_exit_arc(). The default definition is used unless T = DGVertex (see algebra....
const SafePtr< T > & left() const
Returns the left argument.
Definition: algebra.h:93
AlgebraicOperator(const SafePtr< AlgebraicOperator > &A, const SafePtr< T > &left, const SafePtr< T > &right)
Clone A but replace operands with left and right.
Definition: algebra.h:66
ClassID typeid_
typeid stores the ClassID of the concrete type.
Definition: dgvertex.h:68
ArcSetType::const_iterator first_exit_arc() const
returns children::begin()
Definition: dgvertex.h:115
unsigned int num_exit_arcs() const
returns the number of children
Definition: dgvertex.cc:247