Declare RR target pointers LIBINT2_RESTRICT

CHANGES

@@ -4,6 +4,8 @@
 Following is a brief summary of changes made in each release of Libint.
 
 - unreleased:
+  - PR #449: generated and generic RR functions declare their target LIBINT2_RESTRICT (new config.h
+             macro), up to ~1.2x faster high-am and derivative ERIs
   - PR #446: a value of 0 for LIBINT2_OPT_AM (or for the per-class LIBINT2_<class>_MAX_AM /
              LIBINT2_<class>_OPT_AM) is now honored; previously config.h dropped the macro
              (#cmakedefine treats 0 as false), which broke the compiler build for LIBINT2_OPT_AM

include/libint2/config.h.cmake.in

@@ -362,6 +362,17 @@
 #define LIBINT_PRAGMA_GCC(x)
 #endif
 
+/* ----------- restrict ---------------*/
+/* C99 restrict for C++. Do not use __restrict with GCC-compatible compilers:
+   macOS <sys/cdefs.h> defines it to nothing in C++ */
+#if defined(__GNUC__) || defined(__clang__)
+#define LIBINT2_RESTRICT __restrict__
+#elif defined(_MSC_VER)
+#define LIBINT2_RESTRICT __restrict
+#else
+#define LIBINT2_RESTRICT
+#endif
+
 #ifdef __has_cpp_attribute
 #if __has_cpp_attribute(deprecated)
 #define LIBINT_DEPRECATED(msg) [[deprecated(msg)]]

src/bin/libint/context.cc

@@ -578,6 +578,9 @@ std::string CppCodeContext::const_modifier() const { return "const "; }
 std::string CppCodeContext::mutable_modifier() const {
   return "#ifdef __cplusplus\nmutable \n#endif\n";
 }
+std::string CppCodeContext::restrict_modifier() const {
+  return "LIBINT2_RESTRICT ";
+}
 
 std::string CppCodeContext::inteval_type_name(const std::string& tlabel) const {
   if (cparams()->single_evaltype())

src/bin/libint/context.h

@@ -170,6 +170,8 @@ class CodeContext {
   virtual std::string const_modifier() const = 0;
   /// returns the modifier for mutable variables
   virtual std::string mutable_modifier() const = 0;
+  /// returns the modifier for pointers that do not alias any other pointer
+  virtual std::string restrict_modifier() const = 0;
 
  protected:
   /// Lone constructor takes CompilationParams
@@ -321,6 +323,7 @@ class CppCodeContext : public CodeContext,
   std::string ptr_fp_type() const override;
   std::string const_modifier() const override;
   std::string mutable_modifier() const override;
+  std::string restrict_modifier() const override;
 
   std::string start_expr() const;
   std::string end_expr() const;

src/bin/libint/dg.cc

@@ -1012,8 +1012,12 @@ std::string declare_function(const std::shared_ptr<CodeContext>& context,
       << "* inteval";
   const unsigned int nargs = args->n();
   if (nargs > 0) {
-    // first argument is always the target, which is never a const
-    oss << ", " << context->type_name<double*>() << " " << args->symbol(0);
+    // first argument is always the target (the stack, for prerequisite
+    // functions), which is never a const and is not accessed through any
+    // other pointer: RR targets never overlap their children (checked in
+    // DirectedGraph::print_def)
+    oss << ", " << context->type_name<double*>() << " "
+        << context->restrict_modifier() << args->symbol(0);
     for (unsigned int a = 1; a < nargs; a++) {
       oss << ", " << context->type_name<const double*>() << " "
           << args->symbol(a);
@@ -2093,6 +2097,23 @@ void DirectedGraph::print_def(const std::shared_ptr<CodeContext>& context,
                 *(current_vertex->first_exit_arc()));
         if (arc_ptr) {
           std::shared_ptr<RecurrenceRelation> rr = arc_ptr->rr();
+          // RR functions declare their target restrict, hence it must not
+          // overlap any of the children
+          if (current_vertex->address_set()) {
+            const auto t_begin = current_vertex->address();
+            const auto t_end = t_begin + current_vertex->size();
+            for (auto a = current_vertex->first_exit_arc();
+                 a != current_vertex->plast_exit_arc(); ++a) {
+              const std::shared_ptr<DGVertex>& child = (*a)->dest();
+              if (!child->address_set()) continue;
+              const auto c_begin = child->address();
+              const auto c_end = c_begin + child->size();
+              if (c_begin < t_end && t_begin < c_end)
+                throw ProgrammingError(
+                    "DirectedGraph::print_def() -- target of " + rr->label() +
+                    " overlaps its child " + child->label());
+            }
+          }
           os << rr->spfunction_call(context, dims);
           nflops_total += rr->nflops();
 

src/lib/libint/GenericGaussDeriv.h

@@ -39,7 +39,7 @@ template <int L,  //!< the angular momentum of this shell (i.e. the shell being
           bool vectorize>
 struct GenericGaussDeriv {
   static void compute(
-      const Libint_t* inteval, LIBINT2_REALTYPE* target,
+      const Libint_t* inteval, LIBINT2_REALTYPE* LIBINT2_RESTRICT target,
       const LIBINT2_REALTYPE* src0, const LIBINT2_REALTYPE* src1,
       unsigned int
           highdim,  //!< number of functions more significant than this shell

src/lib/libint/GenericGaussDeriv.impl.h

@@ -33,7 +33,7 @@ template <int L,  //!< the angular momentum of this shell (i.e. the shell being
                   //!< differentiated)
           bool vectorize>
 void GenericGaussDeriv<L, vectorize>::compute(
-    const Libint_t* inteval, LIBINT2_REALTYPE* target,
+    const Libint_t* inteval, LIBINT2_REALTYPE* LIBINT2_RESTRICT target,
     const LIBINT2_REALTYPE* src0, const LIBINT2_REALTYPE* src1,
     unsigned int
         highdim,  //!< number of functions more significant than this shell

src/lib/libint/ITR_xs_xs.h

@@ -31,7 +31,8 @@ namespace libint2 {
 
 template <int part, int La, int Lc, bool InBra, bool vectorize>
 struct ITR_xs_xs {
-  static void compute(const Libint_t* inteval, LIBINT2_REALTYPE* target,
+  static void compute(const Libint_t* inteval,
+                      LIBINT2_REALTYPE* LIBINT2_RESTRICT target,
                       const LIBINT2_REALTYPE* src0,
                       const LIBINT2_REALTYPE* src1,
                       const LIBINT2_REALTYPE* src2,
@@ -46,7 +47,8 @@ struct ITR_xs_xs {
  **/
 template <int La, int Lc, bool InBra, bool vectorize>
 struct ITR_xs_xs<0, La, Lc, InBra, vectorize> {
-  static void compute(const Libint_t* inteval, LIBINT2_REALTYPE* target,
+  static void compute(const Libint_t* inteval,
+                      LIBINT2_REALTYPE* LIBINT2_RESTRICT target,
                       const LIBINT2_REALTYPE* src0,
                       const LIBINT2_REALTYPE* src1,
                       const LIBINT2_REALTYPE* src2,
@@ -238,7 +240,7 @@ struct ITR_xs_xs<0, La, Lc, InBra, vectorize> {
   template <int La, int Lc, bool InBra, bool vectorize> struct ITR_xs_xs<1,La,Lc,InBra,vectorize> {
 
     static void compute(const Libint_t* inteval,
-        LIBINT2_REALTYPE* target,
+        LIBINT2_REALTYPE* LIBINT2_RESTRICT target,
         const LIBINT2_REALTYPE* src0,
         const LIBINT2_REALTYPE* src1,
         const LIBINT2_REALTYPE* src2,

src/lib/libint/OSVRR_sx_sx.h

@@ -35,7 +35,8 @@ namespace libint2 {
 
 template <int part, int Lb, int Ld, bool unit_a, bool vectorize>
 struct OSVRR_sx_sx {
-  static void compute(const Libint_t* inteval, LIBINT2_REALTYPE* target,
+  static void compute(const Libint_t* inteval,
+                      LIBINT2_REALTYPE* LIBINT2_RESTRICT target,
                       const LIBINT2_REALTYPE* src1,
                       const LIBINT2_REALTYPE* src0,
                       const LIBINT2_REALTYPE* src2,
@@ -52,7 +53,8 @@ struct OSVRR_sx_sx {
  **/
 template <int Lb, int Ld, bool unit_a, bool vectorize>
 struct OSVRR_sx_sx<0, Lb, Ld, unit_a, vectorize> {
-  static void compute(const Libint_t* inteval, LIBINT2_REALTYPE* target,
+  static void compute(const Libint_t* inteval,
+                      LIBINT2_REALTYPE* LIBINT2_RESTRICT target,
                       const LIBINT2_REALTYPE* src0,
                       const LIBINT2_REALTYPE* src1,
                       const LIBINT2_REALTYPE* src2,
@@ -196,7 +198,7 @@ struct OSVRR_sx_sx<0, Lb, Ld, unit_a, vectorize> {
   template <int Lb, int Ld, bool vectorize> struct OSVRR_sx_sx<1,Lb,Ld,vectorize> {
 
     static void compute(const Libint_t* inteval,
-        LIBINT2_REALTYPE* target,
+        LIBINT2_REALTYPE* LIBINT2_RESTRICT target,
         const LIBINT2_REALTYPE* src0,
         const LIBINT2_REALTYPE* src1,
         const LIBINT2_REALTYPE* src2,
@@ -343,7 +345,8 @@ struct OSVRR_sx_sx<0, Lb, Ld, unit_a, vectorize> {
 
 template <int part, int Lb, int Ld, bool vectorize>
 struct OSAVRR_sx_sx {
-  static void compute(const Libint_t* inteval, LIBINT2_REALTYPE* target,
+  static void compute(const Libint_t* inteval,
+                      LIBINT2_REALTYPE* LIBINT2_RESTRICT target,
                       const LIBINT2_REALTYPE* src1,
                       const LIBINT2_REALTYPE* src4);
 };
@@ -354,7 +357,8 @@ struct OSAVRR_sx_sx {
  **/
 template <int Lb, int Ld, bool vectorize>
 struct OSAVRR_sx_sx<0, Lb, Ld, vectorize> {
-  static void compute(const Libint_t* inteval, LIBINT2_REALTYPE* target,
+  static void compute(const Libint_t* inteval,
+                      LIBINT2_REALTYPE* LIBINT2_RESTRICT target,
                       const LIBINT2_REALTYPE* src1,
                       const LIBINT2_REALTYPE* src4) {
     // works for (sd|sp) and higher

src/lib/libint/OSVRR_sx_sx_deriv.h

@@ -34,7 +34,7 @@ template <int part, int Lb, int Ld, int Da_x, int Da_y, int Da_z, int Db_x,
           int Dd_z, bool unit_a, bool vectorize>
 struct OSVRR_sx_sx_deriv {
   static void compute(
-      const Libint_t* inteval, LIBINT2_REALTYPE* target,
+      const Libint_t* inteval, LIBINT2_REALTYPE* LIBINT2_RESTRICT target,
       const LIBINT2_REALTYPE* src0, const LIBINT2_REALTYPE* src1,
       const LIBINT2_REALTYPE* src2, const LIBINT2_REALTYPE* src3,
       const LIBINT2_REALTYPE* src4, const LIBINT2_REALTYPE* src5,
@@ -84,7 +84,7 @@ template <int Lb, int Ld, int Da_x, int Da_y, int Da_z, int Db_x, int Db_y,
 struct OSVRR_sx_sx_deriv<0, Lb, Ld, Da_x, Da_y, Da_z, Db_x, Db_y, Db_z, Dc_x,
                          Dc_y, Dc_z, Dd_x, Dd_y, Dd_z, unit_a, vectorize> {
   static void compute(
-      const Libint_t* inteval, LIBINT2_REALTYPE* target,
+      const Libint_t* inteval, LIBINT2_REALTYPE* LIBINT2_RESTRICT target,
       const LIBINT2_REALTYPE* src0, const LIBINT2_REALTYPE* src1,
       const LIBINT2_REALTYPE* src2, const LIBINT2_REALTYPE* src3,
       const LIBINT2_REALTYPE* src4, const LIBINT2_REALTYPE* src5,
@@ -389,7 +389,7 @@ template <int part, int Lb, int Ld, int Da_x, int Da_y, int Da_z, int Db_x,
           int Dd_z, bool vectorize>
 struct OSAVRR_sx_sx_deriv {
   static void compute(
-      const Libint_t* inteval, LIBINT2_REALTYPE* target,
+      const Libint_t* inteval, LIBINT2_REALTYPE* LIBINT2_RESTRICT target,
       const LIBINT2_REALTYPE* src1, const LIBINT2_REALTYPE* src4,
       const LIBINT2_REALTYPE* src5, const LIBINT2_REALTYPE* src6,
       const LIBINT2_REALTYPE* src7, const LIBINT2_REALTYPE* src8,
@@ -434,7 +434,7 @@ template <int Lb, int Ld, int Da_x, int Da_y, int Da_z, int Db_x, int Db_y,
 struct OSAVRR_sx_sx_deriv<0, Lb, Ld, Da_x, Da_y, Da_z, Db_x, Db_y, Db_z, Dc_x,
                           Dc_y, Dc_z, Dd_x, Dd_y, Dd_z, vectorize> {
   static void compute(
-      const Libint_t* inteval, LIBINT2_REALTYPE* target,
+      const Libint_t* inteval, LIBINT2_REALTYPE* LIBINT2_RESTRICT target,
       const LIBINT2_REALTYPE* src1, const LIBINT2_REALTYPE* src4,
       const LIBINT2_REALTYPE* src5, const LIBINT2_REALTYPE* src6,
       const LIBINT2_REALTYPE* src7, const LIBINT2_REALTYPE* src8,

src/lib/libint/OSVRR_xs_xs.h

@@ -31,7 +31,8 @@ namespace libint2 {
 
 template <int part, int La, int Lc, bool unit_b, bool vectorize>
 struct OSVRR_xs_xs {
-  static void compute(const Libint_t* inteval, LIBINT2_REALTYPE* target,
+  static void compute(const Libint_t* inteval,
+                      LIBINT2_REALTYPE* LIBINT2_RESTRICT target,
                       const LIBINT2_REALTYPE* src0,
                       const LIBINT2_REALTYPE* src1,
                       const LIBINT2_REALTYPE* src2,
@@ -48,7 +49,8 @@ struct OSVRR_xs_xs {
  **/
 template <int La, int Lc, bool unit_b, bool vectorize>
 struct OSVRR_xs_xs<0, La, Lc, unit_b, vectorize> {
-  static void compute(const Libint_t* inteval, LIBINT2_REALTYPE* target,
+  static void compute(const Libint_t* inteval,
+                      LIBINT2_REALTYPE* LIBINT2_RESTRICT target,
                       const LIBINT2_REALTYPE* src0,
                       const LIBINT2_REALTYPE* src1,
                       const LIBINT2_REALTYPE* src2,
@@ -187,7 +189,8 @@ struct OSVRR_xs_xs<0, La, Lc, unit_b, vectorize> {
 // Ahlrichs' extension of OS VRR
 template <int part, int La, int Lc, bool vectorize>
 struct OSAVRR_xs_xs {
-  static void compute(const Libint_t* inteval, LIBINT2_REALTYPE* target,
+  static void compute(const Libint_t* inteval,
+                      LIBINT2_REALTYPE* LIBINT2_RESTRICT target,
                       const LIBINT2_REALTYPE* src1,
                       const LIBINT2_REALTYPE* src4);
 };
@@ -198,7 +201,8 @@ struct OSAVRR_xs_xs {
  **/
 template <int La, int Lc, bool vectorize>
 struct OSAVRR_xs_xs<0, La, Lc, vectorize> {
-  static void compute(const Libint_t* inteval, LIBINT2_REALTYPE* target,
+  static void compute(const Libint_t* inteval,
+                      LIBINT2_REALTYPE* LIBINT2_RESTRICT target,
                       const LIBINT2_REALTYPE* src1,
                       const LIBINT2_REALTYPE* src4) {
     // works for (ps|ps) and higher

src/lib/libint/OSVRR_xs_xs_deriv.h

@@ -35,7 +35,7 @@ template <int part, int La, int Lc, int Da_x, int Da_y, int Da_z, int Db_x,
           int Dd_z, bool unit_b, bool vectorize>
 struct OSVRR_xs_xs_deriv {
   static void compute(
-      const Libint_t* inteval, LIBINT2_REALTYPE* target,
+      const Libint_t* inteval, LIBINT2_REALTYPE* LIBINT2_RESTRICT target,
       const LIBINT2_REALTYPE* src0, const LIBINT2_REALTYPE* src1,
       const LIBINT2_REALTYPE* src2, const LIBINT2_REALTYPE* src3,
       const LIBINT2_REALTYPE* src4, const LIBINT2_REALTYPE* src5,
@@ -85,7 +85,7 @@ template <int La, int Lc, int Da_x, int Da_y, int Da_z, int Db_x, int Db_y,
 struct OSVRR_xs_xs_deriv<0, La, Lc, Da_x, Da_y, Da_z, Db_x, Db_y, Db_z, Dc_x,
                          Dc_y, Dc_z, Dd_x, Dd_y, Dd_z, unit_b, vectorize> {
   static void compute(
-      const Libint_t* inteval, LIBINT2_REALTYPE* target,
+      const Libint_t* inteval, LIBINT2_REALTYPE* LIBINT2_RESTRICT target,
       const LIBINT2_REALTYPE* src0, const LIBINT2_REALTYPE* src1,
       const LIBINT2_REALTYPE* src2, const LIBINT2_REALTYPE* src3,
       const LIBINT2_REALTYPE* src4, const LIBINT2_REALTYPE* src5,
@@ -394,7 +394,7 @@ template <int part, int La, int Lc, int Da_x, int Da_y, int Da_z, int Db_x,
           bool vectorize>
 struct OSAVRR_xs_xs_deriv {
   static void compute(
-      const Libint_t* inteval, LIBINT2_REALTYPE* target,
+      const Libint_t* inteval, LIBINT2_REALTYPE* LIBINT2_RESTRICT target,
       const LIBINT2_REALTYPE* src1, const LIBINT2_REALTYPE* src4,
       const LIBINT2_REALTYPE* src5, const LIBINT2_REALTYPE* src6,
       const LIBINT2_REALTYPE* src7, const LIBINT2_REALTYPE* src8,
@@ -439,7 +439,7 @@ template <int La, int Lc, int Da_x, int Da_y, int Da_z, int Db_x, int Db_y,
 struct OSAVRR_xs_xs_deriv<0, La, Lc, Da_x, Da_y, Da_z, Db_x, Db_y, Db_z, Dc_x,
                           Dc_y, Dc_z, Dd_x, Dd_y, Dd_z, vectorize> {
   static void compute(
-      const Libint_t* inteval, LIBINT2_REALTYPE* target,
+      const Libint_t* inteval, LIBINT2_REALTYPE* LIBINT2_RESTRICT target,
       const LIBINT2_REALTYPE* src1, const LIBINT2_REALTYPE* src4,
       const LIBINT2_REALTYPE* src5, const LIBINT2_REALTYPE* src6,
       const LIBINT2_REALTYPE* src7, const LIBINT2_REALTYPE* src8,

src/lib/libint/VRR_GTG_1d_xx_xx.h

@@ -42,7 +42,8 @@ namespace libint2 {
 template <unsigned int CartesianAxis, int La, int Lb, int Lc, int Ld,
           bool vectorize>
 struct VRR_GTG_1d_xx_xx {
-  static void compute(const Libint_t* inteval, LIBINT2_REALTYPE* target,
+  static void compute(const Libint_t* inteval,
+                      LIBINT2_REALTYPE* LIBINT2_RESTRICT target,
                       const LIBINT2_REALTYPE* src0) {
     enum XYZ { x = 0, y = 1, z = 2 };
     assert(CartesianAxis == x || CartesianAxis == y || CartesianAxis == z);

src/lib/libint/VRR_GTG_1d_xx_xx_vec.h

@@ -42,7 +42,8 @@ namespace libint2 {
 template <unsigned int CartesianAxis, int La, int Lb, int Lc, int Ld,
           bool vectorize>
 struct VRR_GTG_1d_xx_xx {
-  static void compute(const Libint_t* inteval, VectorSIMD<double, npts>* target,
+  static void compute(const Libint_t* inteval,
+                      VectorSIMD<double, npts>* LIBINT2_RESTRICT target,
                       VectorSIMD<double, npts>* src0) {
     enum XYZ { x = 0, y = 1, z = 2 };
     assert(CartesianAxis == x || CartesianAxis == y || CartesianAxis == z);

src/lib/libint/VRR_r12kg12_xs_xs.h

@@ -31,13 +31,11 @@ namespace libint2 {
 
 template <int part, int La, int Lc, int K, bool vectorize>
 struct VRR_r12kg12_xs_xs {
-  static void compute(const Libint_t* inteval, LIBINT2_REALTYPE* target,
-                      const LIBINT2_REALTYPE* src0,
-                      const LIBINT2_REALTYPE* src1,
-                      const LIBINT2_REALTYPE* src2,
-                      const LIBINT2_REALTYPE* src3,
-                      const LIBINT2_REALTYPE* src4,
-                      const LIBINT2_REALTYPE* src5);
+  static void compute(
+      const Libint_t* inteval, LIBINT2_REALTYPE* LIBINT2_RESTRICT target,
+      const LIBINT2_REALTYPE* src0, const LIBINT2_REALTYPE* src1,
+      const LIBINT2_REALTYPE* src2, const LIBINT2_REALTYPE* src3,
+      const LIBINT2_REALTYPE* src4, const LIBINT2_REALTYPE* src5);
 };
 
 /** builds (a0| G_K |c0), where G_K = r12^K * G12, for K >= 0. For K == -1 use
@@ -51,13 +49,11 @@ struct VRR_r12kg12_xs_xs {
  **/
 template <int La, int Lc, int K, bool vectorize>
 struct VRR_r12kg12_xs_xs<0, La, Lc, K, vectorize> {
-  static void compute(const Libint_t* inteval, LIBINT2_REALTYPE* target,
-                      const LIBINT2_REALTYPE* src0,
-                      const LIBINT2_REALTYPE* src1,
-                      const LIBINT2_REALTYPE* src2,
-                      const LIBINT2_REALTYPE* src3,
-                      const LIBINT2_REALTYPE* src4,
-                      const LIBINT2_REALTYPE* src5) {
+  static void compute(
+      const Libint_t* inteval, LIBINT2_REALTYPE* LIBINT2_RESTRICT target,
+      const LIBINT2_REALTYPE* src0, const LIBINT2_REALTYPE* src1,
+      const LIBINT2_REALTYPE* src2, const LIBINT2_REALTYPE* src3,
+      const LIBINT2_REALTYPE* src4, const LIBINT2_REALTYPE* src5) {
     // works for (ds|ps) and higher, K >= 0
     if ((La < 2 && Lc < 1) || K < 0) abort();