vector.h

/*
 *  Copyright (C) 2004-2021 Edward F. Valeev
 *
 *  This file is part of Libint.
 *
 *  Libint is free software: you can redistribute it and/or modify
 *  it under the terms of the GNU Lesser General Public License as published by
 *  the Free Software Foundation, either version 3 of the License, or
 *  (at your option) any later version.
 *
 *  Libint is distributed in the hope that it will be useful,
 *  but WITHOUT ANY WARRANTY; without even the implied warranty of
 *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 *  GNU Lesser General Public License for more details.
 *
 *  You should have received a copy of the GNU Lesser General Public License
 *  along with Libint.  If not, see <http://www.gnu.org/licenses/>.
 *
 */
 
#ifndef _libint2_src_lib_libint_vector_h_
#define _libint2_src_lib_libint_vector_h_
 
#if defined(__cplusplus)
 
#include <algorithm>
#include <cstddef>
 
#include <libint2/util/type_traits.h>
 
namespace libint2 {
 
  namespace simd {
 
  // add __declspec(align(N*sizeof(T))) ?
 
  template <std::size_t N, typename T>
    struct Vector {
 
      T d[N];
 
      Vector() {}
 
      Vector(T a) {
        std::fill_n(&(d[0]), N, a);
      }
 
      Vector(T (&a)[N]) {
        std::copy(&a[0], &a[0]+N, &d[0]);
      }
 
      Vector& operator=(T a) {
        for(std::size_t i=0; i<N; ++i)
          d[i] = a;
        return *this;
      }
 
      Vector& operator+=(Vector a) {
        for(std::size_t i=0; i<N; ++i)
          d[i] += a.d[i];
        return *this;
      }
 
      Vector& operator-=(Vector a) {
        for(std::size_t i=0; i<N; ++i)
          d[i] -= a.d[i];
        return *this;
      }
 
      operator double() const {
        return d[0];
      }
 
  };
 
  template <std::size_t N, typename T>
  Vector<N,T> operator*(T a, Vector<N,T> b) {
    Vector<N,T> c;
    for(std::size_t i=0; i<N; ++i)
      c.d[i] = a * b.d[i];
    return c;
  }
 
  template <std::size_t N, typename T>
  Vector<N,T> operator*(Vector<N,T> a, T b) {
    Vector<N,T> c;
    for(std::size_t i=0; i<N; ++i)
      c.d[i] = b * a.d[i];
    return c;
  }
 
  template <std::size_t N, typename T>
  Vector<N,T> operator*(int a, Vector<N,T> b) {
    if (a == 1)
      return b;
    else {
      Vector<N, T> c;
      for (std::size_t i = 0; i < N; ++i)
        c.d[i] = a * b.d[i];
      return c;
    }
  }
 
  template <std::size_t N, typename T>
  Vector<N,T> operator*(Vector<N,T> a, int b) {
    if (b == 1)
      return a;
    else {
      Vector<N, T> c;
      for (std::size_t i = 0; i < N; ++i)
        c.d[i] = b * a.d[i];
      return c;
    }
  }
 
  template <std::size_t N, typename T>
  Vector<N,T> operator*(Vector<N,T> a, Vector<N,T> b) {
    Vector<N,T> c;
    for(std::size_t i=0; i<N; ++i)
      c.d[i] = a.d[i] * b.d[i];
    return c;
  }
 
  template <std::size_t N, typename T>
  Vector<N,T> operator+(Vector<N,T> a, Vector<N,T> b) {
    Vector<N,T> c;
    for(std::size_t i=0; i<N; ++i)
      c.d[i] = a.d[i] + b.d[i];
    return c;
  }
 
  template <std::size_t N, typename T>
  Vector<N,T> operator-(Vector<N,T> a, Vector<N,T> b) {
    Vector<N,T> c;
    for(std::size_t i=0; i<N; ++i)
      c.d[i] = a.d[i] - b.d[i];
    return c;
  }
 
  template <std::size_t N, typename T>
  Vector<N,T> operator/(Vector<N,T> a, Vector<N,T> b) {
    Vector<N,T> c;
    for(std::size_t i=0; i<N; ++i)
      c.d[i] = a.d[i] / b.d[i];
    return c;
  }
 
 
 
};}; // namespace libint2::simd
 
namespace libint2 {
 
  template <std::size_t N, typename T>
  struct is_vector<simd::Vector<N,T> > {
      static const bool value = true;
  };
 
  template <std::size_t N, typename T>
  struct vector_traits<simd::Vector<N,T> > {
      typedef T value_type;
      static const std::size_t extent = N;
  };
 
} // namespace libint2
 
#include "vector_x86.h"
#include "vector_ppc.h"
 
#endif // C++ only
 
#endif