1 /*
2  *  Copyright (C) 2004-2021 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_rr_h_
22 #define _libint2_src_bin_libint_rr_h_
23 
24 #include <iostream>
25 #include <sstream>
26 #include <string>
27 #include <vector>
28 #include <stdexcept>
29 #include <exception.h>
30 #include <bfset.h>
31 #include <smart_ptr.h>
32 #include <polyconstr.h>
33 #include <singl_stack.h>
34 #include <code.h>
35 #include <default_params.h>
36 #include <util_types.h>
37 #include <global_macros.h>
38 
39 namespace libint2 {
40 
41   class DGVertex;
42   class CodeContext;
43   class ImplicitDimensions;
44   class DirectedGraph;
45   template <typename V> class AlgebraicOperator;
46 
47   /** RRStack implements a stack of RecurrenceRelation's which can only hold
48       one instance of a given RR. RecurrenceRelation::label() is used for hashing
49     */
50   template <typename RR>
51   class RRStackBase : public SingletonStack<RR,std::string>
52   {
53   private:
54 
55     static SafePtr<RRStackBase<RR> > rrstack_;
56 
57     // private constructor because it's a Singleton
RRStackBase()58     RRStackBase() : parent_type(& RR::label) {}
59 
60   public:
61     typedef SingletonStack<RR,std::string> parent_type;
62     typedef typename parent_type::data_type data_type;
63     typedef typename parent_type::value_type value_type;
64     typedef typename parent_type::iter_type iter_type;
65     typedef typename parent_type::citer_type citer_type;
66     typedef typename parent_type::InstanceID InstanceID;
67 
68     /// Obtain the unique Instance of RRStack
Instance()69     static SafePtr<RRStackBase<RR> >& Instance() {
70       if (!rrstack_) {
71         SafePtr<RRStackBase<RR> > tmpstack(new RRStackBase<RR>);
72         rrstack_ = tmpstack;
73       }
74       return rrstack_;
75     }
76 
~RRStackBase()77     virtual ~RRStackBase() {}
78 
79     /// adds content of rrs to this stack
add(const SafePtr<RRStackBase<RR>> & rrs)80     void add(const SafePtr<RRStackBase<RR> >& rrs) {
81       for(citer_type it=rrs->begin(); it != rrs->end(); it++) {
82         find((*it).second.second);
83       }
84     }
85 
86     /// removes rr from the stack
remove(const data_type & rr)87     void remove(const data_type& rr) {
88       parent_type::remove(rr);
89     }
90   };
91 
92   template <typename RR>
93   SafePtr<RRStackBase<RR> > RRStackBase<RR>::rrstack_;
94 
95 
96   /**
97      RecurrenceRelation describes all recurrence relations
98   */
99   class RecurrenceRelation : public EnableSafePtrFromThis<RecurrenceRelation> {
100   public:
101     typedef RecurrenceRelation this_type;
102 
103     RecurrenceRelation();
104     virtual ~RecurrenceRelation();
105 
106     /** num_children() returns the actual number of children.
107         For example, VRR for ERIs has 5 children on the right-hand side,
108         however, for some ERI classes (ss|ps) the actual number may be
109         smaller.
110     */
111     virtual unsigned int num_children() const =0;
112     /// Returns i-th child
113     virtual SafePtr<DGVertex> rr_child(unsigned int i) const =0;
114     /// Returns the target
115     virtual SafePtr<DGVertex> rr_target() const =0;
116 
117     /** Numerical expression of a recurrence relation is always expressed as
118         an AlgebraicOperator<DGVertex> */
119     typedef AlgebraicOperator<DGVertex> ExprType;
120     /// Returns the expression
rr_expr()121     const SafePtr<ExprType>& rr_expr() const { return expr_; }
122 
123     /**
124        Returns true is this recurrence relation is simple enough to optimize away.
125        As a result of such optimization, standalone function will NOT be
126        generated for this recurrence relation. Instead, it's source will be
127        inlined and optimized.
128     */
129     virtual bool is_simple() const =0;
130     /**
131        Returns true is the type of target and all children are exactly the same
132     */
133     virtual bool invariant_type() const;
134     /**
135      * @return 1 if recurrence relation transfers quanta from particle \c from to particle \c to  where \c from < \c to , -1 if \c from > \c to , and 0 if neither
136      */
partindex_direction()137     virtual int partindex_direction() const { return 0; }
138     /**
139      * @return BraketDirection::BraToKet if recurrence relation transfers quanta from function in bra to function in ket,
140      *         BraketDirection::KetToBra if the transfer is from ket to bra, and BraketDirection::None if neither.
141      */
braket_direction()142     virtual BraketDirection braket_direction() const { return BraketDirection::None; }
143     /**
144      * @return the total size of the children of this RR
145      */
146     size_t size_of_children() const;
147 
148     /**
149       label() returns a unique, short, descriptive label of this RR
150       (e.g. "VRR A (p s | 1/r_{12} | d s )" for Obara-Saika recurrence relation
151       applied to center A to compute (ps|ds) ERI)
152     */
153     const std::string& label() const;
154 
155     /**
156       description() returns a verbose description of this RR
157     */
158     virtual std::string description() const;
159 
160     /// Generate declaration and definition for the recurrence relation
161     virtual void generate_code(const SafePtr<CodeContext>& context,
162                                const SafePtr<ImplicitDimensions>& dims,
163                                const std::string& funcname,
164                                std::ostream& decl, std::ostream& def);
165 
166     /// Generate declaration and definition for the recurrence relation
167     /// using generic code (typically, a manually written code)
168     virtual void generate_generic_code(const SafePtr<CodeContext>& context,
169                                        const SafePtr<ImplicitDimensions>& dims,
170                                        const std::string& funcname,
171                                        std::ostream& decl, std::ostream& def);
172 
173     /// Generate a callback for this recurrence relation
174     virtual std::string spfunction_call(const SafePtr<CodeContext>& context,
175                                         const SafePtr<ImplicitDimensions>& dims) const;
176 
177     /// Return the number of FLOPs per this recurrence relation
nflops()178     unsigned int nflops() const { return nflops_; }
179 
180     /// RecurrenceRelation is managed by SingletonStack but doesn't need to keep track of instance ID
inst_id(const SingletonStack<RecurrenceRelation,std::string>::InstanceID & i)181     void inst_id(const SingletonStack<RecurrenceRelation, std::string>::InstanceID& i) {}
182 
183     protected:
184     unsigned int nflops_;
185     mutable std::string label_;
186     SafePtr<ExprType> expr_;
187     /// Adds a (or -a, if minus = -1) to expr_.
188     void add_expr(const SafePtr<ExprType>& a, int minus=1);
189     /// Generates the label
190     virtual std::string generate_label() const =0;
191     /// Registers with the stack
register_with_rrstack()192     template <class RR> bool register_with_rrstack() {
193       // only register RRs with for shell sets
194       if (TrivialBFSet<typename RR::BasisFunctionType>::result)
195         return false;
196       SafePtr<RRStackBase<RecurrenceRelation> > rrstack = RRStackBase<RecurrenceRelation>::Instance();
197       SafePtr<RR> this_ptr =
198         const_pointer_cast<RR,const RR>(
199           static_pointer_cast<const RR, const RecurrenceRelation>(
200             EnableSafePtrFromThis<RecurrenceRelation>::SafePtr_from_this()
201           )
202         );
203       rrstack->find(this_ptr);
204 #if DEBUG || DEBUG_CONSTRUCTION
205       std::cout << "register_with_rrstack: registered " << this_ptr->label() << std::endl;
206 #endif
207       return true;
208     }
209     private:
210 
211     /** used by generate_code to initialize the computation graph that computes sets of integrals using the RR
212     */
213     SafePtr<DirectedGraph> generate_graph_(const SafePtr<DirectedGraph>& dg);
214     /** assigns "target" symbol to the target vertex and "src<i>" to the i-th child vertex. Also
215         appends these symbols to S. */
216     void assign_symbols_(SafePtr<CodeSymbols>& S);
217     /** given an ImplicitDimension for the computation, adapt it for this recurrence
218         relation. Default version does not do anything. */
219     virtual SafePtr<ImplicitDimensions> adapt_dims_(const SafePtr<ImplicitDimensions>& dims) const;
220 
221     /// does this recurrent relation have a generic equivalent? Default is no.
222     virtual bool has_generic(const SafePtr<CompilationParameters>& cparams) const;
223     /// return the name of a header file with the declaration of the generic code
224     virtual std::string generic_header() const;
225     /// return the implementation of this recurrence relation in terms of generic code
226     virtual std::string generic_instance(const SafePtr<CodeContext>& context, const SafePtr<CodeSymbols>& args) const;
227 
228   };
229 
230   namespace algebra {
231     /// these operators are extremely useful to write compact expressions
232     SafePtr<RecurrenceRelation::ExprType> operator+(const SafePtr<DGVertex>& A,
233                                                     const SafePtr<DGVertex>& B);
234     SafePtr<RecurrenceRelation::ExprType> operator-(const SafePtr<DGVertex>& A,
235                                                     const SafePtr<DGVertex>& B);
236     SafePtr<RecurrenceRelation::ExprType> operator*(const SafePtr<DGVertex>& A,
237                                                     const SafePtr<DGVertex>& B);
238     SafePtr<RecurrenceRelation::ExprType> operator/(const SafePtr<DGVertex>& A,
239                                                     const SafePtr<DGVertex>& B);
240     const SafePtr<RecurrenceRelation::ExprType>& operator+=(SafePtr<RecurrenceRelation::ExprType>& A,
241                                                             const SafePtr<DGVertex>& B);
242     const SafePtr<RecurrenceRelation::ExprType>& operator-=(SafePtr<RecurrenceRelation::ExprType>& A,
243                                                             const SafePtr<DGVertex>& B);
244     const SafePtr<RecurrenceRelation::ExprType>& operator*=(SafePtr<RecurrenceRelation::ExprType>& A,
245                                                             const SafePtr<DGVertex>& B);
246     const SafePtr<RecurrenceRelation::ExprType>& operator/=(SafePtr<RecurrenceRelation::ExprType>& A,
247                                                             const SafePtr<DGVertex>& B);
248 
249     class Entity;
250     template <class T> class RTimeEntity;
251     template <class T> class CTimeEntity;
252     // seems to confound Intel compiler on Linux?
253     //SafePtr<RecurrenceRelation::ExprType> operator*(const SafePtr<Entity>& A,
254     //                                                const SafePtr<DGVertex>& B);
255     template<typename T> SafePtr<RecurrenceRelation::ExprType> operator*(const SafePtr<RTimeEntity<T> >& A,
256                                                                          const SafePtr<DGVertex>& B);
257     template<typename T> SafePtr<RecurrenceRelation::ExprType> operator*(const SafePtr<CTimeEntity<T> >& A,
258                                                                          const SafePtr<DGVertex>& B);
259   };
260 
261   // Instantiate the RRStack
262   typedef RRStackBase<RecurrenceRelation> RRStack;
263 
264 };
265 
266 //#include <vrr_11_twoprep_11.h>
267 //#include <hrr.h>
268 //#include <shell_to_ints.h>
269 //#include <iter.h>
270 
271 #endif
272 
273