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