Tensor.hh

Go to the documentation of this file.

#ifndef TENSOR_HH
#define TENSOR_HH
 
#include <cmath>
#include <iostream>
 
#include "casm/container/Counter.hh"
#include "casm/container/Permutation.hh"
#include "casm/container/Template_Algorithms.hh"
#include "casm/symmetry/SymGroup.hh"
#include "casm/symmetry/SymGroupRep.hh"
 
namespace CASM {
 
//~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
template <class T>
class DynamicMatrix;
 
template <class T>
class ReturnTensor;
 
//~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
template <class T>
class Tensor : public Array<T> {
  Index Nrank;
  Array<Index> Ndim, idim, max_ind;
 
 public:
  using Array<T>::resize;
  using Array<T>::size;
  using Array<T>::at;
 
  // These constructors take two arguments even though they really just need
  // one, since rank==dim.size() if you're doing it right.
  Tensor(Index tNrank, Array<Index> tNdim) : Nrank(tNrank), Ndim(tNdim) {
    Index N = 1;
    idim.resize(Nrank);
    max_ind.resize(Nrank);
    for (Index i = 0; i < Nrank; i++) {
      max_ind[i] = Ndim[i] - 1;
      idim[i] = N;
      N *= Ndim[i];
    }
 
    resize(N);
  };
 
  Tensor(Index tNrank, Array<Index> tNdim, T init_val)
      : Nrank(tNrank), Ndim(tNdim) {
    Index N = 1;
    idim.resize(Nrank);
    max_ind.resize(Nrank);
    for (Index i = 0; i < Nrank; i++) {
      max_ind[i] = Ndim[i] - 1;
      idim[i] = N;
      N *= Ndim[i];
    }
 
    resize(N, init_val);
  };
 
  Tensor(Index tNrank = 0) : Nrank(tNrank), Ndim(tNrank, 3) {
    Index N = 1;
    idim.resize(Nrank);
    max_ind.resize(Nrank);
    for (Index i = 0; i < Nrank; i++) {
      max_ind[i] = Ndim[i] - 1;
      idim[i] = N;
      N *= Ndim[i];
    }
    resize(N);
  };
 
  Tensor(const Eigen::Matrix<T, 3, 3> &mat_init) : Nrank(2), Ndim(2, 3) {
    Index N = 1;
    idim.resize(Nrank);
    max_ind.resize(Nrank);
    for (Index i = 0; i < Nrank; i++) {
      max_ind[i] = Ndim[i] - 1;
      idim[i] = N;
      N *= Ndim[i];
    }
    resize(N);
    for (int i = 0; i < 3; i++)
      for (int j = 0; j < 3; j++) (*this)(i, j) = mat_init(i, j);
  }
 
  Tensor(ReturnTensor<T> &init_cont) { swap(init_cont); };
 
  Tensor &operator=(ReturnTensor<T> &RHS) {
    swap(RHS);
    return *this;
  };
 
  Tensor &operator=(const T &RHS) {
    for (Index i = 0; i < size(); i++) at(i) = RHS;
    return *this;
  }
 
  template <int dim1, int dim2, int flag1>
  Tensor &operator=(const Eigen::Matrix<T, dim1, dim2, flag1> &LHS);
 
  static ReturnTensor<T> identity(Index trank, Index tdim);
 
  void swap(Tensor &RHS);
  void redefine(Array<Index> _Ndim);
  void redefine(Array<Index> _Ndim, T fill_val);
 
  Tensor &dim_permute(const Array<Index> &iperm);
  Tensor &dim_permute(const Permutation &perm);
 
  Tensor &dim_ipermute(const Array<Index> &perm);
  Tensor &dim_ipermute(const Permutation &perm);
 
  Tensor &dim_unpermute();
  bool next_permute();
  void permute_symmetrize();
  void permute_symmetrize(const Array<Array<Index> > &unique_dim);
  void reset_idim();
 
  Counter<Array<Index> > element_counter() const;
 
  Index rank() const;
  const Array<Index> &ind_max() const;
  const Array<Index> &dim() const;
  Index dim(Index i) const;
 
  const Array<Index> &mult_array() const;
 
  const T &get(const Array<Index> &inds) const;
  T &at(const Array<Index> &inds);
  const T &at(const Array<Index> &inds) const;
  T &operator()(const Array<Index> &inds);
  const T &operator()(const Array<Index> &inds) const;
  T &operator()(Index i, Index j);
  const T &operator()(Index i, Index j) const;
 
  Tensor &operator+=(const Tensor &RHS);
  Tensor &operator-=(const Tensor &RHS);
  Tensor &operator*=(T RHS);
  Tensor &operator/=(T RHS);
 
  ReturnTensor<T> operator-();
  ReturnTensor<T> operator+(const Tensor &RHS);
  ReturnTensor<T> operator-(const Tensor &RHS);
 
  // We should find a way to do tensor contraction/multiplication
 
  T scalar_prod(const Tensor &RHS) const;
  ReturnTensor<T> tensor_prod(const Tensor &RHS);
  Tensor &normalize();
 
  Tensor &transform(const Eigen::MatrixXd &op);
 
  // rep_IDs specify the representations that transform each dimension of tensor
  Tensor &transform(const SymOp &op, Array<SymGroupRepID> rep_IDs);
 
  Tensor &apply_sym(const SymOp &op);
 
  Tensor &symmetrize_index(const SymOp &op, Array<SymGroupRepID> rep_IDs,
                           Array<Index> inner_ind);
 
  Tensor slice(Array<Index> slice_Ndim, Array<Index> slice_ind);
 
  void read(std::istream &stream);  // Added by Ivy -- TODO: Make more general
                                    // for not just rank 2 tensors
 
  // Tensor<T> matrix_multiply(Vector3<T> &tvec);
 
  Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic> convert_to_Eigen() const;
};
 
template <class T>
Tensor<T> operator*(const T &LHS, const Tensor<T> &RHS);
 
template <class T>
Tensor<T> operator*(const SymOp &LHS, const Tensor<T> &RHS);
 
template <class T>
std::ostream &operator<<(std::ostream &stream, const Tensor<T> &RHS);
 
template <class T>
std::istream &operator>>(std::istream &stream, Tensor<T> &RHS);
 
//~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
template <class T>
class ReturnTensor : public Tensor<T> {
 public:
  using Tensor<T>::swap;
 
  ReturnTensor(Tensor<T> &init_tens) { swap(init_tens); }
 
  ReturnTensor &operator=(Tensor<T> &RHS) {
    swap(RHS);
    return *this;
  };
};
 
//~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
template <class T>
class TensorBasis : public Array<Tensor<T> > {
 public:
  using Array<Tensor<T> >::size;
  using Array<Tensor<T> >::at;
  using Array<Tensor<T> >::clear;
 
  Array<double> coeffs;
 
  double eci(Index i) const { return coeffs[i]; };
  double &eci(Index i) { return coeffs[i]; };  // Added by Ivy 09/25/2012
 
  void generate_basis(Index Nrank, const SymGroup &sym_group);
  void generate_basis(Index Nrank, const SymGroup &sym_group, Index Rep_ID);
  void generate_basis(Index Nrank, const SymGroup &sym_group,
                      const SymGroupRep &perm_group);
 
  void make_orthogonal_to(const TensorBasis &ortho_basis);
  TensorBasis &apply_sym(const SymOp &op);
  void normalize();
  void idealize();
  bool read(std::istream &stream);  // Added by Ivy
};
 
//************************************************************
 
template <class T>
void Tensor<T>::swap(Tensor<T> &RHS) {
  Index trank(Nrank);
  Nrank = RHS.Nrank;
  RHS.Nrank = trank;
  Ndim.swap(RHS.Ndim);
  idim.swap(RHS.idim);
  max_ind.swap(RHS.max_ind);
  Array<T>::swap(RHS);
  return;
}
 
//************************************************************
template <class T>
template <int dim1, int dim2, int flag1>
Tensor<T> &Tensor<T>::operator=(
    const Eigen::Matrix<T, dim1, dim2, flag1> &LHS) {
  Array<Index> tdim(2);
  if (rank() != 2 || dim(0) != LHS.rows() || dim(1) != LHS.cols()) {
    tdim[0] = LHS.rows();
    tdim[1] = LHS.cols();
    redefine(tdim);
  }
  for (int i = 0; i < LHS.rows(); i++) {
    tdim[0] = i;
    for (int j = 0; j < LHS.rows(); j++) {
      tdim[1] = j;
      at(tdim) = LHS(i, j);
    }
  }
  return (*this);
}
 
//************************************************************
 
template <class T>
void Tensor<T>::redefine(Array<Index> _Ndim) {
  Nrank = _Ndim.size();
  Ndim = _Ndim;
  Index N = 1;
  idim.resize(Nrank);
  max_ind.resize(Nrank);
  for (Index i = 0; i < Nrank; i++) {
    max_ind[i] = Ndim[i] - 1;
    idim[i] = N;
    N *= Ndim[i];
  }
 
  resize(N);
}
 
//************************************************************
 
template <class T>
void Tensor<T>::redefine(Array<Index> _Ndim, T fill_val) {
  Nrank = _Ndim.size();
  Ndim = _Ndim;
  Index N = 1;
  idim.resize(Nrank);
  max_ind.resize(Nrank);
  for (Index i = 0; i < Nrank; i++) {
    max_ind[i] = Ndim[i] - 1;
    idim[i] = N;
    N *= Ndim[i];
  }
 
  resize(N, fill_val);
}
 
//************************************************************
 
template <class T>
Tensor<T> &Tensor<T>::dim_permute(const Array<Index> &iperm) {
  if (iperm.size() != Nrank) {
    // temporarily changed this to cout
    throw std::runtime_error(
        "Tensor<T>::dim_permute, permutation array is incompatible with tensor "
        "rank!");
    return *this;
  }
 
  idim.permute(iperm);
  Ndim.permute(iperm);
  max_ind.permute(iperm);
 
  return *this;
}
 
//************************************************************
 
template <class T>
Tensor<T> &Tensor<T>::dim_permute(const Permutation &perm) {
  return dim_permute(perm.perm_array());
}
//************************************************************
 
template <class T>
Tensor<T> &Tensor<T>::dim_ipermute(const Array<Index> &iperm) {
  if (iperm.size() != Nrank) {
    std::cerr << "WARNING: In Tensor<T>::dim_permute, permutation array is "
                 "incompatible with tensor rank! Continuing without permuting "
                 "indeces...\n";
    return *this;
  }
 
  idim.ipermute(iperm);
  Ndim.ipermute(iperm);
  max_ind.ipermute(iperm);
 
  return *this;
}
 
//************************************************************
 
template <class T>
Tensor<T> &Tensor<T>::dim_ipermute(const Permutation &perm) {
  return dim_ipermute(perm.perm_array());
}
//************************************************************
 
template <class T>
Tensor<T> &Tensor<T>::dim_unpermute() {
  Array<Index> tperm;
 
  idim.sort(tperm);
  Ndim.permute(tperm);
  max_ind.permute(tperm);
 
  return *this;
}
 
//************************************************************
 
template <class T>
bool Tensor<T>::next_permute() {
  Index i = Nrank - 2;
  bool is_valid = true;
  Array<Index> iperm;
  for (Index k = 0; k < Nrank; k++) iperm[k] = k;
 
  while (valid_index(i) && idim[i + 1] <= idim[i]) i--;
 
  Index j = Nrank - 1;
 
  if (valid_index(i)) {
    while (idim[j] <= idim[i]) j--;
 
    idim.swap_elem(i, j);
    Ndim.swap_elem(i, j);
    max_ind.swap_elem(i, j);
 
    i++;
    j = Nrank - 1;
  }
 
  else {
    is_valid = false;
    i = 0;
  }
 
  while (i < j && valid_index(j)) {
    idim.swap_elem(i, j);
    Ndim.swap_elem(i, j);
    max_ind.swap_elem(i, j);
 
    i++;
    j--;
  }
 
  return is_valid;
}
 
//************************************************************
 
template <class T>
void Tensor<T>::reset_idim() {
  if (idim.is_ascending()) return;
 
  Tensor ttens(Nrank, Ndim);
  Counter<Array<Index> > icount(ttens.element_counter());
 
  do {
    ttens(icount()) = at(icount());
  } while ((++icount).valid());
 
  swap(ttens);
 
  return;
}
 
//************************************************************
 
template <class T>
void Tensor<T>::permute_symmetrize() {
  reset_idim();
  Tensor<T> ttens(*this);
  Index nperm = 1;
  while (ttens.next_permute()) {
    (*this) += ttens;
    nperm++;
  }
  (*this) /= nperm;
  return;
}
 
//************************************************************
 
template <class T>
void Tensor<T>::permute_symmetrize(const Array<Array<Index> > &unique_dim) {
  reset_idim();
  Array<Array<Index> > tunique(unique_dim);
  for (Index i = 0; i < tunique.size(); i++) tunique[i].sort();
  Array<Index> tperm = Ndim;
  Tensor<T> ttens(*this);
  Index iset(0), i(0), j(0), nperm(1);
  (*this) = 0;
  while (iset < tunique.size()) {
    for (i = 0; i < tunique.size(); i++)
      for (j = 0; j < tunique[i].size(); j++)
        tperm[unique_dim[i][j]] = tunique[i][j];
    (*this) += ttens.dim_permute(tperm);
    ttens.dim_unpermute();
    nperm++;
    iset = 0;
    while (iset < tunique.size() && !tunique[iset].next_permute()) iset++;
  }
  (*this) /= double(nperm);
}
 
//************************************************************
 
template <class T>
Counter<Array<Index> > Tensor<T>::element_counter() const {
  return Counter<Array<Index> >(Array<Index>(Nrank, 0), max_ind,
                                Array<Index>(Nrank, 1));
}
 
//************************************************************
 
template <class T>
const Array<Index> &Tensor<T>::ind_max() const {
  return max_ind;
}
 
//************************************************************
 
template <class T>
const Array<Index> &Tensor<T>::dim() const {
  return Ndim;
}
 
//************************************************************
 
template <class T>
Index Tensor<T>::dim(Index i) const {
  return Ndim[i];
}
 
//************************************************************
 
template <class T>
const Array<Index> &Tensor<T>::mult_array() const {
  return idim;
}
 
//************************************************************
 
template <class T>
Index Tensor<T>::rank() const {
  return Nrank;
}
 
//************************************************************
 
template <class T>
const T &Tensor<T>::get(const Array<Index> &inds) const {
  Index lin_index = 0;
  for (Index i = 0; i < Nrank; i++) lin_index += inds[i] * idim[i];
 
  return at(lin_index);
};
 
//************************************************************
 
template <class T>
T &Tensor<T>::at(const Array<Index> &inds) {
  Index lin_index = 0;
  for (Index i = 0; i < Nrank; i++) lin_index += inds[i] * idim[i];
 
  return at(lin_index);
};
 
//************************************************************
 
template <class T>
const T &Tensor<T>::at(const Array<Index> &inds) const {
  Index lin_index = 0;
  for (Index i = 0; i < Nrank; i++) lin_index += inds[i] * idim[i];
 
  return at(lin_index);
};
 
//************************************************************
 
template <class T>
T &Tensor<T>::operator()(const Array<Index> &inds) {
  Index lin_index = 0;
  for (Index i = 0; i < Nrank; i++) lin_index += inds[i] * idim[i];
 
  return at(lin_index);
};
 
//************************************************************
 
template <class T>
const T &Tensor<T>::operator()(const Array<Index> &inds) const {
  Index lin_index = 0;
  for (Index i = 0; i < Nrank; i++) lin_index += inds[i] * idim[i];
 
  return at(lin_index);
};
//************************************************************
 
template <class T>
T &Tensor<T>::operator()(Index i, Index j) {
  return at(i * idim[0] + j * idim[1]);
};
 
//************************************************************
 
template <class T>
const T &Tensor<T>::operator()(Index i, Index j) const {
  return at(i * idim[0] + j * idim[1]);
};
 
//************************************************************
 
template <class T>
Tensor<T> &Tensor<T>::operator+=(const Tensor &RHS) {
  if (Ndim != RHS.Ndim) {
    std::cerr << "Attempting to add two incompatible tensors!\n";
    exit(1);
  }
 
  if (idim == RHS.idim) {
    for (Index i = 0; i < size(); i++) at(i) += RHS[i];
 
    return *this;
  }
 
  Counter<Array<Index> > icount(Array<Index>(Nrank, 0), max_ind,
                                Array<Index>(Nrank, 1));
  do {
    at(icount()) += RHS(icount());
 
  } while ((++icount).valid());
 
  return *this;
}
 
//************************************************************
 
template <class T>
Tensor<T> &Tensor<T>::operator-=(const Tensor &RHS) {
  if (Ndim != RHS.Ndim) {
    std::cerr << "Attempting to add two incompatible tensors!\n";
    exit(1);
  }
 
  if (idim == RHS.idim) {
    for (Index i = 0; i < size(); i++) at(i) -= RHS[i];
 
    return *this;
  }
 
  Counter<Array<Index> > icount(Array<Index>(Nrank, 0), max_ind,
                                Array<Index>(Nrank, 1));
  do {
    at(icount()) -= RHS(icount());
  } while ((++icount).valid());
 
  return *this;
}
 
//************************************************************
 
template <class T>
Tensor<T> &Tensor<T>::operator*=(T RHS) {
  for (Index i = 0; i < size(); i++) at(i) *= RHS;
 
  return *this;
}
 
//************************************************************
 
template <class T>
Tensor<T> &Tensor<T>::operator/=(T RHS) {
  for (Index i = 0; i < size(); i++) at(i) /= RHS;
 
  return *this;
}
 
//************************************************************
 
template <class T>
ReturnTensor<T> Tensor<T>::operator-() {
  Tensor<T> ttens(*this);
  for (Index i = 0; i < ttens.size(); i++) ttens[i] = -at(i);
  return ReturnTensor<T>(ttens);
}
 
//************************************************************
 
template <class T>
ReturnTensor<T> Tensor<T>::operator+(const Tensor &RHS) {
  return ReturnTensor<T>(Tensor<T>(*this) += RHS);
}
 
//************************************************************
 
template <class T>
ReturnTensor<T> Tensor<T>::operator-(const Tensor &RHS) {
  return ReturnTensor<T>(Tensor<T>(*this) -= RHS);
}
 
//************************************************************
 
template <class T>
T Tensor<T>::scalar_prod(const Tensor &RHS) const {
  if (Ndim != RHS.Ndim) {
    std::cerr
        << "Attempting to take scalar product of two incompatible tensors!\n";
    assert(0);
    // exit(1);
  }
 
  T tprod = 0;
 
  if (idim == RHS.idim) {
    for (Index i = 0; i < size(); i++) tprod += at(i) * RHS[i];
  }
 
  else {
    Counter<Array<Index> > icount(Array<Index>(Nrank, 0), max_ind,
                                  Array<Index>(Nrank, 1));
    do {
      tprod += at(icount()) * RHS(icount());
    } while ((++icount).valid());
  }
 
  return tprod;
}
 
//************************************************************
 
template <class T>
ReturnTensor<T> Tensor<T>::tensor_prod(const Tensor &RHS) {
  Tensor<T> ttens(rank() + RHS.rank(), array_cat(dim(), RHS.dim()), 0);
 
  Counter<Array<Index> > icount(element_counter()),
      jcount(RHS.element_counter());
 
  Index i, j, tot_ind;
  do {
    jcount.reset();
    do {
      tot_ind = 0;
      for (i = 0; i < icount().size(); i++)
        tot_ind += icount[i] * ttens.idim[i];
 
      for (j = 0; j < jcount().size(); j++)
        tot_ind += jcount[j] * ttens.idim[i + j];
 
      ttens[tot_ind] = at(icount()) * RHS(jcount());
 
    } while ((++jcount).valid());
  } while ((++icount).valid());
 
  return ReturnTensor<T>(ttens);
}
 
//************************************************************
 
template <class T>
Tensor<T> &Tensor<T>::normalize() {
  return (*this) /= sqrt(scalar_prod(*this));
}
 
//************************************************************
 
template <class T>
ReturnTensor<T> Tensor<T>::identity(Index trank, Index tdim) {
  Tensor<T> ttens(trank, Array<Index>(trank, tdim), 0);
  for (Index i = 0; i < tdim; i++) ttens(Array<Index>(trank, i)) = 1;
 
  return ReturnTensor<T>(ttens);
}
 
//************************************************************
template <class T>
Tensor<T> &Tensor<T>::transform(const SymOp &op, Array<SymGroupRepID> rep_IDs) {
  if (!rank()) return *this;
  Array<Eigen::MatrixXd const *> rep_mats(op.get_matrix_reps(rep_IDs));
 
  Tensor<T> ttens(Nrank, Ndim);
  Index nd;
 
  T tval;
  Counter<Array<Index> > icount(element_counter());
  Counter<Array<Index> > jcount(icount);
  do {
    ttens(icount()) = 0;
    jcount.reset();
    do {
      tval = at(jcount());
      for (nd = 0; nd < Nrank; nd++)
        tval *= (*rep_mats[nd])(icount[nd], jcount[nd]);
      ttens(icount()) += tval;
 
    } while ((++jcount).valid());
 
  } while ((++icount).valid());
  swap(ttens);
  return *this;
}
 
//************************************************************
template <class T>
Tensor<T> &Tensor<T>::symmetrize_index(const SymOp &op,
                                       Array<SymGroupRepID> rep_IDs,
                                       Array<Index> inner_ind) {
  if (!rank()) return *this;
  Array<Eigen::MatrixXd const *> rep_mats(op.get_matrix_reps(rep_IDs));
 
  Index nd;
 
  T tval;
  Counter<Array<Index> > icount(element_counter());
  do {
    tval = 1.0;
    for (nd = 0; nd < Nrank; nd++)
      tval *= (*rep_mats[nd])(icount[nd], inner_ind[nd]);
    at(icount()) += tval;
 
  } while ((++icount).valid());
  return *this;
}
 
//************************************************************
template <class T>
Tensor<T> &Tensor<T>::apply_sym(const SymOp &op) {
  if (!rank()) return *this;
 
  Tensor<T> ttens(Nrank, Ndim);
  Index nd;
 
  T tval;
  Counter<Array<Index> > icount(element_counter());
  Counter<Array<Index> > jcount(icount);
  do {
    ttens(icount()) = 0;
    jcount.reset();
    do {
      tval = at(jcount());
      for (nd = 0; nd < Nrank; nd++)
        tval *= op.matrix()(icount[nd], jcount[nd]);
      ttens(icount()) += tval;
 
    } while ((++jcount).valid());
 
  } while ((++icount).valid());
  swap(ttens);
  return *this;
}
 
//************************************************************
//************************************************************
template <class T>
void Tensor<T>::read(std::istream &stream) {
  double val;
  Array<Index> min(rank(), 0), max(ind_max()), inc(rank(), 1);
  Counter<Array<Index> > icount(min, max, inc);
 
  do {
    stream >> val;
    at(icount()) = val;
  } while ((++icount).valid());
 
  return;
};
 
//************************************************************
//************************************************************
/*template<class T>
Tensor<T> Tensor<T>::matrix_multiply(Vector3<T> &tvec) {
 
  Index j = 0;
  Array<Index> tensordim(2, 0);
  tensordim[0] = 3;
  tensordim[1] = 1;
  Tensor<T> result(2, tensordim, 0);
 
  for(Index i = 0; i < size(); i++) {
    if(i % 3 == 0) {
      j++;
    }
    result.at(j - 1) += at(i) * tvec.at(i % 3);
  }
 
  return result;
 
};
*/
//************************************************************
 
template <class T>
Tensor<T> &Tensor<T>::transform(const Eigen::MatrixXd &mat) {
  if (!rank()) return *this;
 
  Tensor<T> ttens(Nrank, Ndim);
  Index nd;
 
  T tval;
  Counter<Array<Index> > icount(element_counter());
  Counter<Array<Index> > jcount(icount);
  do {
    ttens(icount()) = 0;
    jcount.reset();
    do {
      tval = at(jcount());
      for (nd = 0; nd < Nrank; nd++) tval *= mat(icount[nd], jcount[nd]);
      ttens(icount()) += tval;
 
    } while ((++jcount).valid());
 
  } while ((++icount).valid());
  swap(ttens);
  return *this;
}
 
//************************************************************
 
template <class T>
Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic> Tensor<T>::convert_to_Eigen()
    const {
  if (rank() != 2) {
    std::cerr << "WARNING: Attempting to assign Eigen::Matrix object using a "
                 "Tensor of rank "
              << rank()
              << "!\nTensor must have rank=2 for assignment to be valid. "
                 "Exiting...\n";
    exit(1);
  }
  Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic> tmat(dim(0), dim(1));
  tmat.setZero();
  for (Index i = 0; i < tmat.rows(); i++) {
    for (Index j = 0; j < tmat.rows(); j++) {
      tmat(i, j) = (*this)(i, j);
    }
  }
  return tmat;
}
 
//************************************************************
 
template <class T>
Tensor<T> operator*(const T &LHS, const Tensor<T> &RHS) {
  return Tensor<T>(RHS) *= LHS;
}
 
//************************************************************
 
template <class T>
Tensor<T> operator*(const SymOp &LHS, const Tensor<T> &RHS) {
  return Tensor<T>(RHS).apply_sym(LHS);
}
 
//************************************************************
 
template <class T>
T dot(const Tensor<T> &LHS, const Tensor<T> &RHS) {
  return LHS.scalar_prod(RHS);
}
 
//************************************************************
 
template <class T>
T norm(const Tensor<T> &ttens) {
  return sqrt(ttens.scalar_prod(ttens));
}
 
//************************************************************
 
template <class T>
std::ostream &operator<<(std::ostream &stream, const Tensor<T> &RHS) {
  std::ios::fmtflags old_flags = stream.flags();
  stream.flags(std::ios::showpoint | std::ios::fixed | std::ios::right);
  int twide = stream.width();
 
  if (!RHS.rank()) {
    stream << RHS[0] << '\n';
    stream.flags(old_flags);
    return stream;
  }
 
  Array<Index> min(RHS.rank(), 0), max(RHS.ind_max()), inc(RHS.rank(), 1);
  Counter<Array<Index> > icount(min, max, inc);
  do {
    if (RHS.rank() > 2 && icount[0] == 0 && icount[1] == 0) {
      stream << "Block (:,:," << icount[2];
      for (Index j = 3; j < RHS.rank(); j++) stream << "," << icount[j];
      stream << ") :\n";
    }
 
    stream << std::setw(0) << "  " << std::setw(twide) << RHS(icount());
    if (icount[0] == RHS.ind_max()[0]) {
      stream << '\n';
 
      // if(RHS.rank() > 1 && icount[1] == RHS.ind_max()[1])
      //   stream << "\n"; //commented out by Ivy 07/01/13 so that an extra new
      //   line isn't printed.
    }
 
  } while ((++icount).valid());
  stream.flags(old_flags);
  return stream;
}
 
//************************************************************
 
template <class T>
std::istream &operator>>(std::istream &stream, Tensor<T> &RHS) {
  std::cerr << "WARNING: Input stream operator has not been implemented for "
               "Tensor class!  Exiting...\n";
  exit(1);
}
 
//************************************************************
 
template <class T>
void TensorBasis<T>::generate_basis(Index Nrank, const SymGroup &sym_group) {
  Index i, ns;
  clear();
 
  Tensor<T> ttens(Nrank, Array<Index>(Nrank, 3), 0);
 
  for (i = 0; i < ttens.size(); i++) {
    Tensor<T> avtens(Nrank, Array<Index>(Nrank, 3), 0);
    ttens[i] = 1;
    for (ns = 0; ns < sym_group.size(); ns++) avtens += sym_group[ns] * ttens;
    if (!this->contains(avtens)) this->push_back(avtens);
    ttens[i] = 0;
  }
  idealize();
 
  coeffs.resize(size(), NAN);  // Changed by Ivy 09/25/2012
}
 
//************************************************************
 
template <class T>
void TensorBasis<T>::generate_basis(Index Nrank, const SymGroup &sym_group,
                                    Index Rep_ID) {
  if (!sym_group.size()) {
    std::cerr << "WARNING: Attempting to generate tensor basis from an invalid "
                 "symmetry group. Exiting...";
    exit(0);
  }
  Index i, ns;
  clear();
 
  Tensor<T> ttens(Nrank, Array<Index>(Nrank, 3), 0);
 
  for (i = 0; i < ttens.size(); i++) {
    Tensor<T> avtens(Nrank, Array<Index>(Nrank, 3), 0);
    ttens[i] = 1;
    for (ns = 0; ns < sym_group.size(); ns++) avtens += sym_group[ns] * ttens;
    if (!contains(avtens)) push_back(avtens);
    ttens[i] = 0;
  }
  idealize();
 
  coeffs.resize(size(), NAN);  // Changed by Ivy 09/25/2012
}
 
//************************************************************
 
template <class T>
void TensorBasis<T>::generate_basis(Index Nrank, const SymGroup &sym_group,
                                    const SymGroupRep &perm_group) {
  Index i, ns;
  clear();
 
  if (sym_group.size() != perm_group.size()) {
    std::cerr << "WARNING: Attempting to generate tensor basis, but symmetry "
                 "group and permutation group are not compatible!\n";
  }
 
  Tensor<T> ttens(Nrank, Array<Index>(Nrank, 3), 0);
 
  for (i = 0; i < ttens.size(); i++) {
    Tensor<T> avtens(Nrank, Array<Index>(Nrank, 3), 0);
    ttens[i] = 1;
    for (ns = 0; ns < sym_group.size(); ns++)
      avtens += (sym_group[ns] * ttens)
                    .dim_permute(*(perm_group[ns]->get_permutation()));
    if (!this->contains(avtens)) this->push_back(avtens);
    ttens[i] = 0;
  }
  idealize();
 
  coeffs.resize(size(), NAN);  // Changed by Ivy 09/25/2012
  return;
}
 
//************************************************************
 
template <class T>
void TensorBasis<T>::make_orthogonal_to(const TensorBasis<T> &ortho_basis) {
  Index i, j;
  for (i = 0; i < size(); i++) {
    for (j = 0; j < ortho_basis.size(); j++) {
      at(i) -= (at(i).scalar_prod(ortho_basis[j]) / norm(ortho_basis[j])) *
               ortho_basis[j];
    }
  }
  return;
}
 
//************************************************************
 
template <class T>
TensorBasis<T> &TensorBasis<T>::apply_sym(const SymOp &op) {
  for (Index i = 0; i < size(); i++) {
    at(i).apply_sym(op);
  }
  return *this;
}
 
//************************************************************
 
template <class T>
void TensorBasis<T>::normalize() {
  for (Index i = 0; i < size(); i++) at(i).normalize();
  return;
}
 
//************************************************************
 
template <class T>
void TensorBasis<T>::idealize() {
  // GaussJordan<Tensor<T>, T>::eliminate(*this);
 
  //    for(Index i=0; i<size(); i++)
  // std::cout << "Basis Tensor " << i << ":  \n" << at(i);
 
  GramSchmidt<Tensor<T>, T>::orthogonalize(*this);
  return;
}
 
//************************************************************
template <class T>
bool TensorBasis<T>::read(std::istream &stream) {
  char ch;
  Index num_elem, rank;
  Array<Index> dim;
  double tcoeff;
 
  ch = stream.peek();
  while (ch < '0' || ch > '9') {
    stream.ignore(256, '\n');
    ch = stream.peek();
  }
  stream >> num_elem;
  if (!num_elem) {
    stream.ignore(256, '\n');
    return false;
  }
  stream >> rank;
 
  dim.resize(rank);
 
  for (Index i = 0; i < rank; i++) {
    stream >> dim[i];
  }
 
  stream.ignore(1000, '\n');
 
  coeffs.resize(num_elem, NAN);
 
  Tensor<T> ttensor(rank, dim);
 
  for (Index j = 0; j < num_elem; j++) {
    ttensor.redefine(dim);
 
    // Check for "<" to find "<ECI>"
    ch = stream.peek();
    while (ch != char(60)) {
      if (((ch < '0') || (ch > '9')) && (ch != '-')) {
        stream.ignore(1000, '\n');
      } else {
        break;
      }
      ch = stream.peek();
    }
 
    // Reading in coefficient
    if (((ch != char(60)) && ((ch >= '0') || (ch <= '9'))) || (ch == '-')) {
      stream >> tcoeff;
      eci(j) = tcoeff;
    }
 
    stream.ignore(1000, '\n');
 
    ttensor.read(stream);
    this->push_back(ttensor);
  }
 
  return true;
}
 
};  // namespace CASM
 
#endif