MatrixBase.h

Go to the documentation of this file.

#ifndef SimTK_SIMMATRIX_MATRIXBASE_H_
#define SimTK_SIMMATRIX_MATRIXBASE_H_

/* -------------------------------------------------------------------------- *
 *                       Simbody(tm): SimTKcommon                             *
 * -------------------------------------------------------------------------- *
 * This is part of the SimTK 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/home/simbody.  *
 *                                                                            *
 * Portions copyright (c) 2005-13 Stanford University and the Authors.        *
 * Authors: Michael Sherman                                                   *
 * Contributors:                                                              *
 *                                                                            *
 * Licensed under the Apache License, Version 2.0 (the "License"); you may    *
 * not use this file except in compliance with the License. You may obtain a  *
 * copy of the License at http://www.apache.org/licenses/LICENSE-2.0.         *
 *                                                                            *
 * Unless required by applicable law or agreed to in writing, software        *
 * distributed under the License is distributed on an "AS IS" BASIS,          *
 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.   *
 * See the License for the specific language governing permissions and        *
 * limitations under the License.                                             *
 * -------------------------------------------------------------------------- */

namespace SimTK {

//==============================================================================
//                               MATRIX BASE
//==============================================================================
//  ----------------------------------------------------------------------------
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<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& elementwiseAssign(int s)
    {   return elementwiseAssign<Real>(Real(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.   
    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();}

    inline VectorView_<ELT> diag() const;
    inline VectorView_<ELT> updDiag();
    VectorView_<ELT> diag() {return 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);
    }

    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());
    }

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

    Vector_<ELT> rowSum() const {
        const int rows = nrow();
        Vector_<ELT> col(rows);
        for (int i = 0; i < rows; ++i)
            helper.rowSum(i, reinterpret_cast<Scalar*>(&col[i]));
        return col;
    }

    //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;

    // ============================= Unimplemented =============================
    // 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;
    }

};

} //namespace SimTK

#endif // SimTK_SIMMATRIX_MATRIXBASE_H_