tensor.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$
*/

#ifndef MADNESS_TENSOR_TENSOR_H__INCLUDED
#define MADNESS_TENSOR_TENSOR_H__INCLUDED

#include <madness/madness_config.h>
#include <madness/misc/ran.h>
#include <madness/world/posixmem.h>
#include <madness/world/shared_ptr.h>

#include <complex>
#include <vector>
#include <cmath>
#include <cstdlib>

#include <madness/world/archive.h>
// #include <madness/world/print.h>
// 
// typedef std::complex<float> float_complex;
// typedef std::complex<double> double_complex;
// 
// // These probably have to be included in this order
// #include <madness/tensor/tensor_macros.h>
// #include <madness/tensor/type_data.h>
// #include <madness/tensor/slice.h>
// #include <madness/tensor/vector_factory.h>
#include <madness/tensor/basetensor.h>
#include <madness/tensor/aligned.h>
#include <madness/tensor/mxm.h>
#include <madness/tensor/mtxmq.h>
#include <madness/tensor/tensorexcept.h>
#include <madness/tensor/tensoriter.h>

#ifdef USE_GENTENSOR
#define HAVE_GENTENSOR 1
#else
#define HAVE_GENTENSOR 0
#endif


#ifndef HAVE_STD_ABS_LONG
#ifndef HAVE_STD_LABS
namespace std {
    static long abs(long a) {
        return a>=0 ? a : -a;
    }
}
#else
namespace std {
    static long abs(long a) {
        return std::labs(a);
    }
}
#endif
#endif


namespace madness {
#define IS_ODD(n) ((n)&0x1)
#define IS_UNALIGNED(p) (((unsigned long)(p))&0x7)


    template <typename Q, bool iscomplex>
    struct conditional_conj_struct {
        static Q op(const Q& coeff) {
            return coeff;
        }
    };

    template <typename Q>
    struct conditional_conj_struct<Q,true> {
        static Q op(const Q& coeff) {
            return conj(coeff);
        }
    };

    template <typename Q>
    Q conditional_conj(const Q& coeff) {
        return conditional_conj_struct<Q,TensorTypeData<Q>::iscomplex>::op(coeff);
    }

    namespace detail {
        template <typename T> T mynorm(T t) {
            return t*t;
        }

        template <typename T> T mynorm(std::complex<T> t) {
            return std::norm(t);
        }
    }

    template <class T> class SliceTensor;


        enum TensorType {TT_NONE, TT_FULL, TT_2D};

    static
    inline
    std::ostream& operator << (std::ostream& s, const TensorType& tt) {
        std::string str="confused tensor type";
        if (tt==TT_FULL) str="full rank tensor";
        if (tt==TT_2D) str="low rank tensor 2-way";
        if (tt==TT_NONE) str="no tensor type specified";
        s << str.c_str();
        return s;
    }




    template <class T> class Tensor : public BaseTensor {
        template <class U> friend class SliceTensor;

    protected:
        T* restrict _p;
        std::shared_ptr<T> _shptr;


        void allocate(long nd, const long d[], bool dozero) {
            _id = TensorTypeData<T>::id;
            if (nd < 0) {
                _p = 0;
                _shptr.reset();
                _size = 0;
                _ndim = -1;
                return;
            }

            TENSOR_ASSERT(nd>0 && nd <= TENSOR_MAXDIM,"invalid ndim in new tensor", nd, 0);
            // sanity check ... 2GB in doubles
            for (int i=0; i<nd; ++i) {
                TENSOR_ASSERT(d[i]>=0 && d[i]<268435456, "invalid dimension size in new tensor",d[i],0);
            }
            set_dims_and_size(nd, d);
            if (_size) {
                TENSOR_ASSERT(_size>=0 && _size<268435456, "invalid size in new tensor",_size,0);
                try {
#if HAVE_IBMBGP
#define TENSOR_ALIGNMENT 16
#elif HAVE_IBMBGQ
#define TENSOR_ALIGNMENT 32
#else
#define TENSOR_ALIGNMENT 16
#endif

#ifdef WORLD_GATHER_MEM_STATS
                    _p = new T[size];
                    _shptr = std::shared_ptr<T>(_p);
#else
                    if (posix_memalign((void **) &_p, TENSOR_ALIGNMENT, sizeof(T)*_size)) throw 1;
                    _shptr.reset(_p, &::madness::detail::checked_free<T>);
#endif
                }
                catch (...) {
                    std::printf("new failed nd=%ld type=%ld size=%ld\n", nd, id(), _size);
                    std::printf("  %ld %ld %ld %ld %ld %ld\n",
                                d[0], d[1], d[2], d[3], d[4], d[5]);
                    TENSOR_EXCEPTION("new failed",_size,this);
                }
                //std::printf("allocated %p [%ld]  %ld\n", _p, size, p.use_count());
                if (dozero) {
                    //T zero = 0; for (long i=0; i<_size; ++i) _p[i] = zero;
                    // or
#ifdef HAVE_MEMSET
                    memset((void *) _p, 0, _size*sizeof(T));
#else
                    aligned_zero(_size, _p);
#endif
                }
            }
            else {
                _p = 0;
                _shptr.reset();
            }
        }

        // Free memory and restore default constructor state
        void deallocate() {
            _p = 0;
            _shptr.reset();
            _size = 0;
            _ndim = -1;
        }

    public:
        typedef T type;

        typedef typename TensorTypeData<T>::scalar_type scalar_type;

        typedef typename TensorTypeData<T>::float_scalar_type float_scalar_type;

        Tensor() : _p(0) {
            _id = TensorTypeData<T>::id;
        }


        Tensor(const Tensor<T>& t) {
            _id = TensorTypeData<T>::id;
            *this = t;
        }


        Tensor<T>& operator=(const Tensor<T>& t) {
            if (this != &t) {
                _p = t._p;
                _shptr = t._shptr;
                _size = t._size;
                _ndim = t._ndim;
                for (int i=0; i<TENSOR_MAXDIM; ++i) {
                    _dim[i] = t._dim[i];
                    _stride[i] = t._stride[i];
                }
            }
            return *this;
        }


        template <class Q> operator Tensor<Q>() const { // type conv => deep copy
            Tensor<Q> result = Tensor<Q>(this->_ndim,this->_dim,false);
            BINARY_OPTIMIZED_ITERATOR(Q, result, const T, (*this), *_p0 = (Q)(*_p1));
            return result;
        }



        explicit Tensor(long d0) : _p(0) {
            _dim[0] = d0;
            allocate(1, _dim, true);
        }


        explicit Tensor(long d0, long d1) : _p(0) {
            _dim[0] = d0; _dim[1] = d1;
            allocate(2, _dim, true);
        }


        explicit Tensor(long d0, long d1, long d2) : _p(0) {
            _dim[0] = d0; _dim[1] = d1; _dim[2] = d2;
            allocate(3, _dim, true);
        }


        explicit Tensor(long d0, long d1, long d2, long d3) : _p(0) {
            _dim[0] = d0; _dim[1] = d1; _dim[2] = d2; _dim[3] = d3;
            allocate(4, _dim, true);
        }


        explicit Tensor(long d0, long d1, long d2, long d3, long d4) : _p(0) {
            _dim[0] = d0; _dim[1] = d1; _dim[2] = d2; _dim[3] = d3; _dim[4] = d4;
            allocate(5, _dim, true);
        }


        explicit Tensor(long d0, long d1, long d2, long d3, long d4, long d5) {
            _dim[0] = d0; _dim[1] = d1; _dim[2] = d2; _dim[3] = d3; _dim[4] = d4; _dim[5] = d5;
            allocate(6, _dim, true);
        }


        explicit Tensor(const std::vector<long>& d, bool dozero=true) : _p(0) {
            allocate(d.size(), d.size() ? &(d[0]) : 0, dozero);
        }


        explicit Tensor(long nd, const long d[], bool dozero=true) : _p(0) {
            allocate(nd,d,dozero);
        }


        Tensor<T>& operator=(T x) {
            UNARY_OPTIMIZED_ITERATOR(T,(*this),*_p0 = x);
            return *this;
        }


        Tensor<T>& fill(T x) {
            *this = x;
            return *this;
        }


        template <typename Q>
        Tensor<T>& operator+=(const Tensor<Q>& t) {
            BINARY_OPTIMIZED_ITERATOR(T, (*this), const T, t, *_p0 += *_p1);
            return *this;
        }


        template <typename Q>
        Tensor<T>& operator-=(const Tensor<Q>& t) {
            BINARY_OPTIMIZED_ITERATOR(T, (*this), const T, t, *_p0 -= *_p1);
            return *this;
        }


        template <typename Q>
        Tensor< TENSOR_RESULT_TYPE(T,Q) > operator+(const Tensor<Q>& t) const {
            typedef TENSOR_RESULT_TYPE(T,Q) resultT;
            Tensor<resultT> result(_ndim,_dim,false);
            TERNARY_OPTIMIZED_ITERATOR(resultT, result, const T, (*this), const Q, t, *_p0 = *_p1 + *_p2);
            return result;
        }


        template <typename Q>
        Tensor< TENSOR_RESULT_TYPE(T,Q) > operator-(const Tensor<Q>& t) const {
            typedef TENSOR_RESULT_TYPE(T,Q) resultT;
            Tensor<resultT> result(_ndim,_dim,false);
            TERNARY_OPTIMIZED_ITERATOR(resultT, result, const T, (*this), const Q, t, *_p0 = *_p1 - *_p2);
            return result;
        }


        template <typename Q>
        typename IsSupported<TensorTypeData<Q>, Tensor<TENSOR_RESULT_TYPE(T,Q)> >::type
        operator*(const Q& x) const {
            typedef TENSOR_RESULT_TYPE(T,Q) resultT;
            Tensor<resultT> result(_ndim,_dim,false);
            BINARY_OPTIMIZED_ITERATOR(resultT, result, const T, (*this), *_p0 = *_p1 * x);
            return result;
        }


        template <typename Q>
        typename IsSupported<TensorTypeData<Q>, Tensor<TENSOR_RESULT_TYPE(T,Q)> >::type
        operator/(const Q& x) const {
            typedef TENSOR_RESULT_TYPE(T,Q) resultT;
            Tensor<resultT> result(_ndim,_dim);
            BINARY_OPTIMIZED_ITERATOR(resultT, result, const T, (*this), *_p0 = *_p1 / x);
            return result;
        }


        template <typename Q>
        typename IsSupported<TensorTypeData<Q>, Tensor<TENSOR_RESULT_TYPE(T,Q)> >::type
        operator+(const Q& x) const {
            typedef TENSOR_RESULT_TYPE(T,Q) resultT;
            Tensor<resultT> result(_ndim,_dim);
            BINARY_OPTIMIZED_ITERATOR(resultT, result, const T, (*this), *_p0 = *_p1 + x);
            return result;
        }


        template <typename Q>
        typename IsSupported<TensorTypeData<Q>, Tensor<TENSOR_RESULT_TYPE(T,Q)> >::type
        operator-(const Q& x) const {
            return (*this) + (-x);
        }


        Tensor<T> operator-() const {
            Tensor<T> result = Tensor<T>(_ndim,_dim,false);
            BINARY_OPTIMIZED_ITERATOR(T, result, const T, (*this), *(_p0) = - (*_p1));
            return result;
        }


        template <typename Q>
        typename IsSupported<TensorTypeData<Q>,Tensor<T>&>::type
        operator*=(const Q& x) {
            UNARY_OPTIMIZED_ITERATOR(T, (*this), *_p0 *= x);
            return *this;
        }


        template <typename Q>
        typename IsSupported<TensorTypeData<Q>,Tensor<T>&>::type
        scale(Q x) {
            return (*this)*=x;
        }


        template <typename Q>
        typename IsSupported<TensorTypeData<Q>,Tensor<T>&>::type
        operator+=(const Q& x) {
            UNARY_OPTIMIZED_ITERATOR(T, (*this), *_p0 += x);
            return *this;
        }


        template <typename Q>
        typename IsSupported<TensorTypeData<Q>,Tensor<T>&>::type
        operator-=(const Q& x) {
            UNARY_OPTIMIZED_ITERATOR(T, (*this), *_p0 -= x);
            return *this;
        }



        Tensor<T>& conj() {
            UNARY_OPTIMIZED_ITERATOR(T, (*this), *_p0 = conditional_conj(*_p0));
            return *this;
        }


        Tensor<T>& fillrandom() {
            if (iscontiguous()) {
                madness::RandomVector<T>(size(), ptr());
            }
            else {
                UNARY_OPTIMIZED_ITERATOR(T,(*this), *_p0 = madness::RandomValue<T>());
            }
            return *this;
        }


        Tensor<T>& fillindex() {
            long count = 0;
            UNARY_UNOPTIMIZED_ITERATOR(T,(*this), *_p0 = count++); // Fusedim would be OK
            return *this;
        }


        Tensor<T>& screen(double x) {
            T zero = 0;
            UNARY_OPTIMIZED_ITERATOR(T,(*this), if (std::abs(*_p0)<x) *_p0=zero);
            return *this;
        }



        static bool bounds_checking() {
#ifdef TENSOR_BOUNDS_CHECKING
            return true;
#else
            return false;
#endif
        }


        T& operator[](long i) {
#ifdef TENSOR_BOUNDS_CHECKING
            TENSOR_ASSERT(i>=0 && i<_dim[0],"1d bounds check failed dim=0",i,this);
#endif
            return _p[i*_stride[0]];
        }


        const T& operator[](long i) const {
#ifdef TENSOR_BOUNDS_CHECKING
            TENSOR_ASSERT(i>=0 && i<_dim[0],"1d bounds check failed dim=0",i,this);
#endif
            return _p[i*_stride[0]];
        }


        T& operator()(long i) {
#ifdef TENSOR_BOUNDS_CHECKING
            TENSOR_ASSERT(i>=0 && i<_dim[0],"1d bounds check failed dim=0",i,this);
#endif
            return _p[i*_stride[0]];
        }


        const T& operator()(long i) const {
#ifdef TENSOR_BOUNDS_CHECKING
            TENSOR_ASSERT(i>=0 && i<_dim[0],"1d bounds check failed dim=0",i,this);
#endif
            return _p[i*_stride[0]];
        }


        T& operator()(long i, long j) {
#ifdef TENSOR_BOUNDS_CHECKING
            TENSOR_ASSERT(i>=0 && i<_dim[0],"2d bounds check failed dim=0",i,this);
            TENSOR_ASSERT(j>=0 && j<_dim[1],"2d bounds check failed dim=1",j,this);
#endif
            return _p[i*_stride[0]+j*_stride[1]];
        }


        const T& operator()(long i, long j) const {
#ifdef TENSOR_BOUNDS_CHECKING
            TENSOR_ASSERT(i>=0 && i<_dim[0],"2d bounds check failed dim=0",i,this);
            TENSOR_ASSERT(j>=0 && j<_dim[1],"2d bounds check failed dim=1",j,this);
#endif
            return _p[i*_stride[0]+j*_stride[1]];
        }


        T& operator()(long i, long j, long k) {
#ifdef TENSOR_BOUNDS_CHECKING
            TENSOR_ASSERT(i>=0 && i<_dim[0],"3d bounds check failed dim=0",i,this);
            TENSOR_ASSERT(j>=0 && j<_dim[1],"3d bounds check failed dim=1",j,this);
            TENSOR_ASSERT(k>=0 && k<_dim[2],"3d bounds check failed dim=2",k,this);
#endif
            return _p[i*_stride[0]+j*_stride[1]+k*_stride[2]];
        }


        const T& operator()(long i, long j, long k) const {
#ifdef TENSOR_BOUNDS_CHECKING
            TENSOR_ASSERT(i>=0 && i<_dim[0],"3d bounds check failed dim=0",i,this);
            TENSOR_ASSERT(j>=0 && j<_dim[1],"3d bounds check failed dim=1",j,this);
            TENSOR_ASSERT(k>=0 && k<_dim[2],"3d bounds check failed dim=2",k,this);
#endif
            return _p[i*_stride[0]+j*_stride[1]+k*_stride[2]];
        }


        T& operator()(long i, long j, long k, long l) {
#ifdef TENSOR_BOUNDS_CHECKING
            TENSOR_ASSERT(i>=0 && i<_dim[0],"4d bounds check failed dim=0",i,this);
            TENSOR_ASSERT(j>=0 && j<_dim[1],"4d bounds check failed dim=1",j,this);
            TENSOR_ASSERT(k>=0 && k<_dim[2],"4d bounds check failed dim=2",k,this);
            TENSOR_ASSERT(l>=0 && l<_dim[3],"4d bounds check failed dim=3",l,this);
#endif
            return _p[i*_stride[0]+j*_stride[1]+k*_stride[2]+
                           l*_stride[3]];
        }


        const T& operator()(long i, long j, long k, long l) const {
#ifdef TENSOR_BOUNDS_CHECKING
            TENSOR_ASSERT(i>=0 && i<_dim[0],"4d bounds check failed dim=0",i,this);
            TENSOR_ASSERT(j>=0 && j<_dim[1],"4d bounds check failed dim=1",j,this);
            TENSOR_ASSERT(k>=0 && k<_dim[2],"4d bounds check failed dim=2",k,this);
            TENSOR_ASSERT(l>=0 && l<_dim[3],"4d bounds check failed dim=3",l,this);
#endif
            return _p[i*_stride[0]+j*_stride[1]+k*_stride[2]+
                           l*_stride[3]];
        }


        T& operator()(long i, long j, long k, long l, long m) {
#ifdef TENSOR_BOUNDS_CHECKING
            TENSOR_ASSERT(i>=0 && i<_dim[0],"5d bounds check failed dim=0",i,this);
            TENSOR_ASSERT(j>=0 && j<_dim[1],"5d bounds check failed dim=1",j,this);
            TENSOR_ASSERT(k>=0 && k<_dim[2],"5d bounds check failed dim=2",k,this);
            TENSOR_ASSERT(l>=0 && l<_dim[3],"5d bounds check failed dim=3",l,this);
            TENSOR_ASSERT(m>=0 && m<_dim[4],"5d bounds check failed dim=4",m,this);
#endif
            return _p[i*_stride[0]+j*_stride[1]+k*_stride[2]+
                           l*_stride[3]+m*_stride[4]];
        }


        const T& operator()(long i, long j, long k, long l, long m) const {
#ifdef TENSOR_BOUNDS_CHECKING
            TENSOR_ASSERT(i>=0 && i<_dim[0],"5d bounds check failed dim=0",i,this);
            TENSOR_ASSERT(j>=0 && j<_dim[1],"5d bounds check failed dim=1",j,this);
            TENSOR_ASSERT(k>=0 && k<_dim[2],"5d bounds check failed dim=2",k,this);
            TENSOR_ASSERT(l>=0 && l<_dim[3],"5d bounds check failed dim=3",l,this);
            TENSOR_ASSERT(m>=0 && m<_dim[4],"5d bounds check failed dim=4",m,this);
#endif
            return _p[i*_stride[0]+j*_stride[1]+k*_stride[2]+
                           l*_stride[3]+m*_stride[4]];
        }


        T& operator()(long i, long j, long k, long l, long m, long n) {
#ifdef TENSOR_BOUNDS_CHECKING
            TENSOR_ASSERT(i>=0 && i<_dim[0],"6d bounds check failed dim=0",i,this);
            TENSOR_ASSERT(j>=0 && j<_dim[1],"6d bounds check failed dim=1",j,this);
            TENSOR_ASSERT(k>=0 && k<_dim[2],"6d bounds check failed dim=2",k,this);
            TENSOR_ASSERT(l>=0 && l<_dim[3],"6d bounds check failed dim=3",l,this);
            TENSOR_ASSERT(m>=0 && m<_dim[4],"6d bounds check failed dim=4",m,this);
            TENSOR_ASSERT(n>=0 && n<_dim[5],"6d bounds check failed dim=5",n,this);
#endif
            return _p[i*_stride[0]+j*_stride[1]+k*_stride[2]+
                           l*_stride[3]+m*_stride[4]+n*_stride[5]];
        }


        const T& operator()(long i, long j, long k, long l, long m, long n) const {
#ifdef TENSOR_BOUNDS_CHECKING
            TENSOR_ASSERT(i>=0 && i<_dim[0],"6d bounds check failed dim=0",i,this);
            TENSOR_ASSERT(j>=0 && j<_dim[1],"6d bounds check failed dim=1",j,this);
            TENSOR_ASSERT(k>=0 && k<_dim[2],"6d bounds check failed dim=2",k,this);
            TENSOR_ASSERT(l>=0 && l<_dim[3],"6d bounds check failed dim=3",l,this);
            TENSOR_ASSERT(m>=0 && m<_dim[4],"6d bounds check failed dim=4",m,this);
            TENSOR_ASSERT(n>=0 && n<_dim[5],"6d bounds check failed dim=5",n,this);
#endif
            return _p[i*_stride[0]+j*_stride[1]+k*_stride[2]+
                           l*_stride[3]+m*_stride[4]+n*_stride[5]];
        }


        T& operator()(const long ind[]) {
            long offset = 0;
            for (int d=0; d<_ndim; ++d) {
                long i = ind[d];
#ifdef TENSOR_BOUNDS_CHECKING
                TENSOR_ASSERT(i>=0 && i<_dim[0],"non-PC general indexing bounds check failed dim=",d,this);
#endif
                offset += i*_stride[d];
            }
            return _p[offset];
        }


        const T& operator()(const long ind[]) const {
            long offset = 0;
            for (int d=0; d<_ndim; ++d) {
                long i = ind[d];
#ifdef TENSOR_BOUNDS_CHECKING
                TENSOR_ASSERT(i>=0 && i<_dim[0],"non-PC general indexing bounds check failed dim=",d,this);
#endif
                offset += i*_stride[d];
            }
            return _p[offset];
        }


        T& operator()(const std::vector<long> ind) {
            TENSOR_ASSERT(ind.size()>=(unsigned int) _ndim,"invalid number of dimensions",ind.size(),this);
            long index=0;
            for (long d=0; d<_ndim; ++d) {
                TENSOR_ASSERT(ind[d]>=0 && ind[d]<_dim[d],"out-of-bounds access",ind[d],this);
                index += ind[d]*_stride[d];
            }
            return _p[index];
        }


        const T& operator()(const std::vector<long> ind) const {
            TENSOR_ASSERT(ind.size()>=(unsigned int) _ndim,"invalid number of dimensions",ind.size(),this);
            long index=0;
            for (long d=0; d<_ndim; ++d) {
                TENSOR_ASSERT(ind[d]>=0 && ind[d]<_dim[d],"out-of-bounds access",ind[d],this);
                index += ind[d]*_stride[d];
            }
            return _p[index];
        }


        SliceTensor<T> operator()(const std::vector<Slice>& s) {
            TENSOR_ASSERT(s.size()>=(unsigned)(this->ndim()), "invalid number of dimensions",
                          this->ndim(),this);
            return SliceTensor<T>(*this,&(s[0]));
        }


        const Tensor<T> operator()(const std::vector<Slice>& s) const {
            TENSOR_ASSERT(s.size()>=(unsigned)(this->ndim()), "invalid number of dimensions",
                          this->ndim(),this);
            return SliceTensor<T>(*this,&(s[0]));
        }


        SliceTensor<T> operator()(const Slice& s0) {
            TENSOR_ASSERT(this->ndim()==1,"invalid number of dimensions",
                          this->ndim(),this);
            Slice s[1] = {s0};
            return SliceTensor<T>(*this,s);
        }


        const Tensor<T> operator()(const Slice& s0) const {
            TENSOR_ASSERT(this->ndim()==1,"invalid number of dimensions",
                          this->ndim(),this);
            Slice s[1] = {s0};
            return SliceTensor<T>(*this,s);
        }


        SliceTensor<T> operator()(long i, const Slice& s1) {
            TENSOR_ASSERT(this->ndim()==2,"invalid number of dimensions",
                          this->ndim(),this);
            Slice s[2] = {Slice(i,i,0),s1};
            return SliceTensor<T>(*this,s);
        }


        const Tensor<T> operator()(long i, const Slice& s1) const {
            TENSOR_ASSERT(this->ndim()==2,"invalid number of dimensions",
                          this->ndim(),this);
            Slice s[2] = {Slice(i,i,0),s1};
            return SliceTensor<T>(*this,s);
        }


        SliceTensor<T> operator()(const Slice& s0, long j) {
            TENSOR_ASSERT(this->ndim()==2,"invalid number of dimensions",
                          this->ndim(),this);
            Slice s[2] = {s0,Slice(j,j,0)};
            return SliceTensor<T>(*this,s);
        }


        const Tensor<T> operator()(const Slice& s0, long j) const {
            TENSOR_ASSERT(this->ndim()==2,"invalid number of dimensions",
                          this->ndim(),this);
            Slice s[2] = {s0,Slice(j,j,0)};
            return SliceTensor<T>(*this,s);
        }


        SliceTensor<T> operator()(const Slice& s0, const Slice& s1) {
            TENSOR_ASSERT(this->ndim()==2,"invalid number of dimensions",
                          this->ndim(),this);
            Slice s[2] = {s0,s1};
            return SliceTensor<T>(*this,s);
        }


        const Tensor<T> operator()(const Slice& s0, const Slice& s1) const {
            TENSOR_ASSERT(this->ndim()==2,"invalid number of dimensions",
                          this->ndim(),this);
            Slice s[2] = {s0,s1};
            return SliceTensor<T>(*this,s);
        }


        SliceTensor<T> operator()(const Slice& s0, const Slice& s1, const Slice& s2) {
            TENSOR_ASSERT(this->ndim()==3,"invalid number of dimensions",
                          this->ndim(),this);
            Slice s[3] = {s0,s1,s2};
            return SliceTensor<T>(*this,s);
        }


        const Tensor<T> operator()(const Slice& s0, const Slice& s1, const Slice& s2) const {
            TENSOR_ASSERT(this->ndim()==3,"invalid number of dimensions",
                          this->ndim(),this);
            Slice s[3] = {s0,s1,s2};
            return SliceTensor<T>(*this,s);
        }


        SliceTensor<T> operator()(long i, const Slice& s1, const Slice& s2) {
            TENSOR_ASSERT(this->ndim()==3,"invalid number of dimensions",
                          this->ndim(),this);
            Slice s[3] = {Slice(i,i,0),s1,s2};
            return SliceTensor<T>(*this,s);
        }


        const Tensor<T> operator()(long i, const Slice& s1, const Slice& s2) const {
            TENSOR_ASSERT(this->ndim()==3,"invalid number of dimensions",
                          this->ndim(),this);
            Slice s[3] = {Slice(i,i,0),s1,s2};
            return SliceTensor<T>(*this,s);
        }


        SliceTensor<T> operator()(const Slice& s0, long j, const Slice& s2) {
            TENSOR_ASSERT(this->ndim()==3,"invalid number of dimensions",
                          this->ndim(),this);
            Slice s[3] = {s0,Slice(j,j,0),s2};
            return SliceTensor<T>(*this,s);
        }


        const Tensor<T> operator()(const Slice& s0, long j, const Slice& s2) const {
            TENSOR_ASSERT(this->ndim()==3,"invalid number of dimensions",
                          this->ndim(),this);
            Slice s[3] = {s0,Slice(j,j,0),s2};
            return SliceTensor<T>(*this,s);
        }


        SliceTensor<T> operator()(const Slice& s0, const Slice& s1, long k) {
            TENSOR_ASSERT(this->ndim()==3,"invalid number of dimensions",
                          this->ndim(),this);
            Slice s[3] = {s0,s1,Slice(k,k,0)};
            return SliceTensor<T>(*this,s);
        }


        const Tensor<T> operator()(const Slice& s0, const Slice& s1, long k) const {
            TENSOR_ASSERT(this->ndim()==3,"invalid number of dimensions",
                          this->ndim(),this);
            Slice s[3] = {s0,s1,Slice(k,k,0)};
            return SliceTensor<T>(*this,s);
        }


        SliceTensor<T> operator()(long i, long j, const Slice& s2) {
            TENSOR_ASSERT(this->ndim()==3,"invalid number of dimensions",
                          this->ndim(),this);
            Slice s[3] = {Slice(i,i,0),Slice(j,j,0),s2};
            return SliceTensor<T>(*this,s);
        }


        const Tensor<T> operator()(long i, long j, const Slice& s2) const {
            TENSOR_ASSERT(this->ndim()==3,"invalid number of dimensions",
                          this->ndim(),this);
            Slice s[3] = {Slice(i,i,0),Slice(j,j,0),s2};
            return SliceTensor<T>(*this,s);
        }


        SliceTensor<T> operator()(long i, const Slice& s1, long k) {
            TENSOR_ASSERT(this->ndim()==3,"invalid number of dimensions",
                          this->ndim(),this);
            Slice s[3] = {Slice(i,i,0),s1,Slice(k,k,0)};
            return SliceTensor<T>(*this,s);
        }


        const Tensor<T> operator()(long i, const Slice& s1, long k) const {
            TENSOR_ASSERT(this->ndim()==3,"invalid number of dimensions",
                          this->ndim(),this);
            Slice s[3] = {Slice(i,i,0),s1,Slice(k,k,0)};
            return SliceTensor<T>(*this,s);
        }


        SliceTensor<T> operator()(const Slice& s0, long j, long k) {
            TENSOR_ASSERT(this->ndim()==3,"invalid number of dimensions",
                          this->ndim(),this);
            Slice s[3] = {s0,Slice(j,j,0),Slice(k,k,0)};
            return SliceTensor<T>(*this,s);
        }


        const Tensor<T> operator()(const Slice& s0, long j, long k) const {
            TENSOR_ASSERT(this->ndim()==3,"invalid number of dimensions",
                          this->ndim(),this);
            Slice s[3] = {s0,Slice(j,j,0),Slice(k,k,0)};
            return SliceTensor<T>(*this,s);
        }


        SliceTensor<T> operator()(const Slice& s0, const Slice& s1, const Slice& s2,
                                  const Slice& s3) {
            TENSOR_ASSERT(this->ndim()==4,"invalid number of dimensions",
                          this->ndim(),this);
            Slice s[4] = {s0,s1,s2,s3};
            return SliceTensor<T>(*this,s);
        }


        const Tensor<T> operator()(const Slice& s0, const Slice& s1, const Slice& s2,
                                  const Slice& s3) const {
            TENSOR_ASSERT(this->ndim()==4,"invalid number of dimensions",
                          this->ndim(),this);
            Slice s[4] = {s0,s1,s2,s3};
            return SliceTensor<T>(*this,s);
        }


        SliceTensor<T> operator()(const Slice& s0, const Slice& s1, const Slice& s2,
                                  const Slice& s3, const Slice& s4) {
            TENSOR_ASSERT(this->ndim()==5,"invalid number of dimensions",
                          this->ndim(),this);
            Slice s[5] = {s0,s1,s2,s3,s4};
            return SliceTensor<T>(*this,s);
        }


        const Tensor<T> operator()(const Slice& s0, const Slice& s1, const Slice& s2,
                                  const Slice& s3, const Slice& s4) const {
            TENSOR_ASSERT(this->ndim()==5,"invalid number of dimensions",
                          this->ndim(),this);
            Slice s[5] = {s0,s1,s2,s3,s4};
            return SliceTensor<T>(*this,s);
        }


        SliceTensor<T> operator()(const Slice& s0, const Slice& s1, const Slice& s2,
                                  const Slice& s3, const Slice& s4, const Slice& s5) {
            TENSOR_ASSERT(this->ndim()==6,"invalid number of dimensions",
                          this->ndim(),this);
            Slice s[6] = {s0,s1,s2,s3,s4,s5};
            return SliceTensor<T>(*this,s);
        }



        const Tensor<T> operator()(const Slice& s0, const Slice& s1, const Slice& s2,
                                  const Slice& s3, const Slice& s4, const Slice& s5) const {
            TENSOR_ASSERT(this->ndim()==6,"invalid number of dimensions",
                          this->ndim(),this);
            Slice s[6] = {s0,s1,s2,s3,s4,s5};
            return SliceTensor<T>(*this,s);
        }


        Tensor<T> reshape(int ndimnew, const long* d) {
            Tensor<T> result(*this);
            result.reshape_inplace(ndimnew,d);
            return result;
        }


        const Tensor<T> reshape(int ndimnew, const long* d) const {
            Tensor<T> result(*const_cast<Tensor<T>*>(this));
            result.reshape_inplace(ndimnew,d);
            return result;
        }


        Tensor<T> reshape(const std::vector<long>& d) {
          return reshape(d.size(), d.size() ? &d[0] : 0);
        }


        const Tensor<T> reshape(const std::vector<long>& d) const {
            return reshape(d.size(), d.size() ? &d[0] : 0);
        }


        Tensor<T> reshape(long dim0) {
            long d[1] = {dim0};
            return reshape(1,d);
        }

        const Tensor<T> reshape(long dim0) const {
            long d[1] = {dim0};
            return reshape(1,d);
        }


        Tensor<T> reshape(long dim0, long dim1) {
            long d[2] = {dim0,dim1};
            return reshape(2,d);
        }


        const Tensor<T> reshape(long dim0, long dim1) const {
            long d[2] = {dim0,dim1};
            return reshape(2,d);
        }


        Tensor<T> reshape(long dim0, long dim1, long dim2) {
            long d[3] = {dim0,dim1,dim2};
            return reshape(3,d);
        }


        const Tensor<T> reshape(long dim0, long dim1, long dim2) const {
            long d[3] = {dim0,dim1,dim2};
            return reshape(3,d);
        }


        Tensor<T> reshape(long dim0, long dim1, long dim2, long dim3) {
            long d[4] = {dim0,dim1,dim2,dim3};
            return reshape(4,d);
        }


        const Tensor<T> reshape(long dim0, long dim1, long dim2, long dim3) const {
            long d[4] = {dim0,dim1,dim2,dim3};
            return reshape(4,d);
        }


        Tensor<T> reshape(long dim0, long dim1, long dim2, long dim3, long dim4) {
            long d[5] = {dim0,dim1,dim2,dim3,dim4};
            return reshape(5,d);
        }


        const Tensor<T> reshape(long dim0, long dim1, long dim2, long dim3, long dim4) const {
            long d[5] = {dim0,dim1,dim2,dim3,dim4};
            return reshape(5,d);
        }


        Tensor<T> reshape(long dim0, long dim1, long dim2, long dim3, long dim4, long dim5) {
            long d[6] = {dim0,dim1,dim2,dim3,dim4,dim5};
            return reshape(6,d);
        }


        const Tensor<T> reshape(long dim0, long dim1, long dim2, long dim3, long dim4, long dim5) const {
            long d[6] = {dim0,dim1,dim2,dim3,dim4,dim5};
            return reshape(6,d);
        }

        Tensor<T> flat() {
            long d[1] = {_size};
            return reshape(1,d);
        }

        const Tensor<T> flat() const {
            long d[1] = {_size};
            return reshape(1,d);
        }


        Tensor<T> splitdim(long i, long dimi0, long dimi1) {
            Tensor<T> result(*this);
            result.splitdim_inplace(i, dimi0, dimi1);
            return result;
        }


        const Tensor<T> splitdim(long i, long dimi0, long dimi1) const {
            Tensor<T> result(*const_cast<Tensor<T>*>(this));
            result.splitdim_inplace(i, dimi0, dimi1);
            return result;
        }


        Tensor<T> fusedim(long i) {
            Tensor<T> result(*this);
            result.fusedim_inplace(i);
            return result;
        }


        const Tensor<T> fusedim(long i) const {
            Tensor<T> result(*const_cast<Tensor<T>*>(this));
            result.fusedim_inplace(i);
            return result;
        }


        Tensor<T> swapdim(long idim, long jdim) {
            Tensor<T> result(*this);
            result.swapdim_inplace(idim, jdim);
            return result;
        }


        const Tensor<T> swapdim(long idim, long jdim) const {
            Tensor<T> result(*const_cast<Tensor<T>*>(this));
            result.swapdim_inplace(idim, jdim);
            return result;
        }


        Tensor<T> mapdim(const std::vector<long>& map) {
            Tensor<T> result(*this);
            result.mapdim_inplace(map);
            return result;
        }


        const Tensor<T> mapdim(const std::vector<long>& map) const {
            Tensor<T> result(*const_cast<Tensor<T>*>(this));
            result.mapdim_inplace(map);
            return result;
        }


        Tensor<T> cycledim(long nshift, long start, long end) {
            Tensor<T> result(*this);
            result.cycledim_inplace(nshift, start, end);
            return result;
        }


        const Tensor<T> cycledim(long nshift, long start, long end) const {
            Tensor<T> result(*const_cast<Tensor<T>*>(this));
            result.cycledim_inplace(nshift, start, end);
            return result;
        }


        template <class Q> bool conforms(const Tensor<Q>& t) const {
            return BaseTensor::conforms(&t);
        }

        T sum() const {
            T result = 0;
            UNARY_OPTIMIZED_ITERATOR(const T,(*this),result += *_p0);
            return result;
        }

        T sumsq() const {
            T result = 0;
            UNARY_OPTIMIZED_ITERATOR(const T,(*this),result += (*_p0) * (*_p0));
            return result;
        }

        T product() const {
            T result = 1;
            UNARY_OPTIMIZED_ITERATOR(const T,(*this),result *= *_p0);
            return result;
        }

        T min(long* ind=0) const {
            T result = *(this->_p);
            if (ind) {
                for (long i=0; i<_ndim; ++i) ind[i]=0;
                long nd = _ndim-1;
                UNARY_UNOPTIMIZED_ITERATOR(const T,(*this),
                                           if (result > *_p0) {
                                               result = *_p0;
                                               for (long i=0; i<nd; ++i) ind[i]=iter.ind[i];
                                               ind[nd] = _j;
                                           }
                                           );
            }
            else {
                UNARY_OPTIMIZED_ITERATOR(const T,(*this),result=std::min<T>(result,*_p0));
            }
            return result;
        }

        T max(long* ind=0) const {
            T result = *(this->_p);
            if (ind) {
                for (long i=0; i<_ndim; ++i) ind[i]=0;
                long nd = _ndim-1;
                UNARY_UNOPTIMIZED_ITERATOR(const T,(*this),
                                           if (result < *_p0) {
                                               result = *_p0;
                                               for (long i=0; i<nd; ++i) ind[i]=iter.ind[i];
                                               ind[nd] = _j;
                                           }
                                           );
            }
            else {
                UNARY_OPTIMIZED_ITERATOR(const T,(*this),result=std::max<T>(result,*_p0));
            }
            return result;
        }

        // For complex types, this next group returns the appropriate real type
        // For real types, the same type as T is returned (type_data.h)

        float_scalar_type normf() const {
            float_scalar_type result = 0;
            UNARY_OPTIMIZED_ITERATOR(const T,(*this),result += ::madness::detail::mynorm(*_p0));
            return (float_scalar_type) std::sqrt(result);
        }

        scalar_type absmin(long *ind = 0) const {
            scalar_type result = std::abs(*(this->_p));
            if (ind) {
                for (long i=0; i<_ndim; ++i) ind[i]=0;
                long nd = _ndim-1;
                UNARY_UNOPTIMIZED_ITERATOR(const T,(*this),
                                           scalar_type absval = std::abs(*_p0);
                                           if (result > absval) {
                                               result = absval;
                                               for (long i=0; i<nd; ++i) ind[i]=iter.ind[i];
                                               ind[nd] = _j;
                                           }
                                           );
            }
            else {
                UNARY_OPTIMIZED_ITERATOR(const T,(*this),result=std::min<scalar_type>(result,std::abs(*_p0)));
            }
            return result;
        }

        scalar_type absmax(long *ind = 0) const {
            scalar_type result = std::abs(*(this->_p));
            if (ind) {
                for (long i=0; i<_ndim; ++i) ind[i]=0;
                long nd = _ndim-1;
                UNARY_UNOPTIMIZED_ITERATOR(T,(*this),
                                           scalar_type absval = std::abs(*_p0);
                                           if (result < absval) {
                                               result = absval;
                                               for (long i=0; i<nd; ++i) ind[i]=iter.ind[i];
                                               ind[nd] = _j;
                                           }
                                           );
            }
            else {
                UNARY_OPTIMIZED_ITERATOR(const T,(*this),result=std::max<scalar_type>(result,std::abs(*_p0)));
            }
            return result;
        }


        T trace(const Tensor<T>& t) const {
            T result = 0;
            BINARY_OPTIMIZED_ITERATOR(const T,(*this),const T,t,result += (*_p0)*(*_p1));
            return result;
        }

        template <class Q>
        TENSOR_RESULT_TYPE(T,Q) trace_conj(const Tensor<Q>& t) const {
            TENSOR_RESULT_TYPE(T,Q) result = 0;
            BINARY_OPTIMIZED_ITERATOR(const T,(*this),const Q,t,result += conditional_conj(*_p0)*(*_p1));
            return result;
        }

        template <typename opT>
        Tensor<T>& unaryop(opT& op) {
            UNARY_OPTIMIZED_ITERATOR(T,(*this),*_p0=op(*_p0));
            return *this;
        }

        Tensor<T>& emul(const Tensor<T>& t) {
            BINARY_OPTIMIZED_ITERATOR(T,(*this),const T,t,*_p0 *= *_p1);
            return *this;
        }

        Tensor<T>& gaxpy(T alpha, const Tensor<T>& t, T beta) {
            if (iscontiguous() && t.iscontiguous()) {
                T* restrict a = ptr();
                const T* restrict b = t.ptr();
                if (alpha == T(1.0)) {
                    for (long i=0; i<_size; ++i) a[i] += b[i]*beta;
                }
                else {
                    for (long i=0; i<_size; ++i) a[i] = a[i]*alpha + b[i]*beta;
                }
            }
            else {
                //BINARYITERATOR(T,(*this),T,t, (*_p0) = alpha*(*_p0) + beta*(*_p1));
                BINARY_OPTIMIZED_ITERATOR(T,(*this),const T,t, (*_p0) = alpha*(*_p0) + beta*(*_p1));
                //ITERATOR((*this),(*this)(IND) = alpha*(*this)(IND) + beta*t(IND));
            }
            return *this;
        }

        T* ptr() {
            return _p;
        }

        const T* ptr() const {
            return _p;
        }

        BaseTensor* base() {
            return static_cast<BaseTensor*>(this);
        }

        const BaseTensor* base() const {
            return static_cast<BaseTensor*>(this);
        }

        TensorIterator<T> unary_iterator(long iterlevel=0,
                                         bool optimize=true,
                                         bool fusedim=true,
                                         long jdim=default_jdim) const {
            return TensorIterator<T>(this,(const Tensor<T>*) 0, (const Tensor<T>*) 0,
                                     iterlevel, optimize, fusedim, jdim);
        }

        template <class Q>
        TensorIterator<T,Q> binary_iterator(const Tensor<Q>& q,
                                            long iterlevel=0,
                                            bool optimize=true,
                                            bool fusedim=true,
                                            long jdim=default_jdim) const {
            return TensorIterator<T,Q>(this,&q,(const Tensor<T>*) 0,
                                       iterlevel, optimize, fusedim, jdim);
        }

        template <class Q, class R>
        TensorIterator<T,Q,R> ternary_iterator(const Tensor<Q>& q,
                                               const Tensor<R>& r,
                                               long iterlevel=0,
                                               bool optimize=true,
                                               bool fusedim=true,
                                               long jdim=default_jdim) const {
            return TensorIterator<T,Q,R>(this,&q,&r,
                                         iterlevel, optimize, fusedim, jdim);
        }

        const TensorIterator<T>& end() const {
            static TensorIterator<T> theend(0,0,0,0,0,0);
            return theend;
        }

        virtual ~Tensor() {}

        void clear() {deallocate();}

        bool has_data() const {return size()!=0;};

    };

    template <class T>
    std::ostream& operator << (std::ostream& out, const Tensor<T>& t);


    namespace archive {
        template <class Archive, typename T>
        struct ArchiveStoreImpl< Archive, Tensor<T> > {
            static void store(const Archive& s, const Tensor<T>& t) {
                if (t.iscontiguous()) {
                    s & t.size() & t.id();
                    if (t.size()) s & t.ndim() & wrap(t.dims(),TENSOR_MAXDIM) & wrap(t.ptr(),t.size());
                }
                else {
                    s & copy(t);
                }
            };
        };


        template <class Archive, typename T>
        struct ArchiveLoadImpl< Archive, Tensor<T> > {
            static void load(const Archive& s, Tensor<T>& t) {
                long sz = 0l, id = 0l;
                s & sz & id;
                if (id != t.id()) throw "type mismatch deserializing a tensor";
                if (sz) {
                    long _ndim = 0l, _dim[TENSOR_MAXDIM];
                    s & _ndim & wrap(_dim,TENSOR_MAXDIM);
                    t = Tensor<T>(_ndim, _dim, false);
                    if (sz != t.size()) throw "size mismatch deserializing a tensor";
                    s & wrap(t.ptr(), t.size());
                }
                else {
                    t = Tensor<T>();
                }
            };
        };

    }


    template <typename T, typename Q>
    typename IsSupported < TensorTypeData<Q>, Tensor<T> >::type
    operator+(Q x, const Tensor<T>& t) {
        return t+x;
    }


    template <typename T, typename Q>
    typename IsSupported < TensorTypeData<Q>, Tensor<T> >::type
    operator*(const Q& x, const Tensor<T>& t) {
        return t*x;
    }


    template <typename T, typename Q>
    typename IsSupported < TensorTypeData<Q>, Tensor<T> >::type
    operator-(Q x, const Tensor<T>& t) {
        return (-t)+=x;
    }


    template <class T> Tensor<T> copy(const Tensor<T>& t) {
        if (t.size()) {
            Tensor<T> result = Tensor<T>(t.ndim(),t.dims(),false);
            BINARY_OPTIMIZED_ITERATOR(T, result, const T, t, *_p0 = *_p1);
            return result;
        }
        else {
            return Tensor<T>();
        }
    }


    template <class T, class Q>
    Tensor<TENSOR_RESULT_TYPE(T,Q)> transform_dir(const Tensor<T>& t, const Tensor<Q>& c, int axis) {
        if (axis == 0) {
            return inner(c,t,0,axis);
        }
        else if (axis == t.ndim()-1) {
            return inner(t,c,axis,0);
        }
        else {
            return copy(inner(t,c,axis,0).cycledim(1,axis, -1)); // Copy to make contiguous
        }
    }


    template <class T>
    Tensor<T> transpose(const Tensor<T>& t) {
        TENSOR_ASSERT(t.ndim() == 2, "transpose requires a matrix", t.ndim(), &t);
        return copy(t.swapdim(0,1));
    }


    template <class T>
    Tensor<T> conj_transpose(const Tensor<T>& t) {
        TENSOR_ASSERT(t.ndim() == 2, "conj_transpose requires a matrix", t.ndim(), &t);
        return conj(t.swapdim(0,1));
    }


    template <class T> class SliceTensor : public Tensor<T> {
    private:
        SliceTensor<T>();

    public:
        SliceTensor(const Tensor<T>& t, const Slice s[])
            : Tensor<T>(const_cast<Tensor<T>&>(t)) 
        {
            // C++ standard says class derived from parameterized base class cannot
            // directly access the base class elements ... must explicitly reference.

            long nd = 0, size=1;
            for (long i=0; i<t._ndim; ++i) {
                long start=s[i].start, end=s[i].end, step=s[i].step;
                //std::printf("%ld input start=%ld end=%ld step=%ld\n",
                //i, start, end, step);
                if (start < 0) start += this->_dim[i];
                if (end < 0) end += this->_dim[i];
                long len = end-start+1;
                if (step) len /= step;  // Rounds len towards zero

                // if input length is not exact multiple of step, round end towards start
                // for the same behaviour of for (i=start; i<=end; i+=step);
                end = start + (len-1)*step;

                //std::printf("%ld munged start=%ld end=%ld step=%ld len=%ld _dim=%ld\n",
                //              i, start, end, step, len, this->_dim[i]);

                TENSOR_ASSERT(start>=0 && start<this->_dim[i],"slice start invalid",start,this);
                TENSOR_ASSERT(end>=0 && end<this->_dim[i],"slice end invalid",end,this);
                TENSOR_ASSERT(len>0,"slice length must be non-zero",len,this);

                this->_p += start * t._stride[i];

                if (step) {
                    size *= len;
                    this->_dim[nd] = len;
                    this->_stride[nd] = step * t._stride[i];
                    ++nd;
                }
            }
            //For Python interface need to be able to return a scalar inside a tensor with nd=0
            //TENSOR_ASSERT(nd>0,"slicing produced a scalar, but cannot return one",nd,this);
            for (long i=nd; i<TENSOR_MAXDIM; ++i) { // So can iterate over missing dimensions
                this->_dim[i] = 1;
                this->_stride[i] = 0;
            }

            this->_ndim = nd;
            this->_size = size;
        }

        SliceTensor<T>& operator=(const SliceTensor<T>& t) {
            BINARY_OPTIMIZED_ITERATOR(T, (*this), const T, t, *_p0 = (T)(*_p1));
            return *this;
        }

        template <class Q>
        SliceTensor<T>& operator=(const SliceTensor<Q>& t) {
            BINARY_OPTIMIZED_ITERATOR(T, (*this), const Q, t, *_p0 = (T)(*_p1));
            return *this;
        }

        SliceTensor<T>& operator=(const Tensor<T>& t) {
            BINARY_OPTIMIZED_ITERATOR(T, (*this), const T, t, *_p0 = (T)(*_p1));
            return *this;
        }

        template <class Q>
        SliceTensor<T>& operator=(const Tensor<Q>& t) {
            BINARY_OPTIMIZED_ITERATOR(T, (*this), const Q, t, *_p0 = (T)(*_p1));
            return *this;
        }

        SliceTensor<T>& operator=(const T& t) {
            UNARY_OPTIMIZED_ITERATOR(T, (*this), *_p0 = t);
            return *this;
        }

        virtual ~SliceTensor() {};              // Tensor<T> destructor does enough
    };


    // Specializations for complex types
    template<> float_complex Tensor<float_complex>::min(long* ind) const ;
    template<> double_complex Tensor<double_complex>::min(long* ind) const ;
    template<> float_complex Tensor<float_complex>::max(long* ind) const ;
    template<> double_complex Tensor<double_complex>::max(long* ind) const ;

    // Stream stuff


    template <class T>
    std::ostream& operator << (std::ostream& s, const Tensor<T>& t) {
        if (t.size() == 0) {
            s << "[empty tensor]\n";
            return s;
        }

        long maxdim = 0;
        long index_width = 0;
        for (int i = 0; i<(t.ndim()-1); ++i) {
            if (maxdim < t.dim(i)) maxdim = t.dim(i);
        }
        if (maxdim < 10)
            index_width = 1;
        else if (maxdim < 100)
            index_width = 2;
        else if (maxdim < 1000)
            index_width = 3;
        else if (maxdim < 10000)
            index_width = 4;
        else
            index_width = 6;

        std::ios::fmtflags oldflags = s.setf(std::ios::scientific);
        long oldprec = s.precision();
        long oldwidth = s.width();

        // C++ formatted IO is worse than Fortran !!
        for (TensorIterator<T> iter=t.unary_iterator(1,false,false); iter!=t.end(); ++iter) {
            const T* p = iter._p0;
            long inc = iter._s0;
            long dimj = iter.dimj;
            s.unsetf(std::ios::scientific);
            s << '[';
            for (long i=0; i<iter.ndim; ++i) {
                s.width(index_width);
                s << iter.ind[i];
                if (i != iter.ndim) s << ",";
            }
            s << "*]";
//flo            s.setf(std::ios::scientific);
            s.setf(std::ios::fixed);
            for (long j=0; j<dimj; ++j, p+=inc) {
//flo                s.precision(4);
                s.precision(8);
                s.width(12);
                s << *p;
            }
            s.unsetf(std::ios::scientific);
            s << std::endl;
        }
        s.setf(oldflags,std::ios::floatfield);
        s.precision(oldprec);
        s.width(oldwidth);

        return s;
    }



    template <class T>
    Tensor<T> outer(const Tensor<T>& left, const Tensor<T>& right) {
        long nd = left.ndim() + right.ndim();
        TENSOR_ASSERT(nd <= TENSOR_MAXDIM,"too many dimensions in result",
                      nd,0);
        long d[TENSOR_MAXDIM];
        for (long i=0; i<left.ndim(); ++i) d[i] = left.dim(i);
        for (long i=0; i<right.ndim(); ++i) d[i+left.ndim()] = right.dim(i);
        Tensor<T> result(nd,d,false);
        T* ptr = result.ptr();

        TensorIterator<T> iter=right.unary_iterator(1,false,true);
        for (TensorIterator<T> p=left.unary_iterator(); p!=left.end(); ++p) {
            T val1 = *p;
            // Cannot reorder dimensions, but can fuse contiguous dimensions
            for (iter.reset(); iter._p0; ++iter) {
                long dimj = iter.dimj;
                T* _p0 = iter._p0;
                long Tstride = iter._s0;
                for (long _j=0; _j<dimj; ++_j, _p0+=Tstride) {
                    *ptr++ = val1 * (*_p0);
                }
            }
        }

        return result;
    }



    template <class T, class Q>
    Tensor<TENSOR_RESULT_TYPE(T,Q)> inner(const Tensor<T>& left, const Tensor<Q>& right,
                                          long k0=-1, long k1=0) {
        if (k0 < 0) k0 += left.ndim();
        if (k1 < 0) k1 += right.ndim();
        long nd = left.ndim() + right.ndim() - 2;
        TENSOR_ASSERT(nd!=0, "result is a scalar but cannot return one ... use dot",
                      nd, &left);
        TENSOR_ASSERT(left.dim(k0) == right.dim(k1),"common index must be same length",
                      right.dim(k1), &left);

        TENSOR_ASSERT(nd > 0 && nd <= TENSOR_MAXDIM,
                      "invalid number of dimensions in the result", nd,0);

        long d[TENSOR_MAXDIM];

        long base=0;
        for (long i=0; i<k0; ++i) d[i] = left.dim(i);
        for (long i=k0+1; i<left.ndim(); ++i) d[i-1] = left.dim(i);
        base = left.ndim()-1;
        for (long i=0; i<k1; ++i) d[i+base] = right.dim(i);
        base--;
        for (long i=k1+1; i<right.ndim(); ++i) d[i+base] = right.dim(i);

        Tensor<TENSOR_RESULT_TYPE(T,Q)> result(nd,d);

        inner_result(left,right,k0,k1,result);

        return result;
    }


    template <class T, class Q>
    void inner_result(const Tensor<T>& left, const Tensor<Q>& right,
                      long k0, long k1, Tensor< TENSOR_RESULT_TYPE(T,Q) >& result) {

        typedef TENSOR_RESULT_TYPE(T,Q) resultT;
        // Need to include explicit optimizations for common special cases
        // E.g., contiguous, matrix-matrix, and 3d-tensor*matrix

        resultT* ptr = result.ptr();

        if (k0 < 0) k0 += left.ndim();
        if (k1 < 0) k1 += right.ndim();

        if (left.iscontiguous() && right.iscontiguous()) {
            if (k0==0 && k1==0) {
                // c[i,j] = a[k,i]*b[k,j] ... collapsing extra indices to i & j
                long dimk = left.dim(k0);
                long dimj = right.stride(0);
                long dimi = left.stride(0);
                ::mTxm(dimi,dimj,dimk,ptr,left.ptr(),right.ptr());
                return;
            }
            else if (k0==(left.ndim()-1) && k1==(right.ndim()-1)) {
                // c[i,j] = a[i,k]*b[j,k] ... collapsing extra indices to i & j
                long dimk = left.dim(k0);
                long dimi = left.size()/dimk;
                long dimj = right.size()/dimk;
                ::mxmT(dimi,dimj,dimk,ptr,left.ptr(),right.ptr());
                return;
            }
            else if (k0==0 && k1==(right.ndim()-1)) {
                // c[i,j] = a[k,i]*b[j,k] ... collapsing extra indices to i & j
                long dimk = left.dim(k0);
                long dimi = left.stride(0);
                long dimj = right.size()/dimk;
                ::mTxmT(dimi,dimj,dimk,ptr,left.ptr(),right.ptr());
                return;
            }
            else if (k0==(left.ndim()-1) && k1==0) {
                // c[i,j] = a[i,k]*b[k,j] ... collapsing extra indices to i & j
                long dimk = left.dim(k0);
                long dimi = left.size()/dimk;
                long dimj = right.stride(0);
                ::mxm(dimi,dimj,dimk,ptr,left.ptr(),right.ptr());
                return;
            }
        }

        long dimj = left.dim(k0);
        TensorIterator<Q> iter1=right.unary_iterator(1,false,false,k1);

        for (TensorIterator<T> iter0=left.unary_iterator(1,false,false,k0);
                iter0._p0; ++iter0) {
            T* restrict xp0 = iter0._p0;
            long s0 = iter0._s0;
            for (iter1.reset(); iter1._p0; ++iter1) {
                T* restrict p0 = xp0;
                Q* restrict p1 = iter1._p0;
                long s1 = iter1._s0;
                resultT sum = 0;
                for (long j=0; j<dimj; ++j,p0+=s0,p1+=s1) {
                    sum += (*p0) * (*p1);
                }
                *ptr++ += sum;
            }
        }
    }


    template <class T, class Q>
    Tensor<TENSOR_RESULT_TYPE(T,Q)> transform(const Tensor<T>& t, const Tensor<Q>& c) {
        typedef TENSOR_RESULT_TYPE(T,Q) resultT;
        TENSOR_ASSERT(c.ndim() == 2,"second argument must be a matrix",c.ndim(),&c);
        if (c.dim(0)==c.dim(1) && t.iscontiguous() && c.iscontiguous()) {
            Tensor<resultT> result(t.ndim(),t.dims(),false);
            Tensor<resultT> work(t.ndim(),t.dims(),false);
            return fast_transform(t, c, result, work);
        }
        else {
            Tensor<resultT> result = t;
            for (long i=0; i<t.ndim(); ++i) {
                result = inner(result,c,0,0);
            }
            return result;
        }
    }


    template <class T, class Q>
    Tensor<TENSOR_RESULT_TYPE(T,Q)> general_transform(const Tensor<T>& t, const Tensor<Q> c[]) {
        typedef TENSOR_RESULT_TYPE(T,Q) resultT;
        Tensor<resultT> result = t;
        for (long i=0; i<t.ndim(); ++i) {
            result = inner(result,c[i],0,0);
        }
        return result;
    }


    template <class T, class Q>
    Tensor< TENSOR_RESULT_TYPE(T,Q) >& fast_transform(const Tensor<T>& t, const Tensor<Q>& c,  Tensor< TENSOR_RESULT_TYPE(T,Q) >& result,
            Tensor< TENSOR_RESULT_TYPE(T,Q) >& workspace) {
        typedef  TENSOR_RESULT_TYPE(T,Q) resultT;
        const Q *pc=c.ptr();
        resultT *t0=workspace.ptr(), *t1=result.ptr();
        if (t.ndim()&1) {
            t0 = result.ptr();
            t1 = workspace.ptr();
        }

        long dimj = c.dim(1);
        long dimi = 1;
        for (int n=1; n<t.ndim(); ++n) dimi *= dimj;

#ifdef AVX_MTXMQ_TEST
        // The new AVX code is smokin' fast and has no restrictions
            mTxmq(dimi, dimj, dimj, t0, t.ptr(), pc);
            for (int n=1; n<t.ndim(); ++n) {
                mTxmq(dimi, dimj, dimj, t1, t0, pc);
                std::swap(t0,t1);
            }
#elif HAVE_IBMBGQ
            long nij = dimi*dimj;
            if (IS_UNALIGNED(pc) || IS_UNALIGNED(t0) || IS_UNALIGNED(t1)) {
                for (long i=0; i<nij; ++i) t0[i] = 0.0;
                mTxm(dimi, dimj, dimj, t0, t.ptr(), pc);
                for (int n=1; n<t.ndim(); ++n) {
                    for (long i=0; i<nij; ++i) t1[i] = 0.0;
                    mTxm(dimi, dimj, dimj, t1, t0, pc);
                    std::swap(t0,t1);
                }
            }
            else {
                mTxmq_padding(dimi, dimj, dimj, dimj, t0, t.ptr(), pc);
                for (int n=1; n<t.ndim(); ++n) {
                    mTxmq_padding(dimi, dimj, dimj, dimj, t1, t0, pc);
                    std::swap(t0,t1);
                }
            }
#else
        long nij = dimi*dimj;
        if (IS_ODD(dimi) || IS_ODD(dimj) ||
                IS_UNALIGNED(pc) || IS_UNALIGNED(t0) || IS_UNALIGNED(t1)) {
            for (long i=0; i<nij; ++i) t0[i] = 0.0;
            mTxm(dimi, dimj, dimj, t0, t.ptr(), pc);
            for (int n=1; n<t.ndim(); ++n) {
                for (long i=0; i<nij; ++i) t1[i] = 0.0;
                mTxm(dimi, dimj, dimj, t1, t0, pc);
                std::swap(t0,t1);
            }
        }
        else {
         mTxmq(dimi, dimj, dimj, t0, t.ptr(), pc);
         for (int n=1; n<t.ndim(); ++n) {
             mTxmq(dimi, dimj, dimj, t1, t0, pc);
             std::swap(t0,t1);
         }
        }
#endif

        return result;
    }


    template <class T>
    Tensor< typename Tensor<T>::scalar_type > abs(const Tensor<T>& t) {
        typedef typename Tensor<T>::scalar_type scalar_type;
        Tensor<scalar_type> result(t.ndim(),t.dims(),false);
        BINARY_OPTIMIZED_ITERATOR(scalar_type,result,const T,t,*_p0 = std::abs(*_p1));
        return result;
    }


    template <class T>
    Tensor< typename Tensor<T>::scalar_type > arg(const Tensor<T>& t) {
        typedef typename Tensor<T>::scalar_type scalar_type;
        Tensor<scalar_type> result(t.ndim(),t.dims(),false);
        BINARY_OPTIMIZED_ITERATOR(scalar_type,result,T,t,*_p0 = std::arg(*_p1));
        return result;
    }


    template <class T>
    Tensor< typename Tensor<T>::scalar_type > real(const Tensor<T>& t) {
        typedef typename Tensor<T>::scalar_type scalar_type;
        Tensor<scalar_type> result(t.ndim(),t.dims(),false);
        BINARY_OPTIMIZED_ITERATOR(scalar_type,result,const T,t,*_p0 = std::real(*_p1));
        return result;
    }


    template <class T>
    Tensor< typename Tensor<T>::scalar_type > imag(const Tensor<T>& t) {
        typedef typename Tensor<T>::scalar_type scalar_type;
        Tensor<scalar_type> result(t.ndim(),t.dims(),false);
        BINARY_OPTIMIZED_ITERATOR(scalar_type,result,const T,t,*_p0 = std::imag(*_p1));
        return result;
    }


    template <class T>
    Tensor<T> conj(const Tensor<T>& t) {
        Tensor<T> result(t.ndim(),t.dims(),false);
        BINARY_OPTIMIZED_ITERATOR(T,result,const T,t,*_p0 = std::conj(*_p1));
        return result;
    }
}

#endif // MADNESS_TENSOR_TENSOR_H__INCLUDED