BigMatrix.h

Go to the documentation of this file.
00001 #ifndef SimTK_SIMMATRIX_BIGMATRIX_H_
00002 #define SimTK_SIMMATRIX_BIGMATRIX_H_
00003 
00004 /* -------------------------------------------------------------------------- *
00005  *                      SimTK Core: SimTK Simmatrix(tm)                       *
00006  * -------------------------------------------------------------------------- *
00007  * This is part of the SimTK Core 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.               *
00011  *                                                                            *
00012  * Portions copyright (c) 2005-8 Stanford University and the Authors.         *
00013  * Authors: Michael Sherman                                                   *
00014  * Contributors:                                                              *
00015  *                                                                            *
00016  * Permission is hereby granted, free of charge, to any person obtaining a    *
00017  * copy of this software and associated documentation files (the "Software"), *
00018  * to deal in the Software without restriction, including without limitation  *
00019  * the rights to use, copy, modify, merge, publish, distribute, sublicense,   *
00020  * and/or sell copies of the Software, and to permit persons to whom the      *
00021  * Software is furnished to do so, subject to the following conditions:       *
00022  *                                                                            *
00023  * The above copyright notice and this permission notice shall be included in *
00024  * all copies or substantial portions of the Software.                        *
00025  *                                                                            *
00026  * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR *
00027  * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,   *
00028  * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL    *
00029  * THE AUTHORS, CONTRIBUTORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM,    *
00030  * DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR      *
00031  * OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE  *
00032  * USE OR OTHER DEALINGS IN THE SOFTWARE.                                     *
00033  * -------------------------------------------------------------------------- */
00034 
00108 #include "SimTKcommon/Scalar.h"
00109 #include "SimTKcommon/SmallMatrix.h"
00110 
00111 #include "SimTKcommon/internal/MatrixHelper.h"
00112 
00113 #include <iostream>
00114 #include <cassert>
00115 #include <complex>
00116 #include <cstddef>
00117 #include <limits>
00118 
00119 namespace SimTK {
00120 
00121 template <class ELT>    class MatrixBase;
00122 template <class ELT>    class VectorBase;
00123 template <class ELT>    class RowVectorBase;
00124 
00125 template <class T=Real> class MatrixView_;
00126 template <class T=Real> class TmpMatrixView_;
00127 template <class T=Real> class Matrix_;
00128 template <class T=Real> class DeadMatrix_;
00129 
00130 template <class T=Real> class VectorView_;
00131 template <class T=Real> class TmpVectorView_;
00132 template <class T=Real> class Vector_;
00133 template <class T=Real> class DeadVector_;
00134 
00135 template <class T=Real> class RowVectorView_;
00136 template <class T=Real> class TmpRowVectorView_;
00137 template <class T=Real> class RowVector_;
00138 template <class T=Real> class DeadRowVector_;
00139 
00140 template <class ELT, class VECTOR_CLASS> class VectorIterator;
00141 
00164 template <class ELT> class MatrixBase {  
00165 public:
00166     // These typedefs are handy, but despite appearances you cannot 
00167     // treat a MatrixBase as a composite numerical type. That is,
00168     // CNT<MatrixBase> will not compile, or if it does it won't be
00169     // meaningful.
00170 
00171     typedef ELT                                 E;
00172     typedef typename CNT<E>::TNeg               ENeg;
00173     typedef typename CNT<E>::TWithoutNegator    EWithoutNegator;
00174     typedef typename CNT<E>::TReal              EReal;
00175     typedef typename CNT<E>::TImag              EImag;
00176     typedef typename CNT<E>::TComplex           EComplex;
00177     typedef typename CNT<E>::THerm              EHerm;       
00178     typedef typename CNT<E>::TPosTrans          EPosTrans;
00179 
00180     typedef typename CNT<E>::TAbs               EAbs;
00181     typedef typename CNT<E>::TStandard          EStandard;
00182     typedef typename CNT<E>::TInvert            EInvert;
00183     typedef typename CNT<E>::TNormalize         ENormalize;
00184     typedef typename CNT<E>::TSqHermT           ESqHermT;
00185     typedef typename CNT<E>::TSqTHerm           ESqTHerm;
00186 
00187     typedef typename CNT<E>::Scalar             EScalar;
00188     typedef typename CNT<E>::Number             ENumber;
00189     typedef typename CNT<E>::StdNumber          EStdNumber;
00190     typedef typename CNT<E>::Precision          EPrecision;
00191     typedef typename CNT<E>::ScalarSq           EScalarSq;
00192 
00193     typedef EScalar    Scalar;        // the underlying Scalar type
00194     typedef ENumber    Number;        // negator removed from Scalar
00195     typedef EStdNumber StdNumber;     // conjugate goes to complex
00196     typedef EPrecision Precision;     // complex removed from StdNumber
00197     typedef EScalarSq  ScalarSq;
00198 
00199     typedef MatrixBase<E>                T;
00200     typedef MatrixBase<ENeg>             TNeg;
00201     typedef MatrixBase<EWithoutNegator>  TWithoutNegator;
00202     typedef MatrixBase<EReal>            TReal;
00203     typedef MatrixBase<EImag>            TImag;
00204     typedef MatrixBase<EComplex>         TComplex;
00205     typedef MatrixBase<EHerm>            THerm; 
00206     typedef MatrixBase<E>                TPosTrans;
00207 
00208     typedef MatrixBase<EAbs>             TAbs;
00209     typedef MatrixBase<EStandard>        TStandard;
00210     typedef MatrixBase<EInvert>          TInvert;
00211     typedef MatrixBase<ENormalize>       TNormalize;
00212     typedef MatrixBase<ESqHermT>         TSqHermT;  // ~Mat*Mat
00213     typedef MatrixBase<ESqTHerm>         TSqTHerm;  // Mat*~Mat
00214 
00215 
00216     void setMatrixStructure(MatrixStructures::Structure structure) {
00217         helper.setMatrixStructure( structure);  // default Uncommitted
00218     }
00219     MatrixStructures::Structure getMatrixStructure() const {
00220         return helper.getMatrixStructure();
00221     }
00222     void setMatrixShape(MatrixShapes::Shape shape) {
00223         helper.setMatrixShape( shape);          // default Uncommitted
00224     }
00225     MatrixShapes::Shape getMatrixShape() const  {
00226         return helper.getMatrixShape();
00227     }
00228     void setMatrixSparsity(MatrixSparseFormats::Sparsity sparsity) {
00229         helper.setMatrixSparsity( sparsity);    // default Uncommitted
00230     }
00231     MatrixSparseFormats::Sparsity getMatrixSparsity() const  {
00232         return helper.getMatrixSparsity();
00233     }
00234     void setMatrixStorage(MatrixStorageFormats::Storage storage) {
00235         helper.setMatrixStorage( storage);      // default Uncommitted
00236     }
00237     MatrixStorageFormats::Storage getMatrixStorage() const  {
00238         return helper.getMatrixStorage();
00239     }
00240 
00241     void setMatrixCondition(MatrixConditions::Condition condition) {
00242         helper.setMatrixCondition(condition);  //  default Uncommitted
00243     }
00244     MatrixConditions::Condition getMatrixCondition() const  {
00245         return helper.getMatrixCondition();
00246     }
00247 
00248 
00249     // This gives the resulting matrix type when (m(i,j) op P) is applied to each element.
00250     // It will have element types which are the regular composite result of E op P.
00251     template <class P> struct EltResult { 
00252         typedef MatrixBase<typename CNT<E>::template Result<P>::Mul> Mul;
00253         typedef MatrixBase<typename CNT<E>::template Result<P>::Dvd> Dvd;
00254         typedef MatrixBase<typename CNT<E>::template Result<P>::Add> Add;
00255         typedef MatrixBase<typename CNT<E>::template Result<P>::Sub> Sub;
00256     };
00257        
00258     // Product of dimensions may be > 32 bits.
00259     long size() const { return helper.size(); }
00260     int  nrow() const { return helper.nrow(); }
00261     int  ncol() const { return helper.ncol(); }
00262 
00263     // Scalar norm square is sum( squares of all scalars )
00264     // TODO: very slow! Should be optimized at least for the case
00265     //       where ELT is a Scalar.
00266     ScalarSq scalarNormSqr() const {
00267         const int nr=nrow(), nc=ncol();
00268         ScalarSq sum(0);
00269         for(int j=0;j<nc;++j) 
00270             for (int i=0; i<nr; ++i)
00271                 sum += CNT<E>::scalarNormSqr((*this)(i,j));
00272         return sum;
00273     }
00274 
00275     // abs() is elementwise absolute value; that is, the return value has the same
00276     // dimension as this Matrix but with each element replaced by whatever it thinks
00277     // its absolute value is.
00278     // TODO: very slow! Should be optimized at least for the case
00279     //       where ELT is a Scalar.
00280     void abs(TAbs& mabs) const {
00281         const int nr=nrow(), nc=ncol();
00282         mabs.resize(nr,nc);
00283         for(int j=0;j<nc;++j) 
00284             for (int i=0; i<nr; ++i)
00285                 mabs(i,j) = CNT<E>::abs((*this)(i,j));
00286     }
00287 
00288     TAbs abs() const { TAbs mabs; abs(mabs); return mabs; }
00289 
00290     // Return a Matrix of the same shape and contents as this one but
00291     // with the element type converted to one based on the standard
00292     // C++ scalar types: float, double, long double or complex<float>,
00293     // complex<double>, complex<long double>. That is, negator<>
00294     // and conjugate<> are eliminated from the element type by 
00295     // performing any needed negations computationally.
00296     TStandard standardize() const {
00297         const int nr=nrow(), nc=ncol();
00298         TStandard mstd(nr, nc);
00299         for(int j=0;j<nc;++j) 
00300             for (int i=0; i<nr; ++i)
00301                 mstd(i,j) = CNT<E>::standardize((*this)(i,j));
00302         return mstd;
00303     }
00304 
00305     enum { 
00306         NScalarsPerElement    = CNT<E>::NActualScalars,
00307         CppNScalarsPerElement = sizeof(E) / sizeof(Scalar)
00308     };
00309   
00310     MatrixBase() : helper(NScalarsPerElement,CppNScalarsPerElement)  { }
00311 
00312     // This restores the MatrixBase to its just-constructed state.
00313     void clear() {
00314         helper.clear();
00315     }
00316 
00317     
00318     // Copy constructor is a deep copy (not appropriate for views!).    
00319     MatrixBase(const MatrixBase& b)
00320       : helper(b.helper, typename MatrixHelper<Scalar>::DeepCopy()) { }
00321     
00322     // Copy assignment is a deep copy but behavior depends on type of lhs: if view, rhs
00323     // must match. If owner, we reallocate and copy rhs.
00324     MatrixBase& copyAssign(const MatrixBase& b) {
00325         helper.copyAssign(b.helper);
00326         return *this;
00327     }
00328     MatrixBase& operator=(const MatrixBase& b) { return copyAssign(b); }
00329 
00330     // View assignment is a shallow copy, meaning that we disconnect the MatrixBase 
00331     // from whatever it used to refer to (destructing as necessary), then make it a new view
00332     // for the data descriptor referenced by the source.
00333     // CAUTION: we always take the source as const, but that is ignored in 
00334     // determining whether the resulting view is writable. Instead, that is
00335     // inherited from the writability status of the source. We have to do this
00336     // in order to allow temporary view objects to be writable -- the compiler
00337     // creates temporaries like m(i,j,m,n) as const.
00338 
00339     MatrixBase& viewAssign(const MatrixBase& src) {
00340         helper.writableViewAssign(const_cast<MatrixHelper<Scalar>&>(src.helper));
00341         return *this;
00342     }
00343 
00344     // default destructor
00345 
00346 
00347     MatrixBase(int m, int n, bool lockNrow=false, bool lockNcol=false) 
00348       : helper(NScalarsPerElement, CppNScalarsPerElement, m, n, lockNrow, lockNcol)  { }
00349         
00350     // Non-resizeable owner sharing pre-allocated raw data
00351     MatrixBase(int m, int n, int leadingDim, const Scalar* s)
00352       : helper(NScalarsPerElement, CppNScalarsPerElement, m, n, leadingDim, s) { }    // read only
00353     MatrixBase(int m, int n, int leadingDim, Scalar* s)
00354       : helper(NScalarsPerElement, CppNScalarsPerElement, m, n, leadingDim, s) { }    // writable
00355             
00356     MatrixBase(int m, int n, const ELT& t) 
00357       : helper(NScalarsPerElement, CppNScalarsPerElement, m, n)
00358       { helper.fillWith(reinterpret_cast<const Scalar*>(&t)); }    
00359     MatrixBase(int m, int n, bool lockNrow, bool lockNcol, const ELT& t) 
00360       : helper(NScalarsPerElement, CppNScalarsPerElement, m, n, lockNrow, lockNcol)
00361       { helper.fillWith(reinterpret_cast<const Scalar*>(&t)); }
00362         
00363     MatrixBase(int m, int n, const ELT* p) 
00364       : helper(NScalarsPerElement, CppNScalarsPerElement, m, n)
00365       { helper.copyInByRows(reinterpret_cast<const Scalar*>(p)); }
00366     MatrixBase(int m, int n,  bool lockNrow, bool lockNcol, const ELT* p) 
00367       : helper(NScalarsPerElement, CppNScalarsPerElement, m, n, lockNrow, lockNcol)
00368       { helper.copyInByRows(reinterpret_cast<const Scalar*>(p)); }
00369         
00370     // Create a new MatrixBase from an existing helper. Both shallow and deep copies are possible.
00371     MatrixBase(MatrixHelper<Scalar>&       h, const typename MatrixHelper<Scalar>::ShallowCopy& s) : helper(h,s) { }
00372     MatrixBase(const MatrixHelper<Scalar>& h, const typename MatrixHelper<Scalar>::ShallowCopy& s) : helper(h,s) { }
00373     MatrixBase(const MatrixHelper<Scalar>& h, const typename MatrixHelper<Scalar>::DeepCopy& d)    : helper(h,d) { }
00374 
00375 
00376     MatrixBase& operator*=(const StdNumber& t)  { helper.scaleBy(t);              return *this; }
00377     MatrixBase& operator/=(const StdNumber& t)  { helper.scaleBy(StdNumber(1)/t); return *this; }
00378     MatrixBase& operator+=(const MatrixBase& r) { helper.addIn(r.helper);         return *this; }
00379     MatrixBase& operator-=(const MatrixBase& r) { helper.subIn(r.helper);         return *this; }  
00380 
00381     template <class EE> MatrixBase(const MatrixBase<EE>& b)
00382       : helper(b.helper, typename MatrixHelper<Scalar>::DeepCopy()) { }
00383 
00384     template <class EE> MatrixBase& operator=(const MatrixBase<EE>& b) 
00385       { helper = b.helper; return *this; }
00386     template <class EE> MatrixBase& operator+=(const MatrixBase<EE>& b) 
00387       { helper.addIn(b.helper); return *this; }
00388     template <class EE> MatrixBase& operator-=(const MatrixBase<EE>& b) 
00389       { helper.subIn(b.helper); return *this; }
00390 
00391     // Matrix assignment to an element sets only the *diagonal* elements to
00392     // the indicated value; everything else is set to zero. This is particularly
00393     // useful for setting a Matrix to zero or to the identity; for other values
00394     // it create as Matrix which acts like the scalar. That is, if the scalar
00395     // is s and we do M=s, then multiplying another Matrix B by the resulting 
00396     // diagonal matrix M gives the same result as multiplying B by s. That is
00397     // (M=s)*B == s*B.
00398     //
00399     // NOTE: this must be overridden for Vector and RowVector since then scalar
00400     // assignment is defined to copy the scalar to every element.
00401     MatrixBase& operator=(const ELT& t) { 
00402         setToZero(); updDiag().setTo(t); 
00403         return *this;
00404     }
00405 
00411     template <class S> inline MatrixBase&
00412     scalarAssign(const S& s) {
00413         setToZero(); updDiag().setTo(s);
00414         return *this;
00415     }
00416 
00420     template <class S> inline MatrixBase&
00421     scalarAddInPlace(const S& s) {
00422         updDiag().elementwiseAddScalarInPlace(s);
00423     }
00424 
00425 
00429     template <class S> inline MatrixBase&
00430     scalarSubtractInPlace(const S& s) {
00431         updDiag().elementwiseSubtractScalarInPlace(s);
00432     }
00433 
00436     template <class S> inline MatrixBase&
00437     scalarSubtractFromLeftInPlace(const S& s) {
00438         negateInPlace();
00439         updDiag().elementwiseAddScalarInPlace(s); // yes, add
00440     }
00441 
00448     template <class S> inline MatrixBase&
00449     scalarMultiplyInPlace(const S&);
00450 
00454     template <class S> inline MatrixBase&
00455     scalarMultiplyFromLeftInPlace(const S&);
00456 
00463     template <class S> inline MatrixBase&
00464     scalarDivideInPlace(const S&);
00465 
00469     template <class S> inline MatrixBase&
00470     scalarDivideFromLeftInPlace(const S&);
00471 
00472 
00473     // M = diag(r) * M; r must have nrow() elements.
00474     // That is, M[i] *= r[i].
00475     template <class EE> inline MatrixBase& 
00476     rowScaleInPlace(const VectorBase<EE>&);
00477 
00478     // Return type is a new matrix which will have the same dimensions as 'this' but
00479     // will have element types appropriate for the elementwise multiply being performed.
00480     template <class EE> inline void 
00481     rowScale(const VectorBase<EE>& r, typename EltResult<EE>::Mul& out) const;
00482 
00483     template <class EE> inline typename EltResult<EE>::Mul 
00484     rowScale(const VectorBase<EE>& r) const {
00485         typename EltResult<EE>::Mul out(nrow(), ncol()); rowScale(r,out); return out;
00486     }
00487 
00488     // M = M * diag(c); c must have ncol() elements
00489     // That is, M(j) *= c[j]
00490     template <class EE> inline MatrixBase& 
00491     colScaleInPlace(const VectorBase<EE>&);
00492 
00493     template <class EE> inline void 
00494     colScale(const VectorBase<EE>& c, typename EltResult<EE>::Mul& out) const;
00495 
00496     template <class EE> inline typename EltResult<EE>::Mul
00497     colScale(const VectorBase<EE>& c) const {
00498         typename EltResult<EE>::Mul out(nrow(), ncol()); colScale(c,out); return out;
00499     }
00500 
00501     // Having a combined row & column scaling operator means we can go through the matrix
00502     // memory once instead of twice.
00503 
00504     // M = diag(r) * M * diag(c); r must have nrow() elements;  must have ncol() elements
00505     // That is, M(i,j) *= r[i]*c[j]
00506     template <class ER, class EC> inline MatrixBase& 
00507     rowAndColScaleInPlace(const VectorBase<ER>& r, const VectorBase<EC>& c);
00508 
00509     template <class ER, class EC> inline void 
00510     rowAndColScale(const VectorBase<ER>& r, const VectorBase<EC>& c, 
00511                    typename EltResult<typename VectorBase<ER>::template EltResult<EC>::Mul>::Mul& out) const;
00512 
00513     template <class ER, class EC> inline typename EltResult<typename VectorBase<ER>::template EltResult<EC>::Mul>::Mul
00514     rowAndColScale(const VectorBase<ER>& r, const VectorBase<EC>& c) const {
00515         typename EltResult<typename VectorBase<ER>::template EltResult<EC>::Mul>::Mul 
00516             out(nrow(), ncol()); 
00517         rowAndColScale(r,c,out); return out;
00518     }
00519 
00527     template <class S> inline MatrixBase&
00528     elementwiseAssign(const S& s);
00529 
00531     MatrixBase& elementwiseInvertInPlace();
00532 
00533     void elementwiseInvert(MatrixBase<typename CNT<E>::TInvert>& out) const;
00534 
00535     MatrixBase<typename CNT<E>::TInvert> elementwiseInvert() const {
00536         MatrixBase<typename CNT<E>::TInvert> out(nrow(), ncol());
00537         elementwiseInvert(out);
00538         return out;
00539     }
00540 
00548     template <class S> inline MatrixBase&
00549     elementwiseAddScalarInPlace(const S& s);
00550 
00551     template <class S> inline void
00552     elementwiseAddScalar(const S& s, typename EltResult<S>::Add&) const;
00553 
00554     template <class S> inline typename EltResult<S>::Add
00555     elementwiseAddScalar(const S& s) const {
00556         typename EltResult<S>::Add out(nrow(), ncol());
00557         elementwiseAddScalar(s,out);
00558         return out;
00559     }
00560 
00568     template <class S> inline MatrixBase&
00569     elementwiseSubtractScalarInPlace(const S& s);
00570 
00571     template <class S> inline void
00572     elementwiseSubtractScalar(const S& s, typename EltResult<S>::Sub&) const;
00573 
00574     template <class S> inline typename EltResult<S>::Sub
00575     elementwiseSubtractScalar(const S& s) const {
00576         typename EltResult<S>::Sub out(nrow(), ncol());
00577         elementwiseSubtractScalar(s,out);
00578         return out;
00579     }
00580 
00589     template <class S> inline MatrixBase&
00590     elementwiseSubtractFromScalarInPlace(const S& s);
00591 
00592     template <class S> inline void
00593     elementwiseSubtractFromScalar(
00594         const S&, 
00595         typename MatrixBase<S>::template EltResult<E>::Sub&) const;
00596 
00597     template <class S> inline typename MatrixBase<S>::template EltResult<E>::Sub
00598     elementwiseSubtractFromScalar(const S& s) const {
00599         typename MatrixBase<S>::template EltResult<E>::Sub out(nrow(), ncol());
00600         elementwiseSubtractFromScalar<S>(s,out);
00601         return out;
00602     }
00603 
00604     // M(i,j) *= R(i,j); R must have same dimensions as this
00605     template <class EE> inline MatrixBase& 
00606     elementwiseMultiplyInPlace(const MatrixBase<EE>&);
00607 
00608     template <class EE> inline void 
00609     elementwiseMultiply(const MatrixBase<EE>&, typename EltResult<EE>::Mul&) const;
00610 
00611     template <class EE> inline typename EltResult<EE>::Mul 
00612     elementwiseMultiply(const MatrixBase<EE>& m) const {
00613         typename EltResult<EE>::Mul out(nrow(), ncol()); 
00614         elementwiseMultiply<EE>(m,out); 
00615         return out;
00616     }
00617 
00618     // M(i,j) = R(i,j) * M(i,j); R must have same dimensions as this
00619     template <class EE> inline MatrixBase& 
00620     elementwiseMultiplyFromLeftInPlace(const MatrixBase<EE>&);
00621 
00622     template <class EE> inline void 
00623     elementwiseMultiplyFromLeft(
00624         const MatrixBase<EE>&, 
00625         typename MatrixBase<EE>::template EltResult<E>::Mul&) const;
00626 
00627     template <class EE> inline typename MatrixBase<EE>::template EltResult<E>::Mul 
00628     elementwiseMultiplyFromLeft(const MatrixBase<EE>& m) const {
00629         typename EltResult<EE>::Mul out(nrow(), ncol()); 
00630         elementwiseMultiplyFromLeft<EE>(m,out); 
00631         return out;
00632     }
00633 
00634     // M(i,j) /= R(i,j); R must have same dimensions as this
00635     template <class EE> inline MatrixBase& 
00636     elementwiseDivideInPlace(const MatrixBase<EE>&);
00637 
00638     template <class EE> inline void 
00639     elementwiseDivide(const MatrixBase<EE>&, typename EltResult<EE>::Dvd&) const;
00640 
00641     template <class EE> inline typename EltResult<EE>::Dvd 
00642     elementwiseDivide(const MatrixBase<EE>& m) const {
00643         typename EltResult<EE>::Dvd out(nrow(), ncol()); 
00644         elementwiseDivide<EE>(m,out); 
00645         return out;
00646     }
00647 
00648     // M(i,j) = R(i,j) / M(i,j); R must have same dimensions as this
00649     template <class EE> inline MatrixBase& 
00650     elementwiseDivideFromLeftInPlace(const MatrixBase<EE>&);
00651 
00652     template <class EE> inline void 
00653     elementwiseDivideFromLeft(
00654         const MatrixBase<EE>&,
00655         typename MatrixBase<EE>::template EltResult<E>::Dvd&) const;
00656 
00657     template <class EE> inline typename MatrixBase<EE>::template EltResult<EE>::Dvd 
00658     elementwiseDivideFromLeft(const MatrixBase<EE>& m) const {
00659         typename MatrixBase<EE>::template EltResult<E>::Dvd out(nrow(), ncol()); 
00660         elementwiseDivideFromLeft<EE>(m,out); 
00661         return out;
00662     }
00663 
00664     // fill every element in current allocation with given element (or NaN or 0)
00665     MatrixBase& setTo(const ELT& t) {helper.fillWith(reinterpret_cast<const Scalar*>(&t)); return *this;}
00666     MatrixBase& setToNaN() {helper.fillWithScalar(CNT<StdNumber>::getNaN()); return *this;}
00667     MatrixBase& setToZero() {helper.fillWithScalar(StdNumber(0)); return *this;}
00668  
00669     // View creating operators. TODO: these should be DeadMatrixViews  
00670     inline RowVectorView_<ELT> row(int i) const;   // select a row
00671     inline RowVectorView_<ELT> updRow(int i);
00672     inline VectorView_<ELT>    col(int j) const;   // select a column
00673     inline VectorView_<ELT>    updCol(int j);
00674 
00675     RowVectorView_<ELT> operator[](int i) const {return row(i);}
00676     RowVectorView_<ELT> operator[](int i)       {return updRow(i);}
00677     VectorView_<ELT>    operator()(int j) const {return col(j);}
00678     VectorView_<ELT>    operator()(int j)       {return updCol(j);}
00679      
00680     // Select a block.
00681     inline MatrixView_<ELT> block(int i, int j, int m, int n) const;
00682     inline MatrixView_<ELT> updBlock(int i, int j, int m, int n);
00683 
00684     MatrixView_<ELT> operator()(int i, int j, int m, int n) const
00685       { return block(i,j,m,n); }
00686     MatrixView_<ELT> operator()(int i, int j, int m, int n)
00687       { return updBlock(i,j,m,n); }
00688  
00689     // Hermitian transpose.
00690     inline MatrixView_<EHerm> transpose() const;
00691     inline MatrixView_<EHerm> updTranspose();
00692 
00693     MatrixView_<EHerm> operator~() const {return transpose();}
00694     MatrixView_<EHerm> operator~()       {return updTranspose();}
00695 
00696     // Select matrix diagonal (of largest leading square if rectangular).
00697     inline VectorView_<ELT> diag() const;
00698     inline VectorView_<ELT> updDiag();
00699 
00700     // Create a view of the real or imaginary elements. TODO
00701     //inline MatrixView_<EReal> real() const;
00702     //inline MatrixView_<EReal> updReal();
00703     //inline MatrixView_<EImag> imag() const;
00704     //inline MatrixView_<EImag> updImag();
00705 
00706     // Overload "real" and "imag" for both read and write as a nod to convention. TODO
00707     //MatrixView_<EReal> real() {return updReal();}
00708     //MatrixView_<EReal> imag() {return updImag();}
00709 
00710     // TODO: this routine seems ill-advised but I need it for the IVM port at the moment
00711     TInvert invert() const {  // return a newly-allocated inverse; dump negator 
00712         TInvert m(*this);
00713         m.helper.invertInPlace();
00714         return m;   // TODO - bad: makes an extra copy
00715     }
00716 
00717     // Matlab-compatible debug output.
00718     void dump(const char* msg=0) const {
00719         helper.dump(msg);
00720     }
00721 
00722 
00723     // This routine is useful for implementing friendlier Matrix expressions and operators.
00724     // It maps closely to the Level-3 BLAS family of pxxmm() routines like sgemm(). The
00725     // operation performed assumes that "this" is the result, and that "this" has 
00726     // already been sized correctly to receive the result. We'll compute
00727     //     this = beta*this + alpha*A*B
00728     // If beta is 0 then "this" can be uninitialized. If alpha is 0 we promise not
00729     // to look at A or B. The routine can work efficiently even if A and/or B are transposed
00730     // by their views, so an expression like
00731     //        C += s * ~A * ~B
00732     // can be performed with the single equivalent call
00733     //        C.matmul(1., s, Tr(A), Tr(B))
00734     // where Tr(A) indicates a transposed view of the original A.
00735     // The ultimate efficiency of this call depends on the data layout and views which
00736     // are used for the three matrices.
00737     // NOTE: neither A nor B can be the same matrix as 'this', nor views of the same data
00738     // which would expose elements of 'this' that will be modified by this operation.
00739     template <class ELT_A, class ELT_B>
00740     MatrixBase& matmul(const StdNumber& beta,   // applied to 'this'
00741                        const StdNumber& alpha, const MatrixBase<ELT_A>& A, const MatrixBase<ELT_B>& B)
00742     {
00743         helper.matmul(beta,alpha,A.helper,B.helper);
00744         return *this;
00745     }
00746 
00747     // Element selection    
00748     const ELT& getElt(int i, int j) const { return *reinterpret_cast<const ELT*>(helper.getElt(i,j)); }
00749     ELT&       updElt(int i, int j)       { return *reinterpret_cast<      ELT*>(helper.updElt(i,j)); }
00750 
00751     const ELT& operator()(int i, int j) const {return getElt(i,j);}
00752     ELT&       operator()(int i, int j)       {return updElt(i,j);}
00753 
00754     // This is the scalar Frobenius norm, and its square. Note: if this is a Matrix then the Frobenius
00755     // norm is NOT the same as the 2-norm, although they are equivalent for Vectors.
00756     ScalarSq normSqr() const { return scalarNormSqr(); }
00757     ScalarSq norm()    const { return std::sqrt(scalarNormSqr()); } // TODO -- not good; unnecessary overflow
00758 
00759     // We only allow RMS norm if the elements are scalars. If there are no elements in this Matrix,
00760     // we'll define its RMS norm to be 0, although NaN might be a better choice.
00761     ScalarSq normRMS() const {
00762         if (!CNT<ELT>::IsScalar)
00763             SimTK_THROW1(Exception::Cant, "normRMS() only defined for scalar elements");
00764         const long nelt = (long)nrow()*(long)ncol();
00765         if (nelt == 0)
00766             return ScalarSq(0);
00767         return std::sqrt(scalarNormSqr()/nelt);
00768     }
00769 
00770 
00771     RowVectorBase<ELT> sum() const {
00772         const int cols = ncol();
00773         RowVectorBase<ELT> row(cols);
00774         for (int i = 0; i < cols; ++i)
00775             helper.colSum(i, reinterpret_cast<Scalar*>(&row[i]));
00776         return row;
00777     }
00778 
00779     //TODO: make unary +/- return a self-view so they won't reallocate?
00780     const MatrixBase& operator+() const {return *this; }
00781     const TNeg&       negate()    const {return *reinterpret_cast<const TNeg*>(this); }
00782     TNeg&             updNegate()       {return *reinterpret_cast<TNeg*>(this); }
00783 
00784     const TNeg&       operator-() const {return negate();}
00785     TNeg&             operator-()       {return updNegate();}
00786 
00787     MatrixBase& negateInPlace() {(*this) *= EPrecision(-1);}
00788  
00789     MatrixBase& resize(int m, int n)     { helper.resize(m,n); return *this; }
00790     MatrixBase& resizeKeep(int m, int n) { helper.resizeKeep(m,n); return *this; }
00791 
00792     // These prevent shape changes in a Matrix that would otherwise allow it. No harm if they
00793     // are called on a Matrix that is locked already; they always succeed.
00794     void lockNRows() {helper.lockNRows();}
00795     void lockNCols() {helper.lockNCols();}
00796     void lockShape() {helper.lockShape();}
00797 
00798     // These allow shape changes again for a Matrix which was constructed to allow them
00799     // but had them locked with the above routines. No harm if these are called on a Matrix
00800     // that is already unlocked, but it is not allowed to call them on a Matrix which
00801     // *never* allowed resizing. An exception will be thrown in that case.
00802     void unlockNRows() {helper.unlockNRows();}
00803     void unlockNCols() {helper.unlockNCols();}
00804     void unlockShape() {helper.unlockShape();}
00805     
00806     // An assortment of handy conversions
00807     const MatrixView_<ELT>& getAsMatrixView() const { return *reinterpret_cast<const MatrixView_<ELT>*>(this); }
00808     MatrixView_<ELT>&       updAsMatrixView()       { return *reinterpret_cast<      MatrixView_<ELT>*>(this); } 
00809     const Matrix_<ELT>&     getAsMatrix()     const { return *reinterpret_cast<const Matrix_<ELT>*>(this); }
00810     Matrix_<ELT>&           updAsMatrix()           { return *reinterpret_cast<      Matrix_<ELT>*>(this); }
00811          
00812     const VectorView_<ELT>& getAsVectorView() const 
00813       { assert(ncol()==1); return *reinterpret_cast<const VectorView_<ELT>*>(this); }
00814     VectorView_<ELT>&       updAsVectorView()       
00815       { assert(ncol()==1); return *reinterpret_cast<      VectorView_<ELT>*>(this); } 
00816     const Vector_<ELT>&     getAsVector()     const 
00817       { assert(ncol()==1); return *reinterpret_cast<const Vector_<ELT>*>(this); }
00818     Vector_<ELT>&           updAsVector()           
00819       { assert(ncol()==1); return *reinterpret_cast<      Vector_<ELT>*>(this); }
00820     const VectorBase<ELT>& getAsVectorBase() const 
00821       { assert(ncol()==1); return *reinterpret_cast<const VectorBase<ELT>*>(this); }
00822     VectorBase<ELT>&       updAsVectorBase()       
00823       { assert(ncol()==1); return *reinterpret_cast<      VectorBase<ELT>*>(this); } 
00824                 
00825     const RowVectorView_<ELT>& getAsRowVectorView() const 
00826       { assert(nrow()==1); return *reinterpret_cast<const RowVectorView_<ELT>*>(this); }
00827     RowVectorView_<ELT>&       updAsRowVectorView()       
00828       { assert(nrow()==1); return *reinterpret_cast<      RowVectorView_<ELT>*>(this); } 
00829     const RowVector_<ELT>&     getAsRowVector()     const 
00830       { assert(nrow()==1); return *reinterpret_cast<const RowVector_<ELT>*>(this); }
00831     RowVector_<ELT>&           updAsRowVector()           
00832       { assert(nrow()==1); return *reinterpret_cast<      RowVector_<ELT>*>(this); }        
00833     const RowVectorBase<ELT>& getAsRowVectorBase() const 
00834       { assert(nrow()==1); return *reinterpret_cast<const RowVectorBase<ELT>*>(this); }
00835     RowVectorBase<ELT>&       updAsRowVectorBase()       
00836       { assert(nrow()==1); return *reinterpret_cast<      RowVectorBase<ELT>*>(this); } 
00837 
00838     // Access to raw data. We have to return the raw data
00839     // pointer as pointer-to-scalar because we may pack the elements tighter
00840     // than a C++ array would.
00841 
00842     // This is the number of consecutive scalars used to represent one
00843     // element of type ELT. This may be fewer than C++ would use for the
00844     // element, since it may introduce some padding.
00845     int getNScalarsPerElement()  const {return NScalarsPerElement;}
00846 
00847     // This is like sizeof(ELT), but returning the number of bytes we use
00848     // to store the element which may be fewer than what C++ would use.
00849     int getPackedSizeofElement() const {return NScalarsPerElement*sizeof(Scalar);}
00850 
00851     bool hasContiguousData() const {return helper.hasContiguousData();}
00852     long getContiguousScalarDataLength() const {
00853         return helper.getContiguousDataLength();
00854     }
00855     const Scalar* getContiguousScalarData() const {
00856         return helper.getContiguousData();
00857     }
00858     Scalar* updContiguousScalarData() {
00859         return helper.updContiguousData();
00860     }
00861     void replaceContiguousScalarData(Scalar* newData, long length, bool takeOwnership) {
00862         helper.replaceContiguousData(newData,length,takeOwnership);
00863     }
00864     void replaceContiguousScalarData(const Scalar* newData, long length) {
00865         helper.replaceContiguousData(newData,length);
00866     }
00867     void swapOwnedContiguousScalarData(Scalar* newData, int length, Scalar*& oldData) {
00868         helper.swapOwnedContiguousData(newData,length,oldData);
00869     }
00870 protected:
00871     MatrixHelper<Scalar> helper; // this is just one pointer
00872 
00873     template <class EE> friend class MatrixBase;
00874 };
00875 
00881 template <class ELT> class VectorBase : public MatrixBase<ELT> {
00882     typedef MatrixBase<ELT>                             Base;
00883     typedef typename CNT<ELT>::Scalar                   Scalar;
00884     typedef typename CNT<ELT>::Number                   Number;
00885     typedef typename CNT<ELT>::StdNumber                StdNumber;
00886     typedef VectorBase<ELT>                             T;
00887     typedef VectorBase<typename CNT<ELT>::TAbs>         TAbs;
00888     typedef VectorBase<typename CNT<ELT>::TNeg>         TNeg;
00889     typedef RowVectorView_<typename CNT<ELT>::THerm>    THerm;
00890 public:     
00891     VectorBase() : Base(0,1,false,true) { } // a 0x1 matrix locked at 1 column
00892 
00893     // Copy constructor is a deep copy (not appropriate for views!).
00894     VectorBase(const VectorBase& b) : Base(b) { }
00895     
00896     // Copy assignment is deep copy but behavior depends on type of lhs: if view, rhs
00897     // must match. If owner, we reallocate and copy rhs.
00898     VectorBase& operator=(const VectorBase& b) {
00899         Base::operator=(b); return *this;
00900     }
00901 
00902     // default destructor
00903 
00904 
00905     explicit VectorBase(int m, bool lockNrow=false)
00906       : Base(m,1,lockNrow,true) { }
00907         
00908     // Non-resizeable owner sharing pre-allocated raw data
00909     VectorBase(int m, int leadingDim, const Scalar* s)
00910       : Base(m,1,leadingDim,s) { }
00911     VectorBase(int m, int leadingDim, Scalar* s)
00912       : Base(m,1,leadingDim,s) { }
00913             
00914     VectorBase(int m, const ELT& t) : Base(m,1,t) { }  
00915     VectorBase(int m, bool lockNrow, const ELT& t) 
00916       : Base(m,1,lockNrow,true,t) { }
00917         
00918     VectorBase(int m, const ELT* p) : Base(m,1,p) { }
00919     VectorBase(int m, bool lockNrow, const ELT* p) 
00920       : Base(m, 1, lockNrow, true, p) { }
00921         
00922     // Create a new VectorBase from an existing helper. Both shallow and deep copies are possible.
00923     VectorBase(MatrixHelper<Scalar>&       h, const typename MatrixHelper<Scalar>::ShallowCopy& s) : Base(h,s) { }
00924     VectorBase(const MatrixHelper<Scalar>& h, const typename MatrixHelper<Scalar>::ShallowCopy& s) : Base(h,s) { }
00925     VectorBase(const MatrixHelper<Scalar>& h, const typename MatrixHelper<Scalar>::DeepCopy& d)    : Base(h,d) { }
00926 
00927     // This gives the resulting vector type when (v[i] op P) is applied to each element.
00928     // It will have element types which are the regular composite result of ELT op P.
00929     template <class P> struct EltResult { 
00930         typedef VectorBase<typename CNT<ELT>::template Result<P>::Mul> Mul;
00931         typedef VectorBase<typename CNT<ELT>::template Result<P>::Dvd> Dvd;
00932         typedef VectorBase<typename CNT<ELT>::template Result<P>::Add> Add;
00933         typedef VectorBase<typename CNT<ELT>::template Result<P>::Sub> Sub;
00934     };
00935 
00936     VectorBase& operator*=(const StdNumber& t)  { Base::operator*=(t); return *this; }
00937     VectorBase& operator/=(const StdNumber& t)  { Base::operator/=(t); return *this; }
00938     VectorBase& operator+=(const VectorBase& r) { Base::operator+=(r); return *this; }
00939     VectorBase& operator-=(const VectorBase& r) { Base::operator-=(r); return *this; }  
00940 
00941 
00942     template <class EE> VectorBase& operator=(const VectorBase<EE>& b) 
00943       { Base::operator=(b);  return *this; } 
00944     template <class EE> VectorBase& operator+=(const VectorBase<EE>& b) 
00945       { Base::operator+=(b); return *this; } 
00946     template <class EE> VectorBase& operator-=(const VectorBase<EE>& b) 
00947       { Base::operator-=(b); return *this; } 
00948 
00949 
00950     // Fill current allocation with copies of element. Note that this is not the 
00951     // same behavior as assignment for Matrices, where only the diagonal is set (and
00952     // everything else is set to zero.)
00953     VectorBase& operator=(const ELT& t) { Base::setTo(t); return *this; }  
00954 
00955     // There's only one column here so it's a bit wierd to use rowScale rather than
00956     // elementwiseMultiply, but there's nothing really wrong with it. Using colScale
00957     // would be really wacky since it is the same as a scalar multiply. We won't support
00958     // colScale here except through inheritance where it will not be much use.
00959     template <class EE> VectorBase& rowScaleInPlace(const VectorBase<EE>& v)
00960       { Base::template rowScaleInPlace<EE>(v); return *this; }
00961     template <class EE> inline void rowScale(const VectorBase<EE>& v, typename EltResult<EE>::Mul& out) const
00962       { Base::rowScale(v,out); }
00963     template <class EE> inline typename EltResult<EE>::Mul rowScale(const VectorBase<EE>& v) const
00964       { typename EltResult<EE>::Mul out(nrow()); Base::rowScale(v,out); return out; }
00965 
00967     VectorBase& elementwiseInvertInPlace() {
00968         Base::elementwiseInvertInPlace();
00969         return *this;
00970     }
00971 
00973     void elementwiseInvert(VectorBase<typename CNT<ELT>::TInvert>& out) const {
00974         Base::elementwiseInvert(out);
00975     }
00976 
00978     VectorBase<typename CNT<ELT>::TInvert> elementwiseInvert() const {
00979         VectorBase<typename CNT<ELT>::TInvert> out(nrow());
00980         Base::elementwiseInvert(out);
00981         return out;
00982     }
00983 
00984         // elementwise multiply
00985     template <class EE> VectorBase& elementwiseMultiplyInPlace(const VectorBase<EE>& r)
00986       { Base::template elementwiseMultiplyInPlace<EE>(r); return *this; }
00987     template <class EE> inline void elementwiseMultiply(const VectorBase<EE>& v, typename EltResult<EE>::Mul& out) const
00988       { Base::template elementwiseMultiply<EE>(v,out); }
00989     template <class EE> inline typename EltResult<EE>::Mul elementwiseMultiply(const VectorBase<EE>& v) const
00990       { typename EltResult<EE>::Mul out(nrow()); Base::template elementwiseMultiply<EE>(v,out); return out; }
00991 
00992         // elementwise multiply from left
00993     template <class EE> VectorBase& elementwiseMultiplyFromLeftInPlace(const VectorBase<EE>& r)
00994       { Base::template elementwiseMultiplyFromLeftInPlace<EE>(r); return *this; }
00995     template <class EE> inline void 
00996     elementwiseMultiplyFromLeft(
00997         const VectorBase<EE>& v, 
00998         typename VectorBase<EE>::template EltResult<ELT>::Mul& out) const
00999     { 
01000         Base::template elementwiseMultiplyFromLeft<EE>(v,out);
01001     }
01002     template <class EE> inline typename VectorBase<EE>::template EltResult<ELT>::Mul 
01003     elementwiseMultiplyFromLeft(const VectorBase<EE>& v) const
01004     { 
01005         typename VectorBase<EE>::template EltResult<ELT>::Mul out(nrow()); 
01006         Base::template elementwiseMultiplyFromLeft<EE>(v,out); 
01007         return out;
01008     }
01009 
01010         // elementwise divide
01011     template <class EE> VectorBase& elementwiseDivideInPlace(const VectorBase<EE>& r)
01012       { Base::template elementwiseDivideInPlace<EE>(r); return *this; }
01013     template <class EE> inline void elementwiseDivide(const VectorBase<EE>& v, typename EltResult<EE>::Dvd& out) const
01014       { Base::template elementwiseDivide<EE>(v,out); }
01015     template <class EE> inline typename EltResult<EE>::Dvd elementwiseDivide(const VectorBase<EE>& v) const
01016       { typename EltResult<EE>::Dvd out(nrow()); Base::template elementwiseDivide<EE>(v,out); return out; }
01017 
01018         // elementwise divide from left
01019     template <class EE> VectorBase& elementwiseDivideFromLeftInPlace(const VectorBase<EE>& r)
01020       { Base::template elementwiseDivideFromLeftInPlace<EE>(r); return *this; }
01021     template <class EE> inline void 
01022     elementwiseDivideFromLeft(
01023         const VectorBase<EE>& v, 
01024         typename VectorBase<EE>::template EltResult<ELT>::Dvd& out) const
01025     { 
01026         Base::template elementwiseDivideFromLeft<EE>(v,out);
01027     }
01028     template <class EE> inline typename VectorBase<EE>::template EltResult<ELT>::Dvd 
01029     elementwiseDivideFromLeft(const VectorBase<EE>& v) const
01030     { 
01031         typename VectorBase<EE>::template EltResult<ELT>::Dvd out(nrow()); 
01032         Base::template elementwiseDivideFromLeft<EE>(v,out); 
01033         return out;
01034     }
01035 
01036     // Implicit conversions are allowed to Vector or Matrix, but not to RowVector.   
01037     operator const Vector_<ELT>&()     const { return *reinterpret_cast<const Vector_<ELT>*>(this); }
01038     operator       Vector_<ELT>&()           { return *reinterpret_cast<      Vector_<ELT>*>(this); }
01039     operator const VectorView_<ELT>&() const { return *reinterpret_cast<const VectorView_<ELT>*>(this); }
01040     operator       VectorView_<ELT>&()       { return *reinterpret_cast<      VectorView_<ELT>*>(this); }
01041     
01042     operator const Matrix_<ELT>&()     const { return *reinterpret_cast<const Matrix_<ELT>*>(this); }
01043     operator       Matrix_<ELT>&()           { return *reinterpret_cast<      Matrix_<ELT>*>(this); } 
01044     operator const MatrixView_<ELT>&() const { return *reinterpret_cast<const MatrixView_<ELT>*>(this); }
01045     operator       MatrixView_<ELT>&()       { return *reinterpret_cast<      MatrixView_<ELT>*>(this); } 
01046 
01047 
01048     // Override MatrixBase size() for Vectors to return int instead of long.
01049     int size() const { 
01050         assert(Base::size() <= (long)std::numeric_limits<int>::max()); 
01051         assert(Base::ncol()==1);
01052         return (int)Base::size();
01053     }
01054     int nrow() const { assert(Base::ncol()==1); return Base::nrow(); }
01055     int ncol() const { assert(Base::ncol()==1); return Base::ncol(); }
01056 
01057     // Override MatrixBase operators to return the right shape
01058     TAbs abs() const {TAbs result; Base::abs(result); return result;}
01059     
01060     // Override MatrixBase indexing operators          
01061     const ELT& operator[](int i) const {return Base::operator()(i,0);}
01062     ELT&       operator[](int i)       {return Base::operator()(i,0);}
01063     const ELT& operator()(int i) const {return Base::operator()(i,0);}
01064     ELT&       operator()(int i)       {return Base::operator()(i,0);}
01065          
01066     // View creation      
01067     VectorView_<ELT> operator()(int i, int m) const {return Base::operator()(i,0,m,1).getAsVectorView();}
01068     VectorView_<ELT> operator()(int i, int m)       {return Base::operator()(i,0,m,1).updAsVectorView();}
01069  
01070     // Hermitian transpose.
01071     THerm transpose() const {return Base::transpose().getAsRowVectorView();}
01072     THerm updTranspose()    {return Base::updTranspose().updAsRowVectorView();}
01073 
01074     THerm operator~() const {return transpose();}
01075     THerm operator~()       {return updTranspose();}
01076 
01077     const VectorBase& operator+() const {return *this; }
01078 
01079     // Negation
01080 
01081     const TNeg& negate()    const {return *reinterpret_cast<const TNeg*>(this); }
01082     TNeg&       updNegate()       {return *reinterpret_cast<TNeg*>(this); }
01083 
01084     const TNeg& operator-() const {return negate();}
01085     TNeg&       operator-()       {return updNegate();}
01086 
01087     VectorBase& resize(int m)     {Base::resize(m,1); return *this;}
01088     VectorBase& resizeKeep(int m) {Base::resizeKeep(m,1); return *this;}
01089 
01090     ELT sum() const {ELT s; Base::helper.sum(reinterpret_cast<Scalar*>(&s)); return s; } // add all the elements        
01091     VectorIterator<ELT, VectorBase<ELT> > begin() {
01092         return VectorIterator<ELT, VectorBase<ELT> >(*this, 0);
01093     }
01094     VectorIterator<ELT, VectorBase<ELT> > end() {
01095         return VectorIterator<ELT, VectorBase<ELT> >(*this, size());
01096     }
01097 private:
01098     // NO DATA MEMBERS ALLOWED
01099 };
01100 
01101 
01107 template <class ELT> class RowVectorBase : public MatrixBase<ELT> {
01108     typedef MatrixBase<ELT>                             Base;
01109     typedef typename CNT<ELT>::Scalar                   Scalar;
01110     typedef typename CNT<ELT>::Number                   Number;
01111     typedef typename CNT<ELT>::StdNumber                StdNumber;
01112     typedef RowVectorBase<ELT>                          T;
01113     typedef RowVectorBase<typename CNT<ELT>::TAbs>      TAbs;
01114     typedef RowVectorBase<typename CNT<ELT>::TNeg>      TNeg;
01115     typedef VectorView_<typename CNT<ELT>::THerm>       THerm;
01116 public:     
01117     RowVectorBase() : Base(1,0,true,false) { } // a 1x0 matrix locked at 1 row
01118     
01119     // Copy constructor is a deep copy (not appropriate for views!).    
01120     RowVectorBase(const RowVectorBase& b) : Base(b) { }
01121 
01122     // Copy assignment is deep copy but behavior depends on type of lhs: if view, rhs
01123     // must match. If owner, we reallocate and copy rhs.
01124     RowVectorBase& operator=(const RowVectorBase& b) {
01125         Base::operator=(b); return *this;
01126     }
01127 
01128     // default destructor
01129 
01130 
01131     explicit RowVectorBase(int n, bool lockNcol=false)
01132       : Base(1,n,true,lockNcol) { }
01133         
01134     // Non-resizeable owner sharing pre-allocated raw data
01135     RowVectorBase(int n, int leadingDim, const Scalar* s)
01136       : Base(1,n,leadingDim,s) { }
01137     RowVectorBase(int n, int leadingDim, Scalar* s)
01138       : Base(1,n,leadingDim,s) { }
01139             
01140     RowVectorBase(int n, const ELT& t) : Base(1,n,t) { }  
01141     RowVectorBase(int n, bool lockNcol, const ELT& t) 
01142       : Base(1,n,true,lockNcol,t) { }
01143         
01144     RowVectorBase(int n, const ELT* p) : Base(1,n,p) { }
01145     RowVectorBase(int n, bool lockNcol, const ELT* p) 
01146       : Base(1, n, true, lockNcol, p) { }
01147         
01148     // Create a new RowVectorBase from an existing helper. Both shallow and deep copies are possible.
01149     RowVectorBase(MatrixHelper<Scalar>&       h, const typename MatrixHelper<Scalar>::ShallowCopy& s) : Base(h,s) { }
01150     RowVectorBase(const MatrixHelper<Scalar>& h, const typename MatrixHelper<Scalar>::ShallowCopy& s) : Base(h,s) { }
01151     RowVectorBase(const MatrixHelper<Scalar>& h, const typename MatrixHelper<Scalar>::DeepCopy& d)    : Base(h,d) { }
01152 
01153     // This gives the resulting rowvector type when (r(i) op P) is applied to each element.
01154     // It will have element types which are the regular composite result of ELT op P.
01155     template <class P> struct EltResult { 
01156         typedef RowVectorBase<typename CNT<ELT>::template Result<P>::Mul> Mul;
01157         typedef RowVectorBase<typename CNT<ELT>::template Result<P>::Dvd> Dvd;
01158         typedef RowVectorBase<typename CNT<ELT>::template Result<P>::Add> Add;
01159         typedef RowVectorBase<typename CNT<ELT>::template Result<P>::Sub> Sub;
01160     };
01161 
01162     RowVectorBase& operator*=(const StdNumber& t)     {Base::operator*=(t); return *this;}
01163     RowVectorBase& operator/=(const StdNumber& t)     {Base::operator/=(t); return *this;}
01164     RowVectorBase& operator+=(const RowVectorBase& r) {Base::operator+=(r); return *this;}
01165     RowVectorBase& operator-=(const RowVectorBase& r) {Base::operator-=(r); return *this;}  
01166 
01167     template <class EE> RowVectorBase& operator=(const RowVectorBase<EE>& b) 
01168       { Base::operator=(b);  return *this; } 
01169     template <class EE> RowVectorBase& operator+=(const RowVectorBase<EE>& b) 
01170       { Base::operator+=(b); return *this; } 
01171     template <class EE> RowVectorBase& operator-=(const RowVectorBase<EE>& b) 
01172       { Base::operator-=(b); return *this; } 
01173 
01174     // default destructor
01175  
01176     // Fill current allocation with copies of element. Note that this is not the 
01177     // same behavior as assignment for Matrices, where only the diagonal is set (and
01178     // everything else is set to zero.)
01179     RowVectorBase& operator=(const ELT& t) { Base::setTo(t); return *this; } 
01180 
01181     // There's only one row here so it's a bit wierd to use colScale rather than
01182     // elementwiseMultiply, but there's nothing really wrong with it. Using rowScale
01183     // would be really wacky since it is the same as a scalar multiply. We won't support
01184     // rowScale here except through inheritance where it will not be much use.
01185 
01186     template <class EE> RowVectorBase& colScaleInPlace(const VectorBase<EE>& v)
01187       { Base::template colScaleInPlace<EE>(v); return *this; }
01188     template <class EE> inline void colScale(const VectorBase<EE>& v, typename EltResult<EE>::Mul& out) const
01189       { return Base::template colScale<EE>(v,out); }
01190     template <class EE> inline typename EltResult<EE>::Mul colScale(const VectorBase<EE>& v) const
01191       { typename EltResult<EE>::Mul out(ncol()); Base::template colScale<EE>(v,out); return out; }
01192 
01193 
01194         // elementwise multiply
01195     template <class EE> RowVectorBase& elementwiseMultiplyInPlace(const RowVectorBase<EE>& r)
01196       { Base::template elementwiseMultiplyInPlace<EE>(r); return *this; }
01197     template <class EE> inline void elementwiseMultiply(const RowVectorBase<EE>& v, typename EltResult<EE>::Mul& out) const
01198       { Base::template elementwiseMultiply<EE>(v,out); }
01199     template <class EE> inline typename EltResult<EE>::Mul elementwiseMultiply(const RowVectorBase<EE>& v) const
01200       { typename EltResult<EE>::Mul out(nrow()); Base::template elementwiseMultiply<EE>(v,out); return out; }
01201 
01202         // elementwise multiply from left
01203     template <class EE> RowVectorBase& elementwiseMultiplyFromLeftInPlace(const RowVectorBase<EE>& r)
01204       { Base::template elementwiseMultiplyFromLeftInPlace<EE>(r); return *this; }
01205     template <class EE> inline void 
01206     elementwiseMultiplyFromLeft(
01207         const RowVectorBase<EE>& v, 
01208         typename RowVectorBase<EE>::template EltResult<ELT>::Mul& out) const
01209     { 
01210         Base::template elementwiseMultiplyFromLeft<EE>(v,out);
01211     }
01212     template <class EE> inline typename RowVectorBase<EE>::template EltResult<ELT>::Mul 
01213     elementwiseMultiplyFromLeft(const RowVectorBase<EE>& v) const
01214     { 
01215         typename RowVectorBase<EE>::template EltResult<ELT>::Mul out(nrow()); 
01216         Base::template elementwiseMultiplyFromLeft<EE>(v,out); 
01217         return out;
01218     }
01219 
01220         // elementwise divide
01221     template <class EE> RowVectorBase& elementwiseDivideInPlace(const RowVectorBase<EE>& r)
01222       { Base::template elementwiseDivideInPlace<EE>(r); return *this; }
01223     template <class EE> inline void elementwiseDivide(const RowVectorBase<EE>& v, typename EltResult<EE>::Dvd& out) const
01224       { Base::template elementwiseDivide<EE>(v,out); }
01225     template <class EE> inline typename EltResult<EE>::Dvd elementwiseDivide(const RowVectorBase<EE>& v) const
01226       { typename EltResult<EE>::Dvd out(nrow()); Base::template elementwiseDivide<EE>(v,out); return out; }
01227 
01228         // elementwise divide from left
01229     template <class EE> RowVectorBase& elementwiseDivideFromLeftInPlace(const RowVectorBase<EE>& r)
01230       { Base::template elementwiseDivideFromLeftInPlace<EE>(r); return *this; }
01231     template <class EE> inline void 
01232     elementwiseDivideFromLeft(
01233         const RowVectorBase<EE>& v, 
01234         typename RowVectorBase<EE>::template EltResult<ELT>::Dvd& out) const
01235     { 
01236         Base::template elementwiseDivideFromLeft<EE>(v,out);
01237     }
01238     template <class EE> inline typename RowVectorBase<EE>::template EltResult<ELT>::Dvd 
01239     elementwiseDivideFromLeft(const RowVectorBase<EE>& v) const
01240     { 
01241         typename RowVectorBase<EE>::template EltResult<ELT>::Dvd out(nrow()); 
01242         Base::template elementwiseDivideFromLeft<EE>(v,out); 
01243         return out;
01244     }
01245 
01246     // Implicit conversions are allowed to RowVector or Matrix, but not to Vector.   
01247     operator const RowVector_<ELT>&()     const {return *reinterpret_cast<const RowVector_<ELT>*>(this);}
01248     operator       RowVector_<ELT>&()           {return *reinterpret_cast<      RowVector_<ELT>*>(this);}
01249     operator const RowVectorView_<ELT>&() const {return *reinterpret_cast<const RowVectorView_<ELT>*>(this);}
01250     operator       RowVectorView_<ELT>&()       {return *reinterpret_cast<      RowVectorView_<ELT>*>(this);}
01251     
01252     operator const Matrix_<ELT>&()     const {return *reinterpret_cast<const Matrix_<ELT>*>(this);}
01253     operator       Matrix_<ELT>&()           {return *reinterpret_cast<      Matrix_<ELT>*>(this);} 
01254     operator const MatrixView_<ELT>&() const {return *reinterpret_cast<const MatrixView_<ELT>*>(this);}
01255     operator       MatrixView_<ELT>&()       {return *reinterpret_cast<      MatrixView_<ELT>*>(this);} 
01256     
01257 
01258     // Override MatrixBase size() for RowVectors to return int instead of long.
01259     int size() const { 
01260         assert(Base::size() <= (long)std::numeric_limits<int>::max()); 
01261         assert(Base::nrow()==1);
01262         return (int)Base::size();
01263     }
01264     int nrow() const { assert(Base::nrow()==1); return Base::nrow(); }
01265     int ncol() const { assert(Base::nrow()==1); return Base::ncol(); }
01266 
01267     // Override MatrixBase operators to return the right shape
01268     TAbs abs() const {
01269         TAbs result; Base::abs(result); return result;
01270     }
01271 
01272     // Override MatrixBase indexing operators          
01273     const ELT& operator[](int j) const {return Base::operator()(0,j);}
01274     ELT&       operator[](int j)       {return Base::operator()(0,j);}
01275     const ELT& operator()(int j) const {return Base::operator()(0,j);}
01276     ELT&       operator()(int j)       {return Base::operator()(0,j);}
01277          
01278     // View creation      
01279     RowVectorView_<ELT> operator()(int j, int n) const {return Base::operator()(0,j,1,n).getAsRowVectorView();}
01280     RowVectorView_<ELT> operator()(int j, int n)       {return Base::operator()(0,j,1,n).updAsRowVectorView();}
01281  
01282     // Hermitian transpose.
01283     THerm transpose() const {return Base::transpose().getAsVectorView();}
01284     THerm updTranspose()    {return Base::updTranspose().updAsVectorView();}
01285 
01286     THerm operator~() const {return transpose();}
01287     THerm operator~()       {return updTranspose();}
01288 
01289     const RowVectorBase& operator+() const {return *this; }
01290 
01291     // Negation
01292 
01293     const TNeg& negate()    const {return *reinterpret_cast<const TNeg*>(this); }
01294     TNeg&       updNegate()       {return *reinterpret_cast<TNeg*>(this); }
01295 
01296     const TNeg& operator-() const {return negate();}
01297     TNeg&       operator-()       {return updNegate();}
01298 
01299     RowVectorBase& resize(int n)     {Base::resize(1,n); return *this;}
01300     RowVectorBase& resizeKeep(int n) {Base::resizeKeep(1,n); return *this;}
01301 
01302     ELT sum() const {ELT s; Base::helper.sum(reinterpret_cast<Scalar*>(&s)); return s; } // add all the elements        
01303     VectorIterator<ELT, RowVectorBase<ELT> > begin() {
01304         return VectorIterator<ELT, RowVectorBase<ELT> >(*this, 0);
01305     }
01306     VectorIterator<ELT, RowVectorBase<ELT> > end() {
01307         return VectorIterator<ELT, RowVectorBase<ELT> >(*this, size());
01308     }
01309 private:
01310     // NO DATA MEMBERS ALLOWED
01311 };
01312 
01316 template <class ELT> class TmpMatrixView_ : public MatrixBase<ELT> {
01317     typedef MatrixBase<ELT> Base;
01318     typedef typename MatrixBase<ELT>::Scalar Scalar;
01319 public:
01320     TmpMatrixView_() : Base() { }
01321     TmpMatrixView_(int m, int n) : Base(m,n) { }
01322     explicit TmpMatrixView_(const MatrixHelper<Scalar>& h) : Base(h) { }
01323     
01324     operator const Matrix_<ELT>&() const { return *reinterpret_cast<const Matrix_<ELT>*>(this); }
01325     operator Matrix_<ELT>&()             { return *reinterpret_cast<Matrix_<ELT>*>(this); }
01326     
01327     TmpMatrixView_* clone() const { return new TmpMatrixView_(*this); } 
01328 private:
01329     // NO DATA MEMBERS ALLOWED
01330 };
01331 
01339 template <class ELT> class MatrixView_ : public MatrixBase<ELT> {
01340     typedef MatrixBase<ELT>                 Base;
01341     typedef typename CNT<ELT>::Scalar       S;
01342     typedef typename CNT<ELT>::StdNumber    StdNumber;
01343 public:
01344     // Default construction is suppressed.
01345     // Uses default destructor.
01346 
01347     // Copy constructor is shallow. CAUTION: despite const argument, this preserves writability
01348     // if it was present in the source. This is necessary to allow temporary views to be
01349     // created and used as lvalues.
01350     MatrixView_(const MatrixView_& m) 
01351       : Base(const_cast<MatrixHelper<S>&>(m.helper), typename MatrixHelper<S>::ShallowCopy()) { }
01352 
01353     // Copy assignment is deep but not reallocating.
01354     MatrixView_& operator=(const MatrixView_& m) {
01355         Base::operator=(m); return *this;
01356     }
01357 
01358     // Ask for shallow copy    
01359     MatrixView_(const MatrixHelper<S>& h) : Base(h, typename MatrixHelper<S>::ShallowCopy()) { }
01360     MatrixView_(MatrixHelper<S>&       h) : Base(h, typename MatrixHelper<S>::ShallowCopy()) { }
01361 
01362     MatrixView_& operator=(const Matrix_<ELT>& v)     { Base::operator=(v); return *this; }
01363     MatrixView_& operator=(const ELT& e)              { Base::operator=(e); return *this; }
01364 
01365     template <class EE> MatrixView_& operator=(const MatrixBase<EE>& m)
01366       { Base::operator=(m); return *this; }
01367     template <class EE> MatrixView_& operator+=(const MatrixBase<EE>& m)
01368       { Base::operator+=(m); return *this; }
01369     template <class EE> MatrixView_& operator-=(const MatrixBase<EE>& m)
01370       { Base::operator-=(m); return *this; }
01371 
01372     MatrixView_& operator*=(const StdNumber& t) { Base::operator*=(t); return *this; }
01373     MatrixView_& operator/=(const StdNumber& t) { Base::operator/=(t); return *this; }
01374     MatrixView_& operator+=(const ELT& r)       { this->updDiag() += r; return *this; }
01375     MatrixView_& operator-=(const ELT& r)       { this->updDiag() -= r; return *this; }  
01376 
01377     operator const Matrix_<ELT>&() const { return *reinterpret_cast<const Matrix_<ELT>*>(this); }
01378     operator Matrix_<ELT>&()             { return *reinterpret_cast<Matrix_<ELT>*>(this); }      
01379 
01380 private:
01381     // NO DATA MEMBERS ALLOWED
01382     MatrixView_(); // default constructor suppressed (what's it a view of?)
01383 };
01384 
01385 template <class ELT> inline MatrixView_<ELT> 
01386 MatrixBase<ELT>::block(int i, int j, int m, int n) const { 
01387     SimTK_INDEXCHECK(0,i,nrow()+1,"MatrixBase::block()");
01388     SimTK_INDEXCHECK(0,j,ncol()+1,"MatrixBase::block()");
01389     SimTK_SIZECHECK(i+m,nrow(),"MatrixBase::block()");
01390     SimTK_SIZECHECK(j+n,ncol(),"MatrixBase::block()");
01391 
01392     MatrixHelper<Scalar> h(helper,i,j,m,n);    
01393     return MatrixView_<ELT>(h); 
01394 }
01395     
01396 template <class ELT> inline MatrixView_<ELT>
01397 MatrixBase<ELT>::updBlock(int i, int j, int m, int n) { 
01398     SimTK_INDEXCHECK(0,i,nrow()+1,"MatrixBase::updBlock()");
01399     SimTK_INDEXCHECK(0,j,ncol()+1,"MatrixBase::updBlock()");
01400     SimTK_SIZECHECK(i+m,nrow(),"MatrixBase::updBlock()");
01401     SimTK_SIZECHECK(j+n,ncol(),"MatrixBase::updBlock()");
01402 
01403     MatrixHelper<Scalar> h(helper,i,j,m,n);        
01404     return MatrixView_<ELT>(h); 
01405 }
01406 
01407 template <class E> inline MatrixView_<typename CNT<E>::THerm>
01408 MatrixBase<E>::transpose() const { 
01409     MatrixHelper<typename CNT<Scalar>::THerm> 
01410         h(helper, typename MatrixHelper<typename CNT<Scalar>::THerm>::TransposeView());
01411     return MatrixView_<typename CNT<E>::THerm>(h); 
01412 }
01413     
01414 template <class E> inline MatrixView_<typename CNT<E>::THerm>
01415 MatrixBase<E>::updTranspose() {     
01416     MatrixHelper<typename CNT<Scalar>::THerm> 
01417         h(helper, typename MatrixHelper<typename CNT<Scalar>::THerm>::TransposeView());
01418     return MatrixView_<typename CNT<E>::THerm>(h); 
01419 }
01420 
01421 template <class E> inline VectorView_<E>
01422 MatrixBase<E>::diag() const { 
01423     MatrixHelper<Scalar> h(helper, typename MatrixHelper<Scalar>::DiagonalView());
01424     return VectorView_<E>(h); 
01425 }
01426     
01427 template <class E> inline VectorView_<E>
01428 MatrixBase<E>::updDiag() {     
01429     MatrixHelper<Scalar> h(helper, typename MatrixHelper<Scalar>::DiagonalView());
01430     return VectorView_<E>(h);
01431 }
01432 
01433 template <class ELT> inline VectorView_<ELT> 
01434 MatrixBase<ELT>::col(int j) const { 
01435     SimTK_INDEXCHECK(0,j,ncol(),"MatrixBase::col()");
01436 
01437     MatrixHelper<Scalar> h(helper,0,j,nrow(),1);    
01438     return VectorView_<ELT>(h); 
01439 }
01440     
01441 template <class ELT> inline VectorView_<ELT>
01442 MatrixBase<ELT>::updCol(int j) {
01443     SimTK_INDEXCHECK(0,j,ncol(),"MatrixBase::updCol()");
01444 
01445     MatrixHelper<Scalar> h(helper,0,j,nrow(),1);        
01446     return VectorView_<ELT>(h); 
01447 }
01448 
01449 template <class ELT> inline RowVectorView_<ELT> 
01450 MatrixBase<ELT>::row(int i) const { 
01451     SimTK_INDEXCHECK(0,i,nrow(),"MatrixBase::row()");
01452 
01453     MatrixHelper<Scalar> h(helper,i,0,1,ncol());    
01454     return RowVectorView_<ELT>(h); 
01455 }
01456     
01457 template <class ELT> inline RowVectorView_<ELT>
01458 MatrixBase<ELT>::updRow(int i) { 
01459     SimTK_INDEXCHECK(0,i,nrow(),"MatrixBase::updRow()");
01460 
01461     MatrixHelper<Scalar> h(helper,i,0,1,ncol());        
01462     return RowVectorView_<ELT>(h); 
01463 }
01464 
01465 // M = diag(v) * M; v must have nrow() elements.
01466 // That is, M[i] *= v[i].
01467 template <class ELT> template <class EE> inline MatrixBase<ELT>& 
01468 MatrixBase<ELT>::rowScaleInPlace(const VectorBase<EE>& v) {
01469     assert(v.nrow() == nrow());
01470     for (int i=0; i < nrow(); ++i)
01471         (*this)[i] *= v[i];
01472     return *this;
01473 }
01474 
01475 template <class ELT> template <class EE> inline void
01476 MatrixBase<ELT>::rowScale(const VectorBase<EE>& v, typename MatrixBase<ELT>::template EltResult<EE>::Mul& out) const {
01477     assert(v.nrow() == nrow());
01478     out.resize(nrow(), ncol());
01479     for (int j=0; j<ncol(); ++j)
01480         for (int i=0; i<nrow(); ++i)
01481            out(i,j) = (*this)(i,j) * v[i];
01482 }
01483 
01484 // M = M * diag(v); v must have ncol() elements
01485 // That is, M(i) *= v[i]
01486 template <class ELT> template <class EE>  inline MatrixBase<ELT>& 
01487 MatrixBase<ELT>::colScaleInPlace(const VectorBase<EE>& v) {
01488     assert(v.nrow() == ncol());
01489     for (int j=0; j < ncol(); ++j)
01490         (*this)(j) *= v[j];
01491     return *this;
01492 }
01493 
01494 template <class ELT> template <class EE> inline void
01495 MatrixBase<ELT>::colScale(const VectorBase<EE>& v, typename MatrixBase<ELT>::template EltResult<EE>::Mul& out) const {
01496     assert(v.nrow() == ncol());
01497     out.resize(nrow(), ncol());
01498     for (int j=0; j<ncol(); ++j)
01499         for (int i=0; i<nrow(); ++i)
01500            out(i,j) = (*this)(i,j) * v[j];
01501 }
01502 
01503 
01504 // M(i,j) *= r[i]*c[j]; r must have nrow() elements; c must have ncol() elements
01505 template <class ELT> template <class ER, class EC> inline MatrixBase<ELT>& 
01506 MatrixBase<ELT>::rowAndColScaleInPlace(const VectorBase<ER>& r, const VectorBase<EC>& c) {
01507     assert(r.nrow()==nrow() && c.nrow()==ncol());
01508     for (int j=0; j<ncol(); ++j)
01509         for (int i=0; i<nrow(); ++i)
01510             (*this)(i,j) *= (r[i]*c[j]);
01511     return *this;
01512 }
01513 
01514 template <class ELT> template <class ER, class EC> inline void
01515 MatrixBase<ELT>::rowAndColScale(
01516     const VectorBase<ER>& r, 
01517     const VectorBase<EC>& c,
01518     typename EltResult<typename VectorBase<ER>::template EltResult<EC>::Mul>::Mul& 
01519                           out) const
01520 {
01521     assert(r.nrow()==nrow() && c.nrow()==ncol());
01522     out.resize(nrow(), ncol());
01523     for (int j=0; j<ncol(); ++j)
01524         for (int i=0; i<nrow(); ++i)
01525             out(i,j) = (*this)(i,j) * (r[i]*c[j]);
01526 }
01527 
01528 
01529 // Set M(i,j) = M(i,j)^-1.
01530 template <class ELT> inline MatrixBase<ELT>& 
01531 MatrixBase<ELT>::elementwiseInvertInPlace() {
01532     const int nr=nrow(), nc=ncol();
01533     for (int j=0; j<nc; ++j)
01534         for (int i=0; i<nr; ++i) {
01535             ELT& e = updElt(i,j);
01536             e = CNT<ELT>::invert(e);
01537         }
01538     return *this;
01539 }
01540 
01541 template <class ELT> inline void
01542 MatrixBase<ELT>::elementwiseInvert(MatrixBase<typename CNT<E>::TInvert>& out) const {
01543     const int nr=nrow(), nc=ncol();
01544     out.resize(nr,nc);
01545     for (int j=0; j<nc; ++j)
01546         for (int i=0; i<nr; ++i)
01547             out(i,j) = CNT<ELT>::invert((*this)(i,j));
01548 }
01549 
01550 // M(i,j) += s
01551 template <class ELT> template <class S> inline MatrixBase<ELT>& 
01552 MatrixBase<ELT>::elementwiseAddScalarInPlace(const S& s) {
01553     for (int j=0; j<ncol(); ++j)
01554         for (int i=0; i<nrow(); ++i)
01555             (*this)(i,j) += s;
01556     return *this;
01557 }
01558 
01559 template <class ELT> template <class S> inline void 
01560 MatrixBase<ELT>::elementwiseAddScalar(
01561     const S& s,
01562     typename MatrixBase<ELT>::template EltResult<S>::Add& out) const
01563 {
01564     const int nr=nrow(), nc=ncol();
01565     out.resize(nr,nc);
01566     for (int j=0; j<nc; ++j)
01567         for (int i=0; i<nr; ++i)
01568             out(i,j) = (*this)(i,j) + s;
01569 }
01570 
01571 // M(i,j) -= s
01572 template <class ELT> template <class S> inline MatrixBase<ELT>& 
01573 MatrixBase<ELT>::elementwiseSubtractScalarInPlace(const S& s) {
01574     for (int j=0; j<ncol(); ++j)
01575         for (int i=0; i<nrow(); ++i)
01576             (*this)(i,j) -= s;
01577     return *this;
01578 }
01579 
01580 template <class ELT> template <class S> inline void 
01581 MatrixBase<ELT>::elementwiseSubtractScalar(
01582     const S& s,
01583     typename MatrixBase<ELT>::template EltResult<S>::Sub& out) const
01584 {
01585     const int nr=nrow(), nc=ncol();
01586     out.resize(nr,nc);
01587     for (int j=0; j<nc; ++j)
01588         for (int i=0; i<nr; ++i)
01589             out(i,j) = (*this)(i,j) - s;
01590 }
01591 
01592 // M(i,j) = s - M(i,j)
01593 template <class ELT> template <class S> inline MatrixBase<ELT>& 
01594 MatrixBase<ELT>::elementwiseSubtractFromScalarInPlace(const S& s) {
01595     const int nr=nrow(), nc=ncol();
01596     for (int j=0; j<nc; ++j)
01597         for (int i=0; i<nr; ++i) {
01598             ELT& e = updElt(i,j);
01599             e = s - e;
01600         }
01601     return *this;
01602 }
01603 
01604 template <class ELT> template <class S> inline void 
01605 MatrixBase<ELT>::elementwiseSubtractFromScalar(
01606     const S& s,
01607     typename MatrixBase<S>::template EltResult<ELT>::Sub& out) const
01608 {
01609     const int nr=nrow(), nc=ncol();
01610     out.resize(nr,nc);
01611     for (int j=0; j<nc; ++j)
01612         for (int i=0; i<nr; ++i)
01613             out(i,j) = s - (*this)(i,j);
01614 }
01615 
01616 // M(i,j) *= R(i,j); R must have same dimensions as this
01617 template <class ELT> template <class EE> inline MatrixBase<ELT>& 
01618 MatrixBase<ELT>::elementwiseMultiplyInPlace(const MatrixBase<EE>& r) {
01619     const int nr=nrow(), nc=ncol();
01620     assert(r.nrow()==nr && r.ncol()==nc);
01621     for (int j=0; j<nc; ++j)
01622         for (int i=0; i<nr; ++i)
01623             (*this)(i,j) *= r(i,j);
01624     return *this;
01625 }
01626 
01627 template <class ELT> template <class EE> inline void 
01628 MatrixBase<ELT>::elementwiseMultiply(
01629     const MatrixBase<EE>& r, 
01630     typename MatrixBase<ELT>::template EltResult<EE>::Mul& out) const
01631 {
01632     const int nr=nrow(), nc=ncol();
01633     assert(r.nrow()==nr && r.ncol()==nc);
01634     out.resize(nr,nc);
01635     for (int j=0; j<nc; ++j)
01636         for (int i=0; i<nr; ++i)
01637             out(i,j) = (*this)(i,j) * r(i,j);
01638 }
01639 
01640 // M(i,j) = R(i,j) * M(i,j); R must have same dimensions as this
01641 template <class ELT> template <class EE> inline MatrixBase<ELT>& 
01642 MatrixBase<ELT>::elementwiseMultiplyFromLeftInPlace(const MatrixBase<EE>& r) {
01643     const int nr=nrow(), nc=ncol();
01644     assert(r.nrow()==nr && r.ncol()==nc);
01645     for (int j=0; j<nc; ++j)
01646         for (int i=0; i<nr; ++i) {
01647             ELT& e = updElt(i,j);
01648             e = r(i,j) * e;
01649         }
01650     return *this;
01651 }
01652 
01653 template <class ELT> template <class EE> inline void 
01654 MatrixBase<ELT>::elementwiseMultiplyFromLeft(
01655     const MatrixBase<EE>& r, 
01656     typename MatrixBase<EE>::template EltResult<ELT>::Mul& out) const
01657 {
01658     const int nr=nrow(), nc=ncol();
01659     assert(r.nrow()==nr && r.ncol()==nc);
01660     out.resize(nr,nc);
01661     for (int j=0; j<nc; ++j)
01662         for (int i=0; i<nr; ++i)
01663             out(i,j) =  r(i,j) * (*this)(i,j);
01664 }
01665 
01666 // M(i,j) /= R(i,j); R must have same dimensions as this
01667 template <class ELT> template <class EE> inline MatrixBase<ELT>& 
01668 MatrixBase<ELT>::elementwiseDivideInPlace(const MatrixBase<EE>& r) {
01669     const int nr=nrow(), nc=ncol();
01670     assert(r.nrow()==nr && r.ncol()==nc);
01671     for (int j=0; j<nc; ++j)
01672         for (int i=0; i<nr; ++i)
01673             (*this)(i,j) /= r(i,j);
01674     return *this;
01675 }
01676 
01677 template <class ELT> template <class EE> inline void 
01678 MatrixBase<ELT>::elementwiseDivide(
01679     const MatrixBase<EE>& r,
01680     typename MatrixBase<ELT>::template EltResult<EE>::Dvd& out) const
01681 {
01682     const int nr=nrow(), nc=ncol();
01683     assert(r.nrow()==nr && r.ncol()==nc);
01684     out.resize(nr,nc);
01685     for (int j=0; j<nc; ++j)
01686         for (int i=0; i<nr; ++i)
01687             out(i,j) = (*this)(i,j) / r(i,j);
01688 }
01689 // M(i,j) = R(i,j) / M(i,j); R must have same dimensions as this
01690 template <class ELT> template <class EE> inline MatrixBase<ELT>& 
01691 MatrixBase<ELT>::elementwiseDivideFromLeftInPlace(const MatrixBase<EE>& r) {
01692     const int nr=nrow(), nc=ncol();
01693     assert(r.nrow()==nr && r.ncol()==nc);
01694     for (int j=0; j<nc; ++j)
01695         for (int i=0; i<nr; ++i) {
01696             ELT& e = updElt(i,j);
01697             e = r(i,j) / e;
01698         }
01699     return *this;
01700 }
01701 
01702 template <class ELT> template <class EE> inline void 
01703 MatrixBase<ELT>::elementwiseDivideFromLeft(
01704     const MatrixBase<EE>& r,
01705     typename MatrixBase<EE>::template EltResult<ELT>::Dvd& out) const
01706 {
01707     const int nr=nrow(), nc=ncol();
01708     assert(r.nrow()==nr && r.ncol()==nc);
01709     out.resize(nr,nc);
01710     for (int j=0; j<nc; ++j)
01711         for (int i=0; i<nr; ++i)
01712             out(i,j) = r(i,j) / (*this)(i,j);
01713 }
01714 
01715 /*
01716 template <class ELT> inline MatrixView_< typename CNT<ELT>::TReal > 
01717 MatrixBase<ELT>::real() const { 
01718     if (!CNT<ELT>::IsComplex) { // known at compile time
01719         return MatrixView_< typename CNT<ELT>::TReal >( // this is just ELT
01720             MatrixHelper(helper,0,0,nrow(),ncol()));    // a view of the whole matrix
01721     }
01722     // Elements are complex -- helper uses underlying precision (real) type.
01723     MatrixHelper<Precision> h(helper,typename MatrixHelper<Precision>::RealView);    
01724     return MatrixView_< typename CNT<ELT>::TReal >(h); 
01725 }
01726 */
01727 
01728 
01733 template <class ELT> class Matrix_ : public MatrixBase<ELT> {
01734     typedef MatrixBase<ELT> Base;
01735     typedef typename CNT<ELT>::Scalar            S;
01736     typedef typename CNT<ELT>::Number            Number;
01737     typedef typename CNT<ELT>::StdNumber         StdNumber;
01738 
01739     typedef Matrix_<ELT>                         T;
01740     typedef MatrixView_<ELT>                     TView;
01741     typedef Matrix_< typename CNT<ELT>::TNeg >   TNeg;
01742 
01743 public:
01744     Matrix_() : Base() { }
01745 
01746     // Copy constructor is deep.
01747     Matrix_(const Matrix_& src) : Base(src) { }
01748 
01749     // Assignment is a deep copy and will also allow reallocation if this Matrix
01750     // doesn't have a view.
01751     Matrix_& operator=(const Matrix_& src) { 
01752         Base::operator=(src); return *this;
01753     }
01754 
01755     // Force a deep copy of the view or whatever this is.    
01756     /*explicit*/ Matrix_(const Base& v) : Base(v) { }   // e.g., MatrixView
01757 
01758     Matrix_(int m, int n) : Base(m,n) { }
01759     Matrix_(int m, int n, const ELT* initValsByRow) : Base(m,n,initValsByRow) { }
01760     Matrix_(int m, int n, const ELT& ival) : Base(m,n,ival) { }
01761     
01762     Matrix_(int m, int n, int leadingDim, const S* s): Base(m,n,leadingDim,s) { }
01763     Matrix_(int m, int n, int leadingDim,       S* s): Base(m,n,leadingDim,s) { }
01764     
01766     template <int M, int N>
01767     explicit Matrix_(const Mat<M,N,ELT>& mat) : Base(M, N) {
01768         for (int i = 0; i < M; ++i)
01769             for (int j = 0; j < N; ++j)
01770                 this->updElt(i, j) = mat(i, j);
01771     }
01772 
01773     Matrix_& operator=(const ELT& v) { Base::operator=(v); return *this; }
01774 
01775     template <class EE> Matrix_& operator=(const MatrixBase<EE>& m)
01776       { Base::operator=(m); return*this; }
01777     template <class EE> Matrix_& operator+=(const MatrixBase<EE>& m)
01778       { Base::operator+=(m); return*this; }
01779     template <class EE> Matrix_& operator-=(const MatrixBase<EE>& m)
01780       { Base::operator-=(m); return*this; }
01781 
01782     Matrix_& operator*=(const StdNumber& t) { Base::operator*=(t); return *this; }
01783     Matrix_& operator/=(const StdNumber& t) { Base::operator/=(t); return *this; }
01784     Matrix_& operator+=(const ELT& r)       { this->updDiag() += r; return *this; }
01785     Matrix_& operator-=(const ELT& r)       { this->updDiag() -= r; return *this; }  
01786 
01787     const TNeg& negate()    const {return *reinterpret_cast<const TNeg*>(this); }
01788     TNeg&       updNegate()       {return *reinterpret_cast<TNeg*>(this); }
01789 
01790     const TNeg& operator-() const {return negate();}
01791     TNeg&       operator-()       {return updNegate();}
01792    
01793 private:
01794     // NO DATA MEMBERS ALLOWED
01795 };
01796 
01797 template <class ELT> class VectorView_ : public VectorBase<ELT> {
01798     typedef VectorBase<ELT>                             Base;
01799     typedef typename CNT<ELT>::Scalar                   S;
01800     typedef typename CNT<ELT>::Number                   Number;
01801     typedef typename CNT<ELT>::StdNumber                StdNumber;
01802     typedef VectorView_<ELT>                            T;
01803     typedef VectorView_< typename CNT<ELT>::TNeg >      TNeg;
01804     typedef RowVectorView_< typename CNT<ELT>::THerm >  THerm;
01805 public:
01806     // Default construction is suppressed.
01807     // Uses default destructor.
01808 
01809     // Copy constructor is shallow. CAUTION: despite const argument, this preserves writability
01810     // if it was present in the source. This is necessary to allow temporary views to be
01811     // created and used as lvalues.
01812     VectorView_(const VectorView_& v) 
01813       : Base(const_cast<MatrixHelper<S>&>(v.helper), typename MatrixHelper<S>::ShallowCopy()) { }
01814 
01815     // Copy assignment is deep but not reallocating.
01816     VectorView_& operator=(const VectorView_& v) {
01817         Base::operator=(v); return *this;
01818     }
01819 
01820     // Ask for shallow copy    
01821     explicit VectorView_(const MatrixHelper<S>& h) : Base(h, typename MatrixHelper<S>::ShallowCopy()) { }
01822     explicit VectorView_(MatrixHelper<S>&       h) : Base(h, typename MatrixHelper<S>::ShallowCopy()) { }
01823     
01824     VectorView_& operator=(const Base& b) { Base::operator=(b); return *this; }
01825 
01826     VectorView_& operator=(const ELT& v) { Base::operator=(v); return *this; } 
01827 
01828     template <class EE> VectorView_& operator=(const VectorBase<EE>& m)
01829       { Base::operator=(m); return*this; }
01830     template <class EE> VectorView_& operator+=(const VectorBase<EE>& m)
01831       { Base::operator+=(m); return*this; }
01832     template <class EE> VectorView_& operator-=(const VectorBase<EE>& m)
01833       { Base::operator-=(m); return*this; }
01834 
01835     VectorView_& operator*=(const StdNumber& t) { Base::operator*=(t); return *this; }
01836     VectorView_& operator/=(const StdNumber& t) { Base::operator/=(t); return *this; }
01837     VectorView_& operator+=(const ELT& b) { elementwiseAddScalarInPlace(b); return *this; } 
01838     VectorView_& operator-=(const ELT& b) { elementwiseSubtractScalarInPlace(b); return *this; } 
01839 
01840 private:
01841     // NO DATA MEMBERS ALLOWED
01842     VectorView_(); // default construction suppressed (what's it a View of?)
01843 };
01844 
01848 template <class ELT> class TmpVectorViewT : public VectorBase<ELT> {
01849     typedef VectorBase<ELT>             Base;
01850     typedef typename CNT<ELT>::Scalar   S;
01851 public:
01852     TmpVectorViewT() : Base() { }
01853     explicit TmpVectorViewT(int m) : Base(m,1,false) { }
01854     explicit TmpVectorViewT(const MatrixHelper<S>& h) : Base(h) { }
01855     
01856     operator const Vector_<ELT>&() const { return *reinterpret_cast<const Vector_<ELT>*>(this); }
01857     operator Vector_<ELT>&()             { return *reinterpret_cast<Vector_<ELT>*>(this); }
01858     
01859     TmpVectorViewT* clone() const { return new TmpVectorViewT(*this); } 
01860 private:
01861     // NO DATA MEMBERS ALLOWED
01862 };
01863 
01864 template <class ELT> class Vector_ : public VectorBase<ELT> {
01865     typedef VectorBase<ELT>                 Base;
01866     typedef typename CNT<ELT>::Scalar       S;
01867     typedef typename CNT<ELT>::Number       Number;
01868     typedef typename CNT<ELT>::StdNumber    StdNumber;
01869 public:
01870     Vector_() : Base() { }  // 0x1 reallocatable
01871     // Uses default destructor.
01872 
01873     // Copy constructor is deep.
01874     Vector_(const Vector_& src) : Base(src) { }
01875 
01876     // Copy assignment is deep and can be reallocating if this Vector
01877     // has no View.
01878     Vector_& operator=(const Vector_& src) {
01879         Base::operator=(src); return*this;
01880     }
01881 
01882     explicit Vector_(const Base& src) : Base(src) { }    // e.g., VectorView
01883 
01884     explicit Vector_(int m) : Base(m,false) { }
01885     Vector_(int m, const ELT* initVals) : Base(m,false,initVals) { }
01886     Vector_(int m, const ELT& ival) : Base(m,false,ival) { }
01887 
01888     // Construct a Vector which uses borrowed space.
01889     // Last parameter is a dummy to avoid overload conflicts when ELT=S.    
01890     Vector_(int m, const S* s, bool): Base(m,m,s) { }
01891     Vector_(int m,       S* s, bool): Base(m,m,s) { }
01892     
01894     template <int M>
01895     explicit Vector_(const Vec<M,ELT>& v) : Base(M) {
01896         for (int i = 0; i < M; ++i)
01897             this->updElt(i, 0) = v(i);
01898     }
01899 
01900     Vector_& operator=(const ELT& v) { Base::operator=(v); return *this; } 
01901 
01902     template <class EE> Vector_& operator=(const VectorBase<EE>& m)
01903       { Base::operator=(m); return*this; }
01904     template <class EE> Vector_& operator+=(const VectorBase<EE>& m)
01905       { Base::operator+=(m); return*this; }
01906     template <class EE> Vector_& operator-=(const VectorBase<EE>& m)
01907       { Base::operator-=(m); return*this; }
01908 
01909     Vector_& operator*=(const StdNumber& t) { Base::operator*=(t); return *this; }
01910     Vector_& operator/=(const StdNumber& t) { Base::operator/=(t); return *this; }
01911     Vector_& operator+=(const ELT& b) { elementwiseAddScalarInPlace(b); return *this; } 
01912     Vector_& operator-=(const ELT& b) { elementwiseSubtractScalarInPlace(b); return *this; } 
01913  
01914 private:
01915     // NO DATA MEMBERS ALLOWED
01916 };
01917 
01918 
01919 template <class ELT> class RowVectorView_ : public RowVectorBase<ELT> {
01920     typedef RowVectorBase<ELT>                              Base;
01921     typedef typename CNT<ELT>::Scalar                       S;
01922     typedef typename CNT<ELT>::Number                       Number;
01923     typedef typename CNT<ELT>::StdNumber                    StdNumber;
01924     typedef RowVectorView_<ELT>                             T;
01925     typedef RowVectorView_< typename CNT<ELT>::TNeg >       TNeg;
01926     typedef VectorView_< typename CNT<ELT>::THerm >         THerm;
01927 public:
01928     // Default construction is suppressed.
01929     // Uses default destructor.
01930 
01931     // Copy constructor is shallow. CAUTION: despite const argument, this preserves writability
01932     // if it was present in the source. This is necessary to allow temporary views to be
01933     // created and used as lvalues.
01934     RowVectorView_(const RowVectorView_& r) 
01935       : Base(const_cast<MatrixHelper<S>&>(r.helper), typename MatrixHelper<S>::ShallowCopy()) { }
01936 
01937     // Copy assignment is deep but not reallocating.
01938     RowVectorView_& operator=(const RowVectorView_& r) {
01939         Base::operator=(r); return *this;
01940     }
01941 
01942     // Ask for shallow copy    
01943     explicit RowVectorView_(const MatrixHelper<S>& h) : Base(h, typename MatrixHelper<S>::ShallowCopy()) { }
01944     explicit RowVectorView_(MatrixHelper<S>&       h) : Base(h, typename MatrixHelper<S>::ShallowCopy()) { }
01945     
01946     RowVectorView_& operator=(const Base& b) { Base::operator=(b); return *this; }
01947 
01948     RowVectorView_& operator=(const ELT& v) { Base::operator=(v); return *this; } 
01949 
01950     template <class EE> RowVectorView_& operator=(const RowVectorBase<EE>& m)
01951       { Base::operator=(m); return*this; }
01952     template <class EE> RowVectorView_& operator+=(const RowVectorBase<EE>& m)
01953       { Base::operator+=(m); return*this; }
01954     template <class EE> RowVectorView_& operator-=(const RowVectorBase<EE>& m)
01955       { Base::operator-=(m); return*this; }
01956 
01957     RowVectorView_& operator*=(const StdNumber& t) { Base::operator*=(t); return *this; }
01958     RowVectorView_& operator/=(const StdNumber& t) { Base::operator/=(t); return *this; }
01959     RowVectorView_& operator+=(const ELT& b) { elementwiseAddScalarInPlace(b); return *this; } 
01960     RowVectorView_& operator-=(const ELT& b) { elementwiseSubtractScalarInPlace(b); return *this; } 
01961 
01962 private:
01963     // NO DATA MEMBERS ALLOWED
01964     RowVectorView_(); // default construction suppressed (what is it a view of?)
01965 };
01966 
01970 template <class ELT> class TmpRowVectorView_ : public RowVectorBase<ELT> {
01971     typedef RowVectorBase<ELT>          Base;
01972     typedef typename CNT<ELT>::Scalar   S;
01973 public:
01974     TmpRowVectorView_() : Base() { }
01975     TmpRowVectorView_(int n) : Base(n,false) { }
01976     TmpRowVectorView_(const MatrixHelper<S>& h) : Base(h) { }
01977     
01978     operator const RowVector_<ELT>&() const { return *reinterpret_cast<const RowVector_<ELT>*>(this); }
01979     operator RowVector_<ELT>&()             { return *reinterpret_cast<RowVector_<ELT>*>(this); }
01980     
01981     TmpRowVectorView_* clone() const { return new TmpRowVectorView_(*this); } 
01982 private:
01983     // NO DATA MEMBERS ALLOWED
01984 };
01985 
01986 
01987 template <class ELT> class RowVector_ : public RowVectorBase<ELT> {
01988     typedef RowVectorBase<ELT>              Base;
01989     typedef typename CNT<ELT>::Scalar       S;
01990     typedef typename CNT<ELT>::Number       Number;
01991     typedef typename CNT<ELT>::StdNumber    StdNumber;
01992 public:
01993     RowVector_() : Base() { }   // 1x0 reallocatable
01994     // Uses default destructor.
01995 
01996     // Copy constructor is deep.
01997     RowVector_(const RowVector_& src) : Base(src) { }
01998 
01999     // Copy assignment is deep and can be reallocating if this RowVector
02000     // has no View.
02001     RowVector_& operator=(const RowVector_& src) {
02002         Base::operator=(src); return*this;
02003     }
02004 
02005     explicit RowVector_(const Base& src) : Base(src) { }    // e.g., RowVectorView
02006 
02007     explicit RowVector_(int n) : Base(n,false) { }
02008     RowVector_(int n, const ELT* initVals) : Base(n,false,initVals) { }
02009     RowVector_(int n, const ELT& ival)     : Base(n,false,ival) { }
02010 
02011     // Construct a RowVector which uses borrowed space.
02012     // Last parameter is a dummy to avoid overload conflicts when ELT=S.    
02013     RowVector_(int n, const S* s, bool): Base(n,1,s) { }
02014     RowVector_(int n,       S* s, bool): Base(n,1,s) { }
02015     
02017     template <int M>
02018     explicit RowVector_(const Row<M,ELT>& v) : Base(M) {
02019         for (int i = 0; i < M; ++i)
02020             this->updElt(0, i) = v(i);
02021     }
02022 
02023     RowVector_& operator=(const ELT& v) { Base::operator=(v); return *this; } 
02024 
02025     template <class EE> RowVector_& operator=(const RowVectorBase<EE>& b)
02026       { Base::operator=(b); return*this; }
02027     template <class EE> RowVector_& operator+=(const RowVectorBase<EE>& b)
02028       { Base::operator+=(b); return*this; }
02029     template <class EE> RowVector_& operator-=(const RowVectorBase<EE>& b)
02030       { Base::operator-=(b); return*this; }
02031 
02032     RowVector_& operator*=(const StdNumber& t) { Base::operator*=(t); return *this; }
02033     RowVector_& operator/=(const StdNumber& t) { Base::operator/=(t); return *this; }
02034     RowVector_& operator+=(const ELT& b) { elementwiseAddScalarInPlace(b); return *this; } 
02035     RowVector_& operator-=(const ELT& b) { elementwiseSubtractScalarInPlace(b); return *this; } 
02036 
02037 private:
02038     // NO DATA MEMBERS ALLOWED
02039 };
02040 
02041     // GLOBAL OPERATORS: Matrix_
02042 
02043 // + and - allow mixed element types, but will fail to compile if the elements aren't
02044 // compatible. At run time these will fail if the dimensions are incompatible.
02045 template <class E1, class E2>
02046 Matrix_<typename CNT<E1>::template Result<E2>::Add>
02047 operator+(const MatrixBase<E1>& l, const MatrixBase<E2>& r) {
02048     return Matrix_<typename CNT<E1>::template Result<E2>::Add>(l) += r;
02049 }
02050 
02051 template <class E>
02052 Matrix_<E> operator+(const MatrixBase<E>& l, const typename CNT<E>::T& r) {
02053     return Matrix_<E>(l) += r;
02054 }
02055 
02056 template <class E>
02057 Matrix_<E> operator+(const typename CNT<E>::T& l, const MatrixBase<E>& r) {
02058     return Matrix_<E>(r) += l;
02059 }
02060 
02061 template <class E1, class E2>
02062 Matrix_<typename CNT<E1>::template Result<E2>::Sub>
02063 operator-(const MatrixBase<E1>& l, const MatrixBase<E2>& r) {
02064     return Matrix_<typename CNT<E1>::template Result<E2>::Sub>(l) -= r;
02065 }
02066 
02067 template <class E>
02068 Matrix_<E> operator-(const MatrixBase<E>& l, const typename CNT<E>::T& r) {
02069     return Matrix_<E>(l) -= r;
02070 }
02071 
02072 template <class E>
02073 Matrix_<E> operator-(const typename CNT<E>::T& l, const MatrixBase<E>& r) {
02074     Matrix_<E> temp(r.nrow(), r.ncol());
02075     temp = l;
02076     return (temp -= r);
02077 }
02078 
02079 // Scalar multiply and divide. You might wish the scalar could be
02080 // a templatized type "E2", but that would create horrible ambiguities since
02081 // E2 would match not only scalar types but everything else including
02082 // matrices.
02083 template <class E> Matrix_<E>
02084 operator*(const MatrixBase<E>& l, const typename CNT<E>::StdNumber& r) 
02085   { return Matrix_<E>(l)*=r; }
02086 
02087 template <class E> Matrix_<E>
02088 operator*(const typename CNT<E>::StdNumber& l, const MatrixBase<E>& r) 
02089   { return Matrix_<E>(r)*=l; }
02090 
02091 template <class E> Matrix_<E>
02092 operator/(const MatrixBase<E>& l, const typename CNT<E>::StdNumber& r) 
02093   { return Matrix_<E>(l)/=r; }
02094 
02095     // GLOBAL OPERATORS: Vector_
02096 
02097 template <class E1, class E2>
02098 Vector_<typename CNT<E1>::template Result<E2>::Add>
02099 operator+(const VectorBase<E1>& l, const VectorBase<E2>& r) {
02100     return Vector_<typename CNT<E1>::template Result<E2>::Add>(l) += r;
02101 }
02102 template <class E>
02103 Vector_<E> operator+(const VectorBase<E>& l, const typename CNT<E>::T& r) {
02104     return Vector_<E>(l) += r;
02105 }
02106 template <class E>
02107 Vector_<E> operator+(const typename CNT<E>::T& l, const VectorBase<E>& r) {
02108     return Vector_<E>(r) += l;
02109 }
02110 template <class E1, class E2>
02111 Vector_<typename CNT<E1>::template Result<E2>::Sub>
02112 operator-(const VectorBase<E1>& l, const VectorBase<E2>& r) {
02113     return Vector_<typename CNT<E1>::template Result<E2>::Sub>(l) -= r;
02114 }
02115 template <class E>
02116 Vector_<E> operator-(const VectorBase<E>& l, const typename CNT<E>::T& r) {
02117     return Vector_<E>(l) -= r;
02118 }
02119 template <class E>
02120 Vector_<E> operator-(const typename CNT<E>::T& l, const VectorBase<E>& r) {
02121     Vector_<E> temp(r.size());
02122     temp = l;
02123     return (temp -= r);
02124 }
02125 
02126 template <class E> Vector_<E>
02127 operator*(const VectorBase<E>& l, const typename CNT<E>::StdNumber& r) 
02128   { return Vector_<E>(l)*=r; }
02129 
02130 template <class E> Vector_<E>
02131 operator*(const typename CNT<E>::StdNumber& l, const VectorBase<E>& r) 
02132   { return Vector_<E>(r)*=l; }
02133 
02134 template <class E> Vector_<E>
02135 operator/(const VectorBase<E>& l, const typename CNT<E>::StdNumber& r) 
02136   { return Vector_<E>(l)/=r; }
02137 
02138     // GLOBAL OPERATORS: RowVector_
02139 
02140 template <class E1, class E2>
02141 RowVector_<typename CNT<E1>::template Result<E2>::Add>
02142 operator+(const RowVectorBase<E1>& l, const RowVectorBase<E2>& r) {
02143     return RowVector_<typename CNT<E1>::template Result<E2>::Add>(l) += r;
02144 }
02145 template <class E>
02146 RowVector_<E> operator+(const RowVectorBase<E>& l, const typename CNT<E>::T& r) {
02147     return RowVector_<E>(l) += r;
02148 }
02149 template <class E>
02150 RowVector_<E> operator+(const typename CNT<E>::T& l, const RowVectorBase<E>& r) {
02151     return RowVector_<E>(r) += l;
02152 }
02153 template <class E1, class E2>
02154 RowVector_<typename CNT<E1>::template Result<E2>::Sub>
02155 operator-(const RowVectorBase<E1>& l, const RowVectorBase<E2>& r) {
02156     return RowVector_<typename CNT<E1>::template Result<E2>::Sub>(l) -= r;
02157 }
02158 template <class E>
02159 RowVector_<E> operator-(const RowVectorBase<E>& l, const typename CNT<E>::T& r) {
02160     return RowVector_<E>(l) -= r;
02161 }
02162 template <class E>
02163 RowVector_<E> operator-(const typename CNT<E>::T& l, const RowVectorBase<E>& r) {
02164     RowVector_<E> temp(r.size());
02165     temp = l;
02166     return (temp -= r);
02167 }
02168 
02169 template <class E> RowVector_<E>
02170 operator*(const RowVectorBase<E>& l, const typename CNT<E>::StdNumber& r) 
02171   { return RowVector_<E>(l)*=r; }
02172 
02173 template <class E> RowVector_<E>
02174 operator*(const typename CNT<E>::StdNumber& l, const RowVectorBase<E>& r) 
02175   { return RowVector_<E>(r)*=l; }
02176 
02177 template <class E> RowVector_<E>
02178 operator/(const RowVectorBase<E>& l, const typename CNT<E>::StdNumber& r) 
02179   { return RowVector_<E>(l)/=r; }
02180 
02181 
02182     // GLOBAL OPERATORS: mixed
02183 
02184     // TODO: these should use LAPACK!
02185 
02186 // Dot product
02187 template <class E1, class E2> 
02188 typename CNT<E1>::template Result<E2>::Mul
02189 operator*(const RowVectorBase<E1>& r, const VectorBase<E2>& v) {
02190     assert(r.ncol() == v.nrow());
02191     typename CNT<E1>::template Result<E2>::Mul sum(0);
02192     for (int j=0; j < r.ncol(); ++j)
02193         sum += r(j) * v[j];
02194     return sum;
02195 }
02196 
02197 template <class E1, class E2> 
02198 Vector_<typename CNT<E1>::template Result<E2>::Mul>
02199 operator*(const MatrixBase<E1>& m, const VectorBase<E2>& v) {
02200     assert(m.ncol() == v.nrow());
02201     Vector_<typename CNT<E1>::template Result<E2>::Mul> res(m.nrow());
02202     for (int i=0; i< m.nrow(); ++i)
02203         res[i] = m[i]*v;
02204     return res;
02205 }
02206 
02207 template <class E1, class E2> 
02208 Matrix_<typename CNT<E1>::template Result<E2>::Mul>
02209 operator*(const MatrixBase<E1>& m1, const MatrixBase<E2>& m2) {
02210     assert(m1.ncol() == m2.nrow());
02211     Matrix_<typename CNT<E1>::template Result<E2>::Mul> 
02212         res(m1.nrow(),m2.ncol());
02213 
02214     for (int j=0; j < res.ncol(); ++j)
02215         for (int i=0; i < res.nrow(); ++i)
02216             res(i,j) = m1[i] * m2(j);
02217 
02218     return res;
02219 }
02220 
02221     // GLOBAL OPERATORS: I/O
02222 
02223 template <class T> inline std::ostream&
02224 operator<<(std::ostream& o, const VectorBase<T>& v)
02225 { o << "~["; for(int i=0;i<v.size();++i) o<<(i>0?" ":"")<<v[i]; 
02226     return o << "]"; }
02227 
02228 template <class T> inline std::ostream&
02229 operator<<(std::ostream& o, const RowVectorBase<T>& v)
02230 { o << "["; for(int i=0;i<v.size();++i) o<<(i>0?" ":"")<<v[i]; 
02231     return o << "]"; }
02232 
02233 template <class T> inline std::ostream&
02234 operator<<(std::ostream& o, const MatrixBase<T>& m) {
02235     for (int i=0;i<m.nrow();++i)
02236         o << std::endl << m[i];
02237     if (m.nrow()) o << std::endl;
02238     return o; 
02239 }
02240 
02241 
02242 // Friendly abbreviations for default precision vectors and matrices.
02243 
02244 typedef Vector_<Real>           Vector;
02245 typedef Vector_<Complex>        ComplexVector;
02246 
02247 typedef VectorView_<Real>       VectorView;
02248 typedef VectorView_<Complex>    ComplexVectorView;
02249 
02250 typedef RowVector_<Real>        RowVector;
02251 typedef RowVector_<Complex>     ComplexRowVector;
02252 
02253 typedef RowVectorView_<Real>    RowVectorView;
02254 typedef RowVectorView_<Complex> ComplexRowVectorView;
02255 
02256 typedef Matrix_<Real>           Matrix;
02257 typedef Matrix_<Complex>        ComplexMatrix;
02258 
02259 typedef MatrixView_<Real>       MatrixView;
02260 typedef MatrixView_<Complex>    ComplexMatrixView;
02261 
02262 
02263     
02264 // Not all combinations of structure/sparsity/storage/condition
02265 // are allowed. 
02266 
02267 class MatrixShape {
02268 public:
02269 
02270     MatrixShape() : shape(MatrixShapes::Uncommitted) { }
02271     // implicit conversion
02272     MatrixShape(MatrixShapes::Shape s) : shape(s) { }
02273     MatrixShapes::Shape shape;
02274 };
02275 
02276 class MatrixSize {
02277 public:
02278     enum Freedom {
02279         Uncommitted = 0x00,    // nrow, ncol variable
02280         FixedNRows  = 0x01,    // can't vary nrows
02281         FixedNCols  = 0x02,    // can't vary ncols
02282         Fixed       = FixedNRows | FixedNCols
02283     };
02284 
02285     MatrixSize() 
02286       : freedom(Uncommitted), nrow(0), ncol(0) { }
02287 
02288     MatrixSize(Freedom f, long nr, long nc) 
02289       : freedom(f), nrow(nr), ncol(nc) { }
02290 
02291     Freedom freedom;
02292     long nrow;
02293     long ncol;
02294 };
02295 
02296 class MatrixStructure {
02297 public:
02298 
02299     MatrixStructure() : structure(MatrixStructures::Uncommitted) { }
02300     
02301     // implicit conversion
02302     MatrixStructure(MatrixStructures::Structure ms) : structure(ms) { }
02303 
02304     MatrixStructures::Structure structure;
02305 };
02306 
02307 class MatrixSparsity {
02308 public:
02309 
02310     // If Banded, how stuck are we on a particular bandwidth?
02311     enum BandwidthFreedom {
02312         Free        = 0x00,    // upper & lower are free
02313         FixedUpper  = 0x01,
02314         FixedLower  = 0x02,
02315         Fixed       = FixedUpper | FixedLower
02316     };
02317 
02318     MatrixSparsity() 
02319       : sparsity(MatrixSparseFormats::Uncommitted), lowerBandwidth(-1), upperBandwidth(-1)
02320     {
02321     }
02322 
02323     MatrixSparsity(int lower, int upper) 
02324       : sparsity(MatrixSparseFormats::Banded), lowerBandwidth(lower), upperBandwidth(upper)
02325     {
02326         assert(lower >= 0 && upper >= 0);
02327     }
02328 
02329     MatrixSparseFormats::Sparsity sparsity;
02330     int lowerBandwidth;
02331     int upperBandwidth;
02332 };
02333 
02334 class MatrixStorage {
02335 public:
02336 
02337     enum Position {
02338         Lower,              // lower is default
02339         Upper
02340     };
02341 
02342     // OR-able
02343     enum Assumptions {
02344         None         = 0x00,
02345         UnitDiagonal = 0x01
02346     };
02347 
02348     MatrixStorage() 
02349       : storage(MatrixStorageFormats::Uncommitted), position(Lower), assumptions(None)
02350     {
02351     }
02352     
02353     // also serves as implicit conversion from Storage type
02354     MatrixStorage(MatrixStorageFormats::Storage s, Position p=Lower, Assumptions a=None) 
02355         : storage(s), position(p), assumptions(a)
02356     {
02357     }
02358 
02359     MatrixStorageFormats::Storage     storage;
02360     Position    position;
02361     Assumptions assumptions;
02362 
02363     // All the 2d formats allow a leading dimension larger
02364     // than the number of rows, producing a gap between
02365     // each column.
02366     int leadingDimension;
02367 
02368     // 1d formats allow spacing between elements. Stride==1
02369     // means packed.
02370     int stride;
02371 };
02372 
02373 class MatrixCondition {
02374 public:
02375 
02376     MatrixCondition() : condition(MatrixConditions::Uncommitted) { }
02377     // implicit conversion
02378     MatrixCondition(MatrixConditions::Condition c) : condition(c) { }
02379 
02380     MatrixConditions::Condition condition;
02381 };
02382 
02387 template <class ELT, class VECTOR_CLASS>
02388 class VectorIterator {
02389 public:
02390     typedef ELT value_type;
02391     typedef int difference_type;
02392     typedef ELT& reference;
02393     typedef ELT* pointer;
02394     typedef std::random_access_iterator_tag iterator_category;
02395     VectorIterator(VECTOR_CLASS& vector, int index) : vector(vector), index(index) {
02396     }
02397     VectorIterator(const VectorIterator& iter) : vector(iter.vector), index(iter.index) {
02398     }
02399     VectorIterator& operator=(const VectorIterator& iter) {
02400         vector = iter.vector;
02401         index = iter.index;
02402         return *this;
02403     }
02404     ELT& operator*() {
02405         assert (index >= 0 && index < vector.size());
02406         return vector[index];
02407     }
02408     ELT& operator[](int i) {
02409         assert (i >= 0 && i < vector.size());
02410         return vector[i];
02411     }
02412     VectorIterator operator++() {
02413         assert (index < vector.size());
02414         ++index;
02415         return *this;
02416     }
02417     VectorIterator operator++(int) {
02418         assert (index < vector.size());
02419         VectorIterator current = *this;
02420         ++index;
02421         return current;
02422     }
02423     VectorIterator operator--() {
02424         assert (index > 0);
02425         --index;
02426         return *this;
02427     }
02428     VectorIterator operator--(int) {
02429         assert (index > 0);
02430         VectorIterator current = *this;
02431         --index;
02432         return current;
02433     }
02434     bool operator<(VectorIterator iter) const {
02435         return (index < iter.index);
02436     }
02437     bool operator>(VectorIterator iter) const {
02438         return (index > iter.index);
02439     }
02440     bool operator<=(VectorIterator iter) const {
02441         return (index <= iter.index);
02442     }
02443     bool operator>=(VectorIterator iter) const {
02444         return (index >= iter.index);
02445     }
02446     int operator-(VectorIterator iter) const {
02447         return (index - iter.index);
02448     }
02449     VectorIterator operator-(int n) const {
02450         return VectorIterator(vector, index-n);
02451     }
02452     VectorIterator operator+(int n) const {
02453         return VectorIterator(vector, index+n);
02454     }
02455     bool operator==(VectorIterator iter) const {
02456         return (index == iter.index);
02457     }
02458     bool operator!=(VectorIterator iter) const {
02459         return (index != iter.index);
02460     }
02461 private:
02462     VECTOR_CLASS& vector;
02463     int index;
02464 };
02465 
02466 } //namespace SimTK
02467 
02468 #endif //SimTK_SIMMATRIX_BIGMATRIX_H_

Generated on Thu Feb 28 01:34:28 2008 for SimTKcommon by  doxygen 1.4.7