BigMatrix.h

Go to the documentation of this file.

#ifndef SimTK_SIMMATRIX_BIGMATRIX_H_
#define SimTK_SIMMATRIX_BIGMATRIX_H_

/* -------------------------------------------------------------------------- *
 *                      SimTK Core: SimTK Simmatrix(tm)                       *
 * -------------------------------------------------------------------------- *
 * This is part of the SimTK Core biosimulation toolkit originating from      *
 * Simbios, the NIH National Center for Physics-Based Simulation of           *
 * Biological Structures at Stanford, funded under the NIH Roadmap for        *
 * Medical Research, grant U54 GM072970. See https://simtk.org.               *
 *                                                                            *
 * Portions copyright (c) 2005-9 Stanford University and the Authors.         *
 * Authors: Michael Sherman                                                   *
 * Contributors:                                                              *
 *                                                                            *
 * Permission is hereby granted, free of charge, to any person obtaining a    *
 * copy of this software and associated documentation files (the "Software"), *
 * to deal in the Software without restriction, including without limitation  *
 * the rights to use, copy, modify, merge, publish, distribute, sublicense,   *
 * and/or sell copies of the Software, and to permit persons to whom the      *
 * Software is furnished to do so, subject to the following conditions:       *
 *                                                                            *
 * The above copyright notice and this permission notice shall be included in *
 * all copies or substantial portions of the Software.                        *
 *                                                                            *
 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR *
 * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,   *
 * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL    *
 * THE AUTHORS, CONTRIBUTORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM,    *
 * DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR      *
 * OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE  *
 * USE OR OTHER DEALINGS IN THE SOFTWARE.                                     *
 * -------------------------------------------------------------------------- */

// TODO: matrix expression templates for delaying operator execution.

#include "SimTKcommon/Scalar.h"
#include "SimTKcommon/SmallMatrix.h"

#include "SimTKcommon/internal/MatrixHelper.h"
#include "SimTKcommon/internal/MatrixCharacteristics.h"

#include <iostream>
#include <cassert>
#include <complex>
#include <cstddef>
#include <limits>

namespace SimTK {

template <class ELT>    class MatrixBase;
template <class ELT>    class VectorBase;
template <class ELT>    class RowVectorBase;

template <class T=Real> class MatrixView_;
template <class T=Real> class DeadMatrixView_;
template <class T=Real> class Matrix_;
template <class T=Real> class DeadMatrix_;

template <class T=Real> class VectorView_;
template <class T=Real> class DeadVectorView_;
template <class T=Real> class Vector_;
template <class T=Real> class DeadVector_;

template <class T=Real> class RowVectorView_;
template <class T=Real> class DeadRowVectorView_;
template <class T=Real> class RowVector_;
template <class T=Real> class DeadRowVector_;

template <class ELT, class VECTOR_CLASS> class VectorIterator;

//  -------------------------------- MatrixBase --------------------------------
//  ----------------------------------------------------------------------------
template <class ELT> class MatrixBase {  
public:
    // These typedefs are handy, but despite appearances you cannot 
    // treat a MatrixBase as a composite numerical type. That is,
    // CNT<MatrixBase> will not compile, or if it does it won't be
    // meaningful.

    typedef ELT                                 E;
    typedef typename CNT<E>::TNeg               ENeg;
    typedef typename CNT<E>::TWithoutNegator    EWithoutNegator;
    typedef typename CNT<E>::TReal              EReal;
    typedef typename CNT<E>::TImag              EImag;
    typedef typename CNT<E>::TComplex           EComplex;
    typedef typename CNT<E>::THerm              EHerm;       
    typedef typename CNT<E>::TPosTrans          EPosTrans;

    typedef typename CNT<E>::TAbs               EAbs;
    typedef typename CNT<E>::TStandard          EStandard;
    typedef typename CNT<E>::TInvert            EInvert;
    typedef typename CNT<E>::TNormalize         ENormalize;
    typedef typename CNT<E>::TSqHermT           ESqHermT;
    typedef typename CNT<E>::TSqTHerm           ESqTHerm;

    typedef typename CNT<E>::Scalar             EScalar;
    typedef typename CNT<E>::Number             ENumber;
    typedef typename CNT<E>::StdNumber          EStdNumber;
    typedef typename CNT<E>::Precision          EPrecision;
    typedef typename CNT<E>::ScalarNormSq       EScalarNormSq;

    typedef EScalar    Scalar;        // the underlying Scalar type
    typedef ENumber    Number;        // negator removed from Scalar
    typedef EStdNumber StdNumber;     // conjugate goes to complex
    typedef EPrecision Precision;     // complex removed from StdNumber
    typedef EScalarNormSq  ScalarNormSq;      // type of scalar^2

    typedef MatrixBase<E>                T;
    typedef MatrixBase<ENeg>             TNeg;
    typedef MatrixBase<EWithoutNegator>  TWithoutNegator;
    typedef MatrixBase<EReal>            TReal;
    typedef MatrixBase<EImag>            TImag;
    typedef MatrixBase<EComplex>         TComplex;
    typedef MatrixBase<EHerm>            THerm; 
    typedef MatrixBase<E>                TPosTrans;

    typedef MatrixBase<EAbs>             TAbs;
    typedef MatrixBase<EStandard>        TStandard;
    typedef MatrixBase<EInvert>          TInvert;
    typedef MatrixBase<ENormalize>       TNormalize;
    typedef MatrixBase<ESqHermT>         TSqHermT;  // ~Mat*Mat
    typedef MatrixBase<ESqTHerm>         TSqTHerm;  // Mat*~Mat

    const MatrixCommitment& getCharacterCommitment() const {return helper.getCharacterCommitment();}
    const MatrixCharacter& getMatrixCharacter()     const {return helper.getMatrixCharacter();}

    void commitTo(const MatrixCommitment& mc)
    {   helper.commitTo(mc); }

    // This gives the resulting matrix type when (m(i,j) op P) is applied to each element.
    // It will have element types which are the regular composite result of E op P.
    template <class P> struct EltResult { 
        typedef MatrixBase<typename CNT<E>::template Result<P>::Mul> Mul;
        typedef MatrixBase<typename CNT<E>::template Result<P>::Dvd> Dvd;
        typedef MatrixBase<typename CNT<E>::template Result<P>::Add> Add;
        typedef MatrixBase<typename CNT<E>::template Result<P>::Sub> Sub;
    };

    int  nrow() const {return helper.nrow();}
    int  ncol() const {return helper.ncol();}

    ptrdiff_t nelt() const {return helper.nelt();}

    bool isResizeable() const {return getCharacterCommitment().isResizeable();}

    enum { 
        NScalarsPerElement    = CNT<E>::NActualScalars,
        CppNScalarsPerElement = sizeof(E) / sizeof(Scalar)
    };
  
    MatrixBase() : helper(NScalarsPerElement,CppNScalarsPerElement) {}

    MatrixBase(int m, int n) 
    :   helper(NScalarsPerElement,CppNScalarsPerElement,MatrixCommitment(),m,n) {}

    explicit MatrixBase(const MatrixCommitment& commitment) 
    :   helper(NScalarsPerElement,CppNScalarsPerElement,commitment) {}


    MatrixBase(const MatrixCommitment& commitment, int m, int n) 
    :   helper(NScalarsPerElement,CppNScalarsPerElement,commitment,m,n) {}

    MatrixBase(const MatrixBase& b)
      : helper(b.helper.getCharacterCommitment(), 
               b.helper, typename MatrixHelper<Scalar>::DeepCopy()) { }

    MatrixBase(const TNeg& b)
      : helper(b.helper.getCharacterCommitment(),
               b.helper, typename MatrixHelper<Scalar>::DeepCopy()) { }
    
    MatrixBase& copyAssign(const MatrixBase& b) {
        helper.copyAssign(b.helper);
        return *this;
    }
    MatrixBase& operator=(const MatrixBase& b) { return copyAssign(b); }


    MatrixBase& viewAssign(const MatrixBase& src) {
        helper.writableViewAssign(const_cast<MatrixHelper<Scalar>&>(src.helper));
        return *this;
    }

    // default destructor

    MatrixBase(const MatrixCommitment& commitment, int m, int n, const ELT& initialValue) 
    :   helper(NScalarsPerElement, CppNScalarsPerElement, commitment, m, n)
    {   helper.fillWith(reinterpret_cast<const Scalar*>(&initialValue)); }  

    MatrixBase(const MatrixCommitment& commitment, int m, int n, 
               const ELT* cppInitialValuesByRow) 
    :   helper(NScalarsPerElement, CppNScalarsPerElement, commitment, m, n)
    {   helper.copyInByRowsFromCpp(reinterpret_cast<const Scalar*>(cppInitialValuesByRow)); }
     

    MatrixBase(const MatrixCommitment& commitment, 
               const MatrixCharacter&  character, 
               int spacing, const Scalar* data) // read only data
    :   helper(NScalarsPerElement, CppNScalarsPerElement, 
               commitment, character, spacing, data) {}  

    MatrixBase(const MatrixCommitment& commitment, 
               const MatrixCharacter&  character, 
               int spacing, Scalar* data) // writable data
    :   helper(NScalarsPerElement, CppNScalarsPerElement, 
               commitment, character, spacing, data) {}  
        
    // Create a new MatrixBase from an existing helper. Both shallow and deep copies are possible.
    MatrixBase(const MatrixCommitment& commitment, 
               MatrixHelper<Scalar>&   source, 
               const typename MatrixHelper<Scalar>::ShallowCopy& shallow) 
    :   helper(commitment, source, shallow) {}
    MatrixBase(const MatrixCommitment&      commitment, 
               const MatrixHelper<Scalar>&  source, 
               const typename MatrixHelper<Scalar>::ShallowCopy& shallow) 
    :   helper(commitment, source, shallow) {}
    MatrixBase(const MatrixCommitment&      commitment, 
               const MatrixHelper<Scalar>&  source, 
               const typename MatrixHelper<Scalar>::DeepCopy& deep)    
    :   helper(commitment, source, deep) {}

    void clear() {helper.clear();}

    MatrixBase& operator*=(const StdNumber& t)  { helper.scaleBy(t);              return *this; }
    MatrixBase& operator/=(const StdNumber& t)  { helper.scaleBy(StdNumber(1)/t); return *this; }
    MatrixBase& operator+=(const MatrixBase& r) { helper.addIn(r.helper);         return *this; }
    MatrixBase& operator-=(const MatrixBase& r) { helper.subIn(r.helper);         return *this; }  

    template <class EE> MatrixBase(const MatrixBase<EE>& b)
      : helper(MatrixCommitment(),b.helper, typename MatrixHelper<Scalar>::DeepCopy()) { }

    template <class EE> MatrixBase& operator=(const MatrixBase<EE>& b) 
      { helper = b.helper; return *this; }
    template <class EE> MatrixBase& operator+=(const MatrixBase<EE>& b) 
      { helper.addIn(b.helper); return *this; }
    template <class EE> MatrixBase& operator-=(const MatrixBase<EE>& b) 
      { helper.subIn(b.helper); return *this; }

    MatrixBase& operator=(const ELT& t) { 
        setToZero(); updDiag().setTo(t); 
        return *this;
    }

    template <class S> inline MatrixBase&
    scalarAssign(const S& s) {
        setToZero(); updDiag().setTo(s);
        return *this;
    }

    template <class S> inline MatrixBase&
    scalarAddInPlace(const S& s) {
        updDiag().elementwiseAddScalarInPlace(s);
    }


    template <class S> inline MatrixBase&
    scalarSubtractInPlace(const S& s) {
        updDiag().elementwiseSubtractScalarInPlace(s);
    }

    template <class S> inline MatrixBase&
    scalarSubtractFromLeftInPlace(const S& s) {
        negateInPlace();
        updDiag().elementwiseAddScalarInPlace(s); // yes, add
    }

    template <class S> inline MatrixBase&
    scalarMultiplyInPlace(const S&);

    template <class S> inline MatrixBase&
    scalarMultiplyFromLeftInPlace(const S&);

    template <class S> inline MatrixBase&
    scalarDivideInPlace(const S&);

    template <class S> inline MatrixBase&
    scalarDivideFromLeftInPlace(const S&);


    template <class EE> inline MatrixBase& 
    rowScaleInPlace(const VectorBase<EE>&);

    template <class EE> inline void 
    rowScale(const VectorBase<EE>& r, typename EltResult<EE>::Mul& out) const;

    template <class EE> inline typename EltResult<EE>::Mul 
    rowScale(const VectorBase<EE>& r) const {
        typename EltResult<EE>::Mul out(nrow(), ncol()); rowScale(r,out); return out;
    }

    template <class EE> inline MatrixBase& 
    colScaleInPlace(const VectorBase<EE>&);

    template <class EE> inline void 
    colScale(const VectorBase<EE>& c, typename EltResult<EE>::Mul& out) const;

    template <class EE> inline typename EltResult<EE>::Mul
    colScale(const VectorBase<EE>& c) const {
        typename EltResult<EE>::Mul out(nrow(), ncol()); colScale(c,out); return out;
    }


    template <class ER, class EC> inline MatrixBase& 
    rowAndColScaleInPlace(const VectorBase<ER>& r, const VectorBase<EC>& c);

    template <class ER, class EC> inline void 
    rowAndColScale(const VectorBase<ER>& r, const VectorBase<EC>& c, 
                   typename EltResult<typename VectorBase<ER>::template EltResult<EC>::Mul>::Mul& out) const;

    template <class ER, class EC> inline typename EltResult<typename VectorBase<ER>::template EltResult<EC>::Mul>::Mul
    rowAndColScale(const VectorBase<ER>& r, const VectorBase<EC>& c) const {
        typename EltResult<typename VectorBase<ER>::template EltResult<EC>::Mul>::Mul 
            out(nrow(), ncol()); 
        rowAndColScale(r,c,out); return out;
    }

    template <class S> inline MatrixBase&
    elementwiseAssign(const S& s);

    MatrixBase& elementwiseInvertInPlace();

    void elementwiseInvert(MatrixBase<typename CNT<E>::TInvert>& out) const;

    MatrixBase<typename CNT<E>::TInvert> elementwiseInvert() const {
        MatrixBase<typename CNT<E>::TInvert> out(nrow(), ncol());
        elementwiseInvert(out);
        return out;
    }

    template <class S> inline MatrixBase&
    elementwiseAddScalarInPlace(const S& s);

    template <class S> inline void
    elementwiseAddScalar(const S& s, typename EltResult<S>::Add&) const;

    template <class S> inline typename EltResult<S>::Add
    elementwiseAddScalar(const S& s) const {
        typename EltResult<S>::Add out(nrow(), ncol());
        elementwiseAddScalar(s,out);
        return out;
    }

    template <class S> inline MatrixBase&
    elementwiseSubtractScalarInPlace(const S& s);

    template <class S> inline void
    elementwiseSubtractScalar(const S& s, typename EltResult<S>::Sub&) const;

    template <class S> inline typename EltResult<S>::Sub
    elementwiseSubtractScalar(const S& s) const {
        typename EltResult<S>::Sub out(nrow(), ncol());
        elementwiseSubtractScalar(s,out);
        return out;
    }

    template <class S> inline MatrixBase&
    elementwiseSubtractFromScalarInPlace(const S& s);

    template <class S> inline void
    elementwiseSubtractFromScalar(
        const S&, 
        typename MatrixBase<S>::template EltResult<E>::Sub&) const;

    template <class S> inline typename MatrixBase<S>::template EltResult<E>::Sub
    elementwiseSubtractFromScalar(const S& s) const {
        typename MatrixBase<S>::template EltResult<E>::Sub out(nrow(), ncol());
        elementwiseSubtractFromScalar<S>(s,out);
        return out;
    }

    template <class EE> inline MatrixBase& 
    elementwiseMultiplyInPlace(const MatrixBase<EE>&);

    template <class EE> inline void 
    elementwiseMultiply(const MatrixBase<EE>&, typename EltResult<EE>::Mul&) const;

    template <class EE> inline typename EltResult<EE>::Mul 
    elementwiseMultiply(const MatrixBase<EE>& m) const {
        typename EltResult<EE>::Mul out(nrow(), ncol()); 
        elementwiseMultiply<EE>(m,out); 
        return out;
    }

    template <class EE> inline MatrixBase& 
    elementwiseMultiplyFromLeftInPlace(const MatrixBase<EE>&);

    template <class EE> inline void 
    elementwiseMultiplyFromLeft(
        const MatrixBase<EE>&, 
        typename MatrixBase<EE>::template EltResult<E>::Mul&) const;

    template <class EE> inline typename MatrixBase<EE>::template EltResult<E>::Mul 
    elementwiseMultiplyFromLeft(const MatrixBase<EE>& m) const {
        typename EltResult<EE>::Mul out(nrow(), ncol()); 
        elementwiseMultiplyFromLeft<EE>(m,out); 
        return out;
    }

    template <class EE> inline MatrixBase& 
    elementwiseDivideInPlace(const MatrixBase<EE>&);

    template <class EE> inline void 
    elementwiseDivide(const MatrixBase<EE>&, typename EltResult<EE>::Dvd&) const;

    template <class EE> inline typename EltResult<EE>::Dvd 
    elementwiseDivide(const MatrixBase<EE>& m) const {
        typename EltResult<EE>::Dvd out(nrow(), ncol()); 
        elementwiseDivide<EE>(m,out); 
        return out;
    }

    template <class EE> inline MatrixBase& 
    elementwiseDivideFromLeftInPlace(const MatrixBase<EE>&);

    template <class EE> inline void 
    elementwiseDivideFromLeft(
        const MatrixBase<EE>&,
        typename MatrixBase<EE>::template EltResult<E>::Dvd&) const;

    template <class EE> inline typename MatrixBase<EE>::template EltResult<EE>::Dvd 
    elementwiseDivideFromLeft(const MatrixBase<EE>& m) const {
        typename MatrixBase<EE>::template EltResult<E>::Dvd out(nrow(), ncol()); 
        elementwiseDivideFromLeft<EE>(m,out); 
        return out;
    }

    MatrixBase& setTo(const ELT& t) {helper.fillWith(reinterpret_cast<const Scalar*>(&t)); return *this;}
    MatrixBase& setToNaN() {helper.fillWithScalar(CNT<StdNumber>::getNaN()); return *this;}
    MatrixBase& setToZero() {helper.fillWithScalar(StdNumber(0)); return *this;}
 
    // View creating operators. TODO: these should be DeadMatrixViews  
    inline RowVectorView_<ELT> row(int i) const;   // select a row
    inline RowVectorView_<ELT> updRow(int i);
    inline VectorView_<ELT>    col(int j) const;   // select a column
    inline VectorView_<ELT>    updCol(int j);

    RowVectorView_<ELT> operator[](int i) const {return row(i);}
    RowVectorView_<ELT> operator[](int i)       {return updRow(i);}
    VectorView_<ELT>    operator()(int j) const {return col(j);}
    VectorView_<ELT>    operator()(int j)       {return updCol(j);}
     
    // Select a block.
    inline MatrixView_<ELT> block(int i, int j, int m, int n) const;
    inline MatrixView_<ELT> updBlock(int i, int j, int m, int n);

    MatrixView_<ELT> operator()(int i, int j, int m, int n) const
      { return block(i,j,m,n); }
    MatrixView_<ELT> operator()(int i, int j, int m, int n)
      { return updBlock(i,j,m,n); }
 
    // Hermitian transpose.
    inline MatrixView_<EHerm> transpose() const;
    inline MatrixView_<EHerm> updTranspose();

    MatrixView_<EHerm> operator~() const {return transpose();}
    MatrixView_<EHerm> operator~()       {return updTranspose();}

    // Select matrix diagonal (of largest leading square if rectangular).
    inline VectorView_<ELT> diag() const;
    inline VectorView_<ELT> updDiag();

    // Create a view of the real or imaginary elements. TODO
    //inline MatrixView_<EReal> real() const;
    //inline MatrixView_<EReal> updReal();
    //inline MatrixView_<EImag> imag() const;
    //inline MatrixView_<EImag> updImag();

    // Overload "real" and "imag" for both read and write as a nod to convention. TODO
    //MatrixView_<EReal> real() {return updReal();}
    //MatrixView_<EReal> imag() {return updImag();}

    // TODO: this routine seems ill-advised but I need it for the IVM port at the moment
    TInvert invert() const {  // return a newly-allocated inverse; dump negator 
        TInvert m(*this);
        m.helper.invertInPlace();
        return m;   // TODO - bad: makes an extra copy
    }

    void invertInPlace() {helper.invertInPlace();}

    void dump(const char* msg=0) const {
        helper.dump(msg);
    }


    // This routine is useful for implementing friendlier Matrix expressions and operators.
    // It maps closely to the Level-3 BLAS family of pxxmm() routines like sgemm(). The
    // operation performed assumes that "this" is the result, and that "this" has 
    // already been sized correctly to receive the result. We'll compute
    //     this = beta*this + alpha*A*B
    // If beta is 0 then "this" can be uninitialized. If alpha is 0 we promise not
    // to look at A or B. The routine can work efficiently even if A and/or B are transposed
    // by their views, so an expression like
    //        C += s * ~A * ~B
    // can be performed with the single equivalent call
    //        C.matmul(1., s, Tr(A), Tr(B))
    // where Tr(A) indicates a transposed view of the original A.
    // The ultimate efficiency of this call depends on the data layout and views which
    // are used for the three matrices.
    // NOTE: neither A nor B can be the same matrix as 'this', nor views of the same data
    // which would expose elements of 'this' that will be modified by this operation.
    template <class ELT_A, class ELT_B>
    MatrixBase& matmul(const StdNumber& beta,   // applied to 'this'
                       const StdNumber& alpha, const MatrixBase<ELT_A>& A, const MatrixBase<ELT_B>& B)
    {
        helper.matmul(beta,alpha,A.helper,B.helper);
        return *this;
    }

    const ELT& getElt(int i, int j) const { return *reinterpret_cast<const ELT*>(helper.getElt(i,j)); }
    ELT&       updElt(int i, int j)       { return *reinterpret_cast<      ELT*>(helper.updElt(i,j)); }

    const ELT& operator()(int i, int j) const {return getElt(i,j);}
    ELT&       operator()(int i, int j)       {return updElt(i,j);}

    void getAnyElt(int i, int j, ELT& value) const
    {   helper.getAnyElt(i,j,reinterpret_cast<Scalar*>(&value)); }
    ELT getAnyElt(int i, int j) const {ELT e; getAnyElt(i,j,e); return e;}

    // TODO: very slow! Should be optimized at least for the case
    //       where ELT is a Scalar.
    ScalarNormSq scalarNormSqr() const {
        const int nr=nrow(), nc=ncol();
        ScalarNormSq sum(0);
        for(int j=0;j<nc;++j) 
            for (int i=0; i<nr; ++i)
                sum += CNT<E>::scalarNormSqr((*this)(i,j));
        return sum;
    }

    // TODO: very slow! Should be optimized at least for the case
    //       where ELT is a Scalar.
    void abs(TAbs& mabs) const {
        const int nr=nrow(), nc=ncol();
        mabs.resize(nr,nc);
        for(int j=0;j<nc;++j) 
            for (int i=0; i<nr; ++i)
                mabs(i,j) = CNT<E>::abs((*this)(i,j));
    }

    TAbs abs() const { TAbs mabs; abs(mabs); return mabs; }

    TStandard standardize() const {
        const int nr=nrow(), nc=ncol();
        TStandard mstd(nr, nc);
        for(int j=0;j<nc;++j) 
            for (int i=0; i<nr; ++i)
                mstd(i,j) = CNT<E>::standardize((*this)(i,j));
        return mstd;
    }

    ScalarNormSq normSqr() const { return scalarNormSqr(); }
    // TODO -- not good; unnecessary overflow
    typename CNT<ScalarNormSq>::TSqrt 
        norm() const { return CNT<ScalarNormSq>::sqrt(scalarNormSqr()); }

    typename CNT<ScalarNormSq>::TSqrt 
    normRMS() const {
        if (!CNT<ELT>::IsScalar)
            SimTK_THROW1(Exception::Cant, "normRMS() only defined for scalar elements");
        if (nelt() == 0)
            return typename CNT<ScalarNormSq>::TSqrt(0);
        return CNT<ScalarNormSq>::sqrt(scalarNormSqr()/nelt());
    }

    RowVectorBase<ELT> sum() const {
        const int cols = ncol();
        RowVectorBase<ELT> row(cols);
        for (int i = 0; i < cols; ++i)
            helper.colSum(i, reinterpret_cast<Scalar*>(&row[i]));
        return row;
    }

    //TODO: make unary +/- return a self-view so they won't reallocate?
    const MatrixBase& operator+() const {return *this; }
    const TNeg&       negate()    const {return *reinterpret_cast<const TNeg*>(this); }
    TNeg&             updNegate()       {return *reinterpret_cast<TNeg*>(this); }

    const TNeg&       operator-() const {return negate();}
    TNeg&             operator-()       {return updNegate();}

    MatrixBase& negateInPlace() {(*this) *= EPrecision(-1); return *this;}
 
    MatrixBase& resize(int m, int n)     { helper.resize(m,n); return *this; }
    MatrixBase& resizeKeep(int m, int n) { helper.resizeKeep(m,n); return *this; }

    // This prevents shape changes in a Matrix that would otherwise allow it. No harm if is
    // are called on a Matrix that is locked already; it always succeeds.
    void lockShape() {helper.lockShape();}

    // This allows shape changes again for a Matrix which was constructed to allow them
    // but had them locked with the above routine. No harm if this is called on a Matrix
    // that is already unlocked, but it is not allowed to call this on a Matrix which
    // *never* allowed resizing. An exception will be thrown in that case.
    void unlockShape() {helper.unlockShape();}
    
    // An assortment of handy conversions
    const MatrixView_<ELT>& getAsMatrixView() const { return *reinterpret_cast<const MatrixView_<ELT>*>(this); }
    MatrixView_<ELT>&       updAsMatrixView()       { return *reinterpret_cast<      MatrixView_<ELT>*>(this); } 
    const Matrix_<ELT>&     getAsMatrix()     const { return *reinterpret_cast<const Matrix_<ELT>*>(this); }
    Matrix_<ELT>&           updAsMatrix()           { return *reinterpret_cast<      Matrix_<ELT>*>(this); }
         
    const VectorView_<ELT>& getAsVectorView() const 
      { assert(ncol()==1); return *reinterpret_cast<const VectorView_<ELT>*>(this); }
    VectorView_<ELT>&       updAsVectorView()       
      { assert(ncol()==1); return *reinterpret_cast<      VectorView_<ELT>*>(this); } 
    const Vector_<ELT>&     getAsVector()     const 
      { assert(ncol()==1); return *reinterpret_cast<const Vector_<ELT>*>(this); }
    Vector_<ELT>&           updAsVector()           
      { assert(ncol()==1); return *reinterpret_cast<      Vector_<ELT>*>(this); }
    const VectorBase<ELT>& getAsVectorBase() const 
      { assert(ncol()==1); return *reinterpret_cast<const VectorBase<ELT>*>(this); }
    VectorBase<ELT>&       updAsVectorBase()       
      { assert(ncol()==1); return *reinterpret_cast<      VectorBase<ELT>*>(this); } 
                
    const RowVectorView_<ELT>& getAsRowVectorView() const 
      { assert(nrow()==1); return *reinterpret_cast<const RowVectorView_<ELT>*>(this); }
    RowVectorView_<ELT>&       updAsRowVectorView()       
      { assert(nrow()==1); return *reinterpret_cast<      RowVectorView_<ELT>*>(this); } 
    const RowVector_<ELT>&     getAsRowVector()     const 
      { assert(nrow()==1); return *reinterpret_cast<const RowVector_<ELT>*>(this); }
    RowVector_<ELT>&           updAsRowVector()           
      { assert(nrow()==1); return *reinterpret_cast<      RowVector_<ELT>*>(this); }        
    const RowVectorBase<ELT>& getAsRowVectorBase() const 
      { assert(nrow()==1); return *reinterpret_cast<const RowVectorBase<ELT>*>(this); }
    RowVectorBase<ELT>&       updAsRowVectorBase()       
      { assert(nrow()==1); return *reinterpret_cast<      RowVectorBase<ELT>*>(this); } 

    // Access to raw data. We have to return the raw data
    // pointer as pointer-to-scalar because we may pack the elements tighter
    // than a C++ array would.

    int getNScalarsPerElement()  const {return NScalarsPerElement;}

    int getPackedSizeofElement() const {return NScalarsPerElement*sizeof(Scalar);}

    bool hasContiguousData() const {return helper.hasContiguousData();}
    ptrdiff_t getContiguousScalarDataLength() const {
        return helper.getContiguousDataLength();
    }
    const Scalar* getContiguousScalarData() const {
        return helper.getContiguousData();
    }
    Scalar* updContiguousScalarData() {
        return helper.updContiguousData();
    }
    void replaceContiguousScalarData(Scalar* newData, ptrdiff_t length, bool takeOwnership) {
        helper.replaceContiguousData(newData,length,takeOwnership);
    }
    void replaceContiguousScalarData(const Scalar* newData, ptrdiff_t length) {
        helper.replaceContiguousData(newData,length);
    }
    void swapOwnedContiguousScalarData(Scalar* newData, ptrdiff_t length, Scalar*& oldData) {
        helper.swapOwnedContiguousData(newData,length,oldData);
    }

    explicit MatrixBase(MatrixHelperRep<Scalar>* hrep) : helper(hrep) {}

protected:
    const MatrixHelper<Scalar>& getHelper() const {return helper;}
    MatrixHelper<Scalar>&       updHelper()       {return helper;}

private:
    MatrixHelper<Scalar> helper; // this is just one pointer

    template <class EE> friend class MatrixBase;
};



//  -------------------------------- VectorBase --------------------------------
//  ----------------------------------------------------------------------------
template <class ELT> class VectorBase : public MatrixBase<ELT> {
    typedef MatrixBase<ELT>                             Base;
    typedef typename CNT<ELT>::Scalar                   Scalar;
    typedef typename CNT<ELT>::Number                   Number;
    typedef typename CNT<ELT>::StdNumber                StdNumber;
    typedef VectorBase<ELT>                             T;
    typedef VectorBase<typename CNT<ELT>::TAbs>         TAbs;
    typedef VectorBase<typename CNT<ELT>::TNeg>         TNeg;
    typedef RowVectorView_<typename CNT<ELT>::THerm>    THerm;
public:  
    //  ------------------------------------------------------------------------

    explicit VectorBase(int m=0) : Base(MatrixCommitment::Vector(), m, 1) {}

    VectorBase(const VectorBase& source) : Base(source) {}

    VectorBase(const TNeg& source) : Base(source) {}

    VectorBase(int m, const ELT& initialValue)
    :   Base(MatrixCommitment::Vector(),m,1,initialValue) {}  

    VectorBase(int m, const ELT* cppInitialValues)
    :   Base(MatrixCommitment::Vector(),m,1,cppInitialValues) {}

    //  ------------------------------------------------------------------------

    VectorBase(int m, int stride, const Scalar* s)
    :   Base(MatrixCommitment::Vector(m), MatrixCharacter::Vector(m),stride,s) { }
    VectorBase(int m, int stride, Scalar* s)
    :   Base(MatrixCommitment::Vector(m), MatrixCharacter::Vector(m),stride,s) { }
        
    //  ------------------------------------------------------------------------

    VectorBase(MatrixHelper<Scalar>& h, const typename MatrixHelper<Scalar>::ShallowCopy& s) 
    :   Base(MatrixCommitment::Vector(), h,s) { }
    VectorBase(const MatrixHelper<Scalar>& h, const typename MatrixHelper<Scalar>::ShallowCopy& s) 
    :   Base(MatrixCommitment::Vector(), h,s) { }
    VectorBase(const MatrixHelper<Scalar>& h, const typename MatrixHelper<Scalar>::DeepCopy& d)    
    :   Base(MatrixCommitment::Vector(), h,d) { }

    // This gives the resulting vector type when (v[i] op P) is applied to each element.
    // It will have element types which are the regular composite result of ELT op P.
    template <class P> struct EltResult { 
        typedef VectorBase<typename CNT<ELT>::template Result<P>::Mul> Mul;
        typedef VectorBase<typename CNT<ELT>::template Result<P>::Dvd> Dvd;
        typedef VectorBase<typename CNT<ELT>::template Result<P>::Add> Add;
        typedef VectorBase<typename CNT<ELT>::template Result<P>::Sub> Sub;
    };

    VectorBase& operator=(const VectorBase& b) {
        Base::operator=(b); return *this;
    }

    // default destructor


    VectorBase& operator*=(const StdNumber& t)  { Base::operator*=(t); return *this; }
    VectorBase& operator/=(const StdNumber& t)  { Base::operator/=(t); return *this; }
    VectorBase& operator+=(const VectorBase& r) { Base::operator+=(r); return *this; }
    VectorBase& operator-=(const VectorBase& r) { Base::operator-=(r); return *this; }  


    template <class EE> VectorBase& operator=(const VectorBase<EE>& b) 
      { Base::operator=(b);  return *this; } 
    template <class EE> VectorBase& operator+=(const VectorBase<EE>& b) 
      { Base::operator+=(b); return *this; } 
    template <class EE> VectorBase& operator-=(const VectorBase<EE>& b) 
      { Base::operator-=(b); return *this; } 


    VectorBase& operator=(const ELT& t) { Base::setTo(t); return *this; }  

    template <class EE> VectorBase& rowScaleInPlace(const VectorBase<EE>& v)
      { Base::template rowScaleInPlace<EE>(v); return *this; }
    template <class EE> inline void rowScale(const VectorBase<EE>& v, typename EltResult<EE>::Mul& out) const
      { Base::rowScale(v,out); }
    template <class EE> inline typename EltResult<EE>::Mul rowScale(const VectorBase<EE>& v) const
      { typename EltResult<EE>::Mul out(nrow()); Base::rowScale(v,out); return out; }

    VectorBase& elementwiseInvertInPlace() {
        Base::elementwiseInvertInPlace();
        return *this;
    }

    void elementwiseInvert(VectorBase<typename CNT<ELT>::TInvert>& out) const {
        Base::elementwiseInvert(out);
    }

    VectorBase<typename CNT<ELT>::TInvert> elementwiseInvert() const {
        VectorBase<typename CNT<ELT>::TInvert> out(nrow());
        Base::elementwiseInvert(out);
        return out;
    }

        // elementwise multiply
    template <class EE> VectorBase& elementwiseMultiplyInPlace(const VectorBase<EE>& r)
      { Base::template elementwiseMultiplyInPlace<EE>(r); return *this; }
    template <class EE> inline void elementwiseMultiply(const VectorBase<EE>& v, typename EltResult<EE>::Mul& out) const
      { Base::template elementwiseMultiply<EE>(v,out); }
    template <class EE> inline typename EltResult<EE>::Mul elementwiseMultiply(const VectorBase<EE>& v) const
      { typename EltResult<EE>::Mul out(nrow()); Base::template elementwiseMultiply<EE>(v,out); return out; }

        // elementwise multiply from left
    template <class EE> VectorBase& elementwiseMultiplyFromLeftInPlace(const VectorBase<EE>& r)
      { Base::template elementwiseMultiplyFromLeftInPlace<EE>(r); return *this; }
    template <class EE> inline void 
    elementwiseMultiplyFromLeft(
        const VectorBase<EE>& v, 
        typename VectorBase<EE>::template EltResult<ELT>::Mul& out) const
    { 
        Base::template elementwiseMultiplyFromLeft<EE>(v,out);
    }
    template <class EE> inline typename VectorBase<EE>::template EltResult<ELT>::Mul 
    elementwiseMultiplyFromLeft(const VectorBase<EE>& v) const
    { 
        typename VectorBase<EE>::template EltResult<ELT>::Mul out(nrow()); 
        Base::template elementwiseMultiplyFromLeft<EE>(v,out); 
        return out;
    }

        // elementwise divide
    template <class EE> VectorBase& elementwiseDivideInPlace(const VectorBase<EE>& r)
      { Base::template elementwiseDivideInPlace<EE>(r); return *this; }
    template <class EE> inline void elementwiseDivide(const VectorBase<EE>& v, typename EltResult<EE>::Dvd& out) const
      { Base::template elementwiseDivide<EE>(v,out); }
    template <class EE> inline typename EltResult<EE>::Dvd elementwiseDivide(const VectorBase<EE>& v) const
      { typename EltResult<EE>::Dvd out(nrow()); Base::template elementwiseDivide<EE>(v,out); return out; }

        // elementwise divide from left
    template <class EE> VectorBase& elementwiseDivideFromLeftInPlace(const VectorBase<EE>& r)
      { Base::template elementwiseDivideFromLeftInPlace<EE>(r); return *this; }
    template <class EE> inline void 
    elementwiseDivideFromLeft(
        const VectorBase<EE>& v, 
        typename VectorBase<EE>::template EltResult<ELT>::Dvd& out) const
    { 
        Base::template elementwiseDivideFromLeft<EE>(v,out);
    }
    template <class EE> inline typename VectorBase<EE>::template EltResult<ELT>::Dvd 
    elementwiseDivideFromLeft(const VectorBase<EE>& v) const
    { 
        typename VectorBase<EE>::template EltResult<ELT>::Dvd out(nrow()); 
        Base::template elementwiseDivideFromLeft<EE>(v,out); 
        return out;
    }

    // Implicit conversions are allowed to Vector or Matrix, but not to RowVector.   
    operator const Vector_<ELT>&()     const { return *reinterpret_cast<const Vector_<ELT>*>(this); }
    operator       Vector_<ELT>&()           { return *reinterpret_cast<      Vector_<ELT>*>(this); }
    operator const VectorView_<ELT>&() const { return *reinterpret_cast<const VectorView_<ELT>*>(this); }
    operator       VectorView_<ELT>&()       { return *reinterpret_cast<      VectorView_<ELT>*>(this); }
    
    operator const Matrix_<ELT>&()     const { return *reinterpret_cast<const Matrix_<ELT>*>(this); }
    operator       Matrix_<ELT>&()           { return *reinterpret_cast<      Matrix_<ELT>*>(this); } 
    operator const MatrixView_<ELT>&() const { return *reinterpret_cast<const MatrixView_<ELT>*>(this); }
    operator       MatrixView_<ELT>&()       { return *reinterpret_cast<      MatrixView_<ELT>*>(this); } 


    // size() for Vectors is Base::nelt() but returns int instead of ptrdiff_t.
    int size() const { 
        assert(Base::nelt() <= (ptrdiff_t)std::numeric_limits<int>::max()); 
        assert(Base::ncol()==1);
        return (int)Base::nelt();
    }
    int       nrow() const {assert(Base::ncol()==1); return Base::nrow();}
    int       ncol() const {assert(Base::ncol()==1); return Base::ncol();}
    ptrdiff_t nelt() const {assert(Base::ncol()==1); return Base::nelt();}

    // Override MatrixBase operators to return the right shape
    TAbs abs() const {TAbs result; Base::abs(result); return result;}
    
    // Override MatrixBase indexing operators          
    const ELT& operator[](int i) const {return *reinterpret_cast<const ELT*>(Base::getHelper().getElt(i));}
    ELT&       operator[](int i)       {return *reinterpret_cast<ELT*>      (Base::updHelper().updElt(i));}
    const ELT& operator()(int i) const {return *reinterpret_cast<const ELT*>(Base::getHelper().getElt(i));}
    ELT&       operator()(int i)       {return *reinterpret_cast<ELT*>      (Base::updHelper().updElt(i));}
         
    // Block (contiguous subvector) view creation      
    VectorView_<ELT> operator()(int i, int m) const {return Base::operator()(i,0,m,1).getAsVectorView();}
    VectorView_<ELT> operator()(int i, int m)       {return Base::operator()(i,0,m,1).updAsVectorView();}

    // Indexed view creation (arbitrary subvector). Indices must be monotonically increasing.
    VectorView_<ELT> index(const Array_<int>& indices) const {
        MatrixHelper<Scalar> h(Base::getHelper().getCharacterCommitment(), Base::getHelper(), indices);
        return VectorView_<ELT>(h);
    }
    VectorView_<ELT> updIndex(const Array_<int>& indices) {
        MatrixHelper<Scalar> h(Base::getHelper().getCharacterCommitment(), Base::updHelper(), indices);
        return VectorView_<ELT>(h);
    }

    VectorView_<ELT> operator()(const Array_<int>& indices) const {return index(indices);}
    VectorView_<ELT> operator()(const Array_<int>& indices)       {return updIndex(indices);}
 
    // Hermitian transpose.
    THerm transpose() const {return Base::transpose().getAsRowVectorView();}
    THerm updTranspose()    {return Base::updTranspose().updAsRowVectorView();}

    THerm operator~() const {return transpose();}
    THerm operator~()       {return updTranspose();}

    const VectorBase& operator+() const {return *this; }

    // Negation

    const TNeg& negate()    const {return *reinterpret_cast<const TNeg*>(this); }
    TNeg&       updNegate()       {return *reinterpret_cast<TNeg*>(this); }

    const TNeg& operator-() const {return negate();}
    TNeg&       operator-()       {return updNegate();}

    VectorBase& resize(int m)     {Base::resize(m,1); return *this;}
    VectorBase& resizeKeep(int m) {Base::resizeKeep(m,1); return *this;}

    //TODO: this is not re-locking the number of columns at 1.
    void clear() {Base::clear(); Base::resize(0,1);}

    ELT sum() const {ELT s; Base::getHelper().sum(reinterpret_cast<Scalar*>(&s)); return s; } // add all the elements        
    VectorIterator<ELT, VectorBase<ELT> > begin() {
        return VectorIterator<ELT, VectorBase<ELT> >(*this, 0);
    }
    VectorIterator<ELT, VectorBase<ELT> > end() {
        return VectorIterator<ELT, VectorBase<ELT> >(*this, size());
    }

protected:
    // Create a VectorBase handle using a given helper rep. 
    explicit VectorBase(MatrixHelperRep<Scalar>* hrep) : Base(hrep) {}

private:
    // NO DATA MEMBERS ALLOWED
};



//  ------------------------------- RowVectorBase ------------------------------
//  ----------------------------------------------------------------------------
template <class ELT> class RowVectorBase : public MatrixBase<ELT> {
    typedef MatrixBase<ELT>                             Base;
    typedef typename CNT<ELT>::Scalar                   Scalar;
    typedef typename CNT<ELT>::Number                   Number;
    typedef typename CNT<ELT>::StdNumber                StdNumber;
    typedef RowVectorBase<ELT>                          T;
    typedef RowVectorBase<typename CNT<ELT>::TAbs>      TAbs;
    typedef RowVectorBase<typename CNT<ELT>::TNeg>      TNeg;
    typedef VectorView_<typename CNT<ELT>::THerm>       THerm;
public: 
    //  ------------------------------------------------------------------------

    explicit RowVectorBase(int n=0) : Base(MatrixCommitment::RowVector(), 1, n) {}
    
    RowVectorBase(const RowVectorBase& source) : Base(source) {}

    RowVectorBase(const TNeg& source) : Base(source) {}

    RowVectorBase(int n, const ELT& initialValue)
    :   Base(MatrixCommitment::RowVector(),1,n,initialValue) {}  

    RowVectorBase(int n, const ELT* cppInitialValues)
    :   Base(MatrixCommitment::RowVector(),1,n,cppInitialValues) {}

    //  ------------------------------------------------------------------------

    RowVectorBase(int n, int stride, const Scalar* s)
    :   Base(MatrixCommitment::RowVector(n), MatrixCharacter::RowVector(n),stride,s) { }
    RowVectorBase(int n, int stride, Scalar* s)
    :   Base(MatrixCommitment::RowVector(n), MatrixCharacter::RowVector(n),stride,s) { }

    //  ------------------------------------------------------------------------

    RowVectorBase(MatrixHelper<Scalar>& h, const typename MatrixHelper<Scalar>::ShallowCopy& s) 
    :   Base(MatrixCommitment::RowVector(), h,s) { }
    RowVectorBase(const MatrixHelper<Scalar>& h, const typename MatrixHelper<Scalar>::ShallowCopy& s) 
    :   Base(MatrixCommitment::RowVector(), h,s) { }
    RowVectorBase(const MatrixHelper<Scalar>& h, const typename MatrixHelper<Scalar>::DeepCopy& d)    
    :   Base(MatrixCommitment::RowVector(), h,d) { }

    // This gives the resulting rowvector type when (r(i) op P) is applied to each element.
    // It will have element types which are the regular composite result of ELT op P.
    template <class P> struct EltResult { 
        typedef RowVectorBase<typename CNT<ELT>::template Result<P>::Mul> Mul;
        typedef RowVectorBase<typename CNT<ELT>::template Result<P>::Dvd> Dvd;
        typedef RowVectorBase<typename CNT<ELT>::template Result<P>::Add> Add;
        typedef RowVectorBase<typename CNT<ELT>::template Result<P>::Sub> Sub;
    };

    RowVectorBase& operator=(const RowVectorBase& b) {
        Base::operator=(b); return *this;
    }

    // default destructor

    RowVectorBase& operator*=(const StdNumber& t)     {Base::operator*=(t); return *this;}
    RowVectorBase& operator/=(const StdNumber& t)     {Base::operator/=(t); return *this;}
    RowVectorBase& operator+=(const RowVectorBase& r) {Base::operator+=(r); return *this;}
    RowVectorBase& operator-=(const RowVectorBase& r) {Base::operator-=(r); return *this;}  

    template <class EE> RowVectorBase& operator=(const RowVectorBase<EE>& b) 
      { Base::operator=(b);  return *this; } 
    template <class EE> RowVectorBase& operator+=(const RowVectorBase<EE>& b) 
      { Base::operator+=(b); return *this; } 
    template <class EE> RowVectorBase& operator-=(const RowVectorBase<EE>& b) 
      { Base::operator-=(b); return *this; } 

    // default destructor
 
    RowVectorBase& operator=(const ELT& t) { Base::setTo(t); return *this; } 

    template <class EE> RowVectorBase& colScaleInPlace(const VectorBase<EE>& v)
      { Base::template colScaleInPlace<EE>(v); return *this; }
    template <class EE> inline void colScale(const VectorBase<EE>& v, typename EltResult<EE>::Mul& out) const
      { return Base::template colScale<EE>(v,out); }
    template <class EE> inline typename EltResult<EE>::Mul colScale(const VectorBase<EE>& v) const
      { typename EltResult<EE>::Mul out(ncol()); Base::template colScale<EE>(v,out); return out; }


        // elementwise multiply
    template <class EE> RowVectorBase& elementwiseMultiplyInPlace(const RowVectorBase<EE>& r)
      { Base::template elementwiseMultiplyInPlace<EE>(r); return *this; }
    template <class EE> inline void elementwiseMultiply(const RowVectorBase<EE>& v, typename EltResult<EE>::Mul& out) const
      { Base::template elementwiseMultiply<EE>(v,out); }
    template <class EE> inline typename EltResult<EE>::Mul elementwiseMultiply(const RowVectorBase<EE>& v) const
      { typename EltResult<EE>::Mul out(nrow()); Base::template elementwiseMultiply<EE>(v,out); return out; }

        // elementwise multiply from left
    template <class EE> RowVectorBase& elementwiseMultiplyFromLeftInPlace(const RowVectorBase<EE>& r)
      { Base::template elementwiseMultiplyFromLeftInPlace<EE>(r); return *this; }
    template <class EE> inline void 
    elementwiseMultiplyFromLeft(
        const RowVectorBase<EE>& v, 
        typename RowVectorBase<EE>::template EltResult<ELT>::Mul& out) const
    { 
        Base::template elementwiseMultiplyFromLeft<EE>(v,out);
    }
    template <class EE> inline typename RowVectorBase<EE>::template EltResult<ELT>::Mul 
    elementwiseMultiplyFromLeft(const RowVectorBase<EE>& v) const
    { 
        typename RowVectorBase<EE>::template EltResult<ELT>::Mul out(nrow()); 
        Base::template elementwiseMultiplyFromLeft<EE>(v,out); 
        return out;
    }

        // elementwise divide
    template <class EE> RowVectorBase& elementwiseDivideInPlace(const RowVectorBase<EE>& r)
      { Base::template elementwiseDivideInPlace<EE>(r); return *this; }
    template <class EE> inline void elementwiseDivide(const RowVectorBase<EE>& v, typename EltResult<EE>::Dvd& out) const
      { Base::template elementwiseDivide<EE>(v,out); }
    template <class EE> inline typename EltResult<EE>::Dvd elementwiseDivide(const RowVectorBase<EE>& v) const
      { typename EltResult<EE>::Dvd out(nrow()); Base::template elementwiseDivide<EE>(v,out); return out; }

        // elementwise divide from left
    template <class EE> RowVectorBase& elementwiseDivideFromLeftInPlace(const RowVectorBase<EE>& r)
      { Base::template elementwiseDivideFromLeftInPlace<EE>(r); return *this; }
    template <class EE> inline void 
    elementwiseDivideFromLeft(
        const RowVectorBase<EE>& v, 
        typename RowVectorBase<EE>::template EltResult<ELT>::Dvd& out) const
    { 
        Base::template elementwiseDivideFromLeft<EE>(v,out);
    }
    template <class EE> inline typename RowVectorBase<EE>::template EltResult<ELT>::Dvd 
    elementwiseDivideFromLeft(const RowVectorBase<EE>& v) const
    { 
        typename RowVectorBase<EE>::template EltResult<ELT>::Dvd out(nrow()); 
        Base::template elementwiseDivideFromLeft<EE>(v,out); 
        return out;
    }

    // Implicit conversions are allowed to RowVector or Matrix, but not to Vector.   
    operator const RowVector_<ELT>&()     const {return *reinterpret_cast<const RowVector_<ELT>*>(this);}
    operator       RowVector_<ELT>&()           {return *reinterpret_cast<      RowVector_<ELT>*>(this);}
    operator const RowVectorView_<ELT>&() const {return *reinterpret_cast<const RowVectorView_<ELT>*>(this);}
    operator       RowVectorView_<ELT>&()       {return *reinterpret_cast<      RowVectorView_<ELT>*>(this);}
    
    operator const Matrix_<ELT>&()     const {return *reinterpret_cast<const Matrix_<ELT>*>(this);}
    operator       Matrix_<ELT>&()           {return *reinterpret_cast<      Matrix_<ELT>*>(this);} 
    operator const MatrixView_<ELT>&() const {return *reinterpret_cast<const MatrixView_<ELT>*>(this);}
    operator       MatrixView_<ELT>&()       {return *reinterpret_cast<      MatrixView_<ELT>*>(this);} 
    

    // size() for RowVectors is Base::nelt() but returns int instead of ptrdiff_t.
    int size() const { 
        assert(Base::nelt() <= (ptrdiff_t)std::numeric_limits<int>::max()); 
        assert(Base::nrow()==1);
        return (int)Base::nelt();
    }
    int       nrow() const {assert(Base::nrow()==1); return Base::nrow();}
    int       ncol() const {assert(Base::nrow()==1); return Base::ncol();}
    ptrdiff_t nelt() const {assert(Base::nrow()==1); return Base::nelt();}

    // Override MatrixBase operators to return the right shape
    TAbs abs() const {
        TAbs result; Base::abs(result); return result;
    }

    // Override MatrixBase indexing operators          
    const ELT& operator[](int j) const {return *reinterpret_cast<const ELT*>(Base::getHelper().getElt(j));}
    ELT&       operator[](int j)       {return *reinterpret_cast<ELT*>      (Base::updHelper().updElt(j));}
    const ELT& operator()(int j) const {return *reinterpret_cast<const ELT*>(Base::getHelper().getElt(j));}
    ELT&       operator()(int j)       {return *reinterpret_cast<ELT*>      (Base::updHelper().updElt(j));}
         
    // Block (contiguous subvector) creation      
    RowVectorView_<ELT> operator()(int j, int n) const {return Base::operator()(0,j,1,n).getAsRowVectorView();}
    RowVectorView_<ELT> operator()(int j, int n)       {return Base::operator()(0,j,1,n).updAsRowVectorView();}

    // Indexed view creation (arbitrary subvector). Indices must be monotonically increasing.
    RowVectorView_<ELT> index(const Array_<int>& indices) const {
        MatrixHelper<Scalar> h(Base::getHelper().getCharacterCommitment(), Base::getHelper(), indices);
        return RowVectorView_<ELT>(h);
    }
    RowVectorView_<ELT> updIndex(const Array_<int>& indices) {
        MatrixHelper<Scalar> h(Base::getHelper().getCharacterCommitment(), Base::updHelper(), indices);
        return RowVectorView_<ELT>(h);
    }

    RowVectorView_<ELT> operator()(const Array_<int>& indices) const {return index(indices);}
    RowVectorView_<ELT> operator()(const Array_<int>& indices)       {return updIndex(indices);}
 
    // Hermitian transpose.
    THerm transpose() const {return Base::transpose().getAsVectorView();}
    THerm updTranspose()    {return Base::updTranspose().updAsVectorView();}

    THerm operator~() const {return transpose();}
    THerm operator~()       {return updTranspose();}

    const RowVectorBase& operator+() const {return *this; }

    // Negation

    const TNeg& negate()    const {return *reinterpret_cast<const TNeg*>(this); }
    TNeg&       updNegate()       {return *reinterpret_cast<TNeg*>(this); }

    const TNeg& operator-() const {return negate();}
    TNeg&       operator-()       {return updNegate();}

    RowVectorBase& resize(int n)     {Base::resize(1,n); return *this;}
    RowVectorBase& resizeKeep(int n) {Base::resizeKeep(1,n); return *this;}

    //TODO: this is not re-locking the number of rows at 1.
    void clear() {Base::clear(); Base::resize(1,0);}

    ELT sum() const {ELT s; Base::getHelper().sum(reinterpret_cast<Scalar*>(&s)); return s; } // add all the elements        
    VectorIterator<ELT, RowVectorBase<ELT> > begin() {
        return VectorIterator<ELT, RowVectorBase<ELT> >(*this, 0);
    }
    VectorIterator<ELT, RowVectorBase<ELT> > end() {
        return VectorIterator<ELT, RowVectorBase<ELT> >(*this, size());
    }

protected:
    // Create a RowVectorBase handle using a given helper rep. 
    explicit RowVectorBase(MatrixHelperRep<Scalar>* hrep) : Base(hrep) {}

private:
    // NO DATA MEMBERS ALLOWED
};



//  ------------------------------- MatrixView_ --------------------------------
//  ----------------------------------------------------------------------------
template <class ELT> class MatrixView_ : public MatrixBase<ELT> {
    typedef MatrixBase<ELT>                 Base;
    typedef typename CNT<ELT>::Scalar       S;
    typedef typename CNT<ELT>::StdNumber    StdNumber;
public:
    // Default construction is suppressed.
    // Uses default destructor.

    // Create a MatrixView_ handle using a given helper rep. 
    explicit MatrixView_(MatrixHelperRep<S>* hrep) : Base(hrep) {}

    // Copy constructor is shallow. CAUTION: despite const argument, this preserves writability
    // if it was present in the source. This is necessary to allow temporary views to be
    // created and used as lvalues.
    MatrixView_(const MatrixView_& m) 
      : Base(MatrixCommitment(),
             const_cast<MatrixHelper<S>&>(m.getHelper()), 
             typename MatrixHelper<S>::ShallowCopy()) {}

    // Copy assignment is deep but not reallocating.
    MatrixView_& operator=(const MatrixView_& m) {
        Base::operator=(m); return *this;
    }

    // Copy construction and copy assignment from a DeadMatrixView steals the helper.
    MatrixView_(DeadMatrixView_<ELT>&);
    MatrixView_& operator=(DeadMatrixView_<ELT>&);

    // Ask for shallow copy    
    MatrixView_(const MatrixHelper<S>& h) : Base(MatrixCommitment(), h, typename MatrixHelper<S>::ShallowCopy()) { }
    MatrixView_(MatrixHelper<S>&       h) : Base(MatrixCommitment(), h, typename MatrixHelper<S>::ShallowCopy()) { }

    MatrixView_& operator=(const Matrix_<ELT>& v)     { Base::operator=(v); return *this; }
    MatrixView_& operator=(const ELT& e)              { Base::operator=(e); return *this; }

    template <class EE> MatrixView_& operator=(const MatrixBase<EE>& m)
      { Base::operator=(m); return *this; }
    template <class EE> MatrixView_& operator+=(const MatrixBase<EE>& m)
      { Base::operator+=(m); return *this; }
    template <class EE> MatrixView_& operator-=(const MatrixBase<EE>& m)
      { Base::operator-=(m); return *this; }

    MatrixView_& operator*=(const StdNumber& t) { Base::operator*=(t); return *this; }
    MatrixView_& operator/=(const StdNumber& t) { Base::operator/=(t); return *this; }
    MatrixView_& operator+=(const ELT& r)       { this->updDiag() += r; return *this; }
    MatrixView_& operator-=(const ELT& r)       { this->updDiag() -= r; return *this; }  

    operator const Matrix_<ELT>&() const { return *reinterpret_cast<const Matrix_<ELT>*>(this); }
    operator Matrix_<ELT>&()             { return *reinterpret_cast<Matrix_<ELT>*>(this); }      

private:
    // NO DATA MEMBERS ALLOWED
    MatrixView_(); // default constructor suppressed (what's it a view of?)
};



//  ----------------------------- DeadMatrixView_ ------------------------------
//  ----------------------------------------------------------------------------
template <class ELT> class DeadMatrixView_ : public MatrixView_<ELT> {
    typedef MatrixView_<ELT>                Base;
    typedef typename CNT<ELT>::Scalar       S;
    typedef typename CNT<ELT>::StdNumber    StdNumber;
public:
    // Default construction is suppressed.
    // Uses default destructor.
    
    // All functionality is passed through to MatrixView_.
    explicit DeadMatrixView_(MatrixHelperRep<S>* hrep) : Base(hrep) {}
    DeadMatrixView_(const Base& m) : Base(m) {}
    DeadMatrixView_& operator=(const Base& m) {
        Base::operator=(m); return *this;
    }

    // Ask for shallow copy    
    DeadMatrixView_(const MatrixHelper<S>& h) : Base(h) {}
    DeadMatrixView_(MatrixHelper<S>&       h) : Base(h) {}

    DeadMatrixView_& operator=(const Matrix_<ELT>& v)     { Base::operator=(v); return *this; }
    DeadMatrixView_& operator=(const ELT& e)              { Base::operator=(e); return *this; }

    template <class EE> DeadMatrixView_& operator=(const MatrixBase<EE>& m)
      { Base::operator=(m); return *this; }
    template <class EE> DeadMatrixView_& operator+=(const MatrixBase<EE>& m)
      { Base::operator+=(m); return *this; }
    template <class EE> DeadMatrixView_& operator-=(const MatrixBase<EE>& m)
      { Base::operator-=(m); return *this; }

    DeadMatrixView_& operator*=(const StdNumber& t) { Base::operator*=(t); return *this; }
    DeadMatrixView_& operator/=(const StdNumber& t) { Base::operator/=(t); return *this; }
    DeadMatrixView_& operator+=(const ELT& r)       { this->updDiag() += r; return *this; }
    DeadMatrixView_& operator-=(const ELT& r)       { this->updDiag() -= r; return *this; }  

private:
    // NO DATA MEMBERS ALLOWED
    DeadMatrixView_(); // default constructor suppressed (what's it a view of?)
};

template <class ELT> inline 
MatrixView_<ELT>::MatrixView_(DeadMatrixView_<ELT>& dead) 
:   Base(dead.updHelper().stealRep()) {}

template <class ELT> inline MatrixView_<ELT>& 
MatrixView_<ELT>::operator=(DeadMatrixView_<ELT>& dead) {
    if (Base::getHelper().getCharacterCommitment().isSatisfiedBy(dead.getMatrixCharacter()))
        Base::updHelper().replaceRep(dead.updHelper().stealRep());
    else
        Base::operator=(dead);
    return *this;
}


//  ---------------------------------- Matrix_ ---------------------------------
//  ----------------------------------------------------------------------------
template <class ELT> class Matrix_ : public MatrixBase<ELT> {
    typedef typename CNT<ELT>::Scalar       S;
    typedef typename CNT<ELT>::Number       Number;
    typedef typename CNT<ELT>::StdNumber    StdNumber;

    typedef typename CNT<ELT>::TNeg         ENeg;
    typedef typename CNT<ELT>::THerm        EHerm;

    typedef MatrixBase<ELT>     Base;
    typedef MatrixBase<ENeg>    BaseNeg;
    typedef MatrixBase<EHerm>   BaseHerm;

    typedef Matrix_<ELT>        T;
    typedef MatrixView_<ELT>    TView;
    typedef Matrix_<ENeg>       TNeg;

public:
    Matrix_() : Base() { }
    Matrix_(const MatrixCommitment& mc) : Base(mc) {}

    // Copy constructor is deep.
    Matrix_(const Matrix_& src) : Base(src) { }

    // Assignment is a deep copy and will also allow reallocation if this Matrix
    // doesn't have a view.
    Matrix_& operator=(const Matrix_& src) { 
        Base::operator=(src); return *this;
    }

    // Force a deep copy of the view or whatever this is.
    // Note that this is an implicit conversion.
    Matrix_(const Base& v) : Base(v) {}   // e.g., MatrixView

    // Allow implicit conversion from a source matrix that
    // has a negated version of ELT.
    Matrix_(const BaseNeg& v) : Base(v) {}

    // TODO: implicit conversion from conjugate. This is trickier
    // since real elements are their own conjugate so you'll get
    // duplicate methods defined from Matrix_(BaseHerm) and Matrix_(Base).

    Matrix_(int m, int n) : Base(MatrixCommitment(), m, n) {}

    Matrix_(int m, int n, const ELT* cppInitialValuesByRow) 
    :   Base(MatrixCommitment(), m, n, cppInitialValuesByRow) {}
    Matrix_(int m, int n, const ELT& initialValue) 
    :   Base(MatrixCommitment(), m, n, initialValue) {}
    
    Matrix_(int m, int n, int leadingDim, const S* data) // read only
    :   Base(MatrixCommitment(), MatrixCharacter::LapackFull(m,n), 
             leadingDim, data) {}
    Matrix_(int m, int n, int leadingDim, S* data) // writable
    :   Base(MatrixCommitment(), MatrixCharacter::LapackFull(m,n), 
             leadingDim, data) {}
    
    template <int M, int N, int CS, int RS>
    explicit Matrix_(const Mat<M,N,ELT,CS,RS>& mat)
    :   Base(MatrixCommitment(), M, N)
    {   for (int i = 0; i < M; ++i)
            for (int j = 0; j < N; ++j)
                this->updElt(i, j) = mat(i, j); }

    Matrix_& operator=(const ELT& v) { Base::operator=(v); return *this; }

    template <class EE> Matrix_& operator=(const MatrixBase<EE>& m)
      { Base::operator=(m); return*this; }
    template <class EE> Matrix_& operator+=(const MatrixBase<EE>& m)
      { Base::operator+=(m); return*this; }
    template <class EE> Matrix_& operator-=(const MatrixBase<EE>& m)
      { Base::operator-=(m); return*this; }

    Matrix_& operator*=(const StdNumber& t) { Base::operator*=(t); return *this; }
    Matrix_& operator/=(const StdNumber& t) { Base::operator/=(t); return *this; }
    Matrix_& operator+=(const ELT& r)       { this->updDiag() += r; return *this; }
    Matrix_& operator-=(const ELT& r)       { this->updDiag() -= r; return *this; }  

    const TNeg& negate()    const {return *reinterpret_cast<const TNeg*>(this); }
    TNeg&       updNegate()       {return *reinterpret_cast<TNeg*>(this); }

    const TNeg& operator-() const {return negate();}
    TNeg&       operator-()       {return updNegate();}
   
private:
    // NO DATA MEMBERS ALLOWED
};



//  -------------------------------- VectorView_ -------------------------------
//  ----------------------------------------------------------------------------
template <class ELT> class VectorView_ : public VectorBase<ELT> {
    typedef VectorBase<ELT>                             Base;
    typedef typename CNT<ELT>::Scalar                   S;
    typedef typename CNT<ELT>::Number                   Number;
    typedef typename CNT<ELT>::StdNumber                StdNumber;
    typedef VectorView_<ELT>                            T;
    typedef VectorView_< typename CNT<ELT>::TNeg >      TNeg;
    typedef RowVectorView_< typename CNT<ELT>::THerm >  THerm;
public:
    // Default construction is suppressed.
    // Uses default destructor.

    // Create a VectorView_ handle using a given helper rep. 
    explicit VectorView_(MatrixHelperRep<S>* hrep) : Base(hrep) {}

    // Copy constructor is shallow. CAUTION: despite const argument, this preserves writability
    // if it was present in the source. This is necessary to allow temporary views to be
    // created and used as lvalues.
    VectorView_(const VectorView_& v) 
      : Base(const_cast<MatrixHelper<S>&>(v.getHelper()), typename MatrixHelper<S>::ShallowCopy()) { }

    // Copy assignment is deep but not reallocating.
    VectorView_& operator=(const VectorView_& v) {
        Base::operator=(v); return *this;
    }

    // Ask for shallow copy    
    explicit VectorView_(const MatrixHelper<S>& h) : Base(h, typename MatrixHelper<S>::ShallowCopy()) { }
    explicit VectorView_(MatrixHelper<S>&       h) : Base(h, typename MatrixHelper<S>::ShallowCopy()) { }
    
    VectorView_& operator=(const Base& b) { Base::operator=(b); return *this; }

    VectorView_& operator=(const ELT& v) { Base::operator=(v); return *this; } 

    template <class EE> VectorView_& operator=(const VectorBase<EE>& m)
      { Base::operator=(m); return*this; }
    template <class EE> VectorView_& operator+=(const VectorBase<EE>& m)
      { Base::operator+=(m); return*this; }
    template <class EE> VectorView_& operator-=(const VectorBase<EE>& m)
      { Base::operator-=(m); return*this; }

    VectorView_& operator*=(const StdNumber& t) { Base::operator*=(t); return *this; }
    VectorView_& operator/=(const StdNumber& t) { Base::operator/=(t); return *this; }
    VectorView_& operator+=(const ELT& b) { elementwiseAddScalarInPlace(b); return *this; } 
    VectorView_& operator-=(const ELT& b) { elementwiseSubtractScalarInPlace(b); return *this; } 

private:
    // NO DATA MEMBERS ALLOWED
    VectorView_(); // default construction suppressed (what's it a View of?)
};



//  ---------------------------------- Vector_ ---------------------------------
//  ----------------------------------------------------------------------------
template <class ELT> class Vector_ : public VectorBase<ELT> {
    typedef typename CNT<ELT>::Scalar       S;
    typedef typename CNT<ELT>::Number       Number;
    typedef typename CNT<ELT>::StdNumber    StdNumber;
    typedef typename CNT<ELT>::TNeg         ENeg;
    typedef VectorBase<ELT>                 Base;
    typedef VectorBase<ENeg>                BaseNeg;
public:
    Vector_() : Base() {}  // 0x1 reallocatable
    // Uses default destructor.

    // Copy constructor is deep.
    Vector_(const Vector_& src) : Base(src) {}

    // Implicit conversions.
    Vector_(const Base& src) : Base(src) {}    // e.g., VectorView
    Vector_(const BaseNeg& src) : Base(src) {}

    // Copy assignment is deep and can be reallocating if this Vector
    // has no View.
    Vector_& operator=(const Vector_& src) {
        Base::operator=(src); return*this;
    }


    explicit Vector_(int m) : Base(m) { }
    Vector_(int m, const ELT* cppInitialValues) : Base(m, cppInitialValues) {}
    Vector_(int m, const ELT& initialValue)     : Base(m, initialValue) {}

    Vector_(int m, const S* cppData, bool): Base(m, Base::CppNScalarsPerElement, cppData) {}
    Vector_(int m,       S* cppData, bool): Base(m, Base::CppNScalarsPerElement, cppData) {}

    Vector_(int m, int stride, const S* data, bool) : Base(m, stride, data) {}
    Vector_(int m, int stride,       S* data, bool) : Base(m, stride, data) {}

    template <int M>
    explicit Vector_(const Vec<M,ELT>& v) : Base(M) {
        for (int i = 0; i < M; ++i)
            this->updElt(i, 0) = v(i);
    }

    Vector_& operator=(const ELT& v) { Base::operator=(v); return *this; } 

    template <class EE> Vector_& operator=(const VectorBase<EE>& m)
      { Base::operator=(m); return*this; }
    template <class EE> Vector_& operator+=(const VectorBase<EE>& m)
      { Base::operator+=(m); return*this; }
    template <class EE> Vector_& operator-=(const VectorBase<EE>& m)
      { Base::operator-=(m); return*this; }

    Vector_& operator*=(const StdNumber& t) { Base::operator*=(t); return *this; }
    Vector_& operator/=(const StdNumber& t) { Base::operator/=(t); return *this; }
    Vector_& operator+=(const ELT& b) { elementwiseAddScalarInPlace(b); return *this; } 
    Vector_& operator-=(const ELT& b) { elementwiseSubtractScalarInPlace(b); return *this; } 
 
private:
    // NO DATA MEMBERS ALLOWED
};



//  ------------------------------ RowVectorView_ ------------------------------
//  ----------------------------------------------------------------------------
template <class ELT> class RowVectorView_ : public RowVectorBase<ELT> {
    typedef RowVectorBase<ELT>                              Base;
    typedef typename CNT<ELT>::Scalar                       S;
    typedef typename CNT<ELT>::Number                       Number;
    typedef typename CNT<ELT>::StdNumber                    StdNumber;
    typedef RowVectorView_<ELT>                             T;
    typedef RowVectorView_< typename CNT<ELT>::TNeg >       TNeg;
    typedef VectorView_< typename CNT<ELT>::THerm >         THerm;
public:
    // Default construction is suppressed.
    // Uses default destructor.

    // Create a RowVectorView_ handle using a given helper rep. 
    explicit RowVectorView_(MatrixHelperRep<S>* hrep) : Base(hrep) {}

    // Copy constructor is shallow. CAUTION: despite const argument, this preserves writability
    // if it was present in the source. This is necessary to allow temporary views to be
    // created and used as lvalues.
    RowVectorView_(const RowVectorView_& r) 
      : Base(const_cast<MatrixHelper<S>&>(r.getHelper()), typename MatrixHelper<S>::ShallowCopy()) { }

    // Copy assignment is deep but not reallocating.
    RowVectorView_& operator=(const RowVectorView_& r) {
        Base::operator=(r); return *this;
    }

    // Ask for shallow copy    
    explicit RowVectorView_(const MatrixHelper<S>& h) : Base(h, typename MatrixHelper<S>::ShallowCopy()) { }
    explicit RowVectorView_(MatrixHelper<S>&       h) : Base(h, typename MatrixHelper<S>::ShallowCopy()) { }
    
    RowVectorView_& operator=(const Base& b) { Base::operator=(b); return *this; }

    RowVectorView_& operator=(const ELT& v) { Base::operator=(v); return *this; } 

    template <class EE> RowVectorView_& operator=(const RowVectorBase<EE>& m)
      { Base::operator=(m); return*this; }
    template <class EE> RowVectorView_& operator+=(const RowVectorBase<EE>& m)
      { Base::operator+=(m); return*this; }
    template <class EE> RowVectorView_& operator-=(const RowVectorBase<EE>& m)
      { Base::operator-=(m); return*this; }

    RowVectorView_& operator*=(const StdNumber& t) { Base::operator*=(t); return *this; }
    RowVectorView_& operator/=(const StdNumber& t) { Base::operator/=(t); return *this; }
    RowVectorView_& operator+=(const ELT& b) { elementwiseAddScalarInPlace(b); return *this; } 
    RowVectorView_& operator-=(const ELT& b) { elementwiseSubtractScalarInPlace(b); return *this; } 

private:
    // NO DATA MEMBERS ALLOWED
    RowVectorView_(); // default construction suppressed (what is it a view of?)
};



//  -------------------------------- RowVector_ --------------------------------
//  ----------------------------------------------------------------------------
template <class ELT> class RowVector_ : public RowVectorBase<ELT> {
    typedef typename CNT<ELT>::Scalar       S;
    typedef typename CNT<ELT>::Number       Number;
    typedef typename CNT<ELT>::StdNumber    StdNumber;
    typedef typename CNT<ELT>::TNeg         ENeg;

    typedef RowVectorBase<ELT>              Base;
    typedef RowVectorBase<ENeg>             BaseNeg;
public:
    RowVector_() : Base() {}   // 1x0 reallocatable
    // Uses default destructor.

    // Copy constructor is deep.
    RowVector_(const RowVector_& src) : Base(src) {}

    // Implicit conversions.
    RowVector_(const Base& src) : Base(src) {}    // e.g., RowVectorView
    RowVector_(const BaseNeg& src) : Base(src) {}  

    // Copy assignment is deep and can be reallocating if this RowVector
    // has no View.
    RowVector_& operator=(const RowVector_& src) {
        Base::operator=(src); return*this;
    }


    explicit RowVector_(int n) : Base(n) { }
    RowVector_(int n, const ELT* cppInitialValues) : Base(n, cppInitialValues) {}
    RowVector_(int n, const ELT& initialValue)     : Base(n, initialValue) {}

    RowVector_(int n, const S* cppData, bool): Base(n, Base::CppNScalarsPerElement, cppData) {}
    RowVector_(int n,       S* cppData, bool): Base(n, Base::CppNScalarsPerElement, cppData) {}

    RowVector_(int n, int stride, const S* data, bool) : Base(n, stride, data) {}
    RowVector_(int n, int stride,       S* data, bool) : Base(n, stride, data) {}
    
    template <int M>
    explicit RowVector_(const Row<M,ELT>& v) : Base(M) {
        for (int i = 0; i < M; ++i)
            this->updElt(0, i) = v(i);
    }

    RowVector_& operator=(const ELT& v) { Base::operator=(v); return *this; } 

    template <class EE> RowVector_& operator=(const RowVectorBase<EE>& b)
      { Base::operator=(b); return*this; }
    template <class EE> RowVector_& operator+=(const RowVectorBase<EE>& b)
      { Base::operator+=(b); return*this; }
    template <class EE> RowVector_& operator-=(const RowVectorBase<EE>& b)
      { Base::operator-=(b); return*this; }

    RowVector_& operator*=(const StdNumber& t) { Base::operator*=(t); return *this; }
    RowVector_& operator/=(const StdNumber& t) { Base::operator/=(t); return *this; }
    RowVector_& operator+=(const ELT& b) { elementwiseAddScalarInPlace(b); return *this; } 
    RowVector_& operator-=(const ELT& b) { elementwiseSubtractScalarInPlace(b); return *this; } 

private:
    // NO DATA MEMBERS ALLOWED
};



//  ------------------------ MatrixBase definitions ----------------------------

template <class ELT> inline MatrixView_<ELT> 
MatrixBase<ELT>::block(int i, int j, int m, int n) const { 
    SimTK_INDEXCHECK(i,nrow()+1,"MatrixBase::block()");
    SimTK_INDEXCHECK(j,ncol()+1,"MatrixBase::block()");
    SimTK_SIZECHECK(i+m,nrow(),"MatrixBase::block()");
    SimTK_SIZECHECK(j+n,ncol(),"MatrixBase::block()");

    MatrixHelper<Scalar> h(MatrixCommitment(),helper,i,j,m,n);    
    return MatrixView_<ELT>(h.stealRep()); 
}
    
template <class ELT> inline MatrixView_<ELT>
MatrixBase<ELT>::updBlock(int i, int j, int m, int n) { 
    SimTK_INDEXCHECK(i,nrow()+1,"MatrixBase::updBlock()");
    SimTK_INDEXCHECK(j,ncol()+1,"MatrixBase::updBlock()");
    SimTK_SIZECHECK(i+m,nrow(),"MatrixBase::updBlock()");
    SimTK_SIZECHECK(j+n,ncol(),"MatrixBase::updBlock()");

    MatrixHelper<Scalar> h(MatrixCommitment(),helper,i,j,m,n);        
    return MatrixView_<ELT>(h.stealRep()); 
}

template <class E> inline MatrixView_<typename CNT<E>::THerm>
MatrixBase<E>::transpose() const { 
    MatrixHelper<typename CNT<Scalar>::THerm> 
        h(MatrixCommitment(),
          helper, typename MatrixHelper<typename CNT<Scalar>::THerm>::TransposeView());
    return MatrixView_<typename CNT<E>::THerm>(h.stealRep()); 
}
    
template <class E> inline MatrixView_<typename CNT<E>::THerm>
MatrixBase<E>::updTranspose() {     
    MatrixHelper<typename CNT<Scalar>::THerm> 
        h(MatrixCommitment(),
          helper, typename MatrixHelper<typename CNT<Scalar>::THerm>::TransposeView());
    return MatrixView_<typename CNT<E>::THerm>(h.stealRep()); 
}

template <class E> inline VectorView_<E>
MatrixBase<E>::diag() const { 
    MatrixHelper<Scalar> h(MatrixCommitment::Vector(),
                           helper, typename MatrixHelper<Scalar>::DiagonalView());
    return VectorView_<E>(h.stealRep()); 
}
    
template <class E> inline VectorView_<E>
MatrixBase<E>::updDiag() {     
    MatrixHelper<Scalar> h(MatrixCommitment::Vector(),
                           helper, typename MatrixHelper<Scalar>::DiagonalView());
    return VectorView_<E>(h.stealRep());
}

template <class ELT> inline VectorView_<ELT> 
MatrixBase<ELT>::col(int j) const { 
    SimTK_INDEXCHECK(j,ncol(),"MatrixBase::col()");

    MatrixHelper<Scalar> h(MatrixCommitment::Vector(),
                           helper,0,j,nrow(),1);    
    return VectorView_<ELT>(h.stealRep()); 
}
    
template <class ELT> inline VectorView_<ELT>
MatrixBase<ELT>::updCol(int j) {
    SimTK_INDEXCHECK(j,ncol(),"MatrixBase::updCol()");

    MatrixHelper<Scalar> h(MatrixCommitment::Vector(),
                           helper,0,j,nrow(),1);        
    return VectorView_<ELT>(h.stealRep()); 
}

template <class ELT> inline RowVectorView_<ELT> 
MatrixBase<ELT>::row(int i) const { 
    SimTK_INDEXCHECK(i,nrow(),"MatrixBase::row()");

    MatrixHelper<Scalar> h(MatrixCommitment::RowVector(),
                           helper,i,0,1,ncol());    
    return RowVectorView_<ELT>(h.stealRep()); 
}
    
template <class ELT> inline RowVectorView_<ELT>
MatrixBase<ELT>::updRow(int i) { 
    SimTK_INDEXCHECK(i,nrow(),"MatrixBase::updRow()");

    MatrixHelper<Scalar> h(MatrixCommitment::RowVector(),
                           helper,i,0,1,ncol());        
    return RowVectorView_<ELT>(h.stealRep()); 
}

// M = diag(v) * M; v must have nrow() elements.
// That is, M[i] *= v[i].
template <class ELT> template <class EE> inline MatrixBase<ELT>& 
MatrixBase<ELT>::rowScaleInPlace(const VectorBase<EE>& v) {
    assert(v.nrow() == nrow());
    for (int i=0; i < nrow(); ++i)
        (*this)[i] *= v[i];
    return *this;
}

template <class ELT> template <class EE> inline void
MatrixBase<ELT>::rowScale(const VectorBase<EE>& v, typename MatrixBase<ELT>::template EltResult<EE>::Mul& out) const {
    assert(v.nrow() == nrow());
    out.resize(nrow(), ncol());
    for (int j=0; j<ncol(); ++j)
        for (int i=0; i<nrow(); ++i)
           out(i,j) = (*this)(i,j) * v[i];
}

// M = M * diag(v); v must have ncol() elements
// That is, M(i) *= v[i]
template <class ELT> template <class EE>  inline MatrixBase<ELT>& 
MatrixBase<ELT>::colScaleInPlace(const VectorBase<EE>& v) {
    assert(v.nrow() == ncol());
    for (int j=0; j < ncol(); ++j)
        (*this)(j) *= v[j];
    return *this;
}

template <class ELT> template <class EE> inline void
MatrixBase<ELT>::colScale(const VectorBase<EE>& v, typename MatrixBase<ELT>::template EltResult<EE>::Mul& out) const {
    assert(v.nrow() == ncol());
    out.resize(nrow(), ncol());
    for (int j=0; j<ncol(); ++j)
        for (int i=0; i<nrow(); ++i)
           out(i,j) = (*this)(i,j) * v[j];
}


// M(i,j) *= r[i]*c[j]; r must have nrow() elements; c must have ncol() elements
template <class ELT> template <class ER, class EC> inline MatrixBase<ELT>& 
MatrixBase<ELT>::rowAndColScaleInPlace(const VectorBase<ER>& r, const VectorBase<EC>& c) {
    assert(r.nrow()==nrow() && c.nrow()==ncol());
    for (int j=0; j<ncol(); ++j)
        for (int i=0; i<nrow(); ++i)
            (*this)(i,j) *= (r[i]*c[j]);
    return *this;
}

template <class ELT> template <class ER, class EC> inline void
MatrixBase<ELT>::rowAndColScale(
    const VectorBase<ER>& r, 
    const VectorBase<EC>& c,
    typename EltResult<typename VectorBase<ER>::template EltResult<EC>::Mul>::Mul& 
                          out) const
{
    assert(r.nrow()==nrow() && c.nrow()==ncol());
    out.resize(nrow(), ncol());
    for (int j=0; j<ncol(); ++j)
        for (int i=0; i<nrow(); ++i)
            out(i,j) = (*this)(i,j) * (r[i]*c[j]);
}


// Set M(i,j) = M(i,j)^-1.
template <class ELT> inline MatrixBase<ELT>& 
MatrixBase<ELT>::elementwiseInvertInPlace() {
    const int nr=nrow(), nc=ncol();
    for (int j=0; j<nc; ++j)
        for (int i=0; i<nr; ++i) {
            ELT& e = updElt(i,j);
            e = CNT<ELT>::invert(e);
        }
    return *this;
}

template <class ELT> inline void
MatrixBase<ELT>::elementwiseInvert(MatrixBase<typename CNT<E>::TInvert>& out) const {
    const int nr=nrow(), nc=ncol();
    out.resize(nr,nc);
    for (int j=0; j<nc; ++j)
        for (int i=0; i<nr; ++i)
            out(i,j) = CNT<ELT>::invert((*this)(i,j));
}

// M(i,j) += s
template <class ELT> template <class S> inline MatrixBase<ELT>& 
MatrixBase<ELT>::elementwiseAddScalarInPlace(const S& s) {
    for (int j=0; j<ncol(); ++j)
        for (int i=0; i<nrow(); ++i)
            (*this)(i,j) += s;
    return *this;
}

template <class ELT> template <class S> inline void 
MatrixBase<ELT>::elementwiseAddScalar(
    const S& s,
    typename MatrixBase<ELT>::template EltResult<S>::Add& out) const
{
    const int nr=nrow(), nc=ncol();
    out.resize(nr,nc);
    for (int j=0; j<nc; ++j)
        for (int i=0; i<nr; ++i)
            out(i,j) = (*this)(i,j) + s;
}

// M(i,j) -= s
template <class ELT> template <class S> inline MatrixBase<ELT>& 
MatrixBase<ELT>::elementwiseSubtractScalarInPlace(const S& s) {
    for (int j=0; j<ncol(); ++j)
        for (int i=0; i<nrow(); ++i)
            (*this)(i,j) -= s;
    return *this;
}

template <class ELT> template <class S> inline void 
MatrixBase<ELT>::elementwiseSubtractScalar(
    const S& s,
    typename MatrixBase<ELT>::template EltResult<S>::Sub& out) const
{
    const int nr=nrow(), nc=ncol();
    out.resize(nr,nc);
    for (int j=0; j<nc; ++j)
        for (int i=0; i<nr; ++i)
            out(i,j) = (*this)(i,j) - s;
}

// M(i,j) = s - M(i,j)
template <class ELT> template <class S> inline MatrixBase<ELT>& 
MatrixBase<ELT>::elementwiseSubtractFromScalarInPlace(const S& s) {
    const int nr=nrow(), nc=ncol();
    for (int j=0; j<nc; ++j)
        for (int i=0; i<nr; ++i) {
            ELT& e = updElt(i,j);
            e = s - e;
        }
    return *this;
}

template <class ELT> template <class S> inline void 
MatrixBase<ELT>::elementwiseSubtractFromScalar(
    const S& s,
    typename MatrixBase<S>::template EltResult<ELT>::Sub& out) const
{
    const int nr=nrow(), nc=ncol();
    out.resize(nr,nc);
    for (int j=0; j<nc; ++j)
        for (int i=0; i<nr; ++i)
            out(i,j) = s - (*this)(i,j);
}

// M(i,j) *= R(i,j); R must have same dimensions as this
template <class ELT> template <class EE> inline MatrixBase<ELT>& 
MatrixBase<ELT>::elementwiseMultiplyInPlace(const MatrixBase<EE>& r) {
    const int nr=nrow(), nc=ncol();
    assert(r.nrow()==nr && r.ncol()==nc);
    for (int j=0; j<nc; ++j)
        for (int i=0; i<nr; ++i)
            (*this)(i,j) *= r(i,j);
    return *this;
}

template <class ELT> template <class EE> inline void 
MatrixBase<ELT>::elementwiseMultiply(
    const MatrixBase<EE>& r, 
    typename MatrixBase<ELT>::template EltResult<EE>::Mul& out) const
{
    const int nr=nrow(), nc=ncol();
    assert(r.nrow()==nr && r.ncol()==nc);
    out.resize(nr,nc);
    for (int j=0; j<nc; ++j)
        for (int i=0; i<nr; ++i)
            out(i,j) = (*this)(i,j) * r(i,j);
}

// M(i,j) = R(i,j) * M(i,j); R must have same dimensions as this
template <class ELT> template <class EE> inline MatrixBase<ELT>& 
MatrixBase<ELT>::elementwiseMultiplyFromLeftInPlace(const MatrixBase<EE>& r) {
    const int nr=nrow(), nc=ncol();
    assert(r.nrow()==nr && r.ncol()==nc);
    for (int j=0; j<nc; ++j)
        for (int i=0; i<nr; ++i) {
            ELT& e = updElt(i,j);
            e = r(i,j) * e;
        }
    return *this;
}

template <class ELT> template <class EE> inline void 
MatrixBase<ELT>::elementwiseMultiplyFromLeft(
    const MatrixBase<EE>& r, 
    typename MatrixBase<EE>::template EltResult<ELT>::Mul& out) const
{
    const int nr=nrow(), nc=ncol();
    assert(r.nrow()==nr && r.ncol()==nc);
    out.resize(nr,nc);
    for (int j=0; j<nc; ++j)
        for (int i=0; i<nr; ++i)
            out(i,j) =  r(i,j) * (*this)(i,j);
}

// M(i,j) /= R(i,j); R must have same dimensions as this
template <class ELT> template <class EE> inline MatrixBase<ELT>& 
MatrixBase<ELT>::elementwiseDivideInPlace(const MatrixBase<EE>& r) {
    const int nr=nrow(), nc=ncol();
    assert(r.nrow()==nr && r.ncol()==nc);
    for (int j=0; j<nc; ++j)
        for (int i=0; i<nr; ++i)
            (*this)(i,j) /= r(i,j);
    return *this;
}

template <class ELT> template <class EE> inline void 
MatrixBase<ELT>::elementwiseDivide(
    const MatrixBase<EE>& r,
    typename MatrixBase<ELT>::template EltResult<EE>::Dvd& out) const
{
    const int nr=nrow(), nc=ncol();
    assert(r.nrow()==nr && r.ncol()==nc);
    out.resize(nr,nc);
    for (int j=0; j<nc; ++j)
        for (int i=0; i<nr; ++i)
            out(i,j) = (*this)(i,j) / r(i,j);
}
// M(i,j) = R(i,j) / M(i,j); R must have same dimensions as this
template <class ELT> template <class EE> inline MatrixBase<ELT>& 
MatrixBase<ELT>::elementwiseDivideFromLeftInPlace(const MatrixBase<EE>& r) {
    const int nr=nrow(), nc=ncol();
    assert(r.nrow()==nr && r.ncol()==nc);
    for (int j=0; j<nc; ++j)
        for (int i=0; i<nr; ++i) {
            ELT& e = updElt(i,j);
            e = r(i,j) / e;
        }
    return *this;
}

template <class ELT> template <class EE> inline void 
MatrixBase<ELT>::elementwiseDivideFromLeft(
    const MatrixBase<EE>& r,
    typename MatrixBase<EE>::template EltResult<ELT>::Dvd& out) const
{
    const int nr=nrow(), nc=ncol();
    assert(r.nrow()==nr && r.ncol()==nc);
    out.resize(nr,nc);
    for (int j=0; j<nc; ++j)
        for (int i=0; i<nr; ++i)
            out(i,j) = r(i,j) / (*this)(i,j);
}

/*
template <class ELT> inline MatrixView_< typename CNT<ELT>::TReal > 
MatrixBase<ELT>::real() const { 
    if (!CNT<ELT>::IsComplex) { // known at compile time
        return MatrixView_< typename CNT<ELT>::TReal >( // this is just ELT
            MatrixHelper(helper,0,0,nrow(),ncol()));    // a view of the whole matrix
    }
    // Elements are complex -- helper uses underlying precision (real) type.
    MatrixHelper<Precision> h(helper,typename MatrixHelper<Precision>::RealView);    
    return MatrixView_< typename CNT<ELT>::TReal >(h); 
}
*/


//  ----------------------------------------------------------------------------

// + and - allow mixed element types, but will fail to compile if the elements aren't
// compatible. At run time these will fail if the dimensions are incompatible.
template <class E1, class E2>
Matrix_<typename CNT<E1>::template Result<E2>::Add>
operator+(const MatrixBase<E1>& l, const MatrixBase<E2>& r) {
    return Matrix_<typename CNT<E1>::template Result<E2>::Add>(l) += r;
}

template <class E>
Matrix_<E> operator+(const MatrixBase<E>& l, const typename CNT<E>::T& r) {
    return Matrix_<E>(l) += r;
}

template <class E>
Matrix_<E> operator+(const typename CNT<E>::T& l, const MatrixBase<E>& r) {
    return Matrix_<E>(r) += l;
}

template <class E1, class E2>
Matrix_<typename CNT<E1>::template Result<E2>::Sub>
operator-(const MatrixBase<E1>& l, const MatrixBase<E2>& r) {
    return Matrix_<typename CNT<E1>::template Result<E2>::Sub>(l) -= r;
}

template <class E>
Matrix_<E> operator-(const MatrixBase<E>& l, const typename CNT<E>::T& r) {
    return Matrix_<E>(l) -= r;
}

template <class E>
Matrix_<E> operator-(const typename CNT<E>::T& l, const MatrixBase<E>& r) {
    Matrix_<E> temp(r.nrow(), r.ncol());
    temp = l;
    return (temp -= r);
}

// Scalar multiply and divide. You might wish the scalar could be
// a templatized type "E2", but that would create horrible ambiguities since
// E2 would match not only scalar types but everything else including
// matrices.
template <class E> Matrix_<E>
operator*(const MatrixBase<E>& l, const typename CNT<E>::StdNumber& r) 
  { return Matrix_<E>(l)*=r; }

template <class E> Matrix_<E>
operator*(const typename CNT<E>::StdNumber& l, const MatrixBase<E>& r) 
  { return Matrix_<E>(r)*=l; }

template <class E> Matrix_<E>
operator/(const MatrixBase<E>& l, const typename CNT<E>::StdNumber& r) 
  { return Matrix_<E>(l)/=r; }

// Handle ints explicitly.
template <class E> Matrix_<E>
operator*(const MatrixBase<E>& l, int r) 
  { return Matrix_<E>(l)*= typename CNT<E>::StdNumber(r); }

template <class E> Matrix_<E>
operator*(int l, const MatrixBase<E>& r) 
  { return Matrix_<E>(r)*= typename CNT<E>::StdNumber(l); }

template <class E> Matrix_<E>
operator/(const MatrixBase<E>& l, int r) 
  { return Matrix_<E>(l)/= typename CNT<E>::StdNumber(r); }



template <class E1, class E2>
Vector_<typename CNT<E1>::template Result<E2>::Add>
operator+(const VectorBase<E1>& l, const VectorBase<E2>& r) {
    return Vector_<typename CNT<E1>::template Result<E2>::Add>(l) += r;
}
template <class E>
Vector_<E> operator+(const VectorBase<E>& l, const typename CNT<E>::T& r) {
    return Vector_<E>(l) += r;
}
template <class E>
Vector_<E> operator+(const typename CNT<E>::T& l, const VectorBase<E>& r) {
    return Vector_<E>(r) += l;
}
template <class E1, class E2>
Vector_<typename CNT<E1>::template Result<E2>::Sub>
operator-(const VectorBase<E1>& l, const VectorBase<E2>& r) {
    return Vector_<typename CNT<E1>::template Result<E2>::Sub>(l) -= r;
}
template <class E>
Vector_<E> operator-(const VectorBase<E>& l, const typename CNT<E>::T& r) {
    return Vector_<E>(l) -= r;
}
template <class E>
Vector_<E> operator-(const typename CNT<E>::T& l, const VectorBase<E>& r) {
    Vector_<E> temp(r.size());
    temp = l;
    return (temp -= r);
}

template <class E> Vector_<E>
operator*(const VectorBase<E>& l, const typename CNT<E>::StdNumber& r) 
  { return Vector_<E>(l)*=r; }

template <class E> Vector_<E>
operator*(const typename CNT<E>::StdNumber& l, const VectorBase<E>& r) 
  { return Vector_<E>(r)*=l; }

template <class E> Vector_<E>
operator/(const VectorBase<E>& l, const typename CNT<E>::StdNumber& r) 
  { return Vector_<E>(l)/=r; }

// Handle ints explicitly
template <class E> Vector_<E>
operator*(const VectorBase<E>& l, int r) 
  { return Vector_<E>(l)*= typename CNT<E>::StdNumber(r); }

template <class E> Vector_<E>
operator*(int l, const VectorBase<E>& r) 
  { return Vector_<E>(r)*= typename CNT<E>::StdNumber(l); }

template <class E> Vector_<E>
operator/(const VectorBase<E>& l, int r) 
  { return Vector_<E>(l)/= typename CNT<E>::StdNumber(r); }



template <class E1, class E2>
RowVector_<typename CNT<E1>::template Result<E2>::Add>
operator+(const RowVectorBase<E1>& l, const RowVectorBase<E2>& r) {
    return RowVector_<typename CNT<E1>::template Result<E2>::Add>(l) += r;
}
template <class E>
RowVector_<E> operator+(const RowVectorBase<E>& l, const typename CNT<E>::T& r) {
    return RowVector_<E>(l) += r;
}
template <class E>
RowVector_<E> operator+(const typename CNT<E>::T& l, const RowVectorBase<E>& r) {
    return RowVector_<E>(r) += l;
}
template <class E1, class E2>
RowVector_<typename CNT<E1>::template Result<E2>::Sub>
operator-(const RowVectorBase<E1>& l, const RowVectorBase<E2>& r) {
    return RowVector_<typename CNT<E1>::template Result<E2>::Sub>(l) -= r;
}
template <class E>
RowVector_<E> operator-(const RowVectorBase<E>& l, const typename CNT<E>::T& r) {
    return RowVector_<E>(l) -= r;
}
template <class E>
RowVector_<E> operator-(const typename CNT<E>::T& l, const RowVectorBase<E>& r) {
    RowVector_<E> temp(r.size());
    temp = l;
    return (temp -= r);
}

template <class E> RowVector_<E>
operator*(const RowVectorBase<E>& l, const typename CNT<E>::StdNumber& r) 
  { return RowVector_<E>(l)*=r; }

template <class E> RowVector_<E>
operator*(const typename CNT<E>::StdNumber& l, const RowVectorBase<E>& r) 
  { return RowVector_<E>(r)*=l; }

template <class E> RowVector_<E>
operator/(const RowVectorBase<E>& l, const typename CNT<E>::StdNumber& r) 
  { return RowVector_<E>(l)/=r; }

// Handle ints explicitly.
template <class E> RowVector_<E>
operator*(const RowVectorBase<E>& l, int r) 
  { return RowVector_<E>(l)*= typename CNT<E>::StdNumber(r); }

template <class E> RowVector_<E>
operator*(int l, const RowVectorBase<E>& r) 
  { return RowVector_<E>(r)*= typename CNT<E>::StdNumber(l); }

template <class E> RowVector_<E>
operator/(const RowVectorBase<E>& l, int r) 
  { return RowVector_<E>(l)/= typename CNT<E>::StdNumber(r); }




    // TODO: these should use LAPACK!

// Dot product
template <class E1, class E2> 
typename CNT<E1>::template Result<E2>::Mul
operator*(const RowVectorBase<E1>& r, const VectorBase<E2>& v) {
    assert(r.ncol() == v.nrow());
    typename CNT<E1>::template Result<E2>::Mul sum(0);
    for (int j=0; j < r.ncol(); ++j)
        sum += r(j) * v[j];
    return sum;
}

template <class E1, class E2> 
Vector_<typename CNT<E1>::template Result<E2>::Mul>
operator*(const MatrixBase<E1>& m, const VectorBase<E2>& v) {
    assert(m.ncol() == v.nrow());
    Vector_<typename CNT<E1>::template Result<E2>::Mul> res(m.nrow());
    for (int i=0; i< m.nrow(); ++i)
        res[i] = m[i]*v;
    return res;
}

template <class E1, class E2> 
Matrix_<typename CNT<E1>::template Result<E2>::Mul>
operator*(const MatrixBase<E1>& m1, const MatrixBase<E2>& m2) {
    assert(m1.ncol() == m2.nrow());
    Matrix_<typename CNT<E1>::template Result<E2>::Mul> 
        res(m1.nrow(),m2.ncol());

    for (int j=0; j < res.ncol(); ++j)
        for (int i=0; i < res.nrow(); ++i)
            res(i,j) = m1[i] * m2(j);

    return res;
}


// This "private" static method is used to implement VectorView's 
// fillVectorViewFromStream() and Vector's readVectorFromStream() 
// namespace-scope static methods, which are in turn used to implement 
// VectorView's and 
// Vector's stream extraction operators ">>". This method has to be in the 
// header file so that we don't need to pass streams through the API, but it 
// is not intended for use by users and has no Doxygen presence, unlike 
// fillArrayFromStream() and readArrayFromStream() and (more commonly)
// the extraction operators.
template <class T> static inline 
std::istream& readVectorFromStreamHelper
   (std::istream& in, bool isFixedSize, Vector_<T>& out)
{
    // If already failed, bad, or eof, set failed bit and return without 
    // touching the Vector.
    if (!in.good()) {in.setstate(std::ios::failbit); return in;}

    // If the passed-in Vector isn't resizeable, then we have to treat it as
    // a fixed size VectorView regardless of the setting of the isFixedSize
    // argument.
    if (!out.isResizeable())
        isFixedSize = true; // might be overriding the argument here

    // numRequired will be ignored unless isFixedSize==true.
    const int numRequired = isFixedSize ? out.size() : 0;

    if (!isFixedSize)
        out.clear(); // We're going to replace the entire contents of the Array.

    // Skip initial whitespace. If that results in eof this may be a successful
    // read of a 0-length, unbracketed Vector. That is OK for either a
    // variable-length Vector or a fixed-length VectorView of length zero.
    std::ws(in); if (in.fail()) return in;
    if (in.eof()) {
        if (isFixedSize && numRequired != 0)
            in.setstate(std::ios_base::failbit); // zero elements not OK
        return in;
    }
    
    // Here the stream is good and the next character is non-white.
    assert(in.good());

    // Use this for raw i/o (peeks and gets).
    typename       std::iostream::int_type ch;
    const typename std::iostream::int_type EOFch = 
        std::iostream::traits_type::eof();

    // First we'll look for the optional "~". If found, the brackets become
    // required.
    bool tildeFound = false;
    ch = in.peek(); if (in.fail()) return in;
    assert(ch != EOFch); // we already checked above
    if ((char)ch == '~') {
        tildeFound = true;
        in.get(); // absorb the tilde
        // Eat whitespace after the tilde to see what's next.
        if (in.good()) std::ws(in);
        // If we hit eof after the tilde we don't like the formatting.
        if (!in.good()) {in.setstate(std::ios_base::failbit); return in;}
    }

    // Here the stream is good, the next character is non-white, and we
    // might have seen a tilde.
    assert(in.good());

    // Now see if the sequence is bare or surrounded by (), or [].
    bool lookForCloser = true;
    char openBracket, closeBracket;
    ch = in.peek(); if (in.fail()) return in;
    assert(ch != EOFch); // we already checked above

    openBracket = (char)ch;
    if      (openBracket=='(') {in.get(); closeBracket = ')';}
    else if (openBracket=='[') {in.get(); closeBracket = ']';}
    else lookForCloser = false;

    // If we found a tilde, the opening bracket was mandatory. If we didn't
    // find one then we reject the formatting.
    if (tildeFound && !lookForCloser)
    {   in.setstate(std::ios_base::failbit); return in;}

    // If lookForCloser is true, then closeBracket contains the terminating
    // delimiter, otherwise we're not going to quit until eof.

    // Eat whitespace after the opening bracket to see what's next.
    if (in.good()) std::ws(in);

    // If we're at eof now it must be because the open bracket was the
    // last non-white character in the stream, which is an error.
    if (!in.good()) {
        if (in.eof()) {
            assert(lookForCloser); // or we haven't read anything that could eof
            in.setstate(std::ios::failbit);
        }
        return in;
    }

    // istream is good and next character is non-white; ready to read first
    // value or terminator.

    // We need to figure out whether the elements are space- or comma-
    // separated and then insist on consistency.
    bool commaOK = true, commaRequired = false;
    bool terminatorSeen = false;
    int nextIndex = 0;
    while (true) {
        char c;

        // Here at the top of this loop, we have already successfully read 
        // n=nextIndex values of type T. For fixed-size reads, it might be
        // the case that n==numRequired already, but we still may need to
        // look for a closing bracket before we can declare victory.
        // The stream is good() (not at eof) but it might be the case that 
        // there is nothing but white space left; we don't know yet because
        // if we have satisfied the fixed-size count and are not expecting
        // a terminator then we should quit without absorbing the trailing
        // white space.
        assert(in.good());

        // Look for closing bracket before trying to read value.
        if (lookForCloser) {
            // Eat white space to find the closing bracket.
            std::ws(in); if (!in.good()) break; // eof?
            ch = in.peek(); assert(ch != EOFch);
            if (!in.good()) break;
            c = (char)ch;
            if (c == closeBracket) {   
                in.get(); // absorb the closing bracket
                terminatorSeen = true; 
                break; 
            }
            // next char not a closing bracket; fall through
        }

        // We didn't look or didn't find a closing bracket. The istream is good 
        // but we might be looking at white space.

        // If we already got all the elements we want, break for final checks.
        if (isFixedSize && (nextIndex == numRequired))
            break; // that's a full count.

        // Look for comma before value, except the first time.
        if (commaOK && nextIndex != 0) {
            // Eat white space to find the comma.
            std::ws(in); if (!in.good()) break; // eof?
            ch = in.peek(); assert(ch != EOFch);
            if (!in.good()) break;
            c = (char)ch;
            if (c == ',') {
                in.get(); // absorb comma
                commaRequired = true; // all commas from now on
            } else { // next char not a comma
                if (commaRequired) // bad, e.g.: v1, v2, v3 v4 
                {   in.setstate(std::ios::failbit); break; }
                else commaOK = false; // saw: v1 v2 (no commas now)
            }
            if (!in.good()) break; // might be eof
        }

        // No closing bracket yet; don't have enough elements; skipped comma 
        // if any; istream is good; might be looking at white space.
        assert(in.good());

        // Now read in an element of type T.
        // The extractor T::operator>>() will ignore leading white space.
        if (!isFixedSize)
            out.resizeKeep(out.size()+1); // grow by one (default consructed)
        in >> out[nextIndex]; if (in.fail()) break;
        ++nextIndex;

        if (!in.good()) break; // might be eof
    }

    // We will get here under a number of circumstances:
    //  - the fail bit is set in the istream, or
    //  - we reached eof
    //  - we saw a closing brace
    //  - we got all the elements we wanted (for a fixed-size read)
    // Note that it is possible that we consumed everything except some
    // trailing white space (meaning we're not technically at eof), but
    // for consistency with built-in operator>>()'s we won't try to absorb
    // that trailing white space.

    if (!in.fail()) {
        if (lookForCloser && !terminatorSeen)
            in.setstate(std::ios::failbit); // missing terminator

        if (isFixedSize && nextIndex != numRequired)
            in.setstate(std::ios::failbit); // wrong number of values
    }

    return in;
}



//------------------------------------------------------------------------------
//                          RELATED GLOBAL OPERATORS
//------------------------------------------------------------------------------
// These are logically part of the Matrix_<T> class but are not actually 
// class members; that is, they are in the SimTK namespace.


template <class T> inline std::ostream&
operator<<(std::ostream& o, const VectorBase<T>& v)
{ o << "~["; for(int i=0;i<v.size();++i) o<<(i>0?" ":"")<<v[i]; 
    return o << "]"; }

template <class T> inline std::ostream&
operator<<(std::ostream& o, const RowVectorBase<T>& v)
{ o << "["; for(int i=0;i<v.size();++i) o<<(i>0?" ":"")<<v[i]; 
    return o << "]"; }

template <class T> inline std::ostream&
operator<<(std::ostream& o, const MatrixBase<T>& m) {
    for (int i=0;i<m.nrow();++i)
        o << std::endl << m[i];
    if (m.nrow()) o << std::endl;
    return o; 
}


template <class T> static inline 
std::istream& readVectorFromStream(std::istream& in, Vector_<T>& out)
{   return readVectorFromStreamHelper<T>(in, false /*variable sizez*/, out); }



template <class T> static inline 
std::istream& fillVectorFromStream(std::istream& in, Vector_<T>& out)
{   return readVectorFromStreamHelper<T>(in, true /*fixed size*/, out); }

template <class T> static inline 
std::istream& fillVectorViewFromStream(std::istream& in, VectorView_<T>& out)
{   return readVectorFromStreamHelper<T>(in, true /*fixed size*/, out); }


template <class T> inline
std::istream& operator>>(std::istream& in, Vector_<T>& out) 
{   return readVectorFromStream<T>(in, out); }

template <class T> inline
std::istream& operator>>(std::istream& in, VectorView_<T>& out) 
{   return fillVectorViewFromStream<T>(in, out); }

// Friendly abbreviations for vectors and matrices with scalar elements.

typedef Vector_<Real>           Vector;
typedef Vector_<float>          fVector;
typedef Vector_<double>         dVector;
typedef Vector_<Complex>        ComplexVector;
typedef Vector_<fComplex>       fComplexVector;
typedef Vector_<dComplex>       dComplexVector;

typedef VectorView_<Real>       VectorView;
typedef VectorView_<float>      fVectorView;
typedef VectorView_<double>     dVectorView;
typedef VectorView_<Complex>    ComplexVectorView;
typedef VectorView_<fComplex>   fComplexVectorView;
typedef VectorView_<dComplex>   dComplexVectorView;

typedef RowVector_<Real>        RowVector;
typedef RowVector_<float>       fRowVector;
typedef RowVector_<double>      dRowVector;
typedef RowVector_<Complex>     ComplexRowVector;
typedef RowVector_<fComplex>    fComplexRowVector;
typedef RowVector_<dComplex>    dComplexRowVector;

typedef RowVectorView_<Real>    RowVectorView;
typedef RowVectorView_<float>   fRowVectorView;
typedef RowVectorView_<double>  dRowVectorView;
typedef RowVectorView_<Complex> ComplexRowVectorView;
typedef RowVectorView_<fComplex> fComplexRowVectorView;
typedef RowVectorView_<dComplex> dComplexRowVectorView;

typedef Matrix_<Real>           Matrix;
typedef Matrix_<float>          fMatrix;
typedef Matrix_<double>         dMatrix;
typedef Matrix_<Complex>        ComplexMatrix;
typedef Matrix_<fComplex>       fComplexMatrix;
typedef Matrix_<dComplex>       dComplexMatrix;

typedef MatrixView_<Real>       MatrixView;
typedef MatrixView_<float>      fMatrixView;
typedef MatrixView_<double>     dMatrixView;
typedef MatrixView_<Complex>    ComplexMatrixView;
typedef MatrixView_<fComplex>   fComplexMatrixView;
typedef MatrixView_<dComplex>   dComplexMatrixView;


template <class ELT, class VECTOR_CLASS>
class VectorIterator {
public:
    typedef ELT value_type;
    typedef int difference_type;
    typedef ELT& reference;
    typedef ELT* pointer;
    typedef std::random_access_iterator_tag iterator_category;
    VectorIterator(VECTOR_CLASS& vector, int index) : vector(vector), index(index) {
    }
    VectorIterator(const VectorIterator& iter) : vector(iter.vector), index(iter.index) {
    }
    VectorIterator& operator=(const VectorIterator& iter) {
        vector = iter.vector;
        index = iter.index;
        return *this;
    }
    ELT& operator*() {
        assert (index >= 0 && index < vector.size());
        return vector[index];
    }
    ELT& operator[](int i) {
        assert (i >= 0 && i < vector.size());
        return vector[i];
    }
    VectorIterator operator++() {
        assert (index < vector.size());
        ++index;
        return *this;
    }
    VectorIterator operator++(int) {
        assert (index < vector.size());
        VectorIterator current = *this;
        ++index;
        return current;
    }
    VectorIterator operator--() {
        assert (index > 0);
        --index;
        return *this;
    }
    VectorIterator operator--(int) {
        assert (index > 0);
        VectorIterator current = *this;
        --index;
        return current;
    }
    bool operator<(VectorIterator iter) const {
        return (index < iter.index);
    }
    bool operator>(VectorIterator iter) const {
        return (index > iter.index);
    }
    bool operator<=(VectorIterator iter) const {
        return (index <= iter.index);
    }
    bool operator>=(VectorIterator iter) const {
        return (index >= iter.index);
    }
    int operator-(VectorIterator iter) const {
        return (index - iter.index);
    }
    VectorIterator operator-(int n) const {
        return VectorIterator(vector, index-n);
    }
    VectorIterator operator+(int n) const {
        return VectorIterator(vector, index+n);
    }
    bool operator==(VectorIterator iter) const {
        return (index == iter.index);
    }
    bool operator!=(VectorIterator iter) const {
        return (index != iter.index);
    }
private:
    VECTOR_CLASS& vector;
    int index;
};

} //namespace SimTK

#endif //SimTK_SIMMATRIX_BIGMATRIX_H_