srconf.h

Go to the documentation of this file.

/*
  This file is part of MADNESS.
 
  Copyright (C) 2007,2010 Oak Ridge National Laboratory
 
  This program 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 2 of the License, or
  (at your option) any later version.
 
  This program 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 this program; if not, write to the Free Software
  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
 
  For more information please contact:
 
  Robert J. Harrison
  Oak Ridge National Laboratory
  One Bethel Valley Road
  P.O. Box 2008, MS-6367
 
  email: harrisonrj@ornl.gov
  tel:   865-241-3937
  fax:   865-572-0680
 
  $Id$
*/
 
/// \file SRConf.h
/// \brief handles the low-level details of a separated representation tensor
 
#ifndef SRCONF_H_
#define SRCONF_H_
 
//#define BENCH 0
 
#include <madness/world/print.h>
#include <madness/tensor/tensor.h>
#include <madness/tensor/clapack.h>
#include <madness/tensor/tensor_lapack.h>
#include <list>
#include <random>
#include<array>
 
namespace madness {
#ifdef BENCH
    static double time_[30];
#endif
 
 
    template <class T>
    class GenTensor;
    template <class T>
    class GenTensor;
    template <class T>
    class SliceLowRankTensor;
 
    /**
     * A SRConf handles all the configurations in a Separated Representation.
     */
 
    template <typename T>
    class SRConf : public BaseTensor {
        friend class GenTensor<T>;
        friend class SliceLowRankTensor<T>;
 
        /// the scalar type of T
        typedef typename Tensor<T>::scalar_type scalar_type;
        typedef typename TensorTypeData<T>::float_scalar_type float_scalar_type;
 
    public:
#ifdef BENCH
        static double& time(int i) { return time_[i]; }
#endif
 
        typedef Tensor<T> tensorT;
 
        /// check orthonormality at low rank additions at every xth addition
        static const long check_orthonormality_frequency = 100000;
 
        /// for each configuration the weight; length should be r
        Tensor<typename Tensor<T>::scalar_type> weights_;
 
        /// for each (physical) dimension one Tensor of (logical) dimension (r,k)
        /// for vectors or (r,kprime,k) for operators
        std::array<Tensor<T>, 2> vector_;
 
        /// separation dimensions: A(n,m) -> A(r,n) B(r,m), with n={k1,k2},m={k3,k4,k5..) multi-indices
        long nci_left = -1;
 
        /// Slice containing the actual data in each vector, ignoring "empty" configurations;
        /// will maintain contiguity of the data.
        std::vector<Slice> s0, s1;
 
    public:
        /// return the index of the last singular vector/value to meet the threshold
        ///        (returns -1 if all meet threshold, i.e. || A ||_2 < threshold)
        /// given a matrix A in SVD form, truncate the singular values such that the
        /// accuracy threshold is still met.
        /// @param[in]  thresh  the threshold eps: || A - A(truncated) || < eps
        /// @param[in]  rank    the number of singular values in w
        /// @param[in]  w       the weights/singular values of A
        /// @return     i       the index of s_max to contribute: w(Slice(0,i)); i.e. inclusive!
        static int max_sigma(const double& thresh, const long& rank, const Tensor<double>& w) {
            if (thresh < 0.0) return rank - 1;
            // find the maximal singular value that's supposed to contribute
            // singular values are ordered (largest first)
            double residual = 0.0;
            long i;
            for (i = rank - 1; i >= 0; i--) {
                residual += w(i) * w(i);
                if (residual > thresh * thresh) break;
            }
            return i;
        }
 
        /// default ctor
        SRConf() : nci_left(-1) {};
 
        SRConf(const long& ndim, const long* dimensions,
               const long nci) : nci_left(nci) {
            BaseTensor::set_dims_and_size(ndim, dimensions);
            if (nci_left < 0) nci_left = ndim / 2; // integer division
            make_structure();
        }
 
        SRConf(const long& ndim, const std::array<long,TENSOR_MAXDIM>& dimensions,
               const long nci) : SRConf(ndim, dimensions.data(), nci) {}
 
        /// copy ctor (tested); shallow copy
        SRConf(const SRConf& rhs) = default;
 
        /// ctor with provided weights and effective vectors; shallow copy
        SRConf(const Tensor<double>& weights, const std::vector<Tensor<T>>& vectors,
               const long& ndim, const long& dims, const long nci)
            : SRConf(ndim, dims, nci) {
            MADNESS_ASSERT(vectors.size()==2);
            set_vectors_and_weights(weights, vectors[0], vectors[1]);
            make_structure();
            MADNESS_ASSERT(has_structure());
            MADNESS_ASSERT(check_dimensions());
        }
 
        /// explicit ctor with two vectors (aka SVD), shallow
        SRConf(const Tensor<double>& weights, const tensorT& vector1, const tensorT& vector2,
               const long& ndim, const long* dims, const long nci)
            : SRConf(ndim, dims, nci) {
            set_vectors_and_weights(weights, vector1, vector2);
            make_structure();
            MADNESS_ASSERT(check_dimensions());
        }
 
        /// assignment operator (tested), shallow copy of vectors
        SRConf& operator=(const SRConf& rhs) {
            // check for self-assignment
            if (&rhs == this) return *this;
            if (rhs.has_no_data()) {
                clear();
                return *this;
            }
 
            // these always hold
            nci_left = rhs.nci_left;
            BaseTensor::set_dims_and_size(rhs.ndim(), rhs.dims());
            s0 = rhs.s0;
            s1 = rhs.s1;
 
            if (rhs.rank() == 0) {
                // construct empty vector
                make_empty_vectors_and_weights(0);
                make_structure();
            }
            else {
                // assign vectors; shallow copy
                for (unsigned int i = 0; i < rhs.vector_.size(); i++) {
                    vector_[i] = rhs.vector_[i];
                }
 
                // shallow copy
                weights_ = (rhs.weights_);
            }
            MADNESS_ASSERT(has_structure());
            return *this;
        }
 
        /// assign a number to this;
        SRConf& operator=(const T& number) {
            // rank will be one
            this->make_empty_vectors_and_weights(1);
            vector_[0] = number;
            vector_[1] = 1.0;
            weights_(0l) = 1.0;
            make_structure();
            normalize();
            return *this;
        }
 
 
        void set_size_and_dim(long ndim, long k) {
            std::array<long,TENSOR_MAXDIM> dims;
            dims.fill(k);
            BaseTensor::set_dims_and_size(ndim, dims.data());
        }
 
        /// deduce the dimensions of the left and right singular vectors from the tensor dimensions
        std::array<std::array<long,TENSOR_MAXDIM>, 2> make_vector_dimensions(const long rank) const {
            std::array<std::array<long,TENSOR_MAXDIM>, 2> dimensions;
            for (int i = 0; i < nci_left; ++i) dimensions[0][i + 1] = this->dim(i);
            for (int i = nci_left; i < ndim(); ++i) dimensions[1][i + 1 - nci_left] = this->dim(i);
 
            dimensions[0][0] = rank;
            dimensions[1][0] = rank;
            return dimensions;
        }
 
        void make_empty_vectors_and_weights(const long rank) {
            auto dimensions = make_vector_dimensions(rank);
            weights_ = Tensor<float_scalar_type>(rank);
            vector_[0] = Tensor<T>(nci_left + 1, dimensions[0].data());
            vector_[1] = Tensor<T>(ndim() - nci_left + 1, dimensions[1].data());
            make_structure();
        }
 
        void set_vectors_and_weights(const Tensor<typename Tensor<T>::scalar_type>& weights,
                                     const Tensor<T>& vector1, const Tensor<T>& vector2) {
            weights_ = weights;
            vector_ = {vector1, vector2};
            make_structure();
        }
 
        void clear() {
            weights_.clear();
            vector_[0].clear();
            vector_[1].clear();
            _size = 0;
            _ndim = -1;
        }
 
        /// return some of the terms of the SRConf (start,..,end), inclusively
        /// shallow copy
        const SRConf get_configs(const int& start, const int& end) const {
            MADNESS_ASSERT((start>=0) and (end<=rank()));
 
            Slice s(start, end);
            tensorT v0, v1;
 
            // slice vectors
            v0 = flat_vector(0)(s, _);
            v1 = flat_vector(1)(s, _);
 
            SRConf<T> result(ndim(), dims(), nci_left);
            result.set_vectors_and_weights(weights_(s), v0, v1);
            MADNESS_ASSERT(result.has_structure());
            return result;
        }
 
        /// dtor
        ~SRConf() {}
 
        template <typename Archive>
        void serialize(Archive& ar) {
            ar & weights_ & vector_[0] & vector_[1] & nci_left & _ndim & _size
                & _id & archive::wrap(_dim,TENSOR_MAXDIM) & s0 & s1;
            //                  make_slices();
            MADNESS_ASSERT(has_structure());
        }
 
        /// does this have any data?
        bool has_data() const {
            return (size() != 0);
        }
 
        /// does this have any data?
        bool has_no_data() const { return !has_data(); }
 
        /// return a Slice that corresponds the that part of vector_ that holds coefficients
        const std::vector<Slice>& c0(const int idim) const {
            if (idim == 0) return s0;
            else if (idim == 1) return s1;
            else {
                MADNESS_EXCEPTION("invalid idim in SRConf::idim", 1);
            }
        }
 
        /// append configurations of rhs to this
 
        /// simplified version of inplace_add for flattened configurations
        /// *this += fac*rhs
        void append(const SRConf<T>& rhs, const double fac = 1.0) {
            // fast return if possible
            if (rhs.has_no_data() or rhs.rank() == 0) return;
            if (this->has_no_data() or rank() == 0) {
                *this = copy(rhs);
                this->scale(fac);
                return;
            }
 
            auto dimensions = make_vector_dimensions(rank() + rhs.rank());
            Tensor<float_scalar_type> weights(rank() + rhs.rank());
            Tensor<T> vector0(nci_left + 1, dimensions[0].data());
            Tensor<T> vector1(ndim() - nci_left + 1, dimensions[1].data());
 
            // assign weights
            weights(Slice(0, rank() - 1)) = weights_(Slice(0, rank() - 1));
            weights(Slice(rank(), rank() + rhs.rank() - 1)) = rhs.weights_(Slice(0, rhs.rank() - 1)) * fac;
 
            vector0(c0(0)) = vector_[0](c0(0));
            vector1(c0(1)) = vector_[1](c0(1));
 
            auto s00 = s0;
            auto s10 = s1;
            s00[0] = Slice(rank(), rank() + rhs.rank() - 1);
            s10[0] = Slice(rank(), rank() + rhs.rank() - 1);
 
            vector0(s00) = rhs.vector_[0](rhs.c0(0));
            vector1(s10) = rhs.vector_[1](rhs.c0(1));
 
            std::swap(weights, weights_);
            std::swap(vector0, vector_[0]);
            std::swap(vector1, vector_[1]);
 
            make_slices();
            MADNESS_ASSERT(has_structure());
        }
 
        void append(const SRConf<T>& rhs, const double_complex fac = 1.0) {
            MADNESS_EXCEPTION("no complex in SRConf", 1);
        }
 
    public:
        /// add two orthonormal configurations, yielding an optimal SVD decomposition
        void add_SVD(const SRConf<T>& rhs, const double& thresh) {
#ifdef BENCH
            double cpu0 = wall_time();
#endif
            if (rhs.has_no_data() or rhs.rank() == 0) return;
            if (has_no_data() or rank() == 0) {
                *this = rhs;
                return;
            }
 
            if (check_orthonormality_frequency > 0) {
                // get random number
                std::random_device rd;
 
                // 2. Initialize the Mersenne Twister engine with the seed
                std::mt19937 gen(rd());
 
                // 3. Define the distribution range [inclusive lower bound, inclusive upper bound]
                std::uniform_int_distribution<int> distrib(1, check_orthonormality_frequency);
 
                // 4. Generate random numbers
                long iran = distrib(gen);
                if (iran==1) {      // any number will do here..
                    print("checking orthonormality in add_SVD");
                    int isuccess=0;
                    if (not check_vector_orthonormality(0)) isuccess++;
                    if (not check_vector_orthonormality(1)) isuccess++;
                    if (not rhs.check_vector_orthonormality(0)) isuccess++;
                    if (not rhs.check_vector_orthonormality(1)) isuccess++;
                    if (isuccess>0) {
                        print("vector in add_svd not orthogonal");
                    }
                }
            }
 
            //            Tensor<T> x1=flat_vector(0);
            //            Tensor<T> x2=flat_vector(1);
            Tensor<T> x1 = vector_[0].reshape(rank(), kVec(0));
            Tensor<T> x2 = vector_[1].reshape(rank(), kVec(1));
            ortho5(x1, x2, weights_,
                   rhs.flat_vector(0), rhs.flat_vector(1), rhs.weights_, thresh);
            std::swap(x1, vector_[0]);
            std::swap(x2, vector_[1]);
 
            make_structure();
            make_slices();
            MADNESS_ASSERT(has_structure());
#ifdef BENCH
            double cpu1 = wall_time();
            time(25) += cpu1 - cpu0;
#endif
        }
 
    protected:
        /// alpha * this(lhs_s) + beta * rhs(rhs_s)
 
        /// bounds checking should have been performed by caller
        /// s denotes where in lhs the new contribution from rhs will be inserted
        void inplace_add(const SRConf<T>& rhs, std::array<Slice,TENSOR_MAXDIM> lhs_s,
                         std::array<Slice,TENSOR_MAXDIM> rhs_s, const double alpha, const double beta) {
            // cannot scale a slice only...
            MADNESS_ASSERT(alpha==1.0);
 
            // fast return if possible; no fast return for this.rank()==0
            // since we might work with slices!
            if (rhs.has_no_data() or rhs.rank() == 0) return;
 
            // prepare the vectors
            SRConf<T> result(ndim(), dims(), nci_left);
            result.make_empty_vectors_and_weights(rank() + rhs.rank());
 
            // insert lhs into result
            if (rank() > 0) {
                result.vector_[0](s0) = vector_[0];
                result.vector_[1](s1) = vector_[1];
                result.weights_(Slice(0, rank() - 1)) = weights_;
            }
            // insert rhs into result
            {
                auto [sr0,sr1] = rhs.make_slices(rhs_s);
                auto [sl0,sl1] = make_slices(lhs_s);
                sl0[0] = Slice(rank(), result.rank() - 1);
                sl1[0] = Slice(rank(), result.rank() - 1);
                result.vector_[0](sl0) = rhs.vector_[0](sr0);
                result.vector_[1](sl1) = rhs.vector_[1](sr1);
                result.weights_(Slice(rank(), result.rank() - 1)) = rhs.weights_ * beta;
            }
            std::swap(*this, result);
            MADNESS_ASSERT(has_structure());
        }
 
        /// deep copy of rhs, shrink
        friend SRConf<T> copy(const SRConf<T>& rhs) {
            if (rhs.has_no_data()) return SRConf<T>();
 
            SRConf<T> result(rhs.ndim(), rhs.dims(), rhs.nci_left);
 
            // if rhs is non-existent simply construct a new SRConf
            if (rhs.has_data() and rhs.rank() > 0) {
                result.set_vectors_and_weights(copy(rhs.weights_(Slice(0, rhs.rank() - 1))),
                                               copy(rhs.vector_[0](rhs.c0(0))), copy(rhs.vector_[1](rhs.c0(1))));
            }
 
            return result;
        }
 
    public:
        /// return a slice of this (deep copy)
        SRConf<T> copy_slice(const std::array<Slice,TENSOR_MAXDIM>& s) const {
            std::array<long,TENSOR_MAXDIM> k;
            for (int i = 0; i < s.size(); ++i) {
                if (s[i].end == -1) k[i] = dim(i);
                else k[i] = s[i].end - s[i].start + 1;
            }
            SRConf<T> result(ndim(), k, nci_left);
 
            // fast return if possible
            if (this->has_no_data() or rank() == 0) return result;
 
            auto [s00,s11] = make_slices(s);
            Tensor<T> vector0 = copy(vector_[0](s00));
            Tensor<T> vector1 = copy(vector_[1](s11));
 
            Tensor<double> weights = copy(this->weights_(Slice(0, rank() - 1)));
            result.set_vectors_and_weights(weights, vector0, vector1);
 
            return result;
        }
 
        /// perform elementwise Hadamard product
        SRConf<T>& emul(const SRConf<T>& other) {
            // consistency check
            MADNESS_ASSERT(compatible(*this,other));
 
            long finalrank = this->rank() * other.rank();
 
            SRConf<T> result(ndim(), dims(), nci_left); // empty tensor
 
            if ((this->rank() == 0) or (other.rank() == 0)) {
                ; // pass
            }
            else {
                result.make_empty_vectors_and_weights(finalrank);
                result.weights_ = outer(weights_, other.weights_).flat();
 
                // left vector
                for (int i = 0; i < 2; ++i) {
                    Tensor<T> a1 = flat_vector(i);
                    Tensor<T> b1 = other.flat_vector(i);
                    Tensor<T> r1 = result.vector_[i].reshape(finalrank, kVec(i));
                    Tensor<T> tmp(finalrank, 1);
                    for (int k = 0; k < a1.dim(1); ++k) {
                        //                      r1(_,Slice(k,k))=outer(a1(_,k),b1(_,k)).reshape(finalrank,1);
                        outer_result(a1(_, k), b1(_, k), tmp);
                        r1(_, Slice(k, k)) = tmp;
                    }
                }
            }
            result.make_structure();
            result.normalize();
 
            *this = result;
            return *this;
        }
 
    protected:
        /// redo the Slices for getting direct access to the configurations
        void make_slices() {
            s0 = std::vector<Slice>(nci_left + 1, _); // first dim is the rank
            s1 = std::vector<Slice>(ndim() - nci_left + 1, _); // first dim is the rank
 
            s0[0] = Slice(0, rank() - 1);
            s1[0] = Slice(0, rank() - 1);
        }
 
        std::array<std::vector<Slice>, 2> make_slices(const std::array<Slice,TENSOR_MAXDIM>& s) const {
            std::vector<Slice> s00(nci_left + 1, _); // first dim is the rank
            std::vector<Slice> s10(ndim() - nci_left + 1, _); // first dim is the rank
 
            s00[0] = Slice(0, rank() - 1);
            s10[0] = Slice(0, rank() - 1);
 
            for (int idim = 0; idim < ndim(); ++idim) {
                if (idim < nci_left) s00[idim + 1] = s[idim];
                if (idim >= nci_left) s10[idim - nci_left + 1] = s[idim];
            }
 
            std::array<std::vector<Slice>, 2> result = {s00, s10};
            return result;
        }
 
        void make_structure(bool force = false) {
            auto dimensions = make_vector_dimensions(rank());
            vector_[0] = vector_[0].reshape(nci_left + 1, &dimensions[0][0]);
            vector_[1] = vector_[1].reshape(ndim() - nci_left + 1, &dimensions[1][0]);
 
            this->make_slices();
        }
 
    public:
        /// return reference to one of the vectors F
        Tensor<T>& ref_vector(const unsigned int& idim) {
            return vector_[idim];
        }
 
        /// return reference to one of the vectors F
        const Tensor<T>& ref_vector(const unsigned int& idim) const {
            return vector_[idim];
        }
 
 
        long kVec(const int idim) const {
            return vector_[idim].size() / vector_[idim].dim(0);
        }
 
        /// return shallow copy of a slice of one of the vectors, flattened to (r,kVec)
        const Tensor<T> flat_vector(const unsigned int& idim) const {
            MADNESS_ASSERT(rank()>0);
            //          const Tensor<T> result=vector_[idim](c0(idim)).reshape(rank(),kVec(idim));
            const Tensor<T> result = vector_[idim].reshape(rank(), kVec(idim));
            return result;
        }
 
        /// return shallow copy of a slice of one of the vectors, flattened to (r,kVec)
        Tensor<T> flat_vector(const unsigned int& idim) {
            MADNESS_ASSERT(rank()>0);
            //          Tensor<T> result=vector_[idim](c0(idim)).reshape(rank(),kVec(idim));
            Tensor<T> result = vector_[idim].reshape(rank(), kVec(idim));
            return result;
        }
 
    protected:
        //      /// fill this SRConf with 1 flattened random configurations (tested)
        //      void fillWithRandom(const long& rank=1) {
        //
        //
        //          // assign; note that Slice(0,_) is inclusive
        //          weights_=Tensor<double>(rank);
        //          weights_=1.0;
        //
        //          for (unsigned int idim=0; idim<this->dim_eff(); idim++) {
        //              vector_[idim]=Tensor<T>(rank,this->kVec());
        //              vector_[idim].fillrandom();
        //          }
        //
        //          this->normalize();
        //          for (unsigned int idim=0; idim<this->dim_eff(); idim++) {
        //              vector_[idim].scale(madness::RandomValue<T>()*scalar_type(10.0));
        //          }
        //          weights_(Slice(0,this->rank()-1)).fillrandom().scale(10.0);
        //          make_slices();
        //            MADNESS_ASSERT(has_structure());
        //      }
 
        /// normalize the vectors (tested)
        void normalize() {
            if (rank() == 0) return;
            MADNESS_ASSERT(has_structure());
 
            // for convenience
            const unsigned int rank = this->rank();
            //          std::vector<Slice> s(2,_);
 
            // we calculate the norm sum_i < F^r_i | F^r_i > for each dimension for each r
 
            // loop over all configurations
            for (unsigned int r = 0; r < rank; r++) {
                // loop over all dimensions
                for (unsigned int idim = 0; idim < 2; idim++) {
                    std::vector<Slice> s(dim_per_vector(idim) + 1, _);
                    s[0] = Slice(r, r);
 
                    //                  const double norm=this->ref_vector(idim)(s).normf();
                    const double norm = this->vector_[idim](s).normf();
                    const double fac = norm;
                    double oofac = 1.0 / fac;
                    if (fac < 1.e-13) oofac = 0.0;
                    weights_(r) *= fac;
                    vector_[idim](s).scale(oofac);
                }
            }
            MADNESS_ASSERT(has_structure());
        }
 
        bool check_dimensions() const {
            bool correct = true;
            if (vector_[0].dim(0) != rank()) return false;
            if (vector_[0].ndim() + vector_[1].ndim() != ndim() + 2) return false;
            if (vector_[0].ndim() != nci_left + 1) return false;
            if (vector_[1].ndim() != ndim() - nci_left + 1) return false;
 
            return correct;
        }
 
        /// check if the terms are orthogonal
 
        /// @param[in] ivec which vector (0 or 1) to check
        bool check_vector_orthonormality(const int ivec) const {
            // fast return if possible
            if (rank() == 0) return true;
 
            const tensorT t1 = ref_vector(ivec)(c0(ivec)).reshape(rank(), kVec(ivec));
            tensorT S = inner(t1, t1, 1, 1);
            for (int i = 0; i < S.dim(0); i++) S(i, i) -= 1.0;
 
            // error per matrix element
            double norm = S.normf();
            double small = sqrt(norm * norm / S.size());
            return (small < 1.e-10);
        }
 
        /// return if this has only one additional dimension (apart from rank)
        bool is_flat() const {
            return (vector_[0].ndim() == 2);
        }
 
    public:
        /// return if this has a tensor structure (has not been flattened)
        bool has_structure() const {
            if (vector_.size() == 2) {
                auto vector_dimensions = make_vector_dimensions(vector_[0].dim(0));
                for (int i = 0; i < vector_[0].ndim(); ++i) {
                    if (vector_dimensions[0][i] != vector_[0].dim(i)) return false;
                }
                for (int i = 0; i < vector_[1].ndim(); ++i) {
                    if (vector_dimensions[1][i] != vector_[1].dim(i)) return false;
                }
            }
            return true;
        }
 
        /// return the logicalrank
        long rank() const { return weights_.size(); };
 
    public:
        /// return the number of physical dimensions
        int dim_per_vector(int idim) const {
            MADNESS_ASSERT(vector_.size()>size_t(idim));
            return vector_[idim].ndim() - 1; // remove dimension for the rank
        }
 
        /// return the weight
        double weights(const unsigned int& i) const { return weights_(i); };
 
        /// reconstruct this to return a full tensor
        Tensor<T> reconstruct() const {
            // fast return if possible
            if (rank() == 0) return Tensor<T>(ndim(), dims(), true);
 
            // include weights in left vector
            Tensor<T> scr = make_left_vector_with_weights();
 
            Tensor<T> result = inner(conj(scr), flat_vector(1), 0, 0);
            return result.reshape(ndim(), dims());
        }
 
    protected:
        Tensor<T> make_left_vector_with_weights() const {
            return make_vector_with_weights(0);
        }
 
    public:
        Tensor<T> make_vector_with_weights(const int dim) const {
            Tensor<T> v = copy(vector_[dim].reshape(rank(), vector_[dim].size() / rank()));
            for (unsigned int r = 0; r < rank(); r++) v(r, _) *= weights(r);
            v = v.reshape(vector_[dim].ndim(), vector_[dim].dims());
            return v;
        }
 
        /// return flat (r,i) view of the tensor with the weights multiplied in
 
        /// return a(r,i) = vec(dim)(r,i) * w(r)
        Tensor<T> flat_vector_with_weights(const int dim) const {
            return make_vector_with_weights(dim).reshape(rank(), vector_[dim].size() / rank());
        }
 
    protected:
        /// return the number of coefficients
        unsigned int nCoeff() const {
            return vector_[0].size() + vector_[1].size() + weights_.size();
        };
 
        /// return the real size of this
        size_t real_size() const {
            size_t n = 0;
            for (size_t i = 0; i < vector_.size(); ++i) {
                n += vector_[i].size() * sizeof(T) + sizeof(tensorT);
            }
            n += weights_.size() * sizeof(double) + sizeof(Tensor<double>);
            n += sizeof(*this);
            n += 2 * s1.size() * sizeof(Slice);
            return n;
        }
 
 
        template <typename Q>
        typename std::enable_if<(TensorTypeData<T>::iscomplex or TensorTypeData<Q>::iscomplex),
                                TENSOR_RESULT_TYPE(T, Q)>::type
        friend trace(const SRConf<T>& rhs, const SRConf<Q>& lhs) {
            MADNESS_EXCEPTION("no complex trace in srconf.h", 1);
            return T(0.0);
        }
 
        /// calculate the Frobenius inner product (tested)
        template <typename Q>
        typename std::enable_if<!(TensorTypeData<T>::iscomplex or TensorTypeData<Q>::iscomplex),
                                TENSOR_RESULT_TYPE(T, Q)>::type
        friend trace(const SRConf<T>& rhs, const SRConf<Q>& lhs) {
            // fast return if either rank is 0
            if ((lhs.has_no_data()) or (rhs.has_no_data())) return 0.0;
            if ((lhs.rank() == 0) or (rhs.rank() == 0)) return 0.0;
 
            /*
             * the structure of an SRConf is (r,k) or (r,k',k), with
             * r the slowest index; the overlap is therefore simply calculated
             * as the matrix multiplication rhs*lhs^T
             */
 
            // some checks
            MADNESS_ASSERT(rhs.ndim()==lhs.ndim());
            MADNESS_ASSERT(rhs.ndim()>0);
 
            typedef TENSOR_RESULT_TYPE(T, Q) resultT;
 
            const unsigned int dim_eff = 2;
 
            // get the weight matrix
            Tensor<resultT> weightMatrix = outer(lhs.weights_(Slice(0, lhs.rank() - 1)),
                                                 rhs.weights_(Slice(0, rhs.rank() - 1)));
 
            // calculate the overlap matrices for each dimension at a time
            for (unsigned int idim = 0; idim < dim_eff; idim++) {
                const Tensor<T> lhs2 = lhs.flat_vector(idim);
                const Tensor<Q> rhs2 = rhs.flat_vector(idim);
                Tensor<resultT> ovlp(lhs.rank(), rhs.rank());
                inner_result(lhs2, rhs2, -1, -1, ovlp);
 
                // multiply all overlap matrices with the weight matrix
                weightMatrix.emul(ovlp);
            }
 
            //  return weightMatrix;
            const TENSOR_RESULT_TYPE(T, Q) overlap = weightMatrix.sum();
            return overlap;
        }
 
        /// calculate the Frobenius norm, if this is in SVD form
        typename TensorTypeData<T>::float_scalar_type svd_normf() const {
            if (has_no_data() or rank() == 0) return 0.0;
            return weights_(Slice(0, rank() - 1)).normf();
        }
 
        /// calculate the Frobenius norm
        typename TensorTypeData<T>::float_scalar_type normf() const {
            // fast return if possible
            if (has_no_data() or rank() == 0) return 0.0;
 
            // some checks
            MADNESS_ASSERT(ndim()>0);
            MADNESS_ASSERT(not TensorTypeData<T>::iscomplex);
 
            // get the weight matrix
            Tensor<T> weightMatrix = outer(weights_(Slice(0, rank() - 1)), weights_(Slice(0, rank() - 1)));
 
            // calculate the overlap matrices for each dimension at a time
            for (unsigned int idim = 0; idim < 2; idim++) {
                const Tensor<T> vec = flat_vector(idim);
                Tensor<T> ovlp(rank(), rank());
                inner_result(vec, vec, -1, -1, ovlp);
 
                // multiply all overlap matrices with the weight matrix
                weightMatrix.emul(ovlp);
            }
 
            typedef typename TensorTypeData<T>::float_scalar_type resultT;
            const resultT overlap = std::abs(weightMatrix.sum());
            return sqrt(overlap);
        }
 
    protected:
        /// scale this by a number
        void scale(const double& fac) { weights_.scale(fac); };
 
        void scale(const double_complex& fac) {
            MADNESS_EXCEPTION("no complex scaling in SRConf", 1);
        }
 
        /// check compatibility
        friend bool compatible(const SRConf& lhs, const SRConf& rhs) {
            return (lhs.conforms(&rhs) and (lhs.dim_per_vector(0) == rhs.dim_per_vector(0)));
        }
 
        /// \code
        ///     result(i,j,k,...) <-- sum(i',j', k',...) t(i',j',k',...)  c(i',i) c(j',j) c(k',k) ...
        /// \endcode
        ///
        /// The input dimensions of \c t must all be the same .
        SRConf<T> transform(const Tensor<T>& c) const {
            // fast return if possible
            if (this->has_no_data() or rank() == 0) return SRConf<T>(ndim(), dims(), nci_left);
 
            // transpose both singular vectors from U(r,i,j,k) to U(i,j,k,r)
            // and run the contraction as in tensor: U(i,j,k,r) t(i,i') = U(j,k,r,i')
            // and so on, yielding U(r,i',j',k')
            Tensor<T> left = copy(vector_[0].cycledim(-1, 0, -1));
            Tensor<T> right = copy(vector_[1].cycledim(-1, 0, -1));
 
            for (long i = 0; i < nci_left; ++i) left = inner(left, c, 0, 0);
            for (long i = nci_left; i < ndim(); ++i) right = inner(right, c, 0, 0);
 
            SRConf<T> result(ndim(), dims(), this->nci_left);
            result.set_vectors_and_weights(copy(weights_), left, right);
 
            MADNESS_ASSERT(result.has_structure());
            return result;
        }
 
    public:
        /// \code
        ///     result(i,j,k,...) <-- sum(i',j', k',...) t(i',j',k',...)  c(i',i) c(j',j) c(k',k) ...
        /// \endcode
        ///
        /// The input dimensions of \c t must all be the same .
        template <typename Q>
        SRConf<TENSOR_RESULT_TYPE(T, Q)> general_transform(const Tensor<Q> c[]) const {
            // fast return if possible
            if (this->has_no_data() or rank() == 0) return SRConf<T>(ndim(), dims(), nci_left);
 
            // copying shrinks the vectors to (r,k,k,..)
            MADNESS_ASSERT(this->has_structure());
 
            // transpose both singular vectors from U(r,i,j,k) to U(i,j,k,r)
            // and run the contraction as in tensor: U(i,j,k,r) t(i,i') = U(j,k,r,i')
            // and so on, yielding U(r,i',j',k')
            Tensor<T> left = copy(vector_[0].cycledim(-1, 0, -1));
            Tensor<T> right = copy(vector_[1].cycledim(-1, 0, -1));
 
            for (long i = 0; i < nci_left; ++i) left = inner(left, c[i], 0, 0);
            for (long i = nci_left; i < ndim(); ++i) right = inner(right, c[i], 0, 0);
 
            SRConf<T> result(ndim(), dims(), this->nci_left);
            result.set_vectors_and_weights(copy(weights_), left, right);
 
            MADNESS_ASSERT(result.has_structure());
            return result;
        }
 
        SRConf<T> transform_dir(const Tensor<T>& c, const int& axis) const {
            if (this->has_no_data() or rank() == 0) return SRConf<T>(ndim(), dims(), nci_left);
 
            // copying shrinks the vectors to (r,k,k,..)
            SRConf<T> result = copy(*this);
 
            // only a matrix is allowed for c
            MADNESS_ASSERT(c.ndim()==2);
 
            // make sure this is not flattened
            MADNESS_ASSERT(this->has_structure());
 
            // compute idim for accessing the vector_, and the dimension inside vector_
            // the +1 on jdim for the rank
            const long idim = axis / this->dim_per_vector(0);
            const long jdim = axis % this->dim_per_vector(0) + 1;
 
            // note: no slicing necessary, since we've copied this to result (incl shrinking)
            result.ref_vector(idim) = madness::transform_dir(this->ref_vector(idim), c, jdim);
            MADNESS_ASSERT(result.has_structure());
 
            return result;
        }
    };
 
    /// sophisticated version of ortho2
 
    /// after calling this we will have an optimally rank-reduced representation
    /// with the left and right subspaces being bi-orthogonal and normalized;
    /// outline of the algorithm:
    ///  - canonical orthogonalization of the subspaces (screen for small eigenvalues)
    ///  - SVD of the modified overlap (incorporates the roots of eigenvalues)
    /// operation count is O(kr^2 + r^3)
    ///
    /// @param[in,out]  x normalized left subspace
    /// @param[in,out]  y normalize right subspace
    /// @param[in,out]  weights weights
    /// @param[in]      thresh  truncation threshold
    template <typename T>
    void ortho3(Tensor<T>& x, Tensor<T>& y, Tensor<typename Tensor<T>::scalar_type>& weights, const double& thresh) {
#ifdef BENCH
        double cpu0 = wall_time();
#endif
        typedef Tensor<T> tensorT;
        typedef typename Tensor<T>::scalar_type scalar_type;
 
        const long rank = x.dim(0);
        const double w_max = weights.absmax() * rank; // max Frobenius norm
 
 
        // overlap of 1 and 2
        tensorT S1 = inner(x, x, 1, 1);
        tensorT S2 = inner(y, y, 1, 1); // 0.5 / 2.1
#ifdef BENCH
        double cpu1 = wall_time();
        SRConf<T>::time(1) += (cpu1 - cpu0);
#endif
 
        //      print("norm(S1)",S1.normf());
        //      print("norm(S2)",S2.normf());
 
        // diagonalize
        tensorT U1, U2;
        Tensor<scalar_type> e1, e2;
        syev(S1, U1, e1);
        syev(S2, U2, e2); // 2.3 / 4.0
#ifdef BENCH
        double cpu3 = wall_time();
        SRConf<T>::time(3) += cpu3 - cpu1;
#endif
 
        const scalar_type e1_max = e1.absmax();
        const scalar_type e2_max = e2.absmax();
 
        // fast return if possible
        if ((e1_max * w_max < thresh) or (e2_max * w_max < thresh)) {
            x.clear();
            y.clear();
            weights.clear();
            return;
        }
 
        // remove small negative eigenvalues
        e1.screen(1.e-13);
        e2.screen(1.e-13);
        Tensor<scalar_type> sqrt_e1(rank), sqrt_e2(rank);
 
 
        // shrink U1, U2
        int lo1 = 0;
        int lo2 = 0;
        for (unsigned int r = 0; r < rank; r++) {
            if (e1(r) * w_max < thresh) lo1 = r + 1;
            if (e2(r) * w_max < thresh) lo2 = r + 1;
            sqrt_e1(r) = sqrt(std::abs(e1(r)));
            sqrt_e2(r) = sqrt(std::abs(e2(r)));
        }
 
        U1 = U1(Slice(_), Slice(lo1, -1));
        U2 = U2(Slice(_), Slice(lo2, -1));
        sqrt_e1 = sqrt_e1(Slice(lo1, -1));
        sqrt_e2 = sqrt_e2(Slice(lo2, -1));
        unsigned int rank1 = rank - lo1;
        unsigned int rank2 = rank - lo2; // 0.0 / 0.0
 
 
        MADNESS_ASSERT(sqrt_e1.size()==rank1);
        MADNESS_ASSERT(sqrt_e2.size()==rank2);
#ifdef BENCH
        double cpu4 = wall_time();
        SRConf<T>::time(4) += cpu4 - cpu3;
#endif
 
 
        // set up overlap M; include X+
        tensorT M(rank1, rank2);
        for (unsigned int i = 0; i < rank1; i++) {
            for (unsigned int j = 0; j < rank2; j++) {
                for (unsigned int r = 0; r < rank; r++) {
                    M(i, j) += U1(r, i) * sqrt_e1(i) * weights(r) * U2(r, j) * sqrt_e2(j);
                    //                      M(i,j)+=U1(r,i)*weights(r)*U2(r,j);
                }
            }
        }
 
 
        // include X-
        for (unsigned int r = 0; r < rank1; r++) {
            scalar_type fac = 1.0 / sqrt_e1(r);
            for (unsigned int t = 0; t < rank; t++) {
                U1(t, r) *= fac;
                //              if (sqrt_e1(r)<thresh) throw;
            }
        }
 
        for (unsigned int r = 0; r < rank2; r++) {
            scalar_type fac = 1.0 / sqrt_e2(r);
            for (unsigned int t = 0; t < rank; t++) {
                U2(t, r) *= fac;
                //              if (sqrt_e2(r)<thresh) throw;
            }
        } // 0.2 / 0.6
#ifdef BENCH
        double cpu5 = wall_time();
        SRConf<T>::time(5) += cpu5 - cpu4;
#endif
 
        // decompose M
        tensorT Up, VTp;
        Tensor<scalar_type> Sp;
        svd(M, Up, Sp, VTp); // 1.5 / 3.0
#ifdef BENCH
        double cpu6 = wall_time();
        SRConf<T>::time(6) += cpu6 - cpu5;
#endif
 
        // make transformation matrices
        Up = inner(Up, U1, 0, 1);
        VTp = inner(VTp, U2, 1, 1);
 
 
        //      // find the maximal singular value that's supposed to contribute
        //      // singular values are ordered (largest first)
        //      double residual=0.0;
        //      long i;
        //      for (i=Sp.dim(0)-1; i>=0; i--) {
        //          residual+=Sp(i)*Sp(i);
        //          if (residual>thresh*thresh) break;
        //      }
        long i = SRConf<T>::max_sigma(thresh, Sp.dim(0), Sp);
 
#ifdef BENCH
        double cpu7 = wall_time();
        SRConf<T>::time(7) += cpu7 - cpu6;
#endif
 
        //      i=std::min(i,long(0));
 
        Up = Up(Slice(0, i), Slice(_));
        VTp = VTp(Slice(0, i), Slice(_));
 
 
        // convert SVD output to our convention
        if (i >= 0) {
            // transform 1 and 2
            x = inner(Up, x, 1, 0);
            y = inner(VTp, y, 1, 0); // 0.5 / 2.5
            weights = Sp(Slice(0, i));
        }
        else {
            x.clear();
            y.clear();
            weights.clear();
        }
#ifdef BENCH
        double cpu8 = wall_time();
        SRConf<T>::time(8) += cpu8 - cpu7;
        SRConf<T>::time(0) += cpu8 - cpu0;
#endif
        return;
    }
 
    /// specialized version of ortho3
 
    /// does the same as ortho3, but takes two bi-orthonormal configs as input
    /// and saves on the inner product. Result will be written onto the first config
    ///
    /// @param[in,out]  x1  left subspace, will hold the result on exit
    /// @param[in,out]  y1  right subspace, will hold the result on exit
    /// @param[in]      x2  left subspace, will be accumulated onto x1
    /// @param[in]      y2  right subspace, will be accumulated onto y1
    template <typename T>
    void ortho5(Tensor<T>& x1, Tensor<T>& y1, Tensor<typename Tensor<T>::scalar_type>& w1,
                const Tensor<T>& x2, const Tensor<T>& y2, const Tensor<typename Tensor<T>::scalar_type>& w2,
                const double& thresh) {
#ifdef BENCH
        double cpu0 = wall_time();
#endif
        typedef Tensor<T> tensorT;
 
        const long rank1 = x1.dim(0);
        const long rank2 = x2.dim(0);
        const long rank = rank1 + rank2;
 
        // for convenience: blocks of the matrices
        const Slice s0(0, rank1 - 1), s1(rank1, rank - 1);
 
        const double w_max = std::max(w1.absmax(), w2.absmax());
        const double norm_max = w_max * rank; // max Frobenius norm
 
        // the overlap between 1 and 2;
        // the overlap of 1 and 1, and 2 and 2 is assumed to be the identity matrix
        tensorT Sx12 = inner(x1, x2, 1, 1);
        tensorT Sy12 = inner(y1, y2, 1, 1);
#ifdef BENCH
        double cpu1 = wall_time();
        SRConf<T>::time(11) += cpu1 - cpu0;
#endif
 
        tensorT Sx(rank, rank);
        tensorT Sy(rank, rank);
 
        // the identity matrix (half of it)
        for (long i = 0; i < rank; i++) {
            Sx(i, i) = 0.5;
            Sy(i, i) = 0.5;
        }
        Sx(s0, s1) = Sx12;
        Sy(s0, s1) = Sy12;
        Sx += transpose(Sx);
        Sy += transpose(Sy);
 
        // overlap of 1 and 2
        //      tensorT S1=inner(x,x,1,1);
        //      tensorT S2=inner(y,y,1,1);  // 0.5 / 2.1
#ifdef BENCH
        double cpu2 = wall_time();
        SRConf<T>::time(12) += cpu2 - cpu1;
#endif
 
        // diagonalize
        tensorT U1, U2;
        Tensor<typename Tensor<T>::scalar_type> e1, e2;
        syev(Sx, U1, e1);
        syev(Sy, U2, e2); // 2.3 / 4.0
#ifdef BENCH
        double cpu3 = wall_time();
        SRConf<T>::time(13) += cpu3 - cpu2;
#endif
 
        //      print("norm(Sx)",Sx.normf());
        //      print("norm(Sy)",Sy.normf());
 
        const double e1_max = e1.absmax();
        const double e2_max = e2.absmax();
 
        // fast return if possible
        if ((e1_max * norm_max < thresh) or (e2_max * norm_max < thresh)) {
            x1.clear();
            y1.clear();
            w1.clear();
            return;
        }
 
        // remove small negative eigenvalues
        e1.screen(1.e-13);
        e2.screen(1.e-13);
        Tensor<double> sqrt_e1(rank), sqrt_e2(rank);
 
 
        // shrink U1, U2
        int lo1 = 0;
        int lo2 = 0;
        for (unsigned int r = 0; r < rank; r++) {
            if (e1(r) < thresh / norm_max) lo1 = r + 1;
            else sqrt_e1(r) = sqrt(std::abs(e1(r)));
            if (e2(r) < thresh / norm_max) lo2 = r + 1;
            else sqrt_e2(r) = sqrt(std::abs(e2(r)));
        }
 
        U1 = U1(Slice(_), Slice(lo1, -1));
        U2 = U2(Slice(_), Slice(lo2, -1));
        sqrt_e1 = sqrt_e1(Slice(lo1, -1));
        sqrt_e2 = sqrt_e2(Slice(lo2, -1));
        unsigned int rank_x = rank - lo1;
        unsigned int rank_y = rank - lo2; // 0.0 / 0.0
 
 
        //      MADNESS_ASSERT(sqrt_e1.size()==rank_x);
        //      MADNESS_ASSERT(sqrt_e2.size()==rank_y);
 
        // set up overlap M; include X+
 
        //      for (unsigned int i=0; i<rank_x; ++i) U(i,_)*=sqrt_e1(i);
        //      for (unsigned int i=0; i<rank_y; ++i) U(i,_)*=sqrt_e2(i);
 
        tensorT UU1 = copy(U1);
        for (unsigned int i = 0; i < rank1; ++i) UU1(i, _) *= w1(i);
        for (unsigned int i = rank1; i < rank; ++i) UU1(i, _) *= w2(i - rank1);
 
        tensorT M = inner(UU1, U2, 0, 0);
        tensorT ee = outer(sqrt_e1, sqrt_e2);
        M.emul(ee);
 
        // include X-
        for (unsigned int r = 0; r < rank_x; r++) {
            double fac = 1.0 / sqrt_e1(r);
            U1(_, r) *= fac;
        }
 
        for (unsigned int r = 0; r < rank_y; r++) {
            double fac = 1.0 / sqrt_e2(r);
            U2(_, r) *= fac;
        } // 0.2 / 0.6
#ifdef BENCH
        double cpu4 = wall_time();
        SRConf<T>::time(14) += cpu4 - cpu3;
#endif
 
 
        // decompose M
        tensorT Up, VTp;
        Tensor<typename Tensor<T>::scalar_type> Sp;
        svd(M, Up, Sp, VTp); // 1.5 / 3.0
#ifdef BENCH
        double cpu5 = wall_time();
        SRConf<T>::time(15) += cpu5 - cpu4;
#endif
 
        // make transformation matrices
        Up = inner(Up, U1, 0, 1);
        VTp = inner(VTp, U2, 1, 1);
#ifdef BENCH
        double cpu6 = wall_time();
        SRConf<T>::time(16) += cpu6 - cpu5;
#endif
 
        // find the maximal singular value that's supposed to contribute
        // singular values are ordered (largest first)
        double residual = 0.0;
        long i;
        for (i = Sp.dim(0) - 1; i >= 0; i--) {
            residual += Sp(i) * Sp(i);
            if (residual > thresh * thresh) break;
        }
 
        // convert SVD output to our convention
        if (i >= 0) {
            // make it contiguous
            tensorT Up1 = transpose(Up(Slice(0, i), s0));
            tensorT Up2 = transpose(Up(Slice(0, i), s1));
            tensorT VTp1 = transpose(VTp(Slice(0, i), s0));
            tensorT VTp2 = transpose(VTp(Slice(0, i), s1));
 
            // transform 1 and 2
            x1 = inner(Up1, x1, 0, 0);
            inner_result(Up2, x2, 0, 0, x1);
            y1 = inner(VTp1, y1, 0, 0);
            inner_result(VTp2, y2, 0, 0, y1);
            w1 = Sp(Slice(0, i));
        }
        else {
            x1.clear();
            y1.clear();
            w1.clear();
        }
#ifdef BENCH
        double cpu7 = wall_time();
        SRConf<T>::time(17) += cpu7 - cpu6;
        SRConf<T>::time(10) += cpu7 - cpu0;
#endif
        return;
    }
 
    template <typename T>
    static inline
    std::ostream& operator<<(std::ostream& s, const SRConf<T>& sr) {
        s << "dim_          " << sr.ndim() << "\n";;
        s << "rank_         " << sr.rank_ << "\n";;
        s << "vector_.size()" << sr.vector_.size() << "\n";
        s << "has_data()    " << sr.has_data() << "\n";
        s << "TensorType    " << sr.type() << "\n\n";
        return s;
    }
}
 
#endif /* SRCONF_H_ */