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