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-9 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 
00167 // TODO: matrix expression templates for delaying operator execution.
00168 
00169 #include "SimTKcommon/Scalar.h"
00170 #include "SimTKcommon/SmallMatrix.h"
00171 
00172 #include "SimTKcommon/internal/MatrixHelper.h"
00173 #include "SimTKcommon/internal/MatrixCharacteristics.h"
00174 
00175 #include <iostream>
00176 #include <cassert>
00177 #include <complex>
00178 #include <cstddef>
00179 #include <limits>
00180 
00181 namespace SimTK {
00182 
00183 template <class ELT>    class MatrixBase;
00184 template <class ELT>    class VectorBase;
00185 template <class ELT>    class RowVectorBase;
00186 
00187 template <class T=Real> class MatrixView_;
00188 template <class T=Real> class DeadMatrixView_;
00189 template <class T=Real> class Matrix_;
00190 template <class T=Real> class DeadMatrix_;
00191 
00192 template <class T=Real> class VectorView_;
00193 template <class T=Real> class DeadVectorView_;
00194 template <class T=Real> class Vector_;
00195 template <class T=Real> class DeadVector_;
00196 
00197 template <class T=Real> class RowVectorView_;
00198 template <class T=Real> class DeadRowVectorView_;
00199 template <class T=Real> class RowVector_;
00200 template <class T=Real> class DeadRowVector_;
00201 
00202 template <class ELT, class VECTOR_CLASS> class VectorIterator;
00203 
00204 //  -------------------------------- MatrixBase --------------------------------
00226 //  ----------------------------------------------------------------------------
00227 template <class ELT> class MatrixBase {  
00228 public:
00229     // These typedefs are handy, but despite appearances you cannot 
00230     // treat a MatrixBase as a composite numerical type. That is,
00231     // CNT<MatrixBase> will not compile, or if it does it won't be
00232     // meaningful.
00233 
00234     typedef ELT                                 E;
00235     typedef typename CNT<E>::TNeg               ENeg;
00236     typedef typename CNT<E>::TWithoutNegator    EWithoutNegator;
00237     typedef typename CNT<E>::TReal              EReal;
00238     typedef typename CNT<E>::TImag              EImag;
00239     typedef typename CNT<E>::TComplex           EComplex;
00240     typedef typename CNT<E>::THerm              EHerm;       
00241     typedef typename CNT<E>::TPosTrans          EPosTrans;
00242 
00243     typedef typename CNT<E>::TAbs               EAbs;
00244     typedef typename CNT<E>::TStandard          EStandard;
00245     typedef typename CNT<E>::TInvert            EInvert;
00246     typedef typename CNT<E>::TNormalize         ENormalize;
00247     typedef typename CNT<E>::TSqHermT           ESqHermT;
00248     typedef typename CNT<E>::TSqTHerm           ESqTHerm;
00249 
00250     typedef typename CNT<E>::Scalar             EScalar;
00251     typedef typename CNT<E>::Number             ENumber;
00252     typedef typename CNT<E>::StdNumber          EStdNumber;
00253     typedef typename CNT<E>::Precision          EPrecision;
00254     typedef typename CNT<E>::ScalarNormSq       EScalarNormSq;
00255 
00256     typedef EScalar    Scalar;        // the underlying Scalar type
00257     typedef ENumber    Number;        // negator removed from Scalar
00258     typedef EStdNumber StdNumber;     // conjugate goes to complex
00259     typedef EPrecision Precision;     // complex removed from StdNumber
00260     typedef EScalarNormSq  ScalarNormSq;      // type of scalar^2
00261 
00262     typedef MatrixBase<E>                T;
00263     typedef MatrixBase<ENeg>             TNeg;
00264     typedef MatrixBase<EWithoutNegator>  TWithoutNegator;
00265     typedef MatrixBase<EReal>            TReal;
00266     typedef MatrixBase<EImag>            TImag;
00267     typedef MatrixBase<EComplex>         TComplex;
00268     typedef MatrixBase<EHerm>            THerm; 
00269     typedef MatrixBase<E>                TPosTrans;
00270 
00271     typedef MatrixBase<EAbs>             TAbs;
00272     typedef MatrixBase<EStandard>        TStandard;
00273     typedef MatrixBase<EInvert>          TInvert;
00274     typedef MatrixBase<ENormalize>       TNormalize;
00275     typedef MatrixBase<ESqHermT>         TSqHermT;  // ~Mat*Mat
00276     typedef MatrixBase<ESqTHerm>         TSqTHerm;  // Mat*~Mat
00277 
00278     const MatrixCommitment& getCharacterCommitment() const {return helper.getCharacterCommitment();}
00279     const MatrixCharacter& getMatrixCharacter()     const {return helper.getMatrixCharacter();}
00280 
00283     void commitTo(const MatrixCommitment& mc)
00284     {   helper.commitTo(mc); }
00285 
00286     // This gives the resulting matrix type when (m(i,j) op P) is applied to each element.
00287     // It will have element types which are the regular composite result of E op P.
00288     template <class P> struct EltResult { 
00289         typedef MatrixBase<typename CNT<E>::template Result<P>::Mul> Mul;
00290         typedef MatrixBase<typename CNT<E>::template Result<P>::Dvd> Dvd;
00291         typedef MatrixBase<typename CNT<E>::template Result<P>::Add> Add;
00292         typedef MatrixBase<typename CNT<E>::template Result<P>::Sub> Sub;
00293     };
00294 
00296     int  nrow() const {return helper.nrow();}
00298     int  ncol() const {return helper.ncol();}
00299 
00307     ptrdiff_t nelt() const {return helper.nelt();}
00308 
00309     enum { 
00310         NScalarsPerElement    = CNT<E>::NActualScalars,
00311         CppNScalarsPerElement = sizeof(E) / sizeof(Scalar)
00312     };
00313   
00317     MatrixBase() : helper(NScalarsPerElement,CppNScalarsPerElement) {}
00318 
00321     MatrixBase(int m, int n) 
00322     :   helper(NScalarsPerElement,CppNScalarsPerElement,MatrixCommitment(),m,n) {}
00323 
00329     explicit MatrixBase(const MatrixCommitment& commitment) 
00330     :   helper(NScalarsPerElement,CppNScalarsPerElement,commitment) {}
00331 
00332 
00337     MatrixBase(const MatrixCommitment& commitment, int m, int n) 
00338     :   helper(NScalarsPerElement,CppNScalarsPerElement,commitment,m,n) {}
00339 
00341     MatrixBase(const MatrixBase& b)
00342       : helper(b.helper.getCharacterCommitment(), 
00343                b.helper, typename MatrixHelper<Scalar>::DeepCopy()) { }
00344 
00347     MatrixBase(const TNeg& b)
00348       : helper(b.helper.getCharacterCommitment(),
00349                b.helper, typename MatrixHelper<Scalar>::DeepCopy()) { }
00350     
00353     MatrixBase& copyAssign(const MatrixBase& b) {
00354         helper.copyAssign(b.helper);
00355         return *this;
00356     }
00357     MatrixBase& operator=(const MatrixBase& b) { return copyAssign(b); }
00358 
00359 
00368     MatrixBase& viewAssign(const MatrixBase& src) {
00369         helper.writableViewAssign(const_cast<MatrixHelper<Scalar>&>(src.helper));
00370         return *this;
00371     }
00372 
00373     // default destructor
00374 
00381     MatrixBase(const MatrixCommitment& commitment, int m, int n, const ELT& initialValue) 
00382     :   helper(NScalarsPerElement, CppNScalarsPerElement, commitment, m, n)
00383     {   helper.fillWith(reinterpret_cast<const Scalar*>(&initialValue)); }  
00384 
00395     MatrixBase(const MatrixCommitment& commitment, int m, int n, 
00396                const ELT* cppInitialValuesByRow) 
00397     :   helper(NScalarsPerElement, CppNScalarsPerElement, commitment, m, n)
00398     {   helper.copyInByRowsFromCpp(reinterpret_cast<const Scalar*>(cppInitialValuesByRow)); }
00399      
00411 
00413     MatrixBase(const MatrixCommitment& commitment, 
00414                const MatrixCharacter&  character, 
00415                int spacing, const Scalar* data) // read only data
00416     :   helper(NScalarsPerElement, CppNScalarsPerElement, 
00417                commitment, character, spacing, data) {}  
00418 
00420     MatrixBase(const MatrixCommitment& commitment, 
00421                const MatrixCharacter&  character, 
00422                int spacing, Scalar* data) // writable data
00423     :   helper(NScalarsPerElement, CppNScalarsPerElement, 
00424                commitment, character, spacing, data) {}  
00426         
00427     // Create a new MatrixBase from an existing helper. Both shallow and deep copies are possible.
00428     MatrixBase(const MatrixCommitment& commitment, 
00429                MatrixHelper<Scalar>&   source, 
00430                const typename MatrixHelper<Scalar>::ShallowCopy& shallow) 
00431     :   helper(commitment, source, shallow) {}
00432     MatrixBase(const MatrixCommitment&      commitment, 
00433                const MatrixHelper<Scalar>&  source, 
00434                const typename MatrixHelper<Scalar>::ShallowCopy& shallow) 
00435     :   helper(commitment, source, shallow) {}
00436     MatrixBase(const MatrixCommitment&      commitment, 
00437                const MatrixHelper<Scalar>&  source, 
00438                const typename MatrixHelper<Scalar>::DeepCopy& deep)    
00439     :   helper(commitment, source, deep) {}
00440 
00444     void clear() {helper.clear();}
00445 
00446     MatrixBase& operator*=(const StdNumber& t)  { helper.scaleBy(t);              return *this; }
00447     MatrixBase& operator/=(const StdNumber& t)  { helper.scaleBy(StdNumber(1)/t); return *this; }
00448     MatrixBase& operator+=(const MatrixBase& r) { helper.addIn(r.helper);         return *this; }
00449     MatrixBase& operator-=(const MatrixBase& r) { helper.subIn(r.helper);         return *this; }  
00450 
00451     template <class EE> MatrixBase(const MatrixBase<EE>& b)
00452       : helper(MatrixCommitment(),b.helper, typename MatrixHelper<Scalar>::DeepCopy()) { }
00453 
00454     template <class EE> MatrixBase& operator=(const MatrixBase<EE>& b) 
00455       { helper = b.helper; return *this; }
00456     template <class EE> MatrixBase& operator+=(const MatrixBase<EE>& b) 
00457       { helper.addIn(b.helper); return *this; }
00458     template <class EE> MatrixBase& operator-=(const MatrixBase<EE>& b) 
00459       { helper.subIn(b.helper); return *this; }
00460 
00471     MatrixBase& operator=(const ELT& t) { 
00472         setToZero(); updDiag().setTo(t); 
00473         return *this;
00474     }
00475 
00481     template <class S> inline MatrixBase&
00482     scalarAssign(const S& s) {
00483         setToZero(); updDiag().setTo(s);
00484         return *this;
00485     }
00486 
00490     template <class S> inline MatrixBase&
00491     scalarAddInPlace(const S& s) {
00492         updDiag().elementwiseAddScalarInPlace(s);
00493     }
00494 
00495 
00499     template <class S> inline MatrixBase&
00500     scalarSubtractInPlace(const S& s) {
00501         updDiag().elementwiseSubtractScalarInPlace(s);
00502     }
00503 
00506     template <class S> inline MatrixBase&
00507     scalarSubtractFromLeftInPlace(const S& s) {
00508         negateInPlace();
00509         updDiag().elementwiseAddScalarInPlace(s); // yes, add
00510     }
00511 
00518     template <class S> inline MatrixBase&
00519     scalarMultiplyInPlace(const S&);
00520 
00524     template <class S> inline MatrixBase&
00525     scalarMultiplyFromLeftInPlace(const S&);
00526 
00533     template <class S> inline MatrixBase&
00534     scalarDivideInPlace(const S&);
00535 
00539     template <class S> inline MatrixBase&
00540     scalarDivideFromLeftInPlace(const S&);
00541 
00542 
00545     template <class EE> inline MatrixBase& 
00546     rowScaleInPlace(const VectorBase<EE>&);
00547 
00550     template <class EE> inline void 
00551     rowScale(const VectorBase<EE>& r, typename EltResult<EE>::Mul& out) const;
00552 
00553     template <class EE> inline typename EltResult<EE>::Mul 
00554     rowScale(const VectorBase<EE>& r) const {
00555         typename EltResult<EE>::Mul out(nrow(), ncol()); rowScale(r,out); return out;
00556     }
00557 
00560     template <class EE> inline MatrixBase& 
00561     colScaleInPlace(const VectorBase<EE>&);
00562 
00563     template <class EE> inline void 
00564     colScale(const VectorBase<EE>& c, typename EltResult<EE>::Mul& out) const;
00565 
00566     template <class EE> inline typename EltResult<EE>::Mul
00567     colScale(const VectorBase<EE>& c) const {
00568         typename EltResult<EE>::Mul out(nrow(), ncol()); colScale(c,out); return out;
00569     }
00570 
00571 
00576     template <class ER, class EC> inline MatrixBase& 
00577     rowAndColScaleInPlace(const VectorBase<ER>& r, const VectorBase<EC>& c);
00578 
00579     template <class ER, class EC> inline void 
00580     rowAndColScale(const VectorBase<ER>& r, const VectorBase<EC>& c, 
00581                    typename EltResult<typename VectorBase<ER>::template EltResult<EC>::Mul>::Mul& out) const;
00582 
00583     template <class ER, class EC> inline typename EltResult<typename VectorBase<ER>::template EltResult<EC>::Mul>::Mul
00584     rowAndColScale(const VectorBase<ER>& r, const VectorBase<EC>& c) const {
00585         typename EltResult<typename VectorBase<ER>::template EltResult<EC>::Mul>::Mul 
00586             out(nrow(), ncol()); 
00587         rowAndColScale(r,c,out); return out;
00588     }
00589 
00597     template <class S> inline MatrixBase&
00598     elementwiseAssign(const S& s);
00599 
00601     MatrixBase& elementwiseInvertInPlace();
00602 
00603     void elementwiseInvert(MatrixBase<typename CNT<E>::TInvert>& out) const;
00604 
00605     MatrixBase<typename CNT<E>::TInvert> elementwiseInvert() const {
00606         MatrixBase<typename CNT<E>::TInvert> out(nrow(), ncol());
00607         elementwiseInvert(out);
00608         return out;
00609     }
00610 
00618     template <class S> inline MatrixBase&
00619     elementwiseAddScalarInPlace(const S& s);
00620 
00621     template <class S> inline void
00622     elementwiseAddScalar(const S& s, typename EltResult<S>::Add&) const;
00623 
00624     template <class S> inline typename EltResult<S>::Add
00625     elementwiseAddScalar(const S& s) const {
00626         typename EltResult<S>::Add out(nrow(), ncol());
00627         elementwiseAddScalar(s,out);
00628         return out;
00629     }
00630 
00638     template <class S> inline MatrixBase&
00639     elementwiseSubtractScalarInPlace(const S& s);
00640 
00641     template <class S> inline void
00642     elementwiseSubtractScalar(const S& s, typename EltResult<S>::Sub&) const;
00643 
00644     template <class S> inline typename EltResult<S>::Sub
00645     elementwiseSubtractScalar(const S& s) const {
00646         typename EltResult<S>::Sub out(nrow(), ncol());
00647         elementwiseSubtractScalar(s,out);
00648         return out;
00649     }
00650 
00659     template <class S> inline MatrixBase&
00660     elementwiseSubtractFromScalarInPlace(const S& s);
00661 
00662     template <class S> inline void
00663     elementwiseSubtractFromScalar(
00664         const S&, 
00665         typename MatrixBase<S>::template EltResult<E>::Sub&) const;
00666 
00667     template <class S> inline typename MatrixBase<S>::template EltResult<E>::Sub
00668     elementwiseSubtractFromScalar(const S& s) const {
00669         typename MatrixBase<S>::template EltResult<E>::Sub out(nrow(), ncol());
00670         elementwiseSubtractFromScalar<S>(s,out);
00671         return out;
00672     }
00673 
00675     template <class EE> inline MatrixBase& 
00676     elementwiseMultiplyInPlace(const MatrixBase<EE>&);
00677 
00678     template <class EE> inline void 
00679     elementwiseMultiply(const MatrixBase<EE>&, typename EltResult<EE>::Mul&) const;
00680 
00681     template <class EE> inline typename EltResult<EE>::Mul 
00682     elementwiseMultiply(const MatrixBase<EE>& m) const {
00683         typename EltResult<EE>::Mul out(nrow(), ncol()); 
00684         elementwiseMultiply<EE>(m,out); 
00685         return out;
00686     }
00687 
00689     template <class EE> inline MatrixBase& 
00690     elementwiseMultiplyFromLeftInPlace(const MatrixBase<EE>&);
00691 
00692     template <class EE> inline void 
00693     elementwiseMultiplyFromLeft(
00694         const MatrixBase<EE>&, 
00695         typename MatrixBase<EE>::template EltResult<E>::Mul&) const;
00696 
00697     template <class EE> inline typename MatrixBase<EE>::template EltResult<E>::Mul 
00698     elementwiseMultiplyFromLeft(const MatrixBase<EE>& m) const {
00699         typename EltResult<EE>::Mul out(nrow(), ncol()); 
00700         elementwiseMultiplyFromLeft<EE>(m,out); 
00701         return out;
00702     }
00703 
00705     template <class EE> inline MatrixBase& 
00706     elementwiseDivideInPlace(const MatrixBase<EE>&);
00707 
00708     template <class EE> inline void 
00709     elementwiseDivide(const MatrixBase<EE>&, typename EltResult<EE>::Dvd&) const;
00710 
00711     template <class EE> inline typename EltResult<EE>::Dvd 
00712     elementwiseDivide(const MatrixBase<EE>& m) const {
00713         typename EltResult<EE>::Dvd out(nrow(), ncol()); 
00714         elementwiseDivide<EE>(m,out); 
00715         return out;
00716     }
00717 
00719     template <class EE> inline MatrixBase& 
00720     elementwiseDivideFromLeftInPlace(const MatrixBase<EE>&);
00721 
00722     template <class EE> inline void 
00723     elementwiseDivideFromLeft(
00724         const MatrixBase<EE>&,
00725         typename MatrixBase<EE>::template EltResult<E>::Dvd&) const;
00726 
00727     template <class EE> inline typename MatrixBase<EE>::template EltResult<EE>::Dvd 
00728     elementwiseDivideFromLeft(const MatrixBase<EE>& m) const {
00729         typename MatrixBase<EE>::template EltResult<E>::Dvd out(nrow(), ncol()); 
00730         elementwiseDivideFromLeft<EE>(m,out); 
00731         return out;
00732     }
00733 
00735     MatrixBase& setTo(const ELT& t) {helper.fillWith(reinterpret_cast<const Scalar*>(&t)); return *this;}
00736     MatrixBase& setToNaN() {helper.fillWithScalar(CNT<StdNumber>::getNaN()); return *this;}
00737     MatrixBase& setToZero() {helper.fillWithScalar(StdNumber(0)); return *this;}
00738  
00739     // View creating operators. TODO: these should be DeadMatrixViews  
00740     inline RowVectorView_<ELT> row(int i) const;   // select a row
00741     inline RowVectorView_<ELT> updRow(int i);
00742     inline VectorView_<ELT>    col(int j) const;   // select a column
00743     inline VectorView_<ELT>    updCol(int j);
00744 
00745     RowVectorView_<ELT> operator[](int i) const {return row(i);}
00746     RowVectorView_<ELT> operator[](int i)       {return updRow(i);}
00747     VectorView_<ELT>    operator()(int j) const {return col(j);}
00748     VectorView_<ELT>    operator()(int j)       {return updCol(j);}
00749      
00750     // Select a block.
00751     inline MatrixView_<ELT> block(int i, int j, int m, int n) const;
00752     inline MatrixView_<ELT> updBlock(int i, int j, int m, int n);
00753 
00754     MatrixView_<ELT> operator()(int i, int j, int m, int n) const
00755       { return block(i,j,m,n); }
00756     MatrixView_<ELT> operator()(int i, int j, int m, int n)
00757       { return updBlock(i,j,m,n); }
00758  
00759     // Hermitian transpose.
00760     inline MatrixView_<EHerm> transpose() const;
00761     inline MatrixView_<EHerm> updTranspose();
00762 
00763     MatrixView_<EHerm> operator~() const {return transpose();}
00764     MatrixView_<EHerm> operator~()       {return updTranspose();}
00765 
00766     // Select matrix diagonal (of largest leading square if rectangular).
00767     inline VectorView_<ELT> diag() const;
00768     inline VectorView_<ELT> updDiag();
00769 
00770     // Create a view of the real or imaginary elements. TODO
00771     //inline MatrixView_<EReal> real() const;
00772     //inline MatrixView_<EReal> updReal();
00773     //inline MatrixView_<EImag> imag() const;
00774     //inline MatrixView_<EImag> updImag();
00775 
00776     // Overload "real" and "imag" for both read and write as a nod to convention. TODO
00777     //MatrixView_<EReal> real() {return updReal();}
00778     //MatrixView_<EReal> imag() {return updImag();}
00779 
00780     // TODO: this routine seems ill-advised but I need it for the IVM port at the moment
00781     TInvert invert() const {  // return a newly-allocated inverse; dump negator 
00782         TInvert m(*this);
00783         m.helper.invertInPlace();
00784         return m;   // TODO - bad: makes an extra copy
00785     }
00786 
00787     void invertInPlace() {helper.invertInPlace();}
00788 
00790     void dump(const char* msg=0) const {
00791         helper.dump(msg);
00792     }
00793 
00794 
00795     // This routine is useful for implementing friendlier Matrix expressions and operators.
00796     // It maps closely to the Level-3 BLAS family of pxxmm() routines like sgemm(). The
00797     // operation performed assumes that "this" is the result, and that "this" has 
00798     // already been sized correctly to receive the result. We'll compute
00799     //     this = beta*this + alpha*A*B
00800     // If beta is 0 then "this" can be uninitialized. If alpha is 0 we promise not
00801     // to look at A or B. The routine can work efficiently even if A and/or B are transposed
00802     // by their views, so an expression like
00803     //        C += s * ~A * ~B
00804     // can be performed with the single equivalent call
00805     //        C.matmul(1., s, Tr(A), Tr(B))
00806     // where Tr(A) indicates a transposed view of the original A.
00807     // The ultimate efficiency of this call depends on the data layout and views which
00808     // are used for the three matrices.
00809     // NOTE: neither A nor B can be the same matrix as 'this', nor views of the same data
00810     // which would expose elements of 'this' that will be modified by this operation.
00811     template <class ELT_A, class ELT_B>
00812     MatrixBase& matmul(const StdNumber& beta,   // applied to 'this'
00813                        const StdNumber& alpha, const MatrixBase<ELT_A>& A, const MatrixBase<ELT_B>& B)
00814     {
00815         helper.matmul(beta,alpha,A.helper,B.helper);
00816         return *this;
00817     }
00818 
00827     const ELT& getElt(int i, int j) const { return *reinterpret_cast<const ELT*>(helper.getElt(i,j)); }
00828     ELT&       updElt(int i, int j)       { return *reinterpret_cast<      ELT*>(helper.updElt(i,j)); }
00829 
00830     const ELT& operator()(int i, int j) const {return getElt(i,j);}
00831     ELT&       operator()(int i, int j)       {return updElt(i,j);}
00832 
00837     void getAnyElt(int i, int j, ELT& value) const
00838     {   helper.getAnyElt(i,j,reinterpret_cast<Scalar*>(&value)); }
00839     ELT getAnyElt(int i, int j) const {ELT e; getAnyElt(i,j,e); return e;}
00840 
00843     // TODO: very slow! Should be optimized at least for the case
00844     //       where ELT is a Scalar.
00845     ScalarNormSq scalarNormSqr() const {
00846         const int nr=nrow(), nc=ncol();
00847         ScalarNormSq sum(0);
00848         for(int j=0;j<nc;++j) 
00849             for (int i=0; i<nr; ++i)
00850                 sum += CNT<E>::scalarNormSqr((*this)(i,j));
00851         return sum;
00852     }
00853 
00857     // TODO: very slow! Should be optimized at least for the case
00858     //       where ELT is a Scalar.
00859     void abs(TAbs& mabs) const {
00860         const int nr=nrow(), nc=ncol();
00861         mabs.resize(nr,nc);
00862         for(int j=0;j<nc;++j) 
00863             for (int i=0; i<nr; ++i)
00864                 mabs(i,j) = CNT<E>::abs((*this)(i,j));
00865     }
00866 
00869     TAbs abs() const { TAbs mabs; abs(mabs); return mabs; }
00870 
00881     TStandard standardize() const {
00882         const int nr=nrow(), nc=ncol();
00883         TStandard mstd(nr, nc);
00884         for(int j=0;j<nc;++j) 
00885             for (int i=0; i<nr; ++i)
00886                 mstd(i,j) = CNT<E>::standardize((*this)(i,j));
00887         return mstd;
00888     }
00889 
00892     ScalarNormSq normSqr() const { return scalarNormSqr(); }
00893     // TODO -- not good; unnecessary overflow
00894     typename CNT<ScalarNormSq>::TSqrt 
00895         norm() const { return CNT<ScalarNormSq>::sqrt(scalarNormSqr()); }
00896 
00899     typename CNT<ScalarNormSq>::TSqrt 
00900     normRMS() const {
00901         if (!CNT<ELT>::IsScalar)
00902             SimTK_THROW1(Exception::Cant, "normRMS() only defined for scalar elements");
00903         if (nelt() == 0)
00904             return typename CNT<ScalarNormSq>::TSqrt(0);
00905         return CNT<ScalarNormSq>::sqrt(scalarNormSqr()/nelt());
00906     }
00907 
00909     RowVectorBase<ELT> sum() const {
00910         const int cols = ncol();
00911         RowVectorBase<ELT> row(cols);
00912         for (int i = 0; i < cols; ++i)
00913             helper.colSum(i, reinterpret_cast<Scalar*>(&row[i]));
00914         return row;
00915     }
00916 
00917     //TODO: make unary +/- return a self-view so they won't reallocate?
00918     const MatrixBase& operator+() const {return *this; }
00919     const TNeg&       negate()    const {return *reinterpret_cast<const TNeg*>(this); }
00920     TNeg&             updNegate()       {return *reinterpret_cast<TNeg*>(this); }
00921 
00922     const TNeg&       operator-() const {return negate();}
00923     TNeg&             operator-()       {return updNegate();}
00924 
00925     MatrixBase& negateInPlace() {(*this) *= EPrecision(-1);}
00926  
00931     MatrixBase& resize(int m, int n)     { helper.resize(m,n); return *this; }
00937     MatrixBase& resizeKeep(int m, int n) { helper.resizeKeep(m,n); return *this; }
00938 
00939     // These prevent shape changes in a Matrix that would otherwise allow it. No harm if they
00940     // are called on a Matrix that is locked already; they always succeed.
00941     void lockNRows() {helper.lockNRows();}
00942     void lockNCols() {helper.lockNCols();}
00943     void lockShape() {helper.lockShape();}
00944 
00945     // These allow shape changes again for a Matrix which was constructed to allow them
00946     // but had them locked with the above routines. No harm if these are called on a Matrix
00947     // that is already unlocked, but it is not allowed to call them on a Matrix which
00948     // *never* allowed resizing. An exception will be thrown in that case.
00949     void unlockNRows() {helper.unlockNRows();}
00950     void unlockNCols() {helper.unlockNCols();}
00951     void unlockShape() {helper.unlockShape();}
00952     
00953     // An assortment of handy conversions
00954     const MatrixView_<ELT>& getAsMatrixView() const { return *reinterpret_cast<const MatrixView_<ELT>*>(this); }
00955     MatrixView_<ELT>&       updAsMatrixView()       { return *reinterpret_cast<      MatrixView_<ELT>*>(this); } 
00956     const Matrix_<ELT>&     getAsMatrix()     const { return *reinterpret_cast<const Matrix_<ELT>*>(this); }
00957     Matrix_<ELT>&           updAsMatrix()           { return *reinterpret_cast<      Matrix_<ELT>*>(this); }
00958          
00959     const VectorView_<ELT>& getAsVectorView() const 
00960       { assert(ncol()==1); return *reinterpret_cast<const VectorView_<ELT>*>(this); }
00961     VectorView_<ELT>&       updAsVectorView()       
00962       { assert(ncol()==1); return *reinterpret_cast<      VectorView_<ELT>*>(this); } 
00963     const Vector_<ELT>&     getAsVector()     const 
00964       { assert(ncol()==1); return *reinterpret_cast<const Vector_<ELT>*>(this); }
00965     Vector_<ELT>&           updAsVector()           
00966       { assert(ncol()==1); return *reinterpret_cast<      Vector_<ELT>*>(this); }
00967     const VectorBase<ELT>& getAsVectorBase() const 
00968       { assert(ncol()==1); return *reinterpret_cast<const VectorBase<ELT>*>(this); }
00969     VectorBase<ELT>&       updAsVectorBase()       
00970       { assert(ncol()==1); return *reinterpret_cast<      VectorBase<ELT>*>(this); } 
00971                 
00972     const RowVectorView_<ELT>& getAsRowVectorView() const 
00973       { assert(nrow()==1); return *reinterpret_cast<const RowVectorView_<ELT>*>(this); }
00974     RowVectorView_<ELT>&       updAsRowVectorView()       
00975       { assert(nrow()==1); return *reinterpret_cast<      RowVectorView_<ELT>*>(this); } 
00976     const RowVector_<ELT>&     getAsRowVector()     const 
00977       { assert(nrow()==1); return *reinterpret_cast<const RowVector_<ELT>*>(this); }
00978     RowVector_<ELT>&           updAsRowVector()           
00979       { assert(nrow()==1); return *reinterpret_cast<      RowVector_<ELT>*>(this); }        
00980     const RowVectorBase<ELT>& getAsRowVectorBase() const 
00981       { assert(nrow()==1); return *reinterpret_cast<const RowVectorBase<ELT>*>(this); }
00982     RowVectorBase<ELT>&       updAsRowVectorBase()       
00983       { assert(nrow()==1); return *reinterpret_cast<      RowVectorBase<ELT>*>(this); } 
00984 
00985     // Access to raw data. We have to return the raw data
00986     // pointer as pointer-to-scalar because we may pack the elements tighter
00987     // than a C++ array would.
00988 
00992     int getNScalarsPerElement()  const {return NScalarsPerElement;}
00993 
00997     int getPackedSizeofElement() const {return NScalarsPerElement*sizeof(Scalar);}
00998 
00999     bool hasContiguousData() const {return helper.hasContiguousData();}
01000     long getContiguousScalarDataLength() const {
01001         return helper.getContiguousDataLength();
01002     }
01003     const Scalar* getContiguousScalarData() const {
01004         return helper.getContiguousData();
01005     }
01006     Scalar* updContiguousScalarData() {
01007         return helper.updContiguousData();
01008     }
01009     void replaceContiguousScalarData(Scalar* newData, long length, bool takeOwnership) {
01010         helper.replaceContiguousData(newData,length,takeOwnership);
01011     }
01012     void replaceContiguousScalarData(const Scalar* newData, long length) {
01013         helper.replaceContiguousData(newData,length);
01014     }
01015     void swapOwnedContiguousScalarData(Scalar* newData, int length, Scalar*& oldData) {
01016         helper.swapOwnedContiguousData(newData,length,oldData);
01017     }
01018 
01023     explicit MatrixBase(MatrixHelperRep<Scalar>* hrep) : helper(hrep) {}
01024 
01025 protected:
01026     const MatrixHelper<Scalar>& getHelper() const {return helper;}
01027     MatrixHelper<Scalar>&       updHelper()       {return helper;}
01028 
01029 private:
01030     MatrixHelper<Scalar> helper; // this is just one pointer
01031 
01032     template <class EE> friend class MatrixBase;
01033 };
01034 
01035 
01036 
01037 //  -------------------------------- VectorBase --------------------------------
01041 //  ----------------------------------------------------------------------------
01042 template <class ELT> class VectorBase : public MatrixBase<ELT> {
01043     typedef MatrixBase<ELT>                             Base;
01044     typedef typename CNT<ELT>::Scalar                   Scalar;
01045     typedef typename CNT<ELT>::Number                   Number;
01046     typedef typename CNT<ELT>::StdNumber                StdNumber;
01047     typedef VectorBase<ELT>                             T;
01048     typedef VectorBase<typename CNT<ELT>::TAbs>         TAbs;
01049     typedef VectorBase<typename CNT<ELT>::TNeg>         TNeg;
01050     typedef RowVectorView_<typename CNT<ELT>::THerm>    THerm;
01051 public:  
01052     //  ------------------------------------------------------------------------
01062 
01065     explicit VectorBase(int m=0) : Base(MatrixCommitment::Vector(), m, 1) {}
01066 
01070     VectorBase(const VectorBase& source) : Base(source) {}
01071 
01073     VectorBase(const TNeg& source) : Base(source) {}
01074 
01077     VectorBase(int m, const ELT& initialValue)
01078     :   Base(MatrixCommitment::Vector(),m,1,initialValue) {}  
01079 
01084     VectorBase(int m, const ELT* cppInitialValues)
01085     :   Base(MatrixCommitment::Vector(),m,1,cppInitialValues) {}
01087 
01088     //  ------------------------------------------------------------------------
01097 
01099     VectorBase(int m, int stride, const Scalar* s)
01100     :   Base(MatrixCommitment::Vector(m), MatrixCharacter::Vector(m),stride,s) { }
01102     VectorBase(int m, int stride, Scalar* s)
01103     :   Base(MatrixCommitment::Vector(m), MatrixCharacter::Vector(m),stride,s) { }
01105         
01106     //  ------------------------------------------------------------------------
01113 
01115     VectorBase(MatrixHelper<Scalar>& h, const typename MatrixHelper<Scalar>::ShallowCopy& s) 
01116     :   Base(MatrixCommitment::Vector(), h,s) { }
01118     VectorBase(const MatrixHelper<Scalar>& h, const typename MatrixHelper<Scalar>::ShallowCopy& s) 
01119     :   Base(MatrixCommitment::Vector(), h,s) { }
01121     VectorBase(const MatrixHelper<Scalar>& h, const typename MatrixHelper<Scalar>::DeepCopy& d)    
01122     :   Base(MatrixCommitment::Vector(), h,d) { }
01124 
01125     // This gives the resulting vector type when (v[i] op P) is applied to each element.
01126     // It will have element types which are the regular composite result of ELT op P.
01127     template <class P> struct EltResult { 
01128         typedef VectorBase<typename CNT<ELT>::template Result<P>::Mul> Mul;
01129         typedef VectorBase<typename CNT<ELT>::template Result<P>::Dvd> Dvd;
01130         typedef VectorBase<typename CNT<ELT>::template Result<P>::Add> Add;
01131         typedef VectorBase<typename CNT<ELT>::template Result<P>::Sub> Sub;
01132     };
01133 
01136     VectorBase& operator=(const VectorBase& b) {
01137         Base::operator=(b); return *this;
01138     }
01139 
01140     // default destructor
01141 
01142 
01143     VectorBase& operator*=(const StdNumber& t)  { Base::operator*=(t); return *this; }
01144     VectorBase& operator/=(const StdNumber& t)  { Base::operator/=(t); return *this; }
01145     VectorBase& operator+=(const VectorBase& r) { Base::operator+=(r); return *this; }
01146     VectorBase& operator-=(const VectorBase& r) { Base::operator-=(r); return *this; }  
01147 
01148 
01149     template <class EE> VectorBase& operator=(const VectorBase<EE>& b) 
01150       { Base::operator=(b);  return *this; } 
01151     template <class EE> VectorBase& operator+=(const VectorBase<EE>& b) 
01152       { Base::operator+=(b); return *this; } 
01153     template <class EE> VectorBase& operator-=(const VectorBase<EE>& b) 
01154       { Base::operator-=(b); return *this; } 
01155 
01156 
01160     VectorBase& operator=(const ELT& t) { Base::setTo(t); return *this; }  
01161 
01166     template <class EE> VectorBase& rowScaleInPlace(const VectorBase<EE>& v)
01167       { Base::template rowScaleInPlace<EE>(v); return *this; }
01168     template <class EE> inline void rowScale(const VectorBase<EE>& v, typename EltResult<EE>::Mul& out) const
01169       { Base::rowScale(v,out); }
01170     template <class EE> inline typename EltResult<EE>::Mul rowScale(const VectorBase<EE>& v) const
01171       { typename EltResult<EE>::Mul out(nrow()); Base::rowScale(v,out); return out; }
01172 
01174     VectorBase& elementwiseInvertInPlace() {
01175         Base::elementwiseInvertInPlace();
01176         return *this;
01177     }
01178 
01180     void elementwiseInvert(VectorBase<typename CNT<ELT>::TInvert>& out) const {
01181         Base::elementwiseInvert(out);
01182     }
01183 
01185     VectorBase<typename CNT<ELT>::TInvert> elementwiseInvert() const {
01186         VectorBase<typename CNT<ELT>::TInvert> out(nrow());
01187         Base::elementwiseInvert(out);
01188         return out;
01189     }
01190 
01191         // elementwise multiply
01192     template <class EE> VectorBase& elementwiseMultiplyInPlace(const VectorBase<EE>& r)
01193       { Base::template elementwiseMultiplyInPlace<EE>(r); return *this; }
01194     template <class EE> inline void elementwiseMultiply(const VectorBase<EE>& v, typename EltResult<EE>::Mul& out) const
01195       { Base::template elementwiseMultiply<EE>(v,out); }
01196     template <class EE> inline typename EltResult<EE>::Mul elementwiseMultiply(const VectorBase<EE>& v) const
01197       { typename EltResult<EE>::Mul out(nrow()); Base::template elementwiseMultiply<EE>(v,out); return out; }
01198 
01199         // elementwise multiply from left
01200     template <class EE> VectorBase& elementwiseMultiplyFromLeftInPlace(const VectorBase<EE>& r)
01201       { Base::template elementwiseMultiplyFromLeftInPlace<EE>(r); return *this; }
01202     template <class EE> inline void 
01203     elementwiseMultiplyFromLeft(
01204         const VectorBase<EE>& v, 
01205         typename VectorBase<EE>::template EltResult<ELT>::Mul& out) const
01206     { 
01207         Base::template elementwiseMultiplyFromLeft<EE>(v,out);
01208     }
01209     template <class EE> inline typename VectorBase<EE>::template EltResult<ELT>::Mul 
01210     elementwiseMultiplyFromLeft(const VectorBase<EE>& v) const
01211     { 
01212         typename VectorBase<EE>::template EltResult<ELT>::Mul out(nrow()); 
01213         Base::template elementwiseMultiplyFromLeft<EE>(v,out); 
01214         return out;
01215     }
01216 
01217         // elementwise divide
01218     template <class EE> VectorBase& elementwiseDivideInPlace(const VectorBase<EE>& r)
01219       { Base::template elementwiseDivideInPlace<EE>(r); return *this; }
01220     template <class EE> inline void elementwiseDivide(const VectorBase<EE>& v, typename EltResult<EE>::Dvd& out) const
01221       { Base::template elementwiseDivide<EE>(v,out); }
01222     template <class EE> inline typename EltResult<EE>::Dvd elementwiseDivide(const VectorBase<EE>& v) const
01223       { typename EltResult<EE>::Dvd out(nrow()); Base::template elementwiseDivide<EE>(v,out); return out; }
01224 
01225         // elementwise divide from left
01226     template <class EE> VectorBase& elementwiseDivideFromLeftInPlace(const VectorBase<EE>& r)
01227       { Base::template elementwiseDivideFromLeftInPlace<EE>(r); return *this; }
01228     template <class EE> inline void 
01229     elementwiseDivideFromLeft(
01230         const VectorBase<EE>& v, 
01231         typename VectorBase<EE>::template EltResult<ELT>::Dvd& out) const
01232     { 
01233         Base::template elementwiseDivideFromLeft<EE>(v,out);
01234     }
01235     template <class EE> inline typename VectorBase<EE>::template EltResult<ELT>::Dvd 
01236     elementwiseDivideFromLeft(const VectorBase<EE>& v) const
01237     { 
01238         typename VectorBase<EE>::template EltResult<ELT>::Dvd out(nrow()); 
01239         Base::template elementwiseDivideFromLeft<EE>(v,out); 
01240         return out;
01241     }
01242 
01243     // Implicit conversions are allowed to Vector or Matrix, but not to RowVector.   
01244     operator const Vector_<ELT>&()     const { return *reinterpret_cast<const Vector_<ELT>*>(this); }
01245     operator       Vector_<ELT>&()           { return *reinterpret_cast<      Vector_<ELT>*>(this); }
01246     operator const VectorView_<ELT>&() const { return *reinterpret_cast<const VectorView_<ELT>*>(this); }
01247     operator       VectorView_<ELT>&()       { return *reinterpret_cast<      VectorView_<ELT>*>(this); }
01248     
01249     operator const Matrix_<ELT>&()     const { return *reinterpret_cast<const Matrix_<ELT>*>(this); }
01250     operator       Matrix_<ELT>&()           { return *reinterpret_cast<      Matrix_<ELT>*>(this); } 
01251     operator const MatrixView_<ELT>&() const { return *reinterpret_cast<const MatrixView_<ELT>*>(this); }
01252     operator       MatrixView_<ELT>&()       { return *reinterpret_cast<      MatrixView_<ELT>*>(this); } 
01253 
01254 
01255     // size() for Vectors is Base::nelt() but returns int instead of ptrdiff_t.
01256     int size() const { 
01257         assert(Base::nelt() <= (ptrdiff_t)std::numeric_limits<int>::max()); 
01258         assert(Base::ncol()==1);
01259         return (int)Base::nelt();
01260     }
01261     int       nrow() const {assert(Base::ncol()==1); return Base::nrow();}
01262     int       ncol() const {assert(Base::ncol()==1); return Base::ncol();}
01263     ptrdiff_t nelt() const {assert(Base::ncol()==1); return Base::nelt();}
01264 
01265     // Override MatrixBase operators to return the right shape
01266     TAbs abs() const {TAbs result; Base::abs(result); return result;}
01267     
01268     // Override MatrixBase indexing operators          
01269     const ELT& operator[](int i) const {return *reinterpret_cast<const ELT*>(Base::getHelper().getElt(i));}
01270     ELT&       operator[](int i)       {return *reinterpret_cast<ELT*>      (Base::updHelper().updElt(i));}
01271     const ELT& operator()(int i) const {return *reinterpret_cast<const ELT*>(Base::getHelper().getElt(i));}
01272     ELT&       operator()(int i)       {return *reinterpret_cast<ELT*>      (Base::updHelper().updElt(i));}
01273          
01274     // Block (contiguous subvector) view creation      
01275     VectorView_<ELT> operator()(int i, int m) const {return Base::operator()(i,0,m,1).getAsVectorView();}
01276     VectorView_<ELT> operator()(int i, int m)       {return Base::operator()(i,0,m,1).updAsVectorView();}
01277 
01278     // Indexed view creation (arbitrary subvector). Indices must be monotonically increasing.
01279     VectorView_<ELT> index(const std::vector<int>& indices) const {
01280         MatrixHelper<Scalar> h(Base::getHelper().getCharacterCommitment(), Base::getHelper(), indices);
01281         return VectorView_<ELT>(h);
01282     }
01283     VectorView_<ELT> updIndex(const std::vector<int>& indices) {
01284         MatrixHelper<Scalar> h(Base::getHelper().getCharacterCommitment(), Base::updHelper(), indices);
01285         return VectorView_<ELT>(h);
01286     }
01287 
01288     VectorView_<ELT> operator()(const std::vector<int>& indices) const {return index(indices);}
01289     VectorView_<ELT> operator()(const std::vector<int>& indices)       {return updIndex(indices);}
01290  
01291     // Hermitian transpose.
01292     THerm transpose() const {return Base::transpose().getAsRowVectorView();}
01293     THerm updTranspose()    {return Base::updTranspose().updAsRowVectorView();}
01294 
01295     THerm operator~() const {return transpose();}
01296     THerm operator~()       {return updTranspose();}
01297 
01298     const VectorBase& operator+() const {return *this; }
01299 
01300     // Negation
01301 
01302     const TNeg& negate()    const {return *reinterpret_cast<const TNeg*>(this); }
01303     TNeg&       updNegate()       {return *reinterpret_cast<TNeg*>(this); }
01304 
01305     const TNeg& operator-() const {return negate();}
01306     TNeg&       operator-()       {return updNegate();}
01307 
01308     VectorBase& resize(int m)     {Base::resize(m,1); return *this;}
01309     VectorBase& resizeKeep(int m) {Base::resizeKeep(m,1); return *this;}
01310 
01311     //TODO: this is not re-locking the number of columns at 1.
01312     void clear() {Base::clear(); Base::resize(0,1);}
01313 
01314     ELT sum() const {ELT s; Base::getHelper().sum(reinterpret_cast<Scalar*>(&s)); return s; } // add all the elements        
01315     VectorIterator<ELT, VectorBase<ELT> > begin() {
01316         return VectorIterator<ELT, VectorBase<ELT> >(*this, 0);
01317     }
01318     VectorIterator<ELT, VectorBase<ELT> > end() {
01319         return VectorIterator<ELT, VectorBase<ELT> >(*this, size());
01320     }
01321 
01322 protected:
01323     // Create a VectorBase handle using a given helper rep. 
01324     explicit VectorBase(MatrixHelperRep<Scalar>* hrep) : Base(hrep) {}
01325 
01326 private:
01327     // NO DATA MEMBERS ALLOWED
01328 };
01329 
01330 
01331 
01332 //  ------------------------------- RowVectorBase ------------------------------
01336 //  ----------------------------------------------------------------------------
01337 template <class ELT> class RowVectorBase : public MatrixBase<ELT> {
01338     typedef MatrixBase<ELT>                             Base;
01339     typedef typename CNT<ELT>::Scalar                   Scalar;
01340     typedef typename CNT<ELT>::Number                   Number;
01341     typedef typename CNT<ELT>::StdNumber                StdNumber;
01342     typedef RowVectorBase<ELT>                          T;
01343     typedef RowVectorBase<typename CNT<ELT>::TAbs>      TAbs;
01344     typedef RowVectorBase<typename CNT<ELT>::TNeg>      TNeg;
01345     typedef VectorView_<typename CNT<ELT>::THerm>       THerm;
01346 public: 
01347     //  ------------------------------------------------------------------------
01357 
01360     explicit RowVectorBase(int n=0) : Base(MatrixCommitment::RowVector(), 1, n) {}
01361     
01365     RowVectorBase(const RowVectorBase& source) : Base(source) {}
01366 
01368     RowVectorBase(const TNeg& source) : Base(source) {}
01369 
01372     RowVectorBase(int n, const ELT& initialValue)
01373     :   Base(MatrixCommitment::RowVector(),1,n,initialValue) {}  
01374 
01379     RowVectorBase(int n, const ELT* cppInitialValues)
01380     :   Base(MatrixCommitment::RowVector(),1,n,cppInitialValues) {}
01382 
01383     //  ------------------------------------------------------------------------
01392 
01394     RowVectorBase(int n, int stride, const Scalar* s)
01395     :   Base(MatrixCommitment::RowVector(n), MatrixCharacter::RowVector(n),stride,s) { }
01397     RowVectorBase(int n, int stride, Scalar* s)
01398     :   Base(MatrixCommitment::RowVector(n), MatrixCharacter::RowVector(n),stride,s) { }
01400 
01401     //  ------------------------------------------------------------------------
01408 
01410     RowVectorBase(MatrixHelper<Scalar>& h, const typename MatrixHelper<Scalar>::ShallowCopy& s) 
01411     :   Base(MatrixCommitment::RowVector(), h,s) { }
01413     RowVectorBase(const MatrixHelper<Scalar>& h, const typename MatrixHelper<Scalar>::ShallowCopy& s) 
01414     :   Base(MatrixCommitment::RowVector(), h,s) { }
01416     RowVectorBase(const MatrixHelper<Scalar>& h, const typename MatrixHelper<Scalar>::DeepCopy& d)    
01417     :   Base(MatrixCommitment::RowVector(), h,d) { }
01419 
01420     // This gives the resulting rowvector type when (r(i) op P) is applied to each element.
01421     // It will have element types which are the regular composite result of ELT op P.
01422     template <class P> struct EltResult { 
01423         typedef RowVectorBase<typename CNT<ELT>::template Result<P>::Mul> Mul;
01424         typedef RowVectorBase<typename CNT<ELT>::template Result<P>::Dvd> Dvd;
01425         typedef RowVectorBase<typename CNT<ELT>::template Result<P>::Add> Add;
01426         typedef RowVectorBase<typename CNT<ELT>::template Result<P>::Sub> Sub;
01427     };
01428 
01431     RowVectorBase& operator=(const RowVectorBase& b) {
01432         Base::operator=(b); return *this;
01433     }
01434 
01435     // default destructor
01436 
01437     RowVectorBase& operator*=(const StdNumber& t)     {Base::operator*=(t); return *this;}
01438     RowVectorBase& operator/=(const StdNumber& t)     {Base::operator/=(t); return *this;}
01439     RowVectorBase& operator+=(const RowVectorBase& r) {Base::operator+=(r); return *this;}
01440     RowVectorBase& operator-=(const RowVectorBase& r) {Base::operator-=(r); return *this;}  
01441 
01442     template <class EE> RowVectorBase& operator=(const RowVectorBase<EE>& b) 
01443       { Base::operator=(b);  return *this; } 
01444     template <class EE> RowVectorBase& operator+=(const RowVectorBase<EE>& b) 
01445       { Base::operator+=(b); return *this; } 
01446     template <class EE> RowVectorBase& operator-=(const RowVectorBase<EE>& b) 
01447       { Base::operator-=(b); return *this; } 
01448 
01449     // default destructor
01450  
01454     RowVectorBase& operator=(const ELT& t) { Base::setTo(t); return *this; } 
01455 
01460     template <class EE> RowVectorBase& colScaleInPlace(const VectorBase<EE>& v)
01461       { Base::template colScaleInPlace<EE>(v); return *this; }
01462     template <class EE> inline void colScale(const VectorBase<EE>& v, typename EltResult<EE>::Mul& out) const
01463       { return Base::template colScale<EE>(v,out); }
01464     template <class EE> inline typename EltResult<EE>::Mul colScale(const VectorBase<EE>& v) const
01465       { typename EltResult<EE>::Mul out(ncol()); Base::template colScale<EE>(v,out); return out; }
01466 
01467 
01468         // elementwise multiply
01469     template <class EE> RowVectorBase& elementwiseMultiplyInPlace(const RowVectorBase<EE>& r)
01470       { Base::template elementwiseMultiplyInPlace<EE>(r); return *this; }
01471     template <class EE> inline void elementwiseMultiply(const RowVectorBase<EE>& v, typename EltResult<EE>::Mul& out) const
01472       { Base::template elementwiseMultiply<EE>(v,out); }
01473     template <class EE> inline typename EltResult<EE>::Mul elementwiseMultiply(const RowVectorBase<EE>& v) const
01474       { typename EltResult<EE>::Mul out(nrow()); Base::template elementwiseMultiply<EE>(v,out); return out; }
01475 
01476         // elementwise multiply from left
01477     template <class EE> RowVectorBase& elementwiseMultiplyFromLeftInPlace(const RowVectorBase<EE>& r)
01478       { Base::template elementwiseMultiplyFromLeftInPlace<EE>(r); return *this; }
01479     template <class EE> inline void 
01480     elementwiseMultiplyFromLeft(
01481         const RowVectorBase<EE>& v, 
01482         typename RowVectorBase<EE>::template EltResult<ELT>::Mul& out) const
01483     { 
01484         Base::template elementwiseMultiplyFromLeft<EE>(v,out);
01485     }
01486     template <class EE> inline typename RowVectorBase<EE>::template EltResult<ELT>::Mul 
01487     elementwiseMultiplyFromLeft(const RowVectorBase<EE>& v) const
01488     { 
01489         typename RowVectorBase<EE>::template EltResult<ELT>::Mul out(nrow()); 
01490         Base::template elementwiseMultiplyFromLeft<EE>(v,out); 
01491         return out;
01492     }
01493 
01494         // elementwise divide
01495     template <class EE> RowVectorBase& elementwiseDivideInPlace(const RowVectorBase<EE>& r)
01496       { Base::template elementwiseDivideInPlace<EE>(r); return *this; }
01497     template <class EE> inline void elementwiseDivide(const RowVectorBase<EE>& v, typename EltResult<EE>::Dvd& out) const
01498       { Base::template elementwiseDivide<EE>(v,out); }
01499     template <class EE> inline typename EltResult<EE>::Dvd elementwiseDivide(const RowVectorBase<EE>& v) const
01500       { typename EltResult<EE>::Dvd out(nrow()); Base::template elementwiseDivide<EE>(v,out); return out; }
01501 
01502         // elementwise divide from left
01503     template <class EE> RowVectorBase& elementwiseDivideFromLeftInPlace(const RowVectorBase<EE>& r)
01504       { Base::template elementwiseDivideFromLeftInPlace<EE>(r); return *this; }
01505     template <class EE> inline void 
01506     elementwiseDivideFromLeft(
01507         const RowVectorBase<EE>& v, 
01508         typename RowVectorBase<EE>::template EltResult<ELT>::Dvd& out) const
01509     { 
01510         Base::template elementwiseDivideFromLeft<EE>(v,out);
01511     }
01512     template <class EE> inline typename RowVectorBase<EE>::template EltResult<ELT>::Dvd 
01513     elementwiseDivideFromLeft(const RowVectorBase<EE>& v) const
01514     { 
01515         typename RowVectorBase<EE>::template EltResult<ELT>::Dvd out(nrow()); 
01516         Base::template elementwiseDivideFromLeft<EE>(v,out); 
01517         return out;
01518     }
01519 
01520     // Implicit conversions are allowed to RowVector or Matrix, but not to Vector.   
01521     operator const RowVector_<ELT>&()     const {return *reinterpret_cast<const RowVector_<ELT>*>(this);}
01522     operator       RowVector_<ELT>&()           {return *reinterpret_cast<      RowVector_<ELT>*>(this);}
01523     operator const RowVectorView_<ELT>&() const {return *reinterpret_cast<const RowVectorView_<ELT>*>(this);}
01524     operator       RowVectorView_<ELT>&()       {return *reinterpret_cast<      RowVectorView_<ELT>*>(this);}
01525     
01526     operator const Matrix_<ELT>&()     const {return *reinterpret_cast<const Matrix_<ELT>*>(this);}
01527     operator       Matrix_<ELT>&()           {return *reinterpret_cast<      Matrix_<ELT>*>(this);} 
01528     operator const MatrixView_<ELT>&() const {return *reinterpret_cast<const MatrixView_<ELT>*>(this);}
01529     operator       MatrixView_<ELT>&()       {return *reinterpret_cast<      MatrixView_<ELT>*>(this);} 
01530     
01531 
01532     // size() for RowVectors is Base::nelt() but returns int instead of ptrdiff_t.
01533     int size() const { 
01534         assert(Base::nelt() <= (ptrdiff_t)std::numeric_limits<int>::max()); 
01535         assert(Base::nrow()==1);
01536         return (int)Base::nelt();
01537     }
01538     int       nrow() const {assert(Base::nrow()==1); return Base::nrow();}
01539     int       ncol() const {assert(Base::nrow()==1); return Base::ncol();}
01540     ptrdiff_t nelt() const {assert(Base::nrow()==1); return Base::nelt();}
01541 
01542     // Override MatrixBase operators to return the right shape
01543     TAbs abs() const {
01544         TAbs result; Base::abs(result); return result;
01545     }
01546 
01547     // Override MatrixBase indexing operators          
01548     const ELT& operator[](int j) const {return *reinterpret_cast<const ELT*>(Base::getHelper().getElt(j));}
01549     ELT&       operator[](int j)       {return *reinterpret_cast<ELT*>      (Base::updHelper().updElt(j));}
01550     const ELT& operator()(int j) const {return *reinterpret_cast<const ELT*>(Base::getHelper().getElt(j));}
01551     ELT&       operator()(int j)       {return *reinterpret_cast<ELT*>      (Base::updHelper().updElt(j));}
01552          
01553     // Block (contiguous subvector) creation      
01554     RowVectorView_<ELT> operator()(int j, int n) const {return Base::operator()(0,j,1,n).getAsRowVectorView();}
01555     RowVectorView_<ELT> operator()(int j, int n)       {return Base::operator()(0,j,1,n).updAsRowVectorView();}
01556 
01557     // Indexed view creation (arbitrary subvector). Indices must be monotonically increasing.
01558     RowVectorView_<ELT> index(const std::vector<int>& indices) const {
01559         MatrixHelper<Scalar> h(Base::getHelper(), indices);
01560         return RowVectorView_<ELT>(h);
01561     }
01562     RowVectorView_<ELT> updIndex(const std::vector<int>& indices) {
01563         MatrixHelper<Scalar> h(Base::updHelper(), indices);
01564         return RowVectorView_<ELT>(h);
01565     }
01566 
01567     RowVectorView_<ELT> operator()(const std::vector<int>& indices) const {return index(indices);}
01568     RowVectorView_<ELT> operator()(const std::vector<int>& indices)       {return updIndex(indices);}
01569  
01570     // Hermitian transpose.
01571     THerm transpose() const {return Base::transpose().getAsVectorView();}
01572     THerm updTranspose()    {return Base::updTranspose().updAsVectorView();}
01573 
01574     THerm operator~() const {return transpose();}
01575     THerm operator~()       {return updTranspose();}
01576 
01577     const RowVectorBase& operator+() const {return *this; }
01578 
01579     // Negation
01580 
01581     const TNeg& negate()    const {return *reinterpret_cast<const TNeg*>(this); }
01582     TNeg&       updNegate()       {return *reinterpret_cast<TNeg*>(this); }
01583 
01584     const TNeg& operator-() const {return negate();}
01585     TNeg&       operator-()       {return updNegate();}
01586 
01587     RowVectorBase& resize(int n)     {Base::resize(1,n); return *this;}
01588     RowVectorBase& resizeKeep(int n) {Base::resizeKeep(1,n); return *this;}
01589 
01590     //TODO: this is not re-locking the number of rows at 1.
01591     void clear() {Base::clear(); Base::resize(1,0);}
01592 
01593     ELT sum() const {ELT s; Base::getHelper().sum(reinterpret_cast<Scalar*>(&s)); return s; } // add all the elements        
01594     VectorIterator<ELT, RowVectorBase<ELT> > begin() {
01595         return VectorIterator<ELT, RowVectorBase<ELT> >(*this, 0);
01596     }
01597     VectorIterator<ELT, RowVectorBase<ELT> > end() {
01598         return VectorIterator<ELT, RowVectorBase<ELT> >(*this, size());
01599     }
01600 
01601 protected:
01602     // Create a RowVectorBase handle using a given helper rep. 
01603     explicit RowVectorBase(MatrixHelperRep<Scalar>* hrep) : Base(hrep) {}
01604 
01605 private:
01606     // NO DATA MEMBERS ALLOWED
01607 };
01608 
01609 
01610 
01611 //  ------------------------------- MatrixView_ --------------------------------
01617 //  ----------------------------------------------------------------------------
01618 template <class ELT> class MatrixView_ : public MatrixBase<ELT> {
01619     typedef MatrixBase<ELT>                 Base;
01620     typedef typename CNT<ELT>::Scalar       S;
01621     typedef typename CNT<ELT>::StdNumber    StdNumber;
01622 public:
01623     // Default construction is suppressed.
01624     // Uses default destructor.
01625 
01626     // Create a MatrixView_ handle using a given helper rep. 
01627     explicit MatrixView_(MatrixHelperRep<S>* hrep) : Base(hrep) {}
01628 
01629     // Copy constructor is shallow. CAUTION: despite const argument, this preserves writability
01630     // if it was present in the source. This is necessary to allow temporary views to be
01631     // created and used as lvalues.
01632     MatrixView_(const MatrixView_& m) 
01633       : Base(MatrixCommitment(),
01634              const_cast<MatrixHelper<S>&>(m.getHelper()), 
01635              typename MatrixHelper<S>::ShallowCopy()) {}
01636 
01637     // Copy assignment is deep but not reallocating.
01638     MatrixView_& operator=(const MatrixView_& m) {
01639         Base::operator=(m); return *this;
01640     }
01641 
01642     // Copy construction and copy assignment from a DeadMatrixView steals the helper.
01643     MatrixView_(DeadMatrixView_<ELT>&);
01644     MatrixView_& operator=(DeadMatrixView_<ELT>&);
01645 
01646     // Ask for shallow copy    
01647     MatrixView_(const MatrixHelper<S>& h) : Base(MatrixCommitment(), h, typename MatrixHelper<S>::ShallowCopy()) { }
01648     MatrixView_(MatrixHelper<S>&       h) : Base(MatrixCommitment(), h, typename MatrixHelper<S>::ShallowCopy()) { }
01649 
01650     MatrixView_& operator=(const Matrix_<ELT>& v)     { Base::operator=(v); return *this; }
01651     MatrixView_& operator=(const ELT& e)              { Base::operator=(e); return *this; }
01652 
01653     template <class EE> MatrixView_& operator=(const MatrixBase<EE>& m)
01654       { Base::operator=(m); return *this; }
01655     template <class EE> MatrixView_& operator+=(const MatrixBase<EE>& m)
01656       { Base::operator+=(m); return *this; }
01657     template <class EE> MatrixView_& operator-=(const MatrixBase<EE>& m)
01658       { Base::operator-=(m); return *this; }
01659 
01660     MatrixView_& operator*=(const StdNumber& t) { Base::operator*=(t); return *this; }
01661     MatrixView_& operator/=(const StdNumber& t) { Base::operator/=(t); return *this; }
01662     MatrixView_& operator+=(const ELT& r)       { this->updDiag() += r; return *this; }
01663     MatrixView_& operator-=(const ELT& r)       { this->updDiag() -= r; return *this; }  
01664 
01665     operator const Matrix_<ELT>&() const { return *reinterpret_cast<const Matrix_<ELT>*>(this); }
01666     operator Matrix_<ELT>&()             { return *reinterpret_cast<Matrix_<ELT>*>(this); }      
01667 
01668 private:
01669     // NO DATA MEMBERS ALLOWED
01670     MatrixView_(); // default constructor suppressed (what's it a view of?)
01671 };
01672 
01673 
01674 
01675 //  ----------------------------- DeadMatrixView_ ------------------------------
01679 //  ----------------------------------------------------------------------------
01680 template <class ELT> class DeadMatrixView_ : public MatrixView_<ELT> {
01681     typedef MatrixView_<ELT>                Base;
01682     typedef typename CNT<ELT>::Scalar       S;
01683     typedef typename CNT<ELT>::StdNumber    StdNumber;
01684 public:
01685     // Default construction is suppressed.
01686     // Uses default destructor.
01687     
01688     // All functionality is passed through to MatrixView_.
01689     explicit DeadMatrixView_(MatrixHelperRep<S>* hrep) : Base(hrep) {}
01690     DeadMatrixView_(const Base& m) : Base(m) {}
01691     DeadMatrixView_& operator=(const Base& m) {
01692         Base::operator=(m); return *this;
01693     }
01694 
01695     // Ask for shallow copy    
01696     DeadMatrixView_(const MatrixHelper<S>& h) : Base(h) {}
01697     DeadMatrixView_(MatrixHelper<S>&       h) : Base(h) {}
01698 
01699     DeadMatrixView_& operator=(const Matrix_<ELT>& v)     { Base::operator=(v); return *this; }
01700     DeadMatrixView_& operator=(const ELT& e)              { Base::operator=(e); return *this; }
01701 
01702     template <class EE> DeadMatrixView_& operator=(const MatrixBase<EE>& m)
01703       { Base::operator=(m); return *this; }
01704     template <class EE> DeadMatrixView_& operator+=(const MatrixBase<EE>& m)
01705       { Base::operator+=(m); return *this; }
01706     template <class EE> DeadMatrixView_& operator-=(const MatrixBase<EE>& m)
01707       { Base::operator-=(m); return *this; }
01708 
01709     DeadMatrixView_& operator*=(const StdNumber& t) { Base::operator*=(t); return *this; }
01710     DeadMatrixView_& operator/=(const StdNumber& t) { Base::operator/=(t); return *this; }
01711     DeadMatrixView_& operator+=(const ELT& r)       { this->updDiag() += r; return *this; }
01712     DeadMatrixView_& operator-=(const ELT& r)       { this->updDiag() -= r; return *this; }  
01713 
01714 private:
01715     // NO DATA MEMBERS ALLOWED
01716     DeadMatrixView_(); // default constructor suppressed (what's it a view of?)
01717 };
01718 
01719 template <class ELT> inline 
01720 MatrixView_<ELT>::MatrixView_(DeadMatrixView_<ELT>& dead) 
01721 :   Base(dead.updHelper().stealRep()) {}
01722 
01723 template <class ELT> inline MatrixView_<ELT>& 
01724 MatrixView_<ELT>::operator=(DeadMatrixView_<ELT>& dead) {
01725     if (Base::getHelper().getCharacterCommitment().isSatisfiedBy(dead.getMatrixCharacter()))
01726         Base::updHelper().replaceRep(dead.updHelper().stealRep());
01727     else
01728         Base::operator=(dead);
01729     return *this;
01730 }
01731 
01732 
01733 //  ---------------------------------- Matrix_ ---------------------------------
01736 //  ----------------------------------------------------------------------------
01737 template <class ELT> class Matrix_ : public MatrixBase<ELT> {
01738     typedef typename CNT<ELT>::Scalar       S;
01739     typedef typename CNT<ELT>::Number       Number;
01740     typedef typename CNT<ELT>::StdNumber    StdNumber;
01741 
01742     typedef typename CNT<ELT>::TNeg         ENeg;
01743     typedef typename CNT<ELT>::THerm        EHerm;
01744 
01745     typedef MatrixBase<ELT>     Base;
01746     typedef MatrixBase<ENeg>    BaseNeg;
01747     typedef MatrixBase<EHerm>   BaseHerm;
01748 
01749     typedef Matrix_<ELT>        T;
01750     typedef MatrixView_<ELT>    TView;
01751     typedef Matrix_<ENeg>       TNeg;
01752 
01753 public:
01754     Matrix_() : Base() { }
01755     Matrix_(const MatrixCommitment& mc) : Base(mc) {}
01756 
01757     // Copy constructor is deep.
01758     Matrix_(const Matrix_& src) : Base(src) { }
01759 
01760     // Assignment is a deep copy and will also allow reallocation if this Matrix
01761     // doesn't have a view.
01762     Matrix_& operator=(const Matrix_& src) { 
01763         Base::operator=(src); return *this;
01764     }
01765 
01766     // Force a deep copy of the view or whatever this is.
01767     // Note that this is an implicit conversion.
01768     Matrix_(const Base& v) : Base(v) {}   // e.g., MatrixView
01769 
01770     // Allow implicit conversion from a source matrix that
01771     // has a negated version of ELT.
01772     Matrix_(const BaseNeg& v) : Base(v) {}
01773 
01774     // TODO: implicit conversion from conjugate. This is trickier
01775     // since real elements are their own conjugate so you'll get
01776     // duplicate methods defined from Matrix_(BaseHerm) and Matrix_(Base).
01777 
01778     Matrix_(int m, int n) : Base(MatrixCommitment(), m, n) {}
01779 
01780     Matrix_(int m, int n, const ELT* cppInitialValuesByRow) 
01781     :   Base(MatrixCommitment(), m, n, cppInitialValuesByRow) {}
01782     Matrix_(int m, int n, const ELT& initialValue) 
01783     :   Base(MatrixCommitment(), m, n, initialValue) {}
01784     
01785     Matrix_(int m, int n, int leadingDim, const S* data) // read only
01786     :   Base(MatrixCommitment(), MatrixCharacter::LapackFull(m,n), 
01787              leadingDim, data) {}
01788     Matrix_(int m, int n, int leadingDim, S* data) // writable
01789     :   Base(MatrixCommitment(), MatrixCharacter::LapackFull(m,n), 
01790              leadingDim, data) {}
01791     
01793     template <int M, int N, int CS, int RS>
01794     explicit Matrix_(const Mat<M,N,ELT,CS,RS>& mat)
01795     :   Base(MatrixCommitment(), M, N)
01796     {   for (int i = 0; i < M; ++i)
01797             for (int j = 0; j < N; ++j)
01798                 this->updElt(i, j) = mat(i, j); }
01799 
01800     Matrix_& operator=(const ELT& v) { Base::operator=(v); return *this; }
01801 
01802     template <class EE> Matrix_& operator=(const MatrixBase<EE>& m)
01803       { Base::operator=(m); return*this; }
01804     template <class EE> Matrix_& operator+=(const MatrixBase<EE>& m)
01805       { Base::operator+=(m); return*this; }
01806     template <class EE> Matrix_& operator-=(const MatrixBase<EE>& m)
01807       { Base::operator-=(m); return*this; }
01808 
01809     Matrix_& operator*=(const StdNumber& t) { Base::operator*=(t); return *this; }
01810     Matrix_& operator/=(const StdNumber& t) { Base::operator/=(t); return *this; }
01811     Matrix_& operator+=(const ELT& r)       { this->updDiag() += r; return *this; }
01812     Matrix_& operator-=(const ELT& r)       { this->updDiag() -= r; return *this; }  
01813 
01814     const TNeg& negate()    const {return *reinterpret_cast<const TNeg*>(this); }
01815     TNeg&       updNegate()       {return *reinterpret_cast<TNeg*>(this); }
01816 
01817     const TNeg& operator-() const {return negate();}
01818     TNeg&       operator-()       {return updNegate();}
01819    
01820 private:
01821     // NO DATA MEMBERS ALLOWED
01822 };
01823 
01824 
01825 
01826 //  -------------------------------- VectorView_ -------------------------------
01832 //  ----------------------------------------------------------------------------
01833 template <class ELT> class VectorView_ : public VectorBase<ELT> {
01834     typedef VectorBase<ELT>                             Base;
01835     typedef typename CNT<ELT>::Scalar                   S;
01836     typedef typename CNT<ELT>::Number                   Number;
01837     typedef typename CNT<ELT>::StdNumber                StdNumber;
01838     typedef VectorView_<ELT>                            T;
01839     typedef VectorView_< typename CNT<ELT>::TNeg >      TNeg;
01840     typedef RowVectorView_< typename CNT<ELT>::THerm >  THerm;
01841 public:
01842     // Default construction is suppressed.
01843     // Uses default destructor.
01844 
01845     // Create a VectorView_ handle using a given helper rep. 
01846     explicit VectorView_(MatrixHelperRep<S>* hrep) : Base(hrep) {}
01847 
01848     // Copy constructor is shallow. CAUTION: despite const argument, this preserves writability
01849     // if it was present in the source. This is necessary to allow temporary views to be
01850     // created and used as lvalues.
01851     VectorView_(const VectorView_& v) 
01852       : Base(const_cast<MatrixHelper<S>&>(v.getHelper()), typename MatrixHelper<S>::ShallowCopy()) { }
01853 
01854     // Copy assignment is deep but not reallocating.
01855     VectorView_& operator=(const VectorView_& v) {
01856         Base::operator=(v); return *this;
01857     }
01858 
01859     // Ask for shallow copy    
01860     explicit VectorView_(const MatrixHelper<S>& h) : Base(h, typename MatrixHelper<S>::ShallowCopy()) { }
01861     explicit VectorView_(MatrixHelper<S>&       h) : Base(h, typename MatrixHelper<S>::ShallowCopy()) { }
01862     
01863     VectorView_& operator=(const Base& b) { Base::operator=(b); return *this; }
01864 
01865     VectorView_& operator=(const ELT& v) { Base::operator=(v); return *this; } 
01866 
01867     template <class EE> VectorView_& operator=(const VectorBase<EE>& m)
01868       { Base::operator=(m); return*this; }
01869     template <class EE> VectorView_& operator+=(const VectorBase<EE>& m)
01870       { Base::operator+=(m); return*this; }
01871     template <class EE> VectorView_& operator-=(const VectorBase<EE>& m)
01872       { Base::operator-=(m); return*this; }
01873 
01874     VectorView_& operator*=(const StdNumber& t) { Base::operator*=(t); return *this; }
01875     VectorView_& operator/=(const StdNumber& t) { Base::operator/=(t); return *this; }
01876     VectorView_& operator+=(const ELT& b) { elementwiseAddScalarInPlace(b); return *this; } 
01877     VectorView_& operator-=(const ELT& b) { elementwiseSubtractScalarInPlace(b); return *this; } 
01878 
01879 private:
01880     // NO DATA MEMBERS ALLOWED
01881     VectorView_(); // default construction suppressed (what's it a View of?)
01882 };
01883 
01884 
01885 
01886 //  ---------------------------------- Vector_ ---------------------------------
01890 //  ----------------------------------------------------------------------------
01891 template <class ELT> class Vector_ : public VectorBase<ELT> {
01892     typedef typename CNT<ELT>::Scalar       S;
01893     typedef typename CNT<ELT>::Number       Number;
01894     typedef typename CNT<ELT>::StdNumber    StdNumber;
01895     typedef typename CNT<ELT>::TNeg         ENeg;
01896     typedef VectorBase<ELT>                 Base;
01897     typedef VectorBase<ENeg>                BaseNeg;
01898 public:
01899     Vector_() : Base() {}  // 0x1 reallocatable
01900     // Uses default destructor.
01901 
01902     // Copy constructor is deep.
01903     Vector_(const Vector_& src) : Base(src) {}
01904 
01905     // Implicit conversions.
01906     Vector_(const Base& src) : Base(src) {}    // e.g., VectorView
01907     Vector_(const BaseNeg& src) : Base(src) {}
01908 
01909     // Copy assignment is deep and can be reallocating if this Vector
01910     // has no View.
01911     Vector_& operator=(const Vector_& src) {
01912         Base::operator=(src); return*this;
01913     }
01914 
01915 
01916     explicit Vector_(int m) : Base(m) { }
01917     Vector_(int m, const ELT* cppInitialValues) : Base(m, cppInitialValues) {}
01918     Vector_(int m, const ELT& initialValue)     : Base(m, initialValue) {}
01919 
01924     Vector_(int m, const S* cppData, bool): Base(m, Base::CppNScalarsPerElement, cppData) {}
01925     Vector_(int m,       S* cppData, bool): Base(m, Base::CppNScalarsPerElement, cppData) {}
01926 
01930     Vector_(int m, int stride, const S* data, bool) : Base(m, stride, data) {}
01931     Vector_(int m, int stride,       S* data, bool) : Base(m, stride, data) {}
01932 
01934     template <int M>
01935     explicit Vector_(const Vec<M,ELT>& v) : Base(M) {
01936         for (int i = 0; i < M; ++i)
01937             this->updElt(i, 0) = v(i);
01938     }
01939 
01940     Vector_& operator=(const ELT& v) { Base::operator=(v); return *this; } 
01941 
01942     template <class EE> Vector_& operator=(const VectorBase<EE>& m)
01943       { Base::operator=(m); return*this; }
01944     template <class EE> Vector_& operator+=(const VectorBase<EE>& m)
01945       { Base::operator+=(m); return*this; }
01946     template <class EE> Vector_& operator-=(const VectorBase<EE>& m)
01947       { Base::operator-=(m); return*this; }
01948 
01949     Vector_& operator*=(const StdNumber& t) { Base::operator*=(t); return *this; }
01950     Vector_& operator/=(const StdNumber& t) { Base::operator/=(t); return *this; }
01951     Vector_& operator+=(const ELT& b) { elementwiseAddScalarInPlace(b); return *this; } 
01952     Vector_& operator-=(const ELT& b) { elementwiseSubtractScalarInPlace(b); return *this; } 
01953  
01954 private:
01955     // NO DATA MEMBERS ALLOWED
01956 };
01957 
01958 
01959 
01960 //  ------------------------------ RowVectorView_ ------------------------------
01966 //  ----------------------------------------------------------------------------
01967 template <class ELT> class RowVectorView_ : public RowVectorBase<ELT> {
01968     typedef RowVectorBase<ELT>                              Base;
01969     typedef typename CNT<ELT>::Scalar                       S;
01970     typedef typename CNT<ELT>::Number                       Number;
01971     typedef typename CNT<ELT>::StdNumber                    StdNumber;
01972     typedef RowVectorView_<ELT>                             T;
01973     typedef RowVectorView_< typename CNT<ELT>::TNeg >       TNeg;
01974     typedef VectorView_< typename CNT<ELT>::THerm >         THerm;
01975 public:
01976     // Default construction is suppressed.
01977     // Uses default destructor.
01978 
01979     // Create a RowVectorView_ handle using a given helper rep. 
01980     explicit RowVectorView_(MatrixHelperRep<S>* hrep) : Base(hrep) {}
01981 
01982     // Copy constructor is shallow. CAUTION: despite const argument, this preserves writability
01983     // if it was present in the source. This is necessary to allow temporary views to be
01984     // created and used as lvalues.
01985     RowVectorView_(const RowVectorView_& r) 
01986       : Base(const_cast<MatrixHelper<S>&>(r.getHelper()), typename MatrixHelper<S>::ShallowCopy()) { }
01987 
01988     // Copy assignment is deep but not reallocating.
01989     RowVectorView_& operator=(const RowVectorView_& r) {
01990         Base::operator=(r); return *this;
01991     }
01992 
01993     // Ask for shallow copy    
01994     explicit RowVectorView_(const MatrixHelper<S>& h) : Base(h, typename MatrixHelper<S>::ShallowCopy()) { }
01995     explicit RowVectorView_(MatrixHelper<S>&       h) : Base(h, typename MatrixHelper<S>::ShallowCopy()) { }
01996     
01997     RowVectorView_& operator=(const Base& b) { Base::operator=(b); return *this; }
01998 
01999     RowVectorView_& operator=(const ELT& v) { Base::operator=(v); return *this; } 
02000 
02001     template <class EE> RowVectorView_& operator=(const RowVectorBase<EE>& m)
02002       { Base::operator=(m); return*this; }
02003     template <class EE> RowVectorView_& operator+=(const RowVectorBase<EE>& m)
02004       { Base::operator+=(m); return*this; }
02005     template <class EE> RowVectorView_& operator-=(const RowVectorBase<EE>& m)
02006       { Base::operator-=(m); return*this; }
02007 
02008     RowVectorView_& operator*=(const StdNumber& t) { Base::operator*=(t); return *this; }
02009     RowVectorView_& operator/=(const StdNumber& t) { Base::operator/=(t); return *this; }
02010     RowVectorView_& operator+=(const ELT& b) { elementwiseAddScalarInPlace(b); return *this; } 
02011     RowVectorView_& operator-=(const ELT& b) { elementwiseSubtractScalarInPlace(b); return *this; } 
02012 
02013 private:
02014     // NO DATA MEMBERS ALLOWED
02015     RowVectorView_(); // default construction suppressed (what is it a view of?)
02016 };
02017 
02018 
02019 
02020 //  -------------------------------- RowVector_ --------------------------------
02025 //  ----------------------------------------------------------------------------
02026 template <class ELT> class RowVector_ : public RowVectorBase<ELT> {
02027     typedef typename CNT<ELT>::Scalar       S;
02028     typedef typename CNT<ELT>::Number       Number;
02029     typedef typename CNT<ELT>::StdNumber    StdNumber;
02030     typedef typename CNT<ELT>::TNeg         ENeg;
02031 
02032     typedef RowVectorBase<ELT>              Base;
02033     typedef RowVectorBase<ENeg>             BaseNeg;
02034 public:
02035     RowVector_() : Base() {}   // 1x0 reallocatable
02036     // Uses default destructor.
02037 
02038     // Copy constructor is deep.
02039     RowVector_(const RowVector_& src) : Base(src) {}
02040 
02041     // Implicit conversions.
02042     RowVector_(const Base& src) : Base(src) {}    // e.g., RowVectorView
02043     RowVector_(const BaseNeg& src) : Base(src) {}  
02044 
02045     // Copy assignment is deep and can be reallocating if this RowVector
02046     // has no View.
02047     RowVector_& operator=(const RowVector_& src) {
02048         Base::operator=(src); return*this;
02049     }
02050 
02051 
02052     explicit RowVector_(int n) : Base(n) { }
02053     RowVector_(int n, const ELT* cppInitialValues) : Base(n, cppInitialValues) {}
02054     RowVector_(int n, const ELT& initialValue)     : Base(n, initialValue) {}
02055 
02060     RowVector_(int n, const S* cppData, bool): Base(n, Base::CppNScalarsPerElement, cppData) {}
02061     RowVector_(int n,       S* cppData, bool): Base(n, Base::CppNScalarsPerElement, cppData) {}
02062 
02066     RowVector_(int n, int stride, const S* data, bool) : Base(n, stride, data) {}
02067     RowVector_(int n, int stride,       S* data, bool) : Base(n, stride, data) {}
02068     
02070     template <int M>
02071     explicit RowVector_(const Row<M,ELT>& v) : Base(M) {
02072         for (int i = 0; i < M; ++i)
02073             this->updElt(0, i) = v(i);
02074     }
02075 
02076     RowVector_& operator=(const ELT& v) { Base::operator=(v); return *this; } 
02077 
02078     template <class EE> RowVector_& operator=(const RowVectorBase<EE>& b)
02079       { Base::operator=(b); return*this; }
02080     template <class EE> RowVector_& operator+=(const RowVectorBase<EE>& b)
02081       { Base::operator+=(b); return*this; }
02082     template <class EE> RowVector_& operator-=(const RowVectorBase<EE>& b)
02083       { Base::operator-=(b); return*this; }
02084 
02085     RowVector_& operator*=(const StdNumber& t) { Base::operator*=(t); return *this; }
02086     RowVector_& operator/=(const StdNumber& t) { Base::operator/=(t); return *this; }
02087     RowVector_& operator+=(const ELT& b) { elementwiseAddScalarInPlace(b); return *this; } 
02088     RowVector_& operator-=(const ELT& b) { elementwiseSubtractScalarInPlace(b); return *this; } 
02089 
02090 private:
02091     // NO DATA MEMBERS ALLOWED
02092 };
02093 
02094 
02095 
02096 //  ------------------------ MatrixBase definitions ----------------------------
02097 
02098 template <class ELT> inline MatrixView_<ELT> 
02099 MatrixBase<ELT>::block(int i, int j, int m, int n) const { 
02100     SimTK_INDEXCHECK(0,i,nrow()+1,"MatrixBase::block()");
02101     SimTK_INDEXCHECK(0,j,ncol()+1,"MatrixBase::block()");
02102     SimTK_SIZECHECK(i+m,nrow(),"MatrixBase::block()");
02103     SimTK_SIZECHECK(j+n,ncol(),"MatrixBase::block()");
02104 
02105     MatrixHelper<Scalar> h(MatrixCommitment(),helper,i,j,m,n);    
02106     return MatrixView_<ELT>(h.stealRep()); 
02107 }
02108     
02109 template <class ELT> inline MatrixView_<ELT>
02110 MatrixBase<ELT>::updBlock(int i, int j, int m, int n) { 
02111     SimTK_INDEXCHECK(0,i,nrow()+1,"MatrixBase::updBlock()");
02112     SimTK_INDEXCHECK(0,j,ncol()+1,"MatrixBase::updBlock()");
02113     SimTK_SIZECHECK(i+m,nrow(),"MatrixBase::updBlock()");
02114     SimTK_SIZECHECK(j+n,ncol(),"MatrixBase::updBlock()");
02115 
02116     MatrixHelper<Scalar> h(MatrixCommitment(),helper,i,j,m,n);        
02117     return MatrixView_<ELT>(h.stealRep()); 
02118 }
02119 
02120 template <class E> inline MatrixView_<typename CNT<E>::THerm>
02121 MatrixBase<E>::transpose() const { 
02122     MatrixHelper<typename CNT<Scalar>::THerm> 
02123         h(MatrixCommitment(),
02124           helper, typename MatrixHelper<typename CNT<Scalar>::THerm>::TransposeView());
02125     return MatrixView_<typename CNT<E>::THerm>(h.stealRep()); 
02126 }
02127     
02128 template <class E> inline MatrixView_<typename CNT<E>::THerm>
02129 MatrixBase<E>::updTranspose() {     
02130     MatrixHelper<typename CNT<Scalar>::THerm> 
02131         h(MatrixCommitment(),
02132           helper, typename MatrixHelper<typename CNT<Scalar>::THerm>::TransposeView());
02133     return MatrixView_<typename CNT<E>::THerm>(h.stealRep()); 
02134 }
02135 
02136 template <class E> inline VectorView_<E>
02137 MatrixBase<E>::diag() const { 
02138     MatrixHelper<Scalar> h(MatrixCommitment::Vector(),
02139                            helper, typename MatrixHelper<Scalar>::DiagonalView());
02140     return VectorView_<E>(h.stealRep()); 
02141 }
02142     
02143 template <class E> inline VectorView_<E>
02144 MatrixBase<E>::updDiag() {     
02145     MatrixHelper<Scalar> h(MatrixCommitment::Vector(),
02146                            helper, typename MatrixHelper<Scalar>::DiagonalView());
02147     return VectorView_<E>(h.stealRep());
02148 }
02149 
02150 template <class ELT> inline VectorView_<ELT> 
02151 MatrixBase<ELT>::col(int j) const { 
02152     SimTK_INDEXCHECK(0,j,ncol(),"MatrixBase::col()");
02153 
02154     MatrixHelper<Scalar> h(MatrixCommitment::Vector(),
02155                            helper,0,j,nrow(),1);    
02156     return VectorView_<ELT>(h.stealRep()); 
02157 }
02158     
02159 template <class ELT> inline VectorView_<ELT>
02160 MatrixBase<ELT>::updCol(int j) {
02161     SimTK_INDEXCHECK(0,j,ncol(),"MatrixBase::updCol()");
02162 
02163     MatrixHelper<Scalar> h(MatrixCommitment::Vector(),
02164                            helper,0,j,nrow(),1);        
02165     return VectorView_<ELT>(h.stealRep()); 
02166 }
02167 
02168 template <class ELT> inline RowVectorView_<ELT> 
02169 MatrixBase<ELT>::row(int i) const { 
02170     SimTK_INDEXCHECK(0,i,nrow(),"MatrixBase::row()");
02171 
02172     MatrixHelper<Scalar> h(MatrixCommitment::RowVector(),
02173                            helper,i,0,1,ncol());    
02174     return RowVectorView_<ELT>(h.stealRep()); 
02175 }
02176     
02177 template <class ELT> inline RowVectorView_<ELT>
02178 MatrixBase<ELT>::updRow(int i) { 
02179     SimTK_INDEXCHECK(0,i,nrow(),"MatrixBase::updRow()");
02180 
02181     MatrixHelper<Scalar> h(MatrixCommitment::RowVector(),
02182                            helper,i,0,1,ncol());        
02183     return RowVectorView_<ELT>(h.stealRep()); 
02184 }
02185 
02186 // M = diag(v) * M; v must have nrow() elements.
02187 // That is, M[i] *= v[i].
02188 template <class ELT> template <class EE> inline MatrixBase<ELT>& 
02189 MatrixBase<ELT>::rowScaleInPlace(const VectorBase<EE>& v) {
02190     assert(v.nrow() == nrow());
02191     for (int i=0; i < nrow(); ++i)
02192         (*this)[i] *= v[i];
02193     return *this;
02194 }
02195 
02196 template <class ELT> template <class EE> inline void
02197 MatrixBase<ELT>::rowScale(const VectorBase<EE>& v, typename MatrixBase<ELT>::template EltResult<EE>::Mul& out) const {
02198     assert(v.nrow() == nrow());
02199     out.resize(nrow(), ncol());
02200     for (int j=0; j<ncol(); ++j)
02201         for (int i=0; i<nrow(); ++i)
02202            out(i,j) = (*this)(i,j) * v[i];
02203 }
02204 
02205 // M = M * diag(v); v must have ncol() elements
02206 // That is, M(i) *= v[i]
02207 template <class ELT> template <class EE>  inline MatrixBase<ELT>& 
02208 MatrixBase<ELT>::colScaleInPlace(const VectorBase<EE>& v) {
02209     assert(v.nrow() == ncol());
02210     for (int j=0; j < ncol(); ++j)
02211         (*this)(j) *= v[j];
02212     return *this;
02213 }
02214 
02215 template <class ELT> template <class EE> inline void
02216 MatrixBase<ELT>::colScale(const VectorBase<EE>& v, typename MatrixBase<ELT>::template EltResult<EE>::Mul& out) const {
02217     assert(v.nrow() == ncol());
02218     out.resize(nrow(), ncol());
02219     for (int j=0; j<ncol(); ++j)
02220         for (int i=0; i<nrow(); ++i)
02221            out(i,j) = (*this)(i,j) * v[j];
02222 }
02223 
02224 
02225 // M(i,j) *= r[i]*c[j]; r must have nrow() elements; c must have ncol() elements
02226 template <class ELT> template <class ER, class EC> inline MatrixBase<ELT>& 
02227 MatrixBase<ELT>::rowAndColScaleInPlace(const VectorBase<ER>& r, const VectorBase<EC>& c) {
02228     assert(r.nrow()==nrow() && c.nrow()==ncol());
02229     for (int j=0; j<ncol(); ++j)
02230         for (int i=0; i<nrow(); ++i)
02231             (*this)(i,j) *= (r[i]*c[j]);
02232     return *this;
02233 }
02234 
02235 template <class ELT> template <class ER, class EC> inline void
02236 MatrixBase<ELT>::rowAndColScale(
02237     const VectorBase<ER>& r, 
02238     const VectorBase<EC>& c,
02239     typename EltResult<typename VectorBase<ER>::template EltResult<EC>::Mul>::Mul& 
02240                           out) const
02241 {
02242     assert(r.nrow()==nrow() && c.nrow()==ncol());
02243     out.resize(nrow(), ncol());
02244     for (int j=0; j<ncol(); ++j)
02245         for (int i=0; i<nrow(); ++i)
02246             out(i,j) = (*this)(i,j) * (r[i]*c[j]);
02247 }
02248 
02249 
02250 // Set M(i,j) = M(i,j)^-1.
02251 template <class ELT> inline MatrixBase<ELT>& 
02252 MatrixBase<ELT>::elementwiseInvertInPlace() {
02253     const int nr=nrow(), nc=ncol();
02254     for (int j=0; j<nc; ++j)
02255         for (int i=0; i<nr; ++i) {
02256             ELT& e = updElt(i,j);
02257             e = CNT<ELT>::invert(e);
02258         }
02259     return *this;
02260 }
02261 
02262 template <class ELT> inline void
02263 MatrixBase<ELT>::elementwiseInvert(MatrixBase<typename CNT<E>::TInvert>& out) const {
02264     const int nr=nrow(), nc=ncol();
02265     out.resize(nr,nc);
02266     for (int j=0; j<nc; ++j)
02267         for (int i=0; i<nr; ++i)
02268             out(i,j) = CNT<ELT>::invert((*this)(i,j));
02269 }
02270 
02271 // M(i,j) += s
02272 template <class ELT> template <class S> inline MatrixBase<ELT>& 
02273 MatrixBase<ELT>::elementwiseAddScalarInPlace(const S& s) {
02274     for (int j=0; j<ncol(); ++j)
02275         for (int i=0; i<nrow(); ++i)
02276             (*this)(i,j) += s;
02277     return *this;
02278 }
02279 
02280 template <class ELT> template <class S> inline void 
02281 MatrixBase<ELT>::elementwiseAddScalar(
02282     const S& s,
02283     typename MatrixBase<ELT>::template EltResult<S>::Add& out) const
02284 {
02285     const int nr=nrow(), nc=ncol();
02286     out.resize(nr,nc);
02287     for (int j=0; j<nc; ++j)
02288         for (int i=0; i<nr; ++i)
02289             out(i,j) = (*this)(i,j) + s;
02290 }
02291 
02292 // M(i,j) -= s
02293 template <class ELT> template <class S> inline MatrixBase<ELT>& 
02294 MatrixBase<ELT>::elementwiseSubtractScalarInPlace(const S& s) {
02295     for (int j=0; j<ncol(); ++j)
02296         for (int i=0; i<nrow(); ++i)
02297             (*this)(i,j) -= s;
02298     return *this;
02299 }
02300 
02301 template <class ELT> template <class S> inline void 
02302 MatrixBase<ELT>::elementwiseSubtractScalar(
02303     const S& s,
02304     typename MatrixBase<ELT>::template EltResult<S>::Sub& out) const
02305 {
02306     const int nr=nrow(), nc=ncol();
02307     out.resize(nr,nc);
02308     for (int j=0; j<nc; ++j)
02309         for (int i=0; i<nr; ++i)
02310             out(i,j) = (*this)(i,j) - s;
02311 }
02312 
02313 // M(i,j) = s - M(i,j)
02314 template <class ELT> template <class S> inline MatrixBase<ELT>& 
02315 MatrixBase<ELT>::elementwiseSubtractFromScalarInPlace(const S& s) {
02316     const int nr=nrow(), nc=ncol();
02317     for (int j=0; j<nc; ++j)
02318         for (int i=0; i<nr; ++i) {
02319             ELT& e = updElt(i,j);
02320             e = s - e;
02321         }
02322     return *this;
02323 }
02324 
02325 template <class ELT> template <class S> inline void 
02326 MatrixBase<ELT>::elementwiseSubtractFromScalar(
02327     const S& s,
02328     typename MatrixBase<S>::template EltResult<ELT>::Sub& out) const
02329 {
02330     const int nr=nrow(), nc=ncol();
02331     out.resize(nr,nc);
02332     for (int j=0; j<nc; ++j)
02333         for (int i=0; i<nr; ++i)
02334             out(i,j) = s - (*this)(i,j);
02335 }
02336 
02337 // M(i,j) *= R(i,j); R must have same dimensions as this
02338 template <class ELT> template <class EE> inline MatrixBase<ELT>& 
02339 MatrixBase<ELT>::elementwiseMultiplyInPlace(const MatrixBase<EE>& r) {
02340     const int nr=nrow(), nc=ncol();
02341     assert(r.nrow()==nr && r.ncol()==nc);
02342     for (int j=0; j<nc; ++j)
02343         for (int i=0; i<nr; ++i)
02344             (*this)(i,j) *= r(i,j);
02345     return *this;
02346 }
02347 
02348 template <class ELT> template <class EE> inline void 
02349 MatrixBase<ELT>::elementwiseMultiply(
02350     const MatrixBase<EE>& r, 
02351     typename MatrixBase<ELT>::template EltResult<EE>::Mul& out) const
02352 {
02353     const int nr=nrow(), nc=ncol();
02354     assert(r.nrow()==nr && r.ncol()==nc);
02355     out.resize(nr,nc);
02356     for (int j=0; j<nc; ++j)
02357         for (int i=0; i<nr; ++i)
02358             out(i,j) = (*this)(i,j) * r(i,j);
02359 }
02360 
02361 // M(i,j) = R(i,j) * M(i,j); R must have same dimensions as this
02362 template <class ELT> template <class EE> inline MatrixBase<ELT>& 
02363 MatrixBase<ELT>::elementwiseMultiplyFromLeftInPlace(const MatrixBase<EE>& r) {
02364     const int nr=nrow(), nc=ncol();
02365     assert(r.nrow()==nr && r.ncol()==nc);
02366     for (int j=0; j<nc; ++j)
02367         for (int i=0; i<nr; ++i) {
02368             ELT& e = updElt(i,j);
02369             e = r(i,j) * e;
02370         }
02371     return *this;
02372 }
02373 
02374 template <class ELT> template <class EE> inline void 
02375 MatrixBase<ELT>::elementwiseMultiplyFromLeft(
02376     const MatrixBase<EE>& r, 
02377     typename MatrixBase<EE>::template EltResult<ELT>::Mul& out) const
02378 {
02379     const int nr=nrow(), nc=ncol();
02380     assert(r.nrow()==nr && r.ncol()==nc);
02381     out.resize(nr,nc);
02382     for (int j=0; j<nc; ++j)
02383         for (int i=0; i<nr; ++i)
02384             out(i,j) =  r(i,j) * (*this)(i,j);
02385 }
02386 
02387 // M(i,j) /= R(i,j); R must have same dimensions as this
02388 template <class ELT> template <class EE> inline MatrixBase<ELT>& 
02389 MatrixBase<ELT>::elementwiseDivideInPlace(const MatrixBase<EE>& r) {
02390     const int nr=nrow(), nc=ncol();
02391     assert(r.nrow()==nr && r.ncol()==nc);
02392     for (int j=0; j<nc; ++j)
02393         for (int i=0; i<nr; ++i)
02394             (*this)(i,j) /= r(i,j);
02395     return *this;
02396 }
02397 
02398 template <class ELT> template <class EE> inline void 
02399 MatrixBase<ELT>::elementwiseDivide(
02400     const MatrixBase<EE>& r,
02401     typename MatrixBase<ELT>::template EltResult<EE>::Dvd& out) const
02402 {
02403     const int nr=nrow(), nc=ncol();
02404     assert(r.nrow()==nr && r.ncol()==nc);
02405     out.resize(nr,nc);
02406     for (int j=0; j<nc; ++j)
02407         for (int i=0; i<nr; ++i)
02408             out(i,j) = (*this)(i,j) / r(i,j);
02409 }
02410 // M(i,j) = R(i,j) / M(i,j); R must have same dimensions as this
02411 template <class ELT> template <class EE> inline MatrixBase<ELT>& 
02412 MatrixBase<ELT>::elementwiseDivideFromLeftInPlace(const MatrixBase<EE>& r) {
02413     const int nr=nrow(), nc=ncol();
02414     assert(r.nrow()==nr && r.ncol()==nc);
02415     for (int j=0; j<nc; ++j)
02416         for (int i=0; i<nr; ++i) {
02417             ELT& e = updElt(i,j);
02418             e = r(i,j) / e;
02419         }
02420     return *this;
02421 }
02422 
02423 template <class ELT> template <class EE> inline void 
02424 MatrixBase<ELT>::elementwiseDivideFromLeft(
02425     const MatrixBase<EE>& r,
02426     typename MatrixBase<EE>::template EltResult<ELT>::Dvd& out) const
02427 {
02428     const int nr=nrow(), nc=ncol();
02429     assert(r.nrow()==nr && r.ncol()==nc);
02430     out.resize(nr,nc);
02431     for (int j=0; j<nc; ++j)
02432         for (int i=0; i<nr; ++i)
02433             out(i,j) = r(i,j) / (*this)(i,j);
02434 }
02435 
02436 /*
02437 template <class ELT> inline MatrixView_< typename CNT<ELT>::TReal > 
02438 MatrixBase<ELT>::real() const { 
02439     if (!CNT<ELT>::IsComplex) { // known at compile time
02440         return MatrixView_< typename CNT<ELT>::TReal >( // this is just ELT
02441             MatrixHelper(helper,0,0,nrow(),ncol()));    // a view of the whole matrix
02442     }
02443     // Elements are complex -- helper uses underlying precision (real) type.
02444     MatrixHelper<Precision> h(helper,typename MatrixHelper<Precision>::RealView);    
02445     return MatrixView_< typename CNT<ELT>::TReal >(h); 
02446 }
02447 */
02448 
02449 
02450 //  ----------------------------------------------------------------------------
02454 
02455 // + and - allow mixed element types, but will fail to compile if the elements aren't
02456 // compatible. At run time these will fail if the dimensions are incompatible.
02457 template <class E1, class E2>
02458 Matrix_<typename CNT<E1>::template Result<E2>::Add>
02459 operator+(const MatrixBase<E1>& l, const MatrixBase<E2>& r) {
02460     return Matrix_<typename CNT<E1>::template Result<E2>::Add>(l) += r;
02461 }
02462 
02463 template <class E>
02464 Matrix_<E> operator+(const MatrixBase<E>& l, const typename CNT<E>::T& r) {
02465     return Matrix_<E>(l) += r;
02466 }
02467 
02468 template <class E>
02469 Matrix_<E> operator+(const typename CNT<E>::T& l, const MatrixBase<E>& r) {
02470     return Matrix_<E>(r) += l;
02471 }
02472 
02473 template <class E1, class E2>
02474 Matrix_<typename CNT<E1>::template Result<E2>::Sub>
02475 operator-(const MatrixBase<E1>& l, const MatrixBase<E2>& r) {
02476     return Matrix_<typename CNT<E1>::template Result<E2>::Sub>(l) -= r;
02477 }
02478 
02479 template <class E>
02480 Matrix_<E> operator-(const MatrixBase<E>& l, const typename CNT<E>::T& r) {
02481     return Matrix_<E>(l) -= r;
02482 }
02483 
02484 template <class E>
02485 Matrix_<E> operator-(const typename CNT<E>::T& l, const MatrixBase<E>& r) {
02486     Matrix_<E> temp(r.nrow(), r.ncol());
02487     temp = l;
02488     return (temp -= r);
02489 }
02490 
02491 // Scalar multiply and divide. You might wish the scalar could be
02492 // a templatized type "E2", but that would create horrible ambiguities since
02493 // E2 would match not only scalar types but everything else including
02494 // matrices.
02495 template <class E> Matrix_<E>
02496 operator*(const MatrixBase<E>& l, const typename CNT<E>::StdNumber& r) 
02497   { return Matrix_<E>(l)*=r; }
02498 
02499 template <class E> Matrix_<E>
02500 operator*(const typename CNT<E>::StdNumber& l, const MatrixBase<E>& r) 
02501   { return Matrix_<E>(r)*=l; }
02502 
02503 template <class E> Matrix_<E>
02504 operator/(const MatrixBase<E>& l, const typename CNT<E>::StdNumber& r) 
02505   { return Matrix_<E>(l)/=r; }
02506 
02508 
02512 
02513 template <class E1, class E2>
02514 Vector_<typename CNT<E1>::template Result<E2>::Add>
02515 operator+(const VectorBase<E1>& l, const VectorBase<E2>& r) {
02516     return Vector_<typename CNT<E1>::template Result<E2>::Add>(l) += r;
02517 }
02518 template <class E>
02519 Vector_<E> operator+(const VectorBase<E>& l, const typename CNT<E>::T& r) {
02520     return Vector_<E>(l) += r;
02521 }
02522 template <class E>
02523 Vector_<E> operator+(const typename CNT<E>::T& l, const VectorBase<E>& r) {
02524     return Vector_<E>(r) += l;
02525 }
02526 template <class E1, class E2>
02527 Vector_<typename CNT<E1>::template Result<E2>::Sub>
02528 operator-(const VectorBase<E1>& l, const VectorBase<E2>& r) {
02529     return Vector_<typename CNT<E1>::template Result<E2>::Sub>(l) -= r;
02530 }
02531 template <class E>
02532 Vector_<E> operator-(const VectorBase<E>& l, const typename CNT<E>::T& r) {
02533     return Vector_<E>(l) -= r;
02534 }
02535 template <class E>
02536 Vector_<E> operator-(const typename CNT<E>::T& l, const VectorBase<E>& r) {
02537     Vector_<E> temp(r.size());
02538     temp = l;
02539     return (temp -= r);
02540 }
02541 
02542 template <class E> Vector_<E>
02543 operator*(const VectorBase<E>& l, const typename CNT<E>::StdNumber& r) 
02544   { return Vector_<E>(l)*=r; }
02545 
02546 template <class E> Vector_<E>
02547 operator*(const typename CNT<E>::StdNumber& l, const VectorBase<E>& r) 
02548   { return Vector_<E>(r)*=l; }
02549 
02550 template <class E> Vector_<E>
02551 operator/(const VectorBase<E>& l, const typename CNT<E>::StdNumber& r) 
02552   { return Vector_<E>(l)/=r; }
02553 
02555 
02559 
02560 template <class E1, class E2>
02561 RowVector_<typename CNT<E1>::template Result<E2>::Add>
02562 operator+(const RowVectorBase<E1>& l, const RowVectorBase<E2>& r) {
02563     return RowVector_<typename CNT<E1>::template Result<E2>::Add>(l) += r;
02564 }
02565 template <class E>
02566 RowVector_<E> operator+(const RowVectorBase<E>& l, const typename CNT<E>::T& r) {
02567     return RowVector_<E>(l) += r;
02568 }
02569 template <class E>
02570 RowVector_<E> operator+(const typename CNT<E>::T& l, const RowVectorBase<E>& r) {
02571     return RowVector_<E>(r) += l;
02572 }
02573 template <class E1, class E2>
02574 RowVector_<typename CNT<E1>::template Result<E2>::Sub>
02575 operator-(const RowVectorBase<E1>& l, const RowVectorBase<E2>& r) {
02576     return RowVector_<typename CNT<E1>::template Result<E2>::Sub>(l) -= r;
02577 }
02578 template <class E>
02579 RowVector_<E> operator-(const RowVectorBase<E>& l, const typename CNT<E>::T& r) {
02580     return RowVector_<E>(l) -= r;
02581 }
02582 template <class E>
02583 RowVector_<E> operator-(const typename CNT<E>::T& l, const RowVectorBase<E>& r) {
02584     RowVector_<E> temp(r.size());
02585     temp = l;
02586     return (temp -= r);
02587 }
02588 
02589 template <class E> RowVector_<E>
02590 operator*(const RowVectorBase<E>& l, const typename CNT<E>::StdNumber& r) 
02591   { return RowVector_<E>(l)*=r; }
02592 
02593 template <class E> RowVector_<E>
02594 operator*(const typename CNT<E>::StdNumber& l, const RowVectorBase<E>& r) 
02595   { return RowVector_<E>(r)*=l; }
02596 
02597 template <class E> RowVector_<E>
02598 operator/(const RowVectorBase<E>& l, const typename CNT<E>::StdNumber& r) 
02599   { return RowVector_<E>(l)/=r; }
02600 
02602 
02603 
02608 
02609     // TODO: these should use LAPACK!
02610 
02611 // Dot product
02612 template <class E1, class E2> 
02613 typename CNT<E1>::template Result<E2>::Mul
02614 operator*(const RowVectorBase<E1>& r, const VectorBase<E2>& v) {
02615     assert(r.ncol() == v.nrow());
02616     typename CNT<E1>::template Result<E2>::Mul sum(0);
02617     for (int j=0; j < r.ncol(); ++j)
02618         sum += r(j) * v[j];
02619     return sum;
02620 }
02621 
02622 template <class E1, class E2> 
02623 Vector_<typename CNT<E1>::template Result<E2>::Mul>
02624 operator*(const MatrixBase<E1>& m, const VectorBase<E2>& v) {
02625     assert(m.ncol() == v.nrow());
02626     Vector_<typename CNT<E1>::template Result<E2>::Mul> res(m.nrow());
02627     for (int i=0; i< m.nrow(); ++i)
02628         res[i] = m[i]*v;
02629     return res;
02630 }
02631 
02632 template <class E1, class E2> 
02633 Matrix_<typename CNT<E1>::template Result<E2>::Mul>
02634 operator*(const MatrixBase<E1>& m1, const MatrixBase<E2>& m2) {
02635     assert(m1.ncol() == m2.nrow());
02636     Matrix_<typename CNT<E1>::template Result<E2>::Mul> 
02637         res(m1.nrow(),m2.ncol());
02638 
02639     for (int j=0; j < res.ncol(); ++j)
02640         for (int i=0; i < res.nrow(); ++i)
02641             res(i,j) = m1[i] * m2(j);
02642 
02643     return res;
02644 }
02645 
02647 
02648     // GLOBAL OPERATORS: I/O
02649 
02650 template <class T> inline std::ostream&
02651 operator<<(std::ostream& o, const VectorBase<T>& v)
02652 { o << "~["; for(int i=0;i<v.size();++i) o<<(i>0?" ":"")<<v[i]; 
02653     return o << "]"; }
02654 
02655 template <class T> inline std::ostream&
02656 operator<<(std::ostream& o, const RowVectorBase<T>& v)
02657 { o << "["; for(int i=0;i<v.size();++i) o<<(i>0?" ":"")<<v[i]; 
02658     return o << "]"; }
02659 
02660 template <class T> inline std::ostream&
02661 operator<<(std::ostream& o, const MatrixBase<T>& m) {
02662     for (int i=0;i<m.nrow();++i)
02663         o << std::endl << m[i];
02664     if (m.nrow()) o << std::endl;
02665     return o; 
02666 }
02667 
02668 
02669 // Friendly abbreviations for default precision vectors and matrices.
02670 
02671 typedef Vector_<Real>           Vector;
02672 typedef Vector_<Complex>        ComplexVector;
02673 
02674 typedef VectorView_<Real>       VectorView;
02675 typedef VectorView_<Complex>    ComplexVectorView;
02676 
02677 typedef RowVector_<Real>        RowVector;
02678 typedef RowVector_<Complex>     ComplexRowVector;
02679 
02680 typedef RowVectorView_<Real>    RowVectorView;
02681 typedef RowVectorView_<Complex> ComplexRowVectorView;
02682 
02683 typedef Matrix_<Real>           Matrix;
02684 typedef Matrix_<Complex>        ComplexMatrix;
02685 
02686 typedef MatrixView_<Real>       MatrixView;
02687 typedef MatrixView_<Complex>    ComplexMatrixView;
02688 
02689 
02694 template <class ELT, class VECTOR_CLASS>
02695 class VectorIterator {
02696 public:
02697     typedef ELT value_type;
02698     typedef int difference_type;
02699     typedef ELT& reference;
02700     typedef ELT* pointer;
02701     typedef std::random_access_iterator_tag iterator_category;
02702     VectorIterator(VECTOR_CLASS& vector, int index) : vector(vector), index(index) {
02703     }
02704     VectorIterator(const VectorIterator& iter) : vector(iter.vector), index(iter.index) {
02705     }
02706     VectorIterator& operator=(const VectorIterator& iter) {
02707         vector = iter.vector;
02708         index = iter.index;
02709         return *this;
02710     }
02711     ELT& operator*() {
02712         assert (index >= 0 && index < vector.size());
02713         return vector[index];
02714     }
02715     ELT& operator[](int i) {
02716         assert (i >= 0 && i < vector.size());
02717         return vector[i];
02718     }
02719     VectorIterator operator++() {
02720         assert (index < vector.size());
02721         ++index;
02722         return *this;
02723     }
02724     VectorIterator operator++(int) {
02725         assert (index < vector.size());
02726         VectorIterator current = *this;
02727         ++index;
02728         return current;
02729     }
02730     VectorIterator operator--() {
02731         assert (index > 0);
02732         --index;
02733         return *this;
02734     }
02735     VectorIterator operator--(int) {
02736         assert (index > 0);
02737         VectorIterator current = *this;
02738         --index;
02739         return current;
02740     }
02741     bool operator<(VectorIterator iter) const {
02742         return (index < iter.index);
02743     }
02744     bool operator>(VectorIterator iter) const {
02745         return (index > iter.index);
02746     }
02747     bool operator<=(VectorIterator iter) const {
02748         return (index <= iter.index);
02749     }
02750     bool operator>=(VectorIterator iter) const {
02751         return (index >= iter.index);
02752     }
02753     int operator-(VectorIterator iter) const {
02754         return (index - iter.index);
02755     }
02756     VectorIterator operator-(int n) const {
02757         return VectorIterator(vector, index-n);
02758     }
02759     VectorIterator operator+(int n) const {
02760         return VectorIterator(vector, index+n);
02761     }
02762     bool operator==(VectorIterator iter) const {
02763         return (index == iter.index);
02764     }
02765     bool operator!=(VectorIterator iter) const {
02766         return (index != iter.index);
02767     }
02768 private:
02769     VECTOR_CLASS& vector;
02770     int index;
02771 };
02772 
02773 } //namespace SimTK
02774 
02775 #endif //SimTK_SIMMATRIX_BIGMATRIX_H_

Generated by  doxygen 1.6.2