Simbody  3.5
MatrixBase.h
Go to the documentation of this file.
00001 #ifndef SimTK_SIMMATRIX_MATRIXBASE_H_
00002 #define SimTK_SIMMATRIX_MATRIXBASE_H_
00003 
00004 /* -------------------------------------------------------------------------- *
00005  *                       Simbody(tm): SimTKcommon                             *
00006  * -------------------------------------------------------------------------- *
00007  * This is part of the SimTK biosimulation toolkit originating from           *
00008  * Simbios, the NIH National Center for Physics-Based Simulation of           *
00009  * Biological Structures at Stanford, funded under the NIH Roadmap for        *
00010  * Medical Research, grant U54 GM072970. See https://simtk.org/home/simbody.  *
00011  *                                                                            *
00012  * Portions copyright (c) 2005-13 Stanford University and the Authors.        *
00013  * Authors: Michael Sherman                                                   *
00014  * Contributors:                                                              *
00015  *                                                                            *
00016  * Licensed under the Apache License, Version 2.0 (the "License"); you may    *
00017  * not use this file except in compliance with the License. You may obtain a  *
00018  * copy of the License at http://www.apache.org/licenses/LICENSE-2.0.         *
00019  *                                                                            *
00020  * Unless required by applicable law or agreed to in writing, software        *
00021  * distributed under the License is distributed on an "AS IS" BASIS,          *
00022  * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.   *
00023  * See the License for the specific language governing permissions and        *
00024  * limitations under the License.                                             *
00025  * -------------------------------------------------------------------------- */
00026 
00031 namespace SimTK {
00032 
00033 //==============================================================================
00034 //                               MATRIX BASE
00035 //==============================================================================
00067 //  ----------------------------------------------------------------------------
00068 template <class ELT> class MatrixBase {  
00069 public:
00070     // These typedefs are handy, but despite appearances you cannot 
00071     // treat a MatrixBase as a composite numerical type. That is,
00072     // CNT<MatrixBase> will not compile, or if it does it won't be
00073     // meaningful.
00074 
00075     typedef ELT                                 E;
00076     typedef typename CNT<E>::TNeg               ENeg;
00077     typedef typename CNT<E>::TWithoutNegator    EWithoutNegator;
00078     typedef typename CNT<E>::TReal              EReal;
00079     typedef typename CNT<E>::TImag              EImag;
00080     typedef typename CNT<E>::TComplex           EComplex;
00081     typedef typename CNT<E>::THerm              EHerm;       
00082     typedef typename CNT<E>::TPosTrans          EPosTrans;
00083 
00084     typedef typename CNT<E>::TAbs               EAbs;
00085     typedef typename CNT<E>::TStandard          EStandard;
00086     typedef typename CNT<E>::TInvert            EInvert;
00087     typedef typename CNT<E>::TNormalize         ENormalize;
00088     typedef typename CNT<E>::TSqHermT           ESqHermT;
00089     typedef typename CNT<E>::TSqTHerm           ESqTHerm;
00090 
00091     typedef typename CNT<E>::Scalar             EScalar;
00092     typedef typename CNT<E>::Number             ENumber;
00093     typedef typename CNT<E>::StdNumber          EStdNumber;
00094     typedef typename CNT<E>::Precision          EPrecision;
00095     typedef typename CNT<E>::ScalarNormSq       EScalarNormSq;
00096 
00097     typedef EScalar    Scalar;        // the underlying Scalar type
00098     typedef ENumber    Number;        // negator removed from Scalar
00099     typedef EStdNumber StdNumber;     // conjugate goes to complex
00100     typedef EPrecision Precision;     // complex removed from StdNumber
00101     typedef EScalarNormSq  ScalarNormSq;      // type of scalar^2
00102 
00103     typedef MatrixBase<E>                T;
00104     typedef MatrixBase<ENeg>             TNeg;
00105     typedef MatrixBase<EWithoutNegator>  TWithoutNegator;
00106     typedef MatrixBase<EReal>            TReal;
00107     typedef MatrixBase<EImag>            TImag;
00108     typedef MatrixBase<EComplex>         TComplex;
00109     typedef MatrixBase<EHerm>            THerm; 
00110     typedef MatrixBase<E>                TPosTrans;
00111 
00112     typedef MatrixBase<EAbs>             TAbs;
00113     typedef MatrixBase<EStandard>        TStandard;
00114     typedef MatrixBase<EInvert>          TInvert;
00115     typedef MatrixBase<ENormalize>       TNormalize;
00116     typedef MatrixBase<ESqHermT>         TSqHermT;  // ~Mat*Mat
00117     typedef MatrixBase<ESqTHerm>         TSqTHerm;  // Mat*~Mat
00118 
00119     const MatrixCommitment& getCharacterCommitment() const {return helper.getCharacterCommitment();}
00120     const MatrixCharacter& getMatrixCharacter()     const {return helper.getMatrixCharacter();}
00121 
00124     void commitTo(const MatrixCommitment& mc)
00125     {   helper.commitTo(mc); }
00126 
00127     // This gives the resulting matrix type when (m(i,j) op P) is applied to each element.
00128     // It will have element types which are the regular composite result of E op P.
00129     template <class P> struct EltResult { 
00130         typedef MatrixBase<typename CNT<E>::template Result<P>::Mul> Mul;
00131         typedef MatrixBase<typename CNT<E>::template Result<P>::Dvd> Dvd;
00132         typedef MatrixBase<typename CNT<E>::template Result<P>::Add> Add;
00133         typedef MatrixBase<typename CNT<E>::template Result<P>::Sub> Sub;
00134     };
00135 
00137     int  nrow() const {return helper.nrow();}
00139     int  ncol() const {return helper.ncol();}
00140 
00148     ptrdiff_t nelt() const {return helper.nelt();}
00149 
00151     bool isResizeable() const {return getCharacterCommitment().isResizeable();}
00152 
00153     enum { 
00154         NScalarsPerElement    = CNT<E>::NActualScalars,
00155         CppNScalarsPerElement = sizeof(E) / sizeof(Scalar)
00156     };
00157   
00161     MatrixBase() : helper(NScalarsPerElement,CppNScalarsPerElement) {}
00162 
00165     MatrixBase(int m, int n) 
00166     :   helper(NScalarsPerElement,CppNScalarsPerElement,MatrixCommitment(),m,n) {}
00167 
00173     explicit MatrixBase(const MatrixCommitment& commitment) 
00174     :   helper(NScalarsPerElement,CppNScalarsPerElement,commitment) {}
00175 
00176 
00181     MatrixBase(const MatrixCommitment& commitment, int m, int n) 
00182     :   helper(NScalarsPerElement,CppNScalarsPerElement,commitment,m,n) {}
00183 
00185     MatrixBase(const MatrixBase& b)
00186       : helper(b.helper.getCharacterCommitment(), 
00187                b.helper, typename MatrixHelper<Scalar>::DeepCopy()) { }
00188 
00191     MatrixBase(const TNeg& b)
00192       : helper(b.helper.getCharacterCommitment(),
00193                b.helper, typename MatrixHelper<Scalar>::DeepCopy()) { }
00194     
00197     MatrixBase& copyAssign(const MatrixBase& b) {
00198         helper.copyAssign(b.helper);
00199         return *this;
00200     }
00201     MatrixBase& operator=(const MatrixBase& b) { return copyAssign(b); }
00202 
00203 
00212     MatrixBase& viewAssign(const MatrixBase& src) {
00213         helper.writableViewAssign(const_cast<MatrixHelper<Scalar>&>(src.helper));
00214         return *this;
00215     }
00216 
00217     // default destructor
00218 
00225     MatrixBase(const MatrixCommitment& commitment, int m, int n, const ELT& initialValue) 
00226     :   helper(NScalarsPerElement, CppNScalarsPerElement, commitment, m, n)
00227     {   helper.fillWith(reinterpret_cast<const Scalar*>(&initialValue)); }  
00228 
00239     MatrixBase(const MatrixCommitment& commitment, int m, int n, 
00240                const ELT* cppInitialValuesByRow) 
00241     :   helper(NScalarsPerElement, CppNScalarsPerElement, commitment, m, n)
00242     {   helper.copyInByRowsFromCpp(reinterpret_cast<const Scalar*>(cppInitialValuesByRow)); }
00243      
00255 
00257     MatrixBase(const MatrixCommitment& commitment, 
00258                const MatrixCharacter&  character, 
00259                int spacing, const Scalar* data) // read only data
00260     :   helper(NScalarsPerElement, CppNScalarsPerElement, 
00261                commitment, character, spacing, data) {}  
00262 
00264     MatrixBase(const MatrixCommitment& commitment, 
00265                const MatrixCharacter&  character, 
00266                int spacing, Scalar* data) // writable data
00267     :   helper(NScalarsPerElement, CppNScalarsPerElement, 
00268                commitment, character, spacing, data) {}  
00270         
00271     // Create a new MatrixBase from an existing helper. Both shallow and deep copies are possible.
00272     MatrixBase(const MatrixCommitment& commitment, 
00273                MatrixHelper<Scalar>&   source, 
00274                const typename MatrixHelper<Scalar>::ShallowCopy& shallow) 
00275     :   helper(commitment, source, shallow) {}
00276     MatrixBase(const MatrixCommitment&      commitment, 
00277                const MatrixHelper<Scalar>&  source, 
00278                const typename MatrixHelper<Scalar>::ShallowCopy& shallow) 
00279     :   helper(commitment, source, shallow) {}
00280     MatrixBase(const MatrixCommitment&      commitment, 
00281                const MatrixHelper<Scalar>&  source, 
00282                const typename MatrixHelper<Scalar>::DeepCopy& deep)    
00283     :   helper(commitment, source, deep) {}
00284 
00288     void clear() {helper.clear();}
00289 
00290     MatrixBase& operator*=(const StdNumber& t)  { helper.scaleBy(t);              return *this; }
00291     MatrixBase& operator/=(const StdNumber& t)  { helper.scaleBy(StdNumber(1)/t); return *this; }
00292     MatrixBase& operator+=(const MatrixBase& r) { helper.addIn(r.helper);         return *this; }
00293     MatrixBase& operator-=(const MatrixBase& r) { helper.subIn(r.helper);         return *this; }  
00294 
00295     template <class EE> MatrixBase(const MatrixBase<EE>& b)
00296       : helper(MatrixCommitment(),b.helper, typename MatrixHelper<Scalar>::DeepCopy()) { }
00297 
00298     template <class EE> MatrixBase& operator=(const MatrixBase<EE>& b) 
00299       { helper = b.helper; return *this; }
00300     template <class EE> MatrixBase& operator+=(const MatrixBase<EE>& b) 
00301       { helper.addIn(b.helper); return *this; }
00302     template <class EE> MatrixBase& operator-=(const MatrixBase<EE>& b) 
00303       { helper.subIn(b.helper); return *this; }
00304 
00315     MatrixBase& operator=(const ELT& t) { 
00316         setToZero(); updDiag().setTo(t); 
00317         return *this;
00318     }
00319 
00325     template <class S> inline MatrixBase&
00326     scalarAssign(const S& s) {
00327         setToZero(); updDiag().setTo(s);
00328         return *this;
00329     }
00330 
00334     template <class S> inline MatrixBase&
00335     scalarAddInPlace(const S& s) {
00336         updDiag().elementwiseAddScalarInPlace(s);
00337     }
00338 
00339 
00343     template <class S> inline MatrixBase&
00344     scalarSubtractInPlace(const S& s) {
00345         updDiag().elementwiseSubtractScalarInPlace(s);
00346     }
00347 
00350     template <class S> inline MatrixBase&
00351     scalarSubtractFromLeftInPlace(const S& s) {
00352         negateInPlace();
00353         updDiag().elementwiseAddScalarInPlace(s); // yes, add
00354     }
00355 
00362     template <class S> inline MatrixBase&
00363     scalarMultiplyInPlace(const S&);
00364 
00368     template <class S> inline MatrixBase&
00369     scalarMultiplyFromLeftInPlace(const S&);
00370 
00377     template <class S> inline MatrixBase&
00378     scalarDivideInPlace(const S&);
00379 
00383     template <class S> inline MatrixBase&
00384     scalarDivideFromLeftInPlace(const S&);
00385 
00386 
00389     template <class EE> inline MatrixBase& 
00390     rowScaleInPlace(const VectorBase<EE>&);
00391 
00394     template <class EE> inline void 
00395     rowScale(const VectorBase<EE>& r, typename EltResult<EE>::Mul& out) const;
00396 
00397     template <class EE> inline typename EltResult<EE>::Mul 
00398     rowScale(const VectorBase<EE>& r) const {
00399         typename EltResult<EE>::Mul out(nrow(), ncol()); rowScale(r,out); return out;
00400     }
00401 
00404     template <class EE> inline MatrixBase& 
00405     colScaleInPlace(const VectorBase<EE>&);
00406 
00407     template <class EE> inline void 
00408     colScale(const VectorBase<EE>& c, typename EltResult<EE>::Mul& out) const;
00409 
00410     template <class EE> inline typename EltResult<EE>::Mul
00411     colScale(const VectorBase<EE>& c) const {
00412         typename EltResult<EE>::Mul out(nrow(), ncol()); colScale(c,out); return out;
00413     }
00414 
00415 
00420     template <class ER, class EC> inline MatrixBase& 
00421     rowAndColScaleInPlace(const VectorBase<ER>& r, const VectorBase<EC>& c);
00422 
00423     template <class ER, class EC> inline void 
00424     rowAndColScale(const VectorBase<ER>& r, const VectorBase<EC>& c, 
00425                    typename EltResult<typename VectorBase<ER>::template EltResult<EC>::Mul>::Mul& out) const;
00426 
00427     template <class ER, class EC> inline typename EltResult<typename VectorBase<ER>::template EltResult<EC>::Mul>::Mul
00428     rowAndColScale(const VectorBase<ER>& r, const VectorBase<EC>& c) const {
00429         typename EltResult<typename VectorBase<ER>::template EltResult<EC>::Mul>::Mul 
00430             out(nrow(), ncol()); 
00431         rowAndColScale(r,c,out); return out;
00432     }
00433 
00441     template <class S> inline MatrixBase&
00442     elementwiseAssign(const S& s);
00443 
00445     MatrixBase& elementwiseAssign(int s)
00446     {   return elementwiseAssign<Real>(Real(s)); }
00447 
00449     MatrixBase& elementwiseInvertInPlace();
00450 
00451     void elementwiseInvert(MatrixBase<typename CNT<E>::TInvert>& out) const;
00452 
00453     MatrixBase<typename CNT<E>::TInvert> elementwiseInvert() const {
00454         MatrixBase<typename CNT<E>::TInvert> out(nrow(), ncol());
00455         elementwiseInvert(out);
00456         return out;
00457     }
00458 
00466     template <class S> inline MatrixBase&
00467     elementwiseAddScalarInPlace(const S& s);
00468 
00469     template <class S> inline void
00470     elementwiseAddScalar(const S& s, typename EltResult<S>::Add&) const;
00471 
00472     template <class S> inline typename EltResult<S>::Add
00473     elementwiseAddScalar(const S& s) const {
00474         typename EltResult<S>::Add out(nrow(), ncol());
00475         elementwiseAddScalar(s,out);
00476         return out;
00477     }
00478 
00486     template <class S> inline MatrixBase&
00487     elementwiseSubtractScalarInPlace(const S& s);
00488 
00489     template <class S> inline void
00490     elementwiseSubtractScalar(const S& s, typename EltResult<S>::Sub&) const;
00491 
00492     template <class S> inline typename EltResult<S>::Sub
00493     elementwiseSubtractScalar(const S& s) const {
00494         typename EltResult<S>::Sub out(nrow(), ncol());
00495         elementwiseSubtractScalar(s,out);
00496         return out;
00497     }
00498 
00507     template <class S> inline MatrixBase&
00508     elementwiseSubtractFromScalarInPlace(const S& s);
00509 
00510     template <class S> inline void
00511     elementwiseSubtractFromScalar(
00512         const S&, 
00513         typename MatrixBase<S>::template EltResult<E>::Sub&) const;
00514 
00515     template <class S> inline typename MatrixBase<S>::template EltResult<E>::Sub
00516     elementwiseSubtractFromScalar(const S& s) const {
00517         typename MatrixBase<S>::template EltResult<E>::Sub out(nrow(), ncol());
00518         elementwiseSubtractFromScalar<S>(s,out);
00519         return out;
00520     }
00521 
00523     template <class EE> inline MatrixBase& 
00524     elementwiseMultiplyInPlace(const MatrixBase<EE>&);
00525 
00526     template <class EE> inline void 
00527     elementwiseMultiply(const MatrixBase<EE>&, typename EltResult<EE>::Mul&) const;
00528 
00529     template <class EE> inline typename EltResult<EE>::Mul 
00530     elementwiseMultiply(const MatrixBase<EE>& m) const {
00531         typename EltResult<EE>::Mul out(nrow(), ncol()); 
00532         elementwiseMultiply<EE>(m,out); 
00533         return out;
00534     }
00535 
00537     template <class EE> inline MatrixBase& 
00538     elementwiseMultiplyFromLeftInPlace(const MatrixBase<EE>&);
00539 
00540     template <class EE> inline void 
00541     elementwiseMultiplyFromLeft(
00542         const MatrixBase<EE>&, 
00543         typename MatrixBase<EE>::template EltResult<E>::Mul&) const;
00544 
00545     template <class EE> inline typename MatrixBase<EE>::template EltResult<E>::Mul 
00546     elementwiseMultiplyFromLeft(const MatrixBase<EE>& m) const {
00547         typename EltResult<EE>::Mul out(nrow(), ncol()); 
00548         elementwiseMultiplyFromLeft<EE>(m,out); 
00549         return out;
00550     }
00551 
00553     template <class EE> inline MatrixBase& 
00554     elementwiseDivideInPlace(const MatrixBase<EE>&);
00555 
00556     template <class EE> inline void 
00557     elementwiseDivide(const MatrixBase<EE>&, typename EltResult<EE>::Dvd&) const;
00558 
00559     template <class EE> inline typename EltResult<EE>::Dvd 
00560     elementwiseDivide(const MatrixBase<EE>& m) const {
00561         typename EltResult<EE>::Dvd out(nrow(), ncol()); 
00562         elementwiseDivide<EE>(m,out); 
00563         return out;
00564     }
00565 
00567     template <class EE> inline MatrixBase& 
00568     elementwiseDivideFromLeftInPlace(const MatrixBase<EE>&);
00569 
00570     template <class EE> inline void 
00571     elementwiseDivideFromLeft(
00572         const MatrixBase<EE>&,
00573         typename MatrixBase<EE>::template EltResult<E>::Dvd&) const;
00574 
00575     template <class EE> inline typename MatrixBase<EE>::template EltResult<EE>::Dvd 
00576     elementwiseDivideFromLeft(const MatrixBase<EE>& m) const {
00577         typename MatrixBase<EE>::template EltResult<E>::Dvd out(nrow(), ncol()); 
00578         elementwiseDivideFromLeft<EE>(m,out); 
00579         return out;
00580     }
00581 
00583     MatrixBase& setTo(const ELT& t) {helper.fillWith(reinterpret_cast<const Scalar*>(&t)); return *this;}
00584     MatrixBase& setToNaN() {helper.fillWithScalar(CNT<StdNumber>::getNaN()); return *this;}
00585     MatrixBase& setToZero() {helper.fillWithScalar(StdNumber(0)); return *this;}
00586  
00587     // View creating operators.   
00588     inline RowVectorView_<ELT> row(int i) const;   // select a row
00589     inline RowVectorView_<ELT> updRow(int i);
00590     inline VectorView_<ELT>    col(int j) const;   // select a column
00591     inline VectorView_<ELT>    updCol(int j);
00592 
00593     RowVectorView_<ELT> operator[](int i) const {return row(i);}
00594     RowVectorView_<ELT> operator[](int i)       {return updRow(i);}
00595     VectorView_<ELT>    operator()(int j) const {return col(j);}
00596     VectorView_<ELT>    operator()(int j)       {return updCol(j);}
00597      
00598     // Select a block.
00599     inline MatrixView_<ELT> block(int i, int j, int m, int n) const;
00600     inline MatrixView_<ELT> updBlock(int i, int j, int m, int n);
00601 
00602     MatrixView_<ELT> operator()(int i, int j, int m, int n) const
00603       { return block(i,j,m,n); }
00604     MatrixView_<ELT> operator()(int i, int j, int m, int n)
00605       { return updBlock(i,j,m,n); }
00606  
00607     // Hermitian transpose.
00608     inline MatrixView_<EHerm> transpose() const;
00609     inline MatrixView_<EHerm> updTranspose();
00610 
00611     MatrixView_<EHerm> operator~() const {return transpose();}
00612     MatrixView_<EHerm> operator~()       {return updTranspose();}
00613 
00616     inline VectorView_<ELT> diag() const;
00619     inline VectorView_<ELT> updDiag();
00622     VectorView_<ELT> diag() {return updDiag();}
00623 
00624     // Create a view of the real or imaginary elements. TODO
00625     //inline MatrixView_<EReal> real() const;
00626     //inline MatrixView_<EReal> updReal();
00627     //inline MatrixView_<EImag> imag() const;
00628     //inline MatrixView_<EImag> updImag();
00629 
00630     // Overload "real" and "imag" for both read and write as a nod to convention. TODO
00631     //MatrixView_<EReal> real() {return updReal();}
00632     //MatrixView_<EReal> imag() {return updImag();}
00633 
00634     // TODO: this routine seems ill-advised but I need it for the IVM port at the moment
00635     TInvert invert() const {  // return a newly-allocated inverse; dump negator 
00636         TInvert m(*this);
00637         m.helper.invertInPlace();
00638         return m;   // TODO - bad: makes an extra copy
00639     }
00640 
00641     void invertInPlace() {helper.invertInPlace();}
00642 
00644     void dump(const char* msg=0) const {
00645         helper.dump(msg);
00646     }
00647 
00656     const ELT& getElt(int i, int j) const { return *reinterpret_cast<const ELT*>(helper.getElt(i,j)); }
00657     ELT&       updElt(int i, int j)       { return *reinterpret_cast<      ELT*>(helper.updElt(i,j)); }
00658 
00659     const ELT& operator()(int i, int j) const {return getElt(i,j);}
00660     ELT&       operator()(int i, int j)       {return updElt(i,j);}
00661 
00666     void getAnyElt(int i, int j, ELT& value) const
00667     {   helper.getAnyElt(i,j,reinterpret_cast<Scalar*>(&value)); }
00668     ELT getAnyElt(int i, int j) const {ELT e; getAnyElt(i,j,e); return e;}
00669 
00672     // TODO: very slow! Should be optimized at least for the case
00673     //       where ELT is a Scalar.
00674     ScalarNormSq scalarNormSqr() const {
00675         const int nr=nrow(), nc=ncol();
00676         ScalarNormSq sum(0);
00677         for(int j=0;j<nc;++j) 
00678             for (int i=0; i<nr; ++i)
00679                 sum += CNT<E>::scalarNormSqr((*this)(i,j));
00680         return sum;
00681     }
00682 
00686     // TODO: very slow! Should be optimized at least for the case
00687     //       where ELT is a Scalar.
00688     void abs(TAbs& mabs) const {
00689         const int nr=nrow(), nc=ncol();
00690         mabs.resize(nr,nc);
00691         for(int j=0;j<nc;++j) 
00692             for (int i=0; i<nr; ++i)
00693                 mabs(i,j) = CNT<E>::abs((*this)(i,j));
00694     }
00695 
00699     TAbs abs() const { TAbs mabs; abs(mabs); return mabs; }
00700 
00711     TStandard standardize() const {
00712         const int nr=nrow(), nc=ncol();
00713         TStandard mstd(nr, nc);
00714         for(int j=0;j<nc;++j) 
00715             for (int i=0; i<nr; ++i)
00716                 mstd(i,j) = CNT<E>::standardize((*this)(i,j));
00717         return mstd;
00718     }
00719 
00723     ScalarNormSq normSqr() const { return scalarNormSqr(); }
00724     // TODO -- not good; unnecessary overflow
00725     typename CNT<ScalarNormSq>::TSqrt 
00726         norm() const { return CNT<ScalarNormSq>::sqrt(scalarNormSqr()); }
00727 
00731     typename CNT<ScalarNormSq>::TSqrt 
00732     normRMS() const {
00733         if (!CNT<ELT>::IsScalar)
00734             SimTK_THROW1(Exception::Cant, "normRMS() only defined for scalar elements");
00735         if (nelt() == 0)
00736             return typename CNT<ScalarNormSq>::TSqrt(0);
00737         return CNT<ScalarNormSq>::sqrt(scalarNormSqr()/nelt());
00738     }
00739 
00741     RowVector_<ELT> colSum() const {
00742         const int cols = ncol();
00743         RowVector_<ELT> row(cols);
00744         for (int j = 0; j < cols; ++j)
00745             helper.colSum(j, reinterpret_cast<Scalar*>(&row[j]));
00746         return row;
00747     }
00749     RowVector_<ELT> sum() const {return colSum();}
00750 
00752     Vector_<ELT> rowSum() const {
00753         const int rows = nrow();
00754         Vector_<ELT> col(rows);
00755         for (int i = 0; i < rows; ++i)
00756             helper.rowSum(i, reinterpret_cast<Scalar*>(&col[i]));
00757         return col;
00758     }
00759 
00760     //TODO: make unary +/- return a self-view so they won't reallocate?
00761     const MatrixBase& operator+() const {return *this; }
00762     const TNeg&       negate()    const {return *reinterpret_cast<const TNeg*>(this); }
00763     TNeg&             updNegate()       {return *reinterpret_cast<TNeg*>(this); }
00764 
00765     const TNeg&       operator-() const {return negate();}
00766     TNeg&             operator-()       {return updNegate();}
00767 
00768     MatrixBase& negateInPlace() {(*this) *= EPrecision(-1); return *this;}
00769  
00774     MatrixBase& resize(int m, int n)     { helper.resize(m,n); return *this; }
00780     MatrixBase& resizeKeep(int m, int n) { helper.resizeKeep(m,n); return *this; }
00781 
00782     // This prevents shape changes in a Matrix that would otherwise allow it. No harm if is
00783     // are called on a Matrix that is locked already; it always succeeds.
00784     void lockShape() {helper.lockShape();}
00785 
00786     // This allows shape changes again for a Matrix which was constructed to allow them
00787     // but had them locked with the above routine. No harm if this is called on a Matrix
00788     // that is already unlocked, but it is not allowed to call this on a Matrix which
00789     // *never* allowed resizing. An exception will be thrown in that case.
00790     void unlockShape() {helper.unlockShape();}
00791     
00792     // An assortment of handy conversions
00793     const MatrixView_<ELT>& getAsMatrixView() const { return *reinterpret_cast<const MatrixView_<ELT>*>(this); }
00794     MatrixView_<ELT>&       updAsMatrixView()       { return *reinterpret_cast<      MatrixView_<ELT>*>(this); } 
00795     const Matrix_<ELT>&     getAsMatrix()     const { return *reinterpret_cast<const Matrix_<ELT>*>(this); }
00796     Matrix_<ELT>&           updAsMatrix()           { return *reinterpret_cast<      Matrix_<ELT>*>(this); }
00797          
00798     const VectorView_<ELT>& getAsVectorView() const 
00799       { assert(ncol()==1); return *reinterpret_cast<const VectorView_<ELT>*>(this); }
00800     VectorView_<ELT>&       updAsVectorView()       
00801       { assert(ncol()==1); return *reinterpret_cast<      VectorView_<ELT>*>(this); } 
00802     const Vector_<ELT>&     getAsVector()     const 
00803       { assert(ncol()==1); return *reinterpret_cast<const Vector_<ELT>*>(this); }
00804     Vector_<ELT>&           updAsVector()           
00805       { assert(ncol()==1); return *reinterpret_cast<      Vector_<ELT>*>(this); }
00806     const VectorBase<ELT>& getAsVectorBase() const 
00807       { assert(ncol()==1); return *reinterpret_cast<const VectorBase<ELT>*>(this); }
00808     VectorBase<ELT>&       updAsVectorBase()       
00809       { assert(ncol()==1); return *reinterpret_cast<      VectorBase<ELT>*>(this); } 
00810                 
00811     const RowVectorView_<ELT>& getAsRowVectorView() const 
00812       { assert(nrow()==1); return *reinterpret_cast<const RowVectorView_<ELT>*>(this); }
00813     RowVectorView_<ELT>&       updAsRowVectorView()       
00814       { assert(nrow()==1); return *reinterpret_cast<      RowVectorView_<ELT>*>(this); } 
00815     const RowVector_<ELT>&     getAsRowVector()     const 
00816       { assert(nrow()==1); return *reinterpret_cast<const RowVector_<ELT>*>(this); }
00817     RowVector_<ELT>&           updAsRowVector()           
00818       { assert(nrow()==1); return *reinterpret_cast<      RowVector_<ELT>*>(this); }        
00819     const RowVectorBase<ELT>& getAsRowVectorBase() const 
00820       { assert(nrow()==1); return *reinterpret_cast<const RowVectorBase<ELT>*>(this); }
00821     RowVectorBase<ELT>&       updAsRowVectorBase()       
00822       { assert(nrow()==1); return *reinterpret_cast<      RowVectorBase<ELT>*>(this); } 
00823 
00824     // Access to raw data. We have to return the raw data
00825     // pointer as pointer-to-scalar because we may pack the elements tighter
00826     // than a C++ array would.
00827 
00831     int getNScalarsPerElement()  const {return NScalarsPerElement;}
00832 
00836     int getPackedSizeofElement() const {return NScalarsPerElement*sizeof(Scalar);}
00837 
00838     bool hasContiguousData() const {return helper.hasContiguousData();}
00839     ptrdiff_t getContiguousScalarDataLength() const {
00840         return helper.getContiguousDataLength();
00841     }
00842     const Scalar* getContiguousScalarData() const {
00843         return helper.getContiguousData();
00844     }
00845     Scalar* updContiguousScalarData() {
00846         return helper.updContiguousData();
00847     }
00848     void replaceContiguousScalarData(Scalar* newData, ptrdiff_t length, bool takeOwnership) {
00849         helper.replaceContiguousData(newData,length,takeOwnership);
00850     }
00851     void replaceContiguousScalarData(const Scalar* newData, ptrdiff_t length) {
00852         helper.replaceContiguousData(newData,length);
00853     }
00854     void swapOwnedContiguousScalarData(Scalar* newData, ptrdiff_t length, Scalar*& oldData) {
00855         helper.swapOwnedContiguousData(newData,length,oldData);
00856     }
00857 
00862     explicit MatrixBase(MatrixHelperRep<Scalar>* hrep) : helper(hrep) {}
00863 
00864 protected:
00865     const MatrixHelper<Scalar>& getHelper() const {return helper;}
00866     MatrixHelper<Scalar>&       updHelper()       {return helper;}
00867 
00868 private:
00869     MatrixHelper<Scalar> helper; // this is just one pointer
00870 
00871     template <class EE> friend class MatrixBase;
00872 
00873     // ============================= Unimplemented =============================
00874     // This routine is useful for implementing friendlier Matrix expressions and operators.
00875     // It maps closely to the Level-3 BLAS family of pxxmm() routines like sgemm(). The
00876     // operation performed assumes that "this" is the result, and that "this" has 
00877     // already been sized correctly to receive the result. We'll compute
00878     //     this = beta*this + alpha*A*B
00879     // If beta is 0 then "this" can be uninitialized. If alpha is 0 we promise not
00880     // to look at A or B. The routine can work efficiently even if A and/or B are transposed
00881     // by their views, so an expression like
00882     //        C += s * ~A * ~B
00883     // can be performed with the single equivalent call
00884     //        C.matmul(1., s, Tr(A), Tr(B))
00885     // where Tr(A) indicates a transposed view of the original A.
00886     // The ultimate efficiency of this call depends on the data layout and views which
00887     // are used for the three matrices.
00888     // NOTE: neither A nor B can be the same matrix as 'this', nor views of the same data
00889     // which would expose elements of 'this' that will be modified by this operation.
00890     template <class ELT_A, class ELT_B>
00891     MatrixBase& matmul(const StdNumber& beta,   // applied to 'this'
00892                        const StdNumber& alpha, const MatrixBase<ELT_A>& A, const MatrixBase<ELT_B>& B)
00893     {
00894         helper.matmul(beta,alpha,A.helper,B.helper);
00895         return *this;
00896     }
00897 
00898 };
00899 
00900 } //namespace SimTK
00901 
00902 #endif // SimTK_SIMMATRIX_MATRIXBASE_H_
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines