src/bin/libint/braket.h

/*
 *  Copyright (C) 2004-2024 Edward F. Valeev
 *
 *  This file is part of Libint compiler.
 *
 *  Libint compiler is free software: you can redistribute it and/or modify
 *  it under the terms of the GNU General Public License as published by
 *  the Free Software Foundation, either version 3 of the License, or
 *  (at your option) any later version.
 *
 *  Libint compiler 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 General Public License for more details.
 *
 *  You should have received a copy of the GNU General Public License
 *  along with Libint compiler.  If not, see <http://www.gnu.org/licenses/>.
 *
 */
 
#ifndef _libint2_src_bin_libint_braket_h_
#define _libint2_src_bin_libint_braket_h_
 
#include <algebra.h>
#include <bfset.h>
#include <global_macros.h>
#include <libint2/util/intrinsic_types.h>
#include <polyconstr.h>
 
namespace libint2 {
 
 
template <class BFS, unsigned int NP>
class ArrayBraket {
 public:
  typedef ArrayBraket<BFS, NP> this_type;
  typedef struct {
  } parent_type;
  typedef ArrayBraket<typename BFS::iter_type, NP> iter_type;
  static constexpr auto num_particles = NP;
  static constexpr auto num_bf = NP;
 
  typedef BFS bfs_type;
  typedef bfs_type& bfs_ref;
  typedef const BFS& bfs_cref;
 
  ArrayBraket();
  ArrayBraket(const BFS* braket);
  ArrayBraket(const std::vector<std::vector<BFS> >& braket);
  ~ArrayBraket();
 
  bool operator==(const this_type&) const;
  bfs_cref member(unsigned int p, unsigned int i = 0) const;
  bfs_ref member(unsigned int p, unsigned int i = 0);
  SubIterator* member_subiter(unsigned int p, unsigned int i = 0) const;
  void set_member(const BFS&, unsigned int p, unsigned int i = 0);
  void set_member(const ConstructablePolymorphically&, unsigned int p,
                  unsigned int i = 0);
  unsigned int num_members(unsigned int p) const {
    assert(p < NP);
    return 1;
  }
  unsigned int num_part() const { return NP; }
  inline LIBINT2_UINT_LEAST64 key() const;
  LIBINT2_UINT_LEAST64 max_key() const;
#if COMPUTE_SIZE_DIRECTLY
  unsigned int size() const;
#endif
 
 private:
  BFS bfs_[NP];
 
  // this function resets all cached values
  void reset_cache();
#if COMPUTE_SIZE_DIRECTLY
  mutable unsigned int size_;
#endif
};
 
template <class BFS, unsigned int NP>
ArrayBraket<BFS, NP>::ArrayBraket() {
#if COMPUTE_SIZE_DIRECTLY
  size_ = 0;
#endif
}
 
template <class BFS, unsigned int NP>
ArrayBraket<BFS, NP>::ArrayBraket(const BFS* braket) {
  for (int i = 0; i < NP; i++) bfs_[i] = braket[i];
#if COMPUTE_SIZE_DIRECTLY
  size_ = 0;
#endif
}
 
template <class BFS, unsigned int NP>
ArrayBraket<BFS, NP>::ArrayBraket(
    const std::vector<std::vector<BFS> >& braket) {
  assert(braket.size() == NP);
  for (unsigned int i = 0; i < NP; i++) {
    assert(braket[i].size() == 1);
    bfs_[i] = braket[i][0];
  }
#if COMPUTE_SIZE_DIRECTLY
  size_ = 0;
#endif
}
 
template <class BFS, unsigned int NP>
ArrayBraket<BFS, NP>::~ArrayBraket() {}
 
template <class BFS, unsigned int NP>
bool ArrayBraket<BFS, NP>::operator==(const ArrayBraket<BFS, NP>& a) const {
  for (int i = 0; i < NP; i++)
    if (bfs_[i] != a[i]) return false;
  return true;
}
 
template <class BFS, unsigned int NP>
typename ArrayBraket<BFS, NP>::bfs_cref ArrayBraket<BFS, NP>::member(
    unsigned int p, unsigned int i) const {
  assert(i == 0 && p < NP);
  return bfs_[p];
}
 
template <class BFS, unsigned int NP>
typename ArrayBraket<BFS, NP>::bfs_ref ArrayBraket<BFS, NP>::member(
    unsigned int p, unsigned int i) {
  assert(i == 0 && p < NP);
  reset_cache();
  return bfs_[p];
}
 
template <class BFS, unsigned int NP>
SubIterator* ArrayBraket<BFS, NP>::member_subiter(unsigned int p,
                                                  unsigned int i) const {
  assert(i == 0 && p < NP);
  return static_cast<SubIterator*>(new SubIteratorBase<BFS>(member(p, i)));
}
 
template <class BFS, unsigned int NP>
void ArrayBraket<BFS, NP>::set_member(const BFS& f, unsigned int p,
                                      unsigned int i) {
  assert(i == 0 && p < NP);
  reset_cache();
  bfs_[p] = f;
}
 
template <class BFS, unsigned int NP>
void ArrayBraket<BFS, NP>::set_member(const ConstructablePolymorphically& f,
                                      unsigned int p, unsigned int i) {
  assert(i == 0 && p < NP);
  reset_cache();
  bfs_[p] = BFS(f);
}
 
template <class BFS, unsigned int NP>
LIBINT2_UINT_LEAST64 ArrayBraket<BFS, NP>::key() const {
  LIBINT2_UINT_LEAST64 pfac = 1;
  LIBINT2_UINT_LEAST64 key = 0;
  for (int p = NP - 1; p >= 0; p--) {
    key += pfac * bfs_[p].key();
    pfac *= BFS::max_key;
  }
  assert(key < this->max_key());
  return key;
}
 
template <class BFS, unsigned int NP>
LIBINT2_UINT_LEAST64 ArrayBraket<BFS, NP>::max_key() const {
  LIBINT2_UINT_LEAST64 max_key = 1;
  for (int p = NP - 1; p >= 0; p--) {
    max_key *= BFS::max_key;
  }
  return max_key;
}
 
#if COMPUTE_SIZE_DIRECTLY
template <class BFS, unsigned int NP>
unsigned int ArrayBraket<BFS, NP>::size() const {
  if (size_ > 0) return size_;
  size_ = 1;
  if (!TrivialBFSet<BFS>::result) {
    for (int p = NP - 1; p >= 0; p--) {
      size_ *= bfs_[p].num_bf();
    }
  }
  return size_;
}
#endif
 
template <class BFS, unsigned int NP>
void ArrayBraket<BFS, NP>::reset_cache() {
  size_ = 0;
}
 
// really need to have typedef template!
template <typename BFS>
struct DefaultOnePBraket {
  typedef ArrayBraket<BFS, 1> Result;
};
 
// really need to have typedef template!
template <typename BFS>
struct DefaultTwoPBraket {
  typedef ArrayBraket<BFS, 2> Result;
};
 
 
enum BraketType { CBra, CKet, PBra, PKet };
template <class BFS, BraketType BKType>
class BraketPair {
 public:
  typedef BraketPair<BFS, BKType> this_type;
  typedef BFS bfs_type;
 
  BraketPair(const BFS& f1, const BFS& f2) : bfs_(f1, f2) {}
  BraketPair(const BraketPair& source) : bfs_(source.bfs_) {}
  BraketPair& operator=(const BraketPair& rhs) {
    bfs_ = rhs.bfs_;
    return *this;
  }
  const BFS& operator[](unsigned int i) const {
    if (i == 0) return bfs_.first;
    if (i == 1) return bfs_.second;
    throw std::logic_error("BraketPair::operator[] -- argument out of range");
  }
  bool operator==(const this_type& rhs) const { return rhs.bfs_ == bfs_; }
 
 private:
  std::pair<BFS, BFS> bfs_;
};
 
namespace braket {
template <class F>
BraketPair<F, PBra> _pbra(const F& f1, const F& f2) {
  return BraketPair<F, PBra>(f1, f2);
}
template <class F>
BraketPair<F, PKet> _pket(const F& f1, const F& f2) {
  return BraketPair<F, PKet>(f1, f2);
}
template <class F>
BraketPair<F, CBra> _cbra(const F& f1, const F& f2) {
  return BraketPair<F, CBra>(f1, f2);
}
template <class F>
BraketPair<F, CKet> _cket(const F& f1, const F& f2) {
  return BraketPair<F, CKet>(f1, f2);
}
};  // namespace braket
 
template <class BFS, BraketType BKTypeL, BraketType BKTypeR>
algebra::Wedge<BraketPair<BFS, BKTypeL>, BraketPair<BFS, BKTypeR> > operator^(
    const BraketPair<BFS, BKTypeL>& L, const BraketPair<BFS, BKTypeR>& R) {
  return algebra::make_wedge(L, R);
}
 
};  // namespace libint2
 
#endif