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 Lesser 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 Lesser General Public License for more details.
15  *
16  *  You should have received a copy of the GNU Lesser General Public License
17  *  along with Libint.  If not, see <http://www.gnu.org/licenses/>.
18  *
19  */
20 
21 #ifndef _libint2_src_lib_libint_vector_h_
22 #define _libint2_src_lib_libint_vector_h_
23 
24 #if defined(__cplusplus)
25 
26 #include <algorithm>
27 #include <cstddef>
28 
29 #include <libint2/util/type_traits.h>
30 
31 namespace libint2 {
32 
33   /**
34      Contains data types that support SIMD-style computation on vectors of numbers.
35   */
36   namespace simd {
37 
38   // add __declspec(align(N*sizeof(T))) ?
39 
40   /**
41    * Vector<N,T> is used by vectorized Libint library as fixed-length vectors amenable for SIMD-style parallelism
42    * Vectorization via this class should be the last-resort measure if no specialized implementation is available
43    */
44   template <std::size_t N, typename T>
45     struct Vector {
46 
47       T d[N];
48 
49       /**
50        * creates a vector of default-initialized values.
51        */
VectorVector52       Vector() {}
53 
54       /** Initializes all elements to the same value
55        *  @param a the value to which all elements will be set
56        */
VectorVector57       Vector(T a) {
58         std::fill_n(&(d[0]), N, a);
59       }
60 
61       /**
62        * creates a vector of values initialized by an ordinary static-sized array
63        */
VectorVector64       Vector(T (&a)[N]) {
65         std::copy(&a[0], &a[0]+N, &d[0]);
66       }
67 
68       Vector& operator=(T a) {
69         for(std::size_t i=0; i<N; ++i)
70           d[i] = a;
71         return *this;
72       }
73 
74       Vector& operator+=(Vector a) {
75         for(std::size_t i=0; i<N; ++i)
76           d[i] += a.d[i];
77         return *this;
78       }
79 
80       Vector& operator-=(Vector a) {
81         for(std::size_t i=0; i<N; ++i)
82           d[i] -= a.d[i];
83         return *this;
84       }
85 
86       operator double() const {
87         return d[0];
88       }
89 
90   };
91 
92   //@{ arithmetic operators
93   template <std::size_t N, typename T>
94   Vector<N,T> operator*(T a, Vector<N,T> b) {
95     Vector<N,T> c;
96     for(std::size_t i=0; i<N; ++i)
97       c.d[i] = a * b.d[i];
98     return c;
99   }
100 
101   template <std::size_t N, typename T>
102   Vector<N,T> operator*(Vector<N,T> a, T b) {
103     Vector<N,T> c;
104     for(std::size_t i=0; i<N; ++i)
105       c.d[i] = b * a.d[i];
106     return c;
107   }
108 
109   template <std::size_t N, typename T>
110   Vector<N,T> operator*(int a, Vector<N,T> b) {
111     if (a == 1)
112       return b;
113     else {
114       Vector<N, T> c;
115       for (std::size_t i = 0; i < N; ++i)
116         c.d[i] = a * b.d[i];
117       return c;
118     }
119   }
120 
121   template <std::size_t N, typename T>
122   Vector<N,T> operator*(Vector<N,T> a, int b) {
123     if (b == 1)
124       return a;
125     else {
126       Vector<N, T> c;
127       for (std::size_t i = 0; i < N; ++i)
128         c.d[i] = b * a.d[i];
129       return c;
130     }
131   }
132 
133   template <std::size_t N, typename T>
134   Vector<N,T> operator*(Vector<N,T> a, Vector<N,T> b) {
135     Vector<N,T> c;
136     for(std::size_t i=0; i<N; ++i)
137       c.d[i] = a.d[i] * b.d[i];
138     return c;
139   }
140 
141   template <std::size_t N, typename T>
142   Vector<N,T> operator+(Vector<N,T> a, Vector<N,T> b) {
143     Vector<N,T> c;
144     for(std::size_t i=0; i<N; ++i)
145       c.d[i] = a.d[i] + b.d[i];
146     return c;
147   }
148 
149   template <std::size_t N, typename T>
150   Vector<N,T> operator-(Vector<N,T> a, Vector<N,T> b) {
151     Vector<N,T> c;
152     for(std::size_t i=0; i<N; ++i)
153       c.d[i] = a.d[i] - b.d[i];
154     return c;
155   }
156 
157   template <std::size_t N, typename T>
158   Vector<N,T> operator/(Vector<N,T> a, Vector<N,T> b) {
159     Vector<N,T> c;
160     for(std::size_t i=0; i<N; ++i)
161       c.d[i] = a.d[i] / b.d[i];
162     return c;
163   }
164 
165 
166   //@}
167 
168 };}; // namespace libint2::simd
169 
170 namespace libint2 {
171 
172   template <std::size_t N, typename T>
173   struct is_vector<simd::Vector<N,T> > {
174       static const bool value = true;
175   };
176 
177   template <std::size_t N, typename T>
178   struct vector_traits<simd::Vector<N,T> > {
179       typedef T value_type;
180       static const std::size_t extent = N;
181   };
182 
183 } // namespace libint2
184 
185 #include "vector_x86.h"
186 #include "vector_ppc.h"
187 
188 #endif // C++ only
189 
190 #endif
191