00001 #ifndef SimTK_SIMMATRIX_BIGMATRIX_H_
00002 #define SimTK_SIMMATRIX_BIGMATRIX_H_
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00167
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
00226
00227 template <class ELT> class MatrixBase {
00228 public:
00229
00230
00231
00232
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;
00257 typedef ENumber Number;
00258 typedef EStdNumber StdNumber;
00259 typedef EPrecision Precision;
00260 typedef EScalarNormSq ScalarNormSq;
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;
00276 typedef MatrixBase<ESqTHerm> TSqTHerm;
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
00287
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
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)
00416 : helper(NScalarsPerElement, CppNScalarsPerElement,
00417 commitment, character, spacing, data) {}
00418
00420 MatrixBase(const MatrixCommitment& commitment,
00421 const MatrixCharacter& character,
00422 int spacing, Scalar* data)
00423 : helper(NScalarsPerElement, CppNScalarsPerElement,
00424 commitment, character, spacing, data) {}
00426
00427
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);
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
00740 inline RowVectorView_<ELT> row(int i) const;
00741 inline RowVectorView_<ELT> updRow(int i);
00742 inline VectorView_<ELT> col(int j) const;
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
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
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
00767 inline VectorView_<ELT> diag() const;
00768 inline VectorView_<ELT> updDiag();
00769
00770
00771
00772
00773
00774
00775
00776
00777
00778
00779
00780
00781 TInvert invert() const {
00782 TInvert m(*this);
00783 m.helper.invertInPlace();
00784 return m;
00785 }
00786
00787 void invertInPlace() {helper.invertInPlace();}
00788
00790 void dump(const char* msg=0) const {
00791 helper.dump(msg);
00792 }
00793
00794
00795
00796
00797
00798
00799
00800
00801
00802
00803
00804
00805
00806
00807
00808
00809
00810
00811 template <class ELT_A, class ELT_B>
00812 MatrixBase& matmul(const StdNumber& beta,
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
00844
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
00858
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
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
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
00940
00941 void lockNRows() {helper.lockNRows();}
00942 void lockNCols() {helper.lockNCols();}
00943 void lockShape() {helper.lockShape();}
00944
00945
00946
00947
00948
00949 void unlockNRows() {helper.unlockNRows();}
00950 void unlockNCols() {helper.unlockNCols();}
00951 void unlockShape() {helper.unlockShape();}
00952
00953
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
00986
00987
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;
01031
01032 template <class EE> friend class MatrixBase;
01033 };
01034
01035
01036
01037
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
01126
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
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
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
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
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
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
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
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
01266 TAbs abs() const {TAbs result; Base::abs(result); return result;}
01267
01268
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
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
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
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
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
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; }
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
01324 explicit VectorBase(MatrixHelperRep<Scalar>* hrep) : Base(hrep) {}
01325
01326 private:
01327
01328 };
01329
01330
01331
01332
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
01421
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
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
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
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
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
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
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
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
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
01543 TAbs abs() const {
01544 TAbs result; Base::abs(result); return result;
01545 }
01546
01547
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
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
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
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
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
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; }
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
01603 explicit RowVectorBase(MatrixHelperRep<Scalar>* hrep) : Base(hrep) {}
01604
01605 private:
01606
01607 };
01608
01609
01610
01611
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
01624
01625
01626
01627 explicit MatrixView_(MatrixHelperRep<S>* hrep) : Base(hrep) {}
01628
01629
01630
01631
01632 MatrixView_(const MatrixView_& m)
01633 : Base(MatrixCommitment(),
01634 const_cast<MatrixHelper<S>&>(m.getHelper()),
01635 typename MatrixHelper<S>::ShallowCopy()) {}
01636
01637
01638 MatrixView_& operator=(const MatrixView_& m) {
01639 Base::operator=(m); return *this;
01640 }
01641
01642
01643 MatrixView_(DeadMatrixView_<ELT>&);
01644 MatrixView_& operator=(DeadMatrixView_<ELT>&);
01645
01646
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
01670 MatrixView_();
01671 };
01672
01673
01674
01675
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
01686
01687
01688
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
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
01716 DeadMatrixView_();
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
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
01758 Matrix_(const Matrix_& src) : Base(src) { }
01759
01760
01761
01762 Matrix_& operator=(const Matrix_& src) {
01763 Base::operator=(src); return *this;
01764 }
01765
01766
01767
01768 Matrix_(const Base& v) : Base(v) {}
01769
01770
01771
01772 Matrix_(const BaseNeg& v) : Base(v) {}
01773
01774
01775
01776
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)
01786 : Base(MatrixCommitment(), MatrixCharacter::LapackFull(m,n),
01787 leadingDim, data) {}
01788 Matrix_(int m, int n, int leadingDim, S* data)
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
01822 };
01823
01824
01825
01826
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
01843
01844
01845
01846 explicit VectorView_(MatrixHelperRep<S>* hrep) : Base(hrep) {}
01847
01848
01849
01850
01851 VectorView_(const VectorView_& v)
01852 : Base(const_cast<MatrixHelper<S>&>(v.getHelper()), typename MatrixHelper<S>::ShallowCopy()) { }
01853
01854
01855 VectorView_& operator=(const VectorView_& v) {
01856 Base::operator=(v); return *this;
01857 }
01858
01859
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
01881 VectorView_();
01882 };
01883
01884
01885
01886
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() {}
01900
01901
01902
01903 Vector_(const Vector_& src) : Base(src) {}
01904
01905
01906 Vector_(const Base& src) : Base(src) {}
01907 Vector_(const BaseNeg& src) : Base(src) {}
01908
01909
01910
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
01956 };
01957
01958
01959
01960
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
01977
01978
01979
01980 explicit RowVectorView_(MatrixHelperRep<S>* hrep) : Base(hrep) {}
01981
01982
01983
01984
01985 RowVectorView_(const RowVectorView_& r)
01986 : Base(const_cast<MatrixHelper<S>&>(r.getHelper()), typename MatrixHelper<S>::ShallowCopy()) { }
01987
01988
01989 RowVectorView_& operator=(const RowVectorView_& r) {
01990 Base::operator=(r); return *this;
01991 }
01992
01993
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
02015 RowVectorView_();
02016 };
02017
02018
02019
02020
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() {}
02036
02037
02038
02039 RowVector_(const RowVector_& src) : Base(src) {}
02040
02041
02042 RowVector_(const Base& src) : Base(src) {}
02043 RowVector_(const BaseNeg& src) : Base(src) {}
02044
02045
02046
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
02092 };
02093
02094
02095
02096
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
02187
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
02206
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
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
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
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
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
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
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
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
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
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
02438
02439
02440
02441
02442
02443
02444
02445
02446
02447
02448
02449
02450
02454
02455
02456
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
02492
02493
02494
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
02610
02611
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
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
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 }
02774
02775 #endif //SimTK_SIMMATRIX_BIGMATRIX_H_