MultiArrayComplex.h

Go to the documentation of this file.

// Copyright András Vukics 2022–2026. Distributed under the Boost Software License, Version 1.0. (See accompanying file LICENSE.txt)
#pragma once
 
#include "MultiArray.h"
 
#include <Eigen/Dense>
 
#include <boost/serialization/complex.hpp>
 
#include <complex>
 
typedef std::complex<double> dcomp;
 
using std::complex_literals::operator""i;
 
inline bool isNonZero(dcomp c) {return bool(real(c)) || bool(imag(c));}
 
inline bool hasRealPart(dcomp c) {return bool(real(c));}
inline bool hasImagPart(dcomp c) {return bool(imag(c));}
 
inline bool  absCompare(dcomp c1, dcomp c2) {return  abs(c1)< abs(c2);}
inline bool realCompare(dcomp c1, dcomp c2) {return real(c1)<real(c2);}
 
auto sqr(auto&& x) // return type is non-reference or trailing
noexcept(noexcept(x*x)) // propagate noexcept
-> decltype(x*x) // enable SFINAE
{ return x * x; }
 
double relativeDeviation(const auto& a, const auto& b)
{
  const double den = std::abs(a)+std::abs(b);
  return den==0. ? 0. : std::abs(a-b)/den;
}
 
inline double sqrAbs(dcomp v) {return std::norm(v);}
 
inline dcomp textbookProduct(dcomp a, dcomp b) {return {a.real()*b.real()-a.imag()*b.imag(), a.real()*b.imag()+a.imag()*b.real()};}
 
 
namespace boost::json {
    template<class T>
    void tag_invoke(value_from_tag, value& jv, const std::complex<T>& c) {
        jv = {c.real(), c.imag()};
    }
 
    template<class T>
    std::complex<T> tag_invoke(value_to_tag<std::complex<T>>, const value& jv) {
        return { json::value_to<T>(jv.at(0)), json::value_to<T>(jv.at(1)) } ;
    }
} // boost::json
 
using CMatrix=Eigen::MatrixX<dcomp>;
using CVector=Eigen::VectorX<dcomp>;
 
 
namespace cppqedutils {
 
template <typename N> concept scalar = std::convertible_to<N,double> || std::convertible_to<N,dcomp>;
 
template <size_t RANK>
auto vectorize(MultiArray<dcomp,RANK>& ma)
{
  return Eigen::Map<CVector>{ma.dataStorage().data(), Eigen::Index(ma.dataView.size())};
}
 
template <size_t RANK>
auto vectorize(const MultiArray<dcomp,RANK>& ma)
{
  return Eigen::Map<const CVector>{ma.dataStorage().data(), Eigen::Index(ma.dataView.size())};
}
 
 
template <size_t TWO_TIMES_RANK> requires ( !(TWO_TIMES_RANK%2) )
Extents<TWO_TIMES_RANK/2> halveExtents(Extents<TWO_TIMES_RANK> extents)
{
  static constexpr size_t RANK=TWO_TIMES_RANK/2;
  Extents<RANK> res;
  std::ranges::for_each(std::views::iota(0ul,RANK), [&] (size_t i) {
    res[i]=extents[i];
#ifndef   NDEBUG
    if (size_t b=extents[i+RANK]; res[i]!=b)
      throw std::invalid_argument("MultiArray does not have the shape of a matrix: index no. "+std::to_string(i)+": "+std::to_string(res[i])+" vs. "+std::to_string(b));
#endif // NDEBUG
  });
  return res;
}
 
template <size_t TWO_TIMES_RANK>
auto matricize(MultiArray<dcomp,TWO_TIMES_RANK>& ma)
{
  const Eigen::Index matrixDim = Eigen::Index(multiarray::calculateExtent(halveExtents(ma.extents)));
  return Eigen::Map<CMatrix>{ma.dataStorage().data(), matrixDim, matrixDim};
}
 
template <size_t TWO_TIMES_RANK>
auto matricize(const MultiArray<dcomp,TWO_TIMES_RANK>& ma)
{
  const Eigen::Index matrixDim = Eigen::Index(multiarray::calculateExtent(halveExtents(ma.extents)));
  return Eigen::Map<const CMatrix>{ma.dataStorage().data(), matrixDim, matrixDim};
}
 
 
template <size_t RANK1, size_t RANK2>
auto directProduct(const MultiArray<dcomp,RANK1>& m1, const MultiArray<dcomp,RANK2>& m2,
                   std::function<dcomp(dcomp,dcomp)> func=std::multiplies<dcomp>{} )
{
  MultiArray<dcomp,RANK1+RANK2> res{concatenate(m1.extents,m2.extents)};
  for (size_t stride=m1.dataView.size(), i=0; i<stride; ++i) for (size_t j=0; j<m2.dataView.size(); ++j) res.dataStorage()[i+stride*j]=func(m1.dataView[i],m2.dataView[j]);
  return res;
}
 
 
template <size_t RANK>
MultiArray<dcomp,RANK>& conj(MultiArray<dcomp,RANK>& ma)
{
  for (dcomp& v : ma.dataStorage()) v=conj(v);
  return ma;
}
 
 
template <size_t RANK>
double frobeniusNormSqr(const MultiArray<dcomp,RANK>& ma) { return std::ranges::fold_left_first( ma.dataView | std::views::transform(sqrAbs), std::plus{} ).value_or(0.); }
 
template <size_t RANK>
double frobeniusNorm(const MultiArray<dcomp,RANK>& ma) { return sqrt( frobeniusNormSqr(ma) ); }
 
template <size_t TWO_TIMES_RANK>
void hermitianConjugateSelf(MultiArrayView<dcomp,TWO_TIMES_RANK> mav)
{
  checkContiguous(mav,"hermitianConjugateSelf");
  const size_t matrixDim = multiarray::calculateExtent(halveExtents(mav.extents));
  auto _=[&] (size_t i, size_t j) -> dcomp& {return mav.dataView[i+matrixDim*j];};
  for (size_t i=0; i<matrixDim; ++i) {
    _(i,i)=conj(_(i,i));
    for (size_t j=i+1; j<matrixDim; ++j) {
      dcomp &u=_(i,j), &l=_(j,i);
      u=conj(u); l=conj(l);
      std::swap(u, l);
    }
  }
}
 
template <size_t TWO_TIMES_RANK>
void twoTimesRealPartOfSelf(MultiArrayView<dcomp,TWO_TIMES_RANK> mav)
{
  checkContiguous(mav,"twoTimesRealPartOfSelf");
  const size_t matrixDim = multiarray::calculateExtent(halveExtents(mav.extents));
  auto _=[&] (size_t i, size_t j) -> dcomp& {return mav.dataView[i+matrixDim*j];};
  for (size_t i=0; i<matrixDim; ++i) {
    _(i,i)=2.*real(_(i,i));
    for (size_t j=i+1; j<matrixDim; ++j) _(j,i)=conj(_(i,j)=_(i,j)+conj(_(j,i)));
  }
}
 
} // cppqedutils