Simbody  3.5
BigMatrix.h
Go to the documentation of this file.
00001 #ifndef SimTK_SIMMATRIX_BIGMATRIX_H_
00002 #define SimTK_SIMMATRIX_BIGMATRIX_H_
00003 
00004 /* -------------------------------------------------------------------------- *
00005  *                       Simbody(tm): SimTKcommon                             *
00006  * -------------------------------------------------------------------------- *
00007  * This is part of the SimTK biosimulation toolkit originating from           *
00008  * Simbios, the NIH National Center for Physics-Based Simulation of           *
00009  * Biological Structures at Stanford, funded under the NIH Roadmap for        *
00010  * Medical Research, grant U54 GM072970. See https://simtk.org/home/simbody.  *
00011  *                                                                            *
00012  * Portions copyright (c) 2005-13 Stanford University and the Authors.        *
00013  * Authors: Michael Sherman                                                   *
00014  * Contributors:                                                              *
00015  *                                                                            *
00016  * Licensed under the Apache License, Version 2.0 (the "License"); you may    *
00017  * not use this file except in compliance with the License. You may obtain a  *
00018  * copy of the License at http://www.apache.org/licenses/LICENSE-2.0.         *
00019  *                                                                            *
00020  * Unless required by applicable law or agreed to in writing, software        *
00021  * distributed under the License is distributed on an "AS IS" BASIS,          *
00022  * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.   *
00023  * See the License for the specific language governing permissions and        *
00024  * limitations under the License.                                             *
00025  * -------------------------------------------------------------------------- */
00026 
00150 #include "SimTKcommon/Scalar.h"
00151 #include "SimTKcommon/SmallMatrix.h"
00152 
00153 #include "SimTKcommon/internal/MatrixHelper.h"
00154 #include "SimTKcommon/internal/MatrixCharacteristics.h"
00155 
00156 #include <iostream>
00157 #include <cassert>
00158 #include <complex>
00159 #include <cstddef>
00160 #include <limits>
00161 
00162 namespace SimTK {
00163     template <class ELT>    class MatrixBase;
00164     template <class ELT>    class VectorBase;
00165     template <class ELT>    class RowVectorBase;
00166 
00167     template <class ELT = Real> class MatrixView_;
00168     template <class ELT = Real> class Matrix_;
00169 
00170     template <class ELT = Real> class VectorView_;
00171     template <class ELT = Real> class Vector_;
00172 
00173     template <class ELT = Real> class RowVectorView_;
00174     template <class ELT = Real> class RowVector_;
00175 
00176     template <class ELT, class VECTOR_CLASS> class VectorIterator;
00177 }
00178 
00179 #include "SimTKcommon/internal/MatrixBase.h"
00180 #include "SimTKcommon/internal/VectorBase.h"
00181 #include "SimTKcommon/internal/RowVectorBase.h"
00182 
00183 #include "SimTKcommon/internal/MatrixView_.h"
00184 #include "SimTKcommon/internal/Matrix_.h"
00185 
00186 #include "SimTKcommon/internal/VectorView_.h"
00187 #include "SimTKcommon/internal/Vector_.h"
00188 
00189 #include "SimTKcommon/internal/RowVectorView_.h"
00190 #include "SimTKcommon/internal/RowVector_.h"
00191 
00192 #include "SimTKcommon/internal/VectorIterator.h"
00193 
00194 
00195 namespace SimTK {
00196 
00197 //  ------------------------ MatrixBase definitions ----------------------------
00198 
00199 template <class ELT> inline MatrixView_<ELT> 
00200 MatrixBase<ELT>::block(int i, int j, int m, int n) const { 
00201     SimTK_INDEXCHECK(i,nrow()+1,"MatrixBase::block()");
00202     SimTK_INDEXCHECK(j,ncol()+1,"MatrixBase::block()");
00203     SimTK_SIZECHECK(i+m,nrow(),"MatrixBase::block()");
00204     SimTK_SIZECHECK(j+n,ncol(),"MatrixBase::block()");
00205 
00206     MatrixHelper<Scalar> h(MatrixCommitment(),helper,i,j,m,n);    
00207     return MatrixView_<ELT>(h.stealRep()); 
00208 }
00209     
00210 template <class ELT> inline MatrixView_<ELT>
00211 MatrixBase<ELT>::updBlock(int i, int j, int m, int n) { 
00212     SimTK_INDEXCHECK(i,nrow()+1,"MatrixBase::updBlock()");
00213     SimTK_INDEXCHECK(j,ncol()+1,"MatrixBase::updBlock()");
00214     SimTK_SIZECHECK(i+m,nrow(),"MatrixBase::updBlock()");
00215     SimTK_SIZECHECK(j+n,ncol(),"MatrixBase::updBlock()");
00216 
00217     MatrixHelper<Scalar> h(MatrixCommitment(),helper,i,j,m,n);        
00218     return MatrixView_<ELT>(h.stealRep()); 
00219 }
00220 
00221 template <class E> inline MatrixView_<typename CNT<E>::THerm>
00222 MatrixBase<E>::transpose() const { 
00223     MatrixHelper<typename CNT<Scalar>::THerm> 
00224         h(MatrixCommitment(),
00225           helper, typename MatrixHelper<typename CNT<Scalar>::THerm>::TransposeView());
00226     return MatrixView_<typename CNT<E>::THerm>(h.stealRep()); 
00227 }
00228     
00229 template <class E> inline MatrixView_<typename CNT<E>::THerm>
00230 MatrixBase<E>::updTranspose() {     
00231     MatrixHelper<typename CNT<Scalar>::THerm> 
00232         h(MatrixCommitment(),
00233           helper, typename MatrixHelper<typename CNT<Scalar>::THerm>::TransposeView());
00234     return MatrixView_<typename CNT<E>::THerm>(h.stealRep()); 
00235 }
00236 
00237 template <class E> inline VectorView_<E>
00238 MatrixBase<E>::diag() const { 
00239     MatrixHelper<Scalar> h(MatrixCommitment::Vector(),
00240                            helper, typename MatrixHelper<Scalar>::DiagonalView());
00241     return VectorView_<E>(h.stealRep()); 
00242 }
00243     
00244 template <class E> inline VectorView_<E>
00245 MatrixBase<E>::updDiag() {     
00246     MatrixHelper<Scalar> h(MatrixCommitment::Vector(),
00247                            helper, typename MatrixHelper<Scalar>::DiagonalView());
00248     return VectorView_<E>(h.stealRep());
00249 }
00250 
00251 template <class ELT> inline VectorView_<ELT> 
00252 MatrixBase<ELT>::col(int j) const { 
00253     SimTK_INDEXCHECK(j,ncol(),"MatrixBase::col()");
00254 
00255     MatrixHelper<Scalar> h(MatrixCommitment::Vector(),
00256                            helper,0,j,nrow(),1);    
00257     return VectorView_<ELT>(h.stealRep()); 
00258 }
00259     
00260 template <class ELT> inline VectorView_<ELT>
00261 MatrixBase<ELT>::updCol(int j) {
00262     SimTK_INDEXCHECK(j,ncol(),"MatrixBase::updCol()");
00263 
00264     MatrixHelper<Scalar> h(MatrixCommitment::Vector(),
00265                            helper,0,j,nrow(),1);        
00266     return VectorView_<ELT>(h.stealRep()); 
00267 }
00268 
00269 template <class ELT> inline RowVectorView_<ELT> 
00270 MatrixBase<ELT>::row(int i) const { 
00271     SimTK_INDEXCHECK(i,nrow(),"MatrixBase::row()");
00272 
00273     MatrixHelper<Scalar> h(MatrixCommitment::RowVector(),
00274                            helper,i,0,1,ncol());    
00275     return RowVectorView_<ELT>(h.stealRep()); 
00276 }
00277     
00278 template <class ELT> inline RowVectorView_<ELT>
00279 MatrixBase<ELT>::updRow(int i) { 
00280     SimTK_INDEXCHECK(i,nrow(),"MatrixBase::updRow()");
00281 
00282     MatrixHelper<Scalar> h(MatrixCommitment::RowVector(),
00283                            helper,i,0,1,ncol());        
00284     return RowVectorView_<ELT>(h.stealRep()); 
00285 }
00286 
00287 // M = diag(v) * M; v must have nrow() elements.
00288 // That is, M[i] *= v[i].
00289 template <class ELT> template <class EE> inline MatrixBase<ELT>& 
00290 MatrixBase<ELT>::rowScaleInPlace(const VectorBase<EE>& v) {
00291     assert(v.nrow() == nrow());
00292     for (int i=0; i < nrow(); ++i)
00293         (*this)[i] *= v[i];
00294     return *this;
00295 }
00296 
00297 template <class ELT> template <class EE> inline void
00298 MatrixBase<ELT>::rowScale(const VectorBase<EE>& v, typename MatrixBase<ELT>::template EltResult<EE>::Mul& out) const {
00299     assert(v.nrow() == nrow());
00300     out.resize(nrow(), ncol());
00301     for (int j=0; j<ncol(); ++j)
00302         for (int i=0; i<nrow(); ++i)
00303            out(i,j) = (*this)(i,j) * v[i];
00304 }
00305 
00306 // M = M * diag(v); v must have ncol() elements
00307 // That is, M(i) *= v[i]
00308 template <class ELT> template <class EE>  inline MatrixBase<ELT>& 
00309 MatrixBase<ELT>::colScaleInPlace(const VectorBase<EE>& v) {
00310     assert(v.nrow() == ncol());
00311     for (int j=0; j < ncol(); ++j)
00312         (*this)(j) *= v[j];
00313     return *this;
00314 }
00315 
00316 template <class ELT> template <class EE> inline void
00317 MatrixBase<ELT>::colScale(const VectorBase<EE>& v, typename MatrixBase<ELT>::template EltResult<EE>::Mul& out) const {
00318     assert(v.nrow() == ncol());
00319     out.resize(nrow(), ncol());
00320     for (int j=0; j<ncol(); ++j)
00321         for (int i=0; i<nrow(); ++i)
00322            out(i,j) = (*this)(i,j) * v[j];
00323 }
00324 
00325 
00326 // M(i,j) *= r[i]*c[j]; r must have nrow() elements; c must have ncol() elements
00327 template <class ELT> template <class ER, class EC> inline MatrixBase<ELT>& 
00328 MatrixBase<ELT>::rowAndColScaleInPlace(const VectorBase<ER>& r, const VectorBase<EC>& c) {
00329     assert(r.nrow()==nrow() && c.nrow()==ncol());
00330     for (int j=0; j<ncol(); ++j)
00331         for (int i=0; i<nrow(); ++i)
00332             (*this)(i,j) *= (r[i]*c[j]);
00333     return *this;
00334 }
00335 
00336 template <class ELT> template <class ER, class EC> inline void
00337 MatrixBase<ELT>::rowAndColScale(
00338     const VectorBase<ER>& r, 
00339     const VectorBase<EC>& c,
00340     typename EltResult<typename VectorBase<ER>::template EltResult<EC>::Mul>::Mul& 
00341                           out) const
00342 {
00343     assert(r.nrow()==nrow() && c.nrow()==ncol());
00344     out.resize(nrow(), ncol());
00345     for (int j=0; j<ncol(); ++j)
00346         for (int i=0; i<nrow(); ++i)
00347             out(i,j) = (*this)(i,j) * (r[i]*c[j]);
00348 }
00349 
00350 // M(i,j) = s
00351 template <class ELT> template <class S> inline MatrixBase<ELT>& 
00352 MatrixBase<ELT>::elementwiseAssign(const S& s) {
00353     for (int j=0; j<ncol(); ++j)
00354     for (int i=0; i<nrow(); ++i)
00355         (*this)(i,j) = s;
00356     return *this;
00357 }
00358 
00359 // Set M(i,j) = M(i,j)^-1.
00360 template <class ELT> inline MatrixBase<ELT>& 
00361 MatrixBase<ELT>::elementwiseInvertInPlace() {
00362     const int nr=nrow(), nc=ncol();
00363     for (int j=0; j<nc; ++j)
00364         for (int i=0; i<nr; ++i) {
00365             ELT& e = updElt(i,j);
00366             e = CNT<ELT>::invert(e);
00367         }
00368     return *this;
00369 }
00370 
00371 template <class ELT> inline void
00372 MatrixBase<ELT>::elementwiseInvert(MatrixBase<typename CNT<E>::TInvert>& out) const {
00373     const int nr=nrow(), nc=ncol();
00374     out.resize(nr,nc);
00375     for (int j=0; j<nc; ++j)
00376         for (int i=0; i<nr; ++i)
00377             out(i,j) = CNT<ELT>::invert((*this)(i,j));
00378 }
00379 
00380 // M(i,j) += s
00381 template <class ELT> template <class S> inline MatrixBase<ELT>& 
00382 MatrixBase<ELT>::elementwiseAddScalarInPlace(const S& s) {
00383     for (int j=0; j<ncol(); ++j)
00384         for (int i=0; i<nrow(); ++i)
00385             (*this)(i,j) += s;
00386     return *this;
00387 }
00388 
00389 template <class ELT> template <class S> inline void 
00390 MatrixBase<ELT>::elementwiseAddScalar(
00391     const S& s,
00392     typename MatrixBase<ELT>::template EltResult<S>::Add& out) const
00393 {
00394     const int nr=nrow(), nc=ncol();
00395     out.resize(nr,nc);
00396     for (int j=0; j<nc; ++j)
00397         for (int i=0; i<nr; ++i)
00398             out(i,j) = (*this)(i,j) + s;
00399 }
00400 
00401 // M(i,j) -= s
00402 template <class ELT> template <class S> inline MatrixBase<ELT>& 
00403 MatrixBase<ELT>::elementwiseSubtractScalarInPlace(const S& s) {
00404     for (int j=0; j<ncol(); ++j)
00405         for (int i=0; i<nrow(); ++i)
00406             (*this)(i,j) -= s;
00407     return *this;
00408 }
00409 
00410 template <class ELT> template <class S> inline void 
00411 MatrixBase<ELT>::elementwiseSubtractScalar(
00412     const S& s,
00413     typename MatrixBase<ELT>::template EltResult<S>::Sub& out) const
00414 {
00415     const int nr=nrow(), nc=ncol();
00416     out.resize(nr,nc);
00417     for (int j=0; j<nc; ++j)
00418         for (int i=0; i<nr; ++i)
00419             out(i,j) = (*this)(i,j) - s;
00420 }
00421 
00422 // M(i,j) = s - M(i,j)
00423 template <class ELT> template <class S> inline MatrixBase<ELT>& 
00424 MatrixBase<ELT>::elementwiseSubtractFromScalarInPlace(const S& s) {
00425     const int nr=nrow(), nc=ncol();
00426     for (int j=0; j<nc; ++j)
00427         for (int i=0; i<nr; ++i) {
00428             ELT& e = updElt(i,j);
00429             e = s - e;
00430         }
00431     return *this;
00432 }
00433 
00434 template <class ELT> template <class S> inline void 
00435 MatrixBase<ELT>::elementwiseSubtractFromScalar(
00436     const S& s,
00437     typename MatrixBase<S>::template EltResult<ELT>::Sub& out) const
00438 {
00439     const int nr=nrow(), nc=ncol();
00440     out.resize(nr,nc);
00441     for (int j=0; j<nc; ++j)
00442         for (int i=0; i<nr; ++i)
00443             out(i,j) = s - (*this)(i,j);
00444 }
00445 
00446 // M(i,j) *= R(i,j); R must have same dimensions as this
00447 template <class ELT> template <class EE> inline MatrixBase<ELT>& 
00448 MatrixBase<ELT>::elementwiseMultiplyInPlace(const MatrixBase<EE>& r) {
00449     const int nr=nrow(), nc=ncol();
00450     assert(r.nrow()==nr && r.ncol()==nc);
00451     for (int j=0; j<nc; ++j)
00452         for (int i=0; i<nr; ++i)
00453             (*this)(i,j) *= r(i,j);
00454     return *this;
00455 }
00456 
00457 template <class ELT> template <class EE> inline void 
00458 MatrixBase<ELT>::elementwiseMultiply(
00459     const MatrixBase<EE>& r, 
00460     typename MatrixBase<ELT>::template EltResult<EE>::Mul& out) const
00461 {
00462     const int nr=nrow(), nc=ncol();
00463     assert(r.nrow()==nr && r.ncol()==nc);
00464     out.resize(nr,nc);
00465     for (int j=0; j<nc; ++j)
00466         for (int i=0; i<nr; ++i)
00467             out(i,j) = (*this)(i,j) * r(i,j);
00468 }
00469 
00470 // M(i,j) = R(i,j) * M(i,j); R must have same dimensions as this
00471 template <class ELT> template <class EE> inline MatrixBase<ELT>& 
00472 MatrixBase<ELT>::elementwiseMultiplyFromLeftInPlace(const MatrixBase<EE>& r) {
00473     const int nr=nrow(), nc=ncol();
00474     assert(r.nrow()==nr && r.ncol()==nc);
00475     for (int j=0; j<nc; ++j)
00476         for (int i=0; i<nr; ++i) {
00477             ELT& e = updElt(i,j);
00478             e = r(i,j) * e;
00479         }
00480     return *this;
00481 }
00482 
00483 template <class ELT> template <class EE> inline void 
00484 MatrixBase<ELT>::elementwiseMultiplyFromLeft(
00485     const MatrixBase<EE>& r, 
00486     typename MatrixBase<EE>::template EltResult<ELT>::Mul& out) const
00487 {
00488     const int nr=nrow(), nc=ncol();
00489     assert(r.nrow()==nr && r.ncol()==nc);
00490     out.resize(nr,nc);
00491     for (int j=0; j<nc; ++j)
00492         for (int i=0; i<nr; ++i)
00493             out(i,j) =  r(i,j) * (*this)(i,j);
00494 }
00495 
00496 // M(i,j) /= R(i,j); R must have same dimensions as this
00497 template <class ELT> template <class EE> inline MatrixBase<ELT>& 
00498 MatrixBase<ELT>::elementwiseDivideInPlace(const MatrixBase<EE>& r) {
00499     const int nr=nrow(), nc=ncol();
00500     assert(r.nrow()==nr && r.ncol()==nc);
00501     for (int j=0; j<nc; ++j)
00502         for (int i=0; i<nr; ++i)
00503             (*this)(i,j) /= r(i,j);
00504     return *this;
00505 }
00506 
00507 template <class ELT> template <class EE> inline void 
00508 MatrixBase<ELT>::elementwiseDivide(
00509     const MatrixBase<EE>& r,
00510     typename MatrixBase<ELT>::template EltResult<EE>::Dvd& out) const
00511 {
00512     const int nr=nrow(), nc=ncol();
00513     assert(r.nrow()==nr && r.ncol()==nc);
00514     out.resize(nr,nc);
00515     for (int j=0; j<nc; ++j)
00516         for (int i=0; i<nr; ++i)
00517             out(i,j) = (*this)(i,j) / r(i,j);
00518 }
00519 // M(i,j) = R(i,j) / M(i,j); R must have same dimensions as this
00520 template <class ELT> template <class EE> inline MatrixBase<ELT>& 
00521 MatrixBase<ELT>::elementwiseDivideFromLeftInPlace(const MatrixBase<EE>& r) {
00522     const int nr=nrow(), nc=ncol();
00523     assert(r.nrow()==nr && r.ncol()==nc);
00524     for (int j=0; j<nc; ++j)
00525         for (int i=0; i<nr; ++i) {
00526             ELT& e = updElt(i,j);
00527             e = r(i,j) / e;
00528         }
00529     return *this;
00530 }
00531 
00532 template <class ELT> template <class EE> inline void 
00533 MatrixBase<ELT>::elementwiseDivideFromLeft(
00534     const MatrixBase<EE>& r,
00535     typename MatrixBase<EE>::template EltResult<ELT>::Dvd& out) const
00536 {
00537     const int nr=nrow(), nc=ncol();
00538     assert(r.nrow()==nr && r.ncol()==nc);
00539     out.resize(nr,nc);
00540     for (int j=0; j<nc; ++j)
00541         for (int i=0; i<nr; ++i)
00542             out(i,j) = r(i,j) / (*this)(i,j);
00543 }
00544 
00545 /*
00546 template <class ELT> inline MatrixView_< typename CNT<ELT>::TReal > 
00547 MatrixBase<ELT>::real() const { 
00548     if (!CNT<ELT>::IsComplex) { // known at compile time
00549         return MatrixView_< typename CNT<ELT>::TReal >( // this is just ELT
00550             MatrixHelper(helper,0,0,nrow(),ncol()));    // a view of the whole matrix
00551     }
00552     // Elements are complex -- helper uses underlying precision (real) type.
00553     MatrixHelper<Precision> h(helper,typename MatrixHelper<Precision>::RealView);    
00554     return MatrixView_< typename CNT<ELT>::TReal >(h); 
00555 }
00556 */
00557 
00558 
00559 //  ----------------------------------------------------------------------------
00563 
00564 // + and - allow mixed element types, but will fail to compile if the elements aren't
00565 // compatible. At run time these will fail if the dimensions are incompatible.
00566 template <class E1, class E2>
00567 Matrix_<typename CNT<E1>::template Result<E2>::Add>
00568 operator+(const MatrixBase<E1>& l, const MatrixBase<E2>& r) {
00569     return Matrix_<typename CNT<E1>::template Result<E2>::Add>(l) += r;
00570 }
00571 
00572 template <class E>
00573 Matrix_<E> operator+(const MatrixBase<E>& l, const typename CNT<E>::T& r) {
00574     return Matrix_<E>(l) += r;
00575 }
00576 
00577 template <class E>
00578 Matrix_<E> operator+(const typename CNT<E>::T& l, const MatrixBase<E>& r) {
00579     return Matrix_<E>(r) += l;
00580 }
00581 
00582 template <class E1, class E2>
00583 Matrix_<typename CNT<E1>::template Result<E2>::Sub>
00584 operator-(const MatrixBase<E1>& l, const MatrixBase<E2>& r) {
00585     return Matrix_<typename CNT<E1>::template Result<E2>::Sub>(l) -= r;
00586 }
00587 
00588 template <class E>
00589 Matrix_<E> operator-(const MatrixBase<E>& l, const typename CNT<E>::T& r) {
00590     return Matrix_<E>(l) -= r;
00591 }
00592 
00593 template <class E>
00594 Matrix_<E> operator-(const typename CNT<E>::T& l, const MatrixBase<E>& r) {
00595     Matrix_<E> temp(r.nrow(), r.ncol());
00596     temp = l;
00597     return (temp -= r);
00598 }
00599 
00600 // Scalar multiply and divide. You might wish the scalar could be
00601 // a templatized type "E2", but that would create horrible ambiguities since
00602 // E2 would match not only scalar types but everything else including
00603 // matrices.
00604 template <class E> Matrix_<E>
00605 operator*(const MatrixBase<E>& l, const typename CNT<E>::StdNumber& r) 
00606   { return Matrix_<E>(l)*=r; }
00607 
00608 template <class E> Matrix_<E>
00609 operator*(const typename CNT<E>::StdNumber& l, const MatrixBase<E>& r) 
00610   { return Matrix_<E>(r)*=l; }
00611 
00612 template <class E> Matrix_<E>
00613 operator/(const MatrixBase<E>& l, const typename CNT<E>::StdNumber& r) 
00614   { return Matrix_<E>(l)/=r; }
00615 
00616 // Handle ints explicitly.
00617 template <class E> Matrix_<E>
00618 operator*(const MatrixBase<E>& l, int r) 
00619   { return Matrix_<E>(l)*= typename CNT<E>::StdNumber(r); }
00620 
00621 template <class E> Matrix_<E>
00622 operator*(int l, const MatrixBase<E>& r) 
00623   { return Matrix_<E>(r)*= typename CNT<E>::StdNumber(l); }
00624 
00625 template <class E> Matrix_<E>
00626 operator/(const MatrixBase<E>& l, int r) 
00627   { return Matrix_<E>(l)/= typename CNT<E>::StdNumber(r); }
00628 
00630 
00634 
00635 template <class E1, class E2>
00636 Vector_<typename CNT<E1>::template Result<E2>::Add>
00637 operator+(const VectorBase<E1>& l, const VectorBase<E2>& r) {
00638     return Vector_<typename CNT<E1>::template Result<E2>::Add>(l) += r;
00639 }
00640 template <class E>
00641 Vector_<E> operator+(const VectorBase<E>& l, const typename CNT<E>::T& r) {
00642     return Vector_<E>(l) += r;
00643 }
00644 template <class E>
00645 Vector_<E> operator+(const typename CNT<E>::T& l, const VectorBase<E>& r) {
00646     return Vector_<E>(r) += l;
00647 }
00648 template <class E1, class E2>
00649 Vector_<typename CNT<E1>::template Result<E2>::Sub>
00650 operator-(const VectorBase<E1>& l, const VectorBase<E2>& r) {
00651     return Vector_<typename CNT<E1>::template Result<E2>::Sub>(l) -= r;
00652 }
00653 template <class E>
00654 Vector_<E> operator-(const VectorBase<E>& l, const typename CNT<E>::T& r) {
00655     return Vector_<E>(l) -= r;
00656 }
00657 template <class E>
00658 Vector_<E> operator-(const typename CNT<E>::T& l, const VectorBase<E>& r) {
00659     Vector_<E> temp(r.size());
00660     temp = l;
00661     return (temp -= r);
00662 }
00663 
00664 // Scalar multiply and divide.
00665 
00666 template <class E> Vector_<E>
00667 operator*(const VectorBase<E>& l, const typename CNT<E>::StdNumber& r) 
00668   { return Vector_<E>(l)*=r; }
00669 
00670 template <class E> Vector_<E>
00671 operator*(const typename CNT<E>::StdNumber& l, const VectorBase<E>& r) 
00672   { return Vector_<E>(r)*=l; }
00673 
00674 template <class E> Vector_<E>
00675 operator/(const VectorBase<E>& l, const typename CNT<E>::StdNumber& r) 
00676   { return Vector_<E>(l)/=r; }
00677 
00678 // Handle ints explicitly
00679 template <class E> Vector_<E>
00680 operator*(const VectorBase<E>& l, int r) 
00681   { return Vector_<E>(l)*= typename CNT<E>::StdNumber(r); }
00682 
00683 template <class E> Vector_<E>
00684 operator*(int l, const VectorBase<E>& r) 
00685   { return Vector_<E>(r)*= typename CNT<E>::StdNumber(l); }
00686 
00687 template <class E> Vector_<E>
00688 operator/(const VectorBase<E>& l, int r) 
00689   { return Vector_<E>(l)/= typename CNT<E>::StdNumber(r); }
00690 
00691 // These are fancier "scalars"; whether they are allowed depends on
00692 // whether the element type and the CNT are compatible.
00693 
00694 // Vector * Vec
00695 template <class E1, int M, class E2, int S> 
00696 Vector_<typename CNT<E1>::template Result< Vec<M,E2,S> >::Mul>
00697 operator*(const VectorBase<E1>& v, const Vec<M,E2,S>& s) {
00698     Vector_<typename CNT<E1>::template Result< Vec<M,E2,S> >::Mul> res(v.nrow());
00699     for (int i=0; i < v.nrow(); ++i)
00700         res[i] = v[i]*s; 
00701     return res;
00702 }
00703 
00704 // Vec * Vector
00705 template <class E1, int M, class E2, int S> 
00706 Vector_<typename Vec<M,E2,S>::template Result<E1>::Mul>
00707 operator*(const Vec<M,E2,S>& s, const VectorBase<E1>& v) {
00708     Vector_<typename Vec<M,E2,S>::template Result<E1>::Mul> res(v.nrow());
00709     for (int i=0; i < v.nrow(); ++i)
00710         res[i] = s*v[i]; 
00711     return res;
00712 }
00713 
00714 // Vector * Row
00715 template <class E1, int N, class E2, int S> 
00716 Vector_<typename CNT<E1>::template Result< Row<N,E2,S> >::Mul>
00717 operator*(const VectorBase<E1>& v, const Row<N,E2,S>& s) {
00718     Vector_<typename CNT<E1>::template Result< Row<N,E2,S> >::Mul> res(v.nrow());
00719     for (int i=0; i < v.nrow(); ++i)
00720         res[i] = v[i]*s; 
00721     return res;
00722 }
00723 
00724 // Row * Vector
00725 template <class E1, int N, class E2, int S> 
00726 Vector_<typename Row<N,E2,S>::template Result<E1>::Mul>
00727 operator*(const Row<N,E2,S>& s, const VectorBase<E1>& v) {
00728     Vector_<typename Row<N,E2,S>::template Result<E1>::Mul> res(v.nrow());
00729     for (int i=0; i < v.nrow(); ++i)
00730         res[i] = s*v[i]; 
00731     return res;
00732 }
00733 
00734 // Vector * Mat
00735 template <class E1, int M, int N, class E2, int S1, int S2> 
00736 Vector_<typename CNT<E1>::template Result< Mat<M,N,E2,S1,S2> >::Mul>
00737 operator*(const VectorBase<E1>& v, const Mat<M,N,E2,S1,S2>& s) {
00738     Vector_<typename CNT<E1>::template Result< Mat<M,N,E2,S1,S2> >::Mul> res(v.nrow());
00739     for (int i=0; i < v.nrow(); ++i)
00740         res[i] = v[i]*s; 
00741     return res;
00742 }
00743 
00744 // Mat * Vector
00745 template <class E1, int M, int N, class E2, int S1, int S2> 
00746 Vector_<typename Mat<M,N,E2,S1,S2>::template Result<E1>::Mul>
00747 operator*(const Mat<M,N,E2,S1,S2>& s, const VectorBase<E1>& v) {
00748     Vector_<typename Mat<M,N,E2,S1,S2>::template Result<E1>::Mul> res(v.nrow());
00749     for (int i=0; i < v.nrow(); ++i)
00750         res[i] = s*v[i]; 
00751     return res;
00752 }
00753 
00754 // Vector * SymMat
00755 template <class E1, int M, class E2, int S> 
00756 Vector_<typename CNT<E1>::template Result< SymMat<M,E2,S> >::Mul>
00757 operator*(const VectorBase<E1>& v, const SymMat<M,E2,S>& s) {
00758     Vector_<typename CNT<E1>::template Result< SymMat<M,E2,S> >::Mul> res(v.nrow());
00759     for (int i=0; i < v.nrow(); ++i)
00760         res[i] = v[i]*s; 
00761     return res;
00762 }
00763 
00764 // SymMat * Vector
00765 template <class E1, int M, class E2, int S> 
00766 Vector_<typename SymMat<M,E2,S>::template Result<E1>::Mul>
00767 operator*(const SymMat<M,E2,S>& s, const VectorBase<E1>& v) {
00768     Vector_<typename SymMat<M,E2,S>::template Result<E1>::Mul> res(v.nrow());
00769     for (int i=0; i < v.nrow(); ++i)
00770         res[i] = s*v[i]; 
00771     return res;
00772 }
00773 
00775 
00779 
00780 template <class E1, class E2>
00781 RowVector_<typename CNT<E1>::template Result<E2>::Add>
00782 operator+(const RowVectorBase<E1>& l, const RowVectorBase<E2>& r) {
00783     return RowVector_<typename CNT<E1>::template Result<E2>::Add>(l) += r;
00784 }
00785 template <class E>
00786 RowVector_<E> operator+(const RowVectorBase<E>& l, const typename CNT<E>::T& r) {
00787     return RowVector_<E>(l) += r;
00788 }
00789 template <class E>
00790 RowVector_<E> operator+(const typename CNT<E>::T& l, const RowVectorBase<E>& r) {
00791     return RowVector_<E>(r) += l;
00792 }
00793 template <class E1, class E2>
00794 RowVector_<typename CNT<E1>::template Result<E2>::Sub>
00795 operator-(const RowVectorBase<E1>& l, const RowVectorBase<E2>& r) {
00796     return RowVector_<typename CNT<E1>::template Result<E2>::Sub>(l) -= r;
00797 }
00798 template <class E>
00799 RowVector_<E> operator-(const RowVectorBase<E>& l, const typename CNT<E>::T& r) {
00800     return RowVector_<E>(l) -= r;
00801 }
00802 template <class E>
00803 RowVector_<E> operator-(const typename CNT<E>::T& l, const RowVectorBase<E>& r) {
00804     RowVector_<E> temp(r.size());
00805     temp = l;
00806     return (temp -= r);
00807 }
00808 
00809 // Scalar multiply and divide 
00810 
00811 template <class E> RowVector_<E>
00812 operator*(const RowVectorBase<E>& l, const typename CNT<E>::StdNumber& r) 
00813   { return RowVector_<E>(l)*=r; }
00814 
00815 template <class E> RowVector_<E>
00816 operator*(const typename CNT<E>::StdNumber& l, const RowVectorBase<E>& r) 
00817   { return RowVector_<E>(r)*=l; }
00818 
00819 template <class E> RowVector_<E>
00820 operator/(const RowVectorBase<E>& l, const typename CNT<E>::StdNumber& r) 
00821   { return RowVector_<E>(l)/=r; }
00822 
00823 // Handle ints explicitly.
00824 template <class E> RowVector_<E>
00825 operator*(const RowVectorBase<E>& l, int r) 
00826   { return RowVector_<E>(l)*= typename CNT<E>::StdNumber(r); }
00827 
00828 template <class E> RowVector_<E>
00829 operator*(int l, const RowVectorBase<E>& r) 
00830   { return RowVector_<E>(r)*= typename CNT<E>::StdNumber(l); }
00831 
00832 template <class E> RowVector_<E>
00833 operator/(const RowVectorBase<E>& l, int r) 
00834   { return RowVector_<E>(l)/= typename CNT<E>::StdNumber(r); }
00835 
00836 
00837 // These are fancier "scalars"; whether they are allowed depends on
00838 // whether the element type and the CNT are compatible.
00839 
00840 // RowVector * Vec
00841 template <class E1, int M, class E2, int S> 
00842 RowVector_<typename CNT<E1>::template Result< Vec<M,E2,S> >::Mul>
00843 operator*(const RowVectorBase<E1>& v, const Vec<M,E2,S>& s) {
00844     RowVector_<typename CNT<E1>::template Result< Vec<M,E2,S> >::Mul> res(v.ncol());
00845     for (int i=0; i < v.ncol(); ++i)
00846         res[i] = v[i]*s; 
00847     return res;
00848 }
00849 
00850 // Vec * RowVector
00851 template <class E1, int M, class E2, int S> 
00852 RowVector_<typename Vec<M,E2,S>::template Result<E1>::Mul>
00853 operator*(const Vec<M,E2,S>& s, const RowVectorBase<E1>& v) {
00854     RowVector_<typename Vec<M,E2,S>::template Result<E1>::Mul> res(v.ncol());
00855     for (int i=0; i < v.ncol(); ++i)
00856         res[i] = s*v[i]; 
00857     return res;
00858 }
00859 
00860 // RowVector * Row
00861 template <class E1, int N, class E2, int S> 
00862 RowVector_<typename CNT<E1>::template Result< Row<N,E2,S> >::Mul>
00863 operator*(const RowVectorBase<E1>& v, const Row<N,E2,S>& s) {
00864     RowVector_<typename CNT<E1>::template Result< Row<N,E2,S> >::Mul> res(v.ncol());
00865     for (int i=0; i < v.ncol(); ++i)
00866         res[i] = v[i]*s; 
00867     return res;
00868 }
00869 
00870 // Row * RowVector
00871 template <class E1, int N, class E2, int S> 
00872 RowVector_<typename Row<N,E2,S>::template Result<E1>::Mul>
00873 operator*(const Row<N,E2,S>& s, const RowVectorBase<E1>& v) {
00874     RowVector_<typename Row<N,E2,S>::template Result<E1>::Mul> res(v.ncol());
00875     for (int i=0; i < v.ncol(); ++i)
00876         res[i] = s*v[i]; 
00877     return res;
00878 }
00879 
00880 // RowVector * Mat
00881 template <class E1, int M, int N, class E2, int S1, int S2> 
00882 RowVector_<typename CNT<E1>::template Result< Mat<M,N,E2,S1,S2> >::Mul>
00883 operator*(const RowVectorBase<E1>& v, const Mat<M,N,E2,S1,S2>& s) {
00884     RowVector_<typename CNT<E1>::template Result< Mat<M,N,E2,S1,S2> >::Mul> res(v.ncol());
00885     for (int i=0; i < v.ncol(); ++i)
00886         res[i] = v[i]*s; 
00887     return res;
00888 }
00889 
00890 // Mat * RowVector
00891 template <class E1, int M, int N, class E2, int S1, int S2> 
00892 RowVector_<typename Mat<M,N,E2,S1,S2>::template Result<E1>::Mul>
00893 operator*(const Mat<M,N,E2,S1,S2>& s, const RowVectorBase<E1>& v) {
00894     RowVector_<typename Mat<M,N,E2,S1,S2>::template Result<E1>::Mul> res(v.ncol());
00895     for (int i=0; i < v.ncol(); ++i)
00896         res[i] = s*v[i]; 
00897     return res;
00898 }
00899 
00900 // RowVector * SymMat
00901 template <class E1, int M, class E2, int S> 
00902 RowVector_<typename CNT<E1>::template Result< SymMat<M,E2,S> >::Mul>
00903 operator*(const RowVectorBase<E1>& v, const SymMat<M,E2,S>& s) {
00904     RowVector_<typename CNT<E1>::template Result< SymMat<M,E2,S> >::Mul> res(v.ncol());
00905     for (int i=0; i < v.ncol(); ++i)
00906         res[i] = v[i]*s; 
00907     return res;
00908 }
00909 
00910 // SymMat * RowVector
00911 template <class E1, int M, class E2, int S> 
00912 RowVector_<typename SymMat<M,E2,S>::template Result<E1>::Mul>
00913 operator*(const SymMat<M,E2,S>& s, const RowVectorBase<E1>& v) {
00914     RowVector_<typename SymMat<M,E2,S>::template Result<E1>::Mul> res(v.ncol());
00915     for (int i=0; i < v.ncol(); ++i)
00916         res[i] = s*v[i]; 
00917     return res;
00918 }
00919 
00921 
00922 
00927 
00928     // TODO: these should use LAPACK!
00929 
00930 // Dot product
00931 template <class E1, class E2> 
00932 typename CNT<E1>::template Result<E2>::Mul
00933 operator*(const RowVectorBase<E1>& r, const VectorBase<E2>& v) {
00934     assert(r.ncol() == v.nrow());
00935     typename CNT<E1>::template Result<E2>::Mul sum(0);
00936     for (int j=0; j < r.ncol(); ++j)
00937         sum += r(j) * v[j];
00938     return sum;
00939 }
00940 
00941 template <class E1, class E2> 
00942 Vector_<typename CNT<E1>::template Result<E2>::Mul>
00943 operator*(const MatrixBase<E1>& m, const VectorBase<E2>& v) {
00944     assert(m.ncol() == v.nrow());
00945     Vector_<typename CNT<E1>::template Result<E2>::Mul> res(m.nrow());
00946     for (int i=0; i< m.nrow(); ++i)
00947         res[i] = m[i]*v;
00948     return res;
00949 }
00950 
00951 template <class E1, class E2> 
00952 Matrix_<typename CNT<E1>::template Result<E2>::Mul>
00953 operator*(const MatrixBase<E1>& m1, const MatrixBase<E2>& m2) {
00954     assert(m1.ncol() == m2.nrow());
00955     Matrix_<typename CNT<E1>::template Result<E2>::Mul> 
00956         res(m1.nrow(),m2.ncol());
00957 
00958     for (int j=0; j < res.ncol(); ++j)
00959         for (int i=0; i < res.nrow(); ++i)
00960             res(i,j) = m1[i] * m2(j);
00961 
00962     return res;
00963 }
00964 
00966 
00967 // This "private" static method is used to implement VectorView's 
00968 // fillVectorViewFromStream() and Vector's readVectorFromStream() 
00969 // namespace-scope static methods, which are in turn used to implement 
00970 // VectorView's and 
00971 // Vector's stream extraction operators ">>". This method has to be in the 
00972 // header file so that we don't need to pass streams through the API, but it 
00973 // is not intended for use by users and has no Doxygen presence, unlike 
00974 // fillArrayFromStream() and readArrayFromStream() and (more commonly)
00975 // the extraction operators.
00976 template <class T> static inline 
00977 std::istream& readVectorFromStreamHelper
00978    (std::istream& in, bool isFixedSize, Vector_<T>& out)
00979 {
00980     // If already failed, bad, or eof, set failed bit and return without 
00981     // touching the Vector.
00982     if (!in.good()) {in.setstate(std::ios::failbit); return in;}
00983 
00984     // If the passed-in Vector isn't resizeable, then we have to treat it as
00985     // a fixed size VectorView regardless of the setting of the isFixedSize
00986     // argument.
00987     if (!out.isResizeable())
00988         isFixedSize = true; // might be overriding the argument here
00989 
00990     // numRequired will be ignored unless isFixedSize==true.
00991     const int numRequired = isFixedSize ? out.size() : 0;
00992 
00993     if (!isFixedSize)
00994         out.clear(); // We're going to replace the entire contents of the Array.
00995 
00996     // Skip initial whitespace. If that results in eof this may be a successful
00997     // read of a 0-length, unbracketed Vector. That is OK for either a
00998     // variable-length Vector or a fixed-length VectorView of length zero.
00999     std::ws(in); if (in.fail()) return in;
01000     if (in.eof()) {
01001         if (isFixedSize && numRequired != 0)
01002             in.setstate(std::ios_base::failbit); // zero elements not OK
01003         return in;
01004     }
01005     
01006     // Here the stream is good and the next character is non-white.
01007     assert(in.good());
01008 
01009     // Use this for raw i/o (peeks and gets).
01010     typename       std::iostream::int_type ch;
01011 #ifndef NDEBUG
01012     const typename std::iostream::int_type EOFch = 
01013         std::iostream::traits_type::eof();
01014 #endif
01015 
01016     // First we'll look for the optional "~". If found, the brackets become
01017     // required.
01018     bool tildeFound = false;
01019     ch = in.peek(); if (in.fail()) return in;
01020     assert(ch != EOFch); // we already checked above
01021     if ((char)ch == '~') {
01022         tildeFound = true;
01023         in.get(); // absorb the tilde
01024         // Eat whitespace after the tilde to see what's next.
01025         if (in.good()) std::ws(in);
01026         // If we hit eof after the tilde we don't like the formatting.
01027         if (!in.good()) {in.setstate(std::ios_base::failbit); return in;}
01028     }
01029 
01030     // Here the stream is good, the next character is non-white, and we
01031     // might have seen a tilde.
01032     assert(in.good());
01033 
01034     // Now see if the sequence is bare or surrounded by (), or [].
01035     bool lookForCloser = true;
01036     char openBracket, closeBracket;
01037     ch = in.peek(); if (in.fail()) return in;
01038     assert(ch != EOFch); // we already checked above
01039 
01040     openBracket = (char)ch;
01041     if      (openBracket=='(') {in.get(); closeBracket = ')';}
01042     else if (openBracket=='[') {in.get(); closeBracket = ']';}
01043     else lookForCloser = false;
01044 
01045     // If we found a tilde, the opening bracket was mandatory. If we didn't
01046     // find one then we reject the formatting.
01047     if (tildeFound && !lookForCloser)
01048     {   in.setstate(std::ios_base::failbit); return in;}
01049 
01050     // If lookForCloser is true, then closeBracket contains the terminating
01051     // delimiter, otherwise we're not going to quit until eof.
01052 
01053     // Eat whitespace after the opening bracket to see what's next.
01054     if (in.good()) std::ws(in);
01055 
01056     // If we're at eof now it must be because the open bracket was the
01057     // last non-white character in the stream, which is an error.
01058     if (!in.good()) {
01059         if (in.eof()) {
01060             assert(lookForCloser); // or we haven't read anything that could eof
01061             in.setstate(std::ios::failbit);
01062         }
01063         return in;
01064     }
01065 
01066     // istream is good and next character is non-white; ready to read first
01067     // value or terminator.
01068 
01069     // We need to figure out whether the elements are space- or comma-
01070     // separated and then insist on consistency.
01071     bool commaOK = true, commaRequired = false;
01072     bool terminatorSeen = false;
01073     int nextIndex = 0;
01074     while (true) {
01075         char c;
01076 
01077         // Here at the top of this loop, we have already successfully read 
01078         // n=nextIndex values of type T. For fixed-size reads, it might be
01079         // the case that n==numRequired already, but we still may need to
01080         // look for a closing bracket before we can declare victory.
01081         // The stream is good() (not at eof) but it might be the case that 
01082         // there is nothing but white space left; we don't know yet because
01083         // if we have satisfied the fixed-size count and are not expecting
01084         // a terminator then we should quit without absorbing the trailing
01085         // white space.
01086         assert(in.good());
01087 
01088         // Look for closing bracket before trying to read value.
01089         if (lookForCloser) {
01090             // Eat white space to find the closing bracket.
01091             std::ws(in); if (!in.good()) break; // eof?
01092             ch = in.peek(); assert(ch != EOFch);
01093             if (!in.good()) break;
01094             c = (char)ch;
01095             if (c == closeBracket) {   
01096                 in.get(); // absorb the closing bracket
01097                 terminatorSeen = true; 
01098                 break; 
01099             }
01100             // next char not a closing bracket; fall through
01101         }
01102 
01103         // We didn't look or didn't find a closing bracket. The istream is good 
01104         // but we might be looking at white space.
01105 
01106         // If we already got all the elements we want, break for final checks.
01107         if (isFixedSize && (nextIndex == numRequired))
01108             break; // that's a full count.
01109 
01110         // Look for comma before value, except the first time.
01111         if (commaOK && nextIndex != 0) {
01112             // Eat white space to find the comma.
01113             std::ws(in); if (!in.good()) break; // eof?
01114             ch = in.peek(); assert(ch != EOFch);
01115             if (!in.good()) break;
01116             c = (char)ch;
01117             if (c == ',') {
01118                 in.get(); // absorb comma
01119                 commaRequired = true; // all commas from now on
01120             } else { // next char not a comma
01121                 if (commaRequired) // bad, e.g.: v1, v2, v3 v4 
01122                 {   in.setstate(std::ios::failbit); break; }
01123                 else commaOK = false; // saw: v1 v2 (no commas now)
01124             }
01125             if (!in.good()) break; // might be eof
01126         }
01127 
01128         // No closing bracket yet; don't have enough elements; skipped comma 
01129         // if any; istream is good; might be looking at white space.
01130         assert(in.good());
01131 
01132         // Now read in an element of type T.
01133         // The extractor T::operator>>() will ignore leading white space.
01134         if (!isFixedSize)
01135             out.resizeKeep(out.size()+1); // grow by one (default consructed)
01136         in >> out[nextIndex]; if (in.fail()) break;
01137         ++nextIndex;
01138 
01139         if (!in.good()) break; // might be eof
01140     }
01141 
01142     // We will get here under a number of circumstances:
01143     //  - the fail bit is set in the istream, or
01144     //  - we reached eof
01145     //  - we saw a closing brace
01146     //  - we got all the elements we wanted (for a fixed-size read)
01147     // Note that it is possible that we consumed everything except some
01148     // trailing white space (meaning we're not technically at eof), but
01149     // for consistency with built-in operator>>()'s we won't try to absorb
01150     // that trailing white space.
01151 
01152     if (!in.fail()) {
01153         if (lookForCloser && !terminatorSeen)
01154             in.setstate(std::ios::failbit); // missing terminator
01155 
01156         if (isFixedSize && nextIndex != numRequired)
01157             in.setstate(std::ios::failbit); // wrong number of values
01158     }
01159 
01160     return in;
01161 }
01162 
01163 
01164 
01165 //------------------------------------------------------------------------------
01166 //                          RELATED GLOBAL OPERATORS
01167 //------------------------------------------------------------------------------
01168 // These are logically part of the Matrix_<T> class but are not actually 
01169 // class members; that is, they are in the SimTK namespace.
01170 
01182 template <class E> inline void
01183 writeUnformatted(std::ostream& o, const VectorBase<E>& v) {
01184     const int sz = v.size();
01185     for (int i=0; i < sz; ++i) {
01186         if (i != 0) o << " ";
01187         writeUnformatted(o, v[i]);
01188     }
01189 } 
01192 template <class E> inline void
01193 writeUnformatted(std::ostream& o, const VectorView_<E>& v) 
01194 {   writeUnformatted(o, static_cast< const VectorBase<E> >(v)); }
01195 
01198 template <class E> inline void
01199 writeUnformatted(std::ostream& o, const Vector_<E>& v) 
01200 {   writeUnformatted(o, static_cast< const VectorBase<E> >(v)); }
01201 
01205 template <class E> inline void
01206 writeUnformatted(std::ostream& o, const RowVectorBase<E>& v) 
01207 {   writeUnformatted(o, ~v); }
01208 
01211 template <class E> inline void
01212 writeUnformatted(std::ostream& o, const RowVectorView_<E>& v) 
01213 {   writeUnformatted(o, static_cast< const RowVectorBase<E> >(v)); }
01214 
01217 template <class E> inline void
01218 writeUnformatted(std::ostream& o, const RowVector_<E>& v) 
01219 {   writeUnformatted(o, static_cast< const RowVectorBase<E> >(v)); }
01220 
01224 template <class E> inline void
01225 writeUnformatted(std::ostream& o, const MatrixBase<E>& v) {
01226     const int nr = v.nrow();
01227     for (int i=0; i < nr; ++i) {
01228         if (i != 0) o << std::endl;
01229         writeUnformatted(o, v[i]);
01230     }
01231 }
01234 template <class E> inline void
01235 writeUnformatted(std::ostream& o, const MatrixView_<E>& v) 
01236 {   writeUnformatted(o, static_cast< const MatrixBase<E> >(v)); }
01237 
01240 template <class E> inline void
01241 writeUnformatted(std::ostream& o, const Matrix_<E>& v) 
01242 {   writeUnformatted(o, static_cast< const MatrixBase<E> >(v)); }
01243 
01247 template <class E> inline bool
01248 readUnformatted(std::istream& in, VectorView_<E>& v) {
01249     for (int i=0; i < v.size(); ++i)
01250         if (!readUnformatted(in, v[i])) return false;
01251     return true;
01252 }
01253 
01256 template <class E> inline bool
01257 readUnformatted(std::istream& in, Vector_<E>& v) {
01258     if (!v.isResizeable())
01259         return readUnformatted(in, v.updAsVectorView());
01260 
01261     Array_<E,int> a;
01262     if (!readUnformatted(in,a)) return false;
01263     v.resize(a.size());
01264     for (int i=0; i<a.size(); ++i)
01265         v[i] = a[i];
01266     return true;
01267 }
01268 
01272 template <class E> inline bool
01273 readUnformatted(std::istream& in, RowVectorView_<E>& v) 
01274 {   VectorView_<E> vt(~v);
01275     return readUnformatted<E>(in, vt); }
01276 
01280 template <class E> inline bool
01281 readUnformatted(std::istream& in, RowVector_<E>& v) 
01282 {   Vector_<E> vt(~v);
01283     return readUnformatted<E>(in, vt); }
01284 
01290 template <class E> inline bool
01291 readUnformatted(std::istream& in, MatrixView_<E>& v) { 
01292     for (int row=0; row < v.nrow(); ++row) {
01293         RowVectorView_<E> oneRow(v[row]);
01294         if (!readUnformatted<E>(in, oneRow)) return false;
01295     }
01296     return true;
01297 }
01298 
01304 template <class E> inline bool
01305 fillUnformatted(std::istream& in, Matrix_<E>& v) {
01306     return readUnformatted<E>(in, v.updAsMatrixView());
01307 }
01308 
01312 template <class E> inline bool
01313 readUnformatted(std::istream& in, Matrix_<E>& v) {
01314     SimTK_ASSERT_ALWAYS(!"implemented", 
01315         "SimTK::readUnformatted(istream, Matrix) is not implemented; try"
01316         " SimTK::fillUnformatted(istream, Matrix) instead.");
01317     return false;
01318 }
01319   
01326 template <class T> inline std::ostream&
01327 operator<<(std::ostream& o, const VectorBase<T>& v)
01328 {   o << "~["; 
01329     if (v.size()) {
01330         o << v[0];
01331         for (int i=1; i < v.size(); ++i) o << " " << v[i];
01332     }
01333     return o << "]"; 
01334 }
01335 
01342 template <class T> inline std::ostream&
01343 operator<<(std::ostream& o, const RowVectorBase<T>& v)
01344 {   o << "["; 
01345     if (v.size()) {
01346         o << v[0];
01347         for (int i=1; i < v.size(); ++i) o << " " << v[i];
01348     }
01349     return o << "]"; 
01350 }
01351 
01359 template <class T> inline std::ostream&
01360 operator<<(std::ostream& o, const MatrixBase<T>& m) {
01361     for (int i=0;i<m.nrow();++i)
01362         o << std::endl << m[i];
01363     if (m.nrow()) o << std::endl;
01364     return o; 
01365 }
01366 
01367 
01399 template <class T> static inline 
01400 std::istream& readVectorFromStream(std::istream& in, Vector_<T>& out)
01401 {   return readVectorFromStreamHelper<T>(in, false /*variable sizez*/, out); }
01402 
01403 
01404 
01429 template <class T> static inline 
01430 std::istream& fillVectorFromStream(std::istream& in, Vector_<T>& out)
01431 {   return readVectorFromStreamHelper<T>(in, true /*fixed size*/, out); }
01432 
01437 template <class T> static inline 
01438 std::istream& fillVectorViewFromStream(std::istream& in, VectorView_<T>& out)
01439 {   return readVectorFromStreamHelper<T>(in, true /*fixed size*/, out); }
01440 
01441 
01451 template <class T> inline
01452 std::istream& operator>>(std::istream& in, Vector_<T>& out) 
01453 {   return readVectorFromStream<T>(in, out); }
01454 
01462 template <class T> inline
01463 std::istream& operator>>(std::istream& in, VectorView_<T>& out) 
01464 {   return fillVectorViewFromStream<T>(in, out); }  // End of Matrix serialization.
01466 
01467 // Friendly abbreviations for vectors and matrices with scalar elements.
01473 typedef Vector_<Real>           Vector;
01474 
01478 typedef Matrix_<Real>           Matrix;
01479 
01483 typedef RowVector_<Real>        RowVector;
01489 typedef Vector_<Complex>        ComplexVector;
01491 typedef Matrix_<Complex>        ComplexMatrix;
01493 typedef RowVector_<Complex>     ComplexRowVector;
01494 
01496 typedef VectorView_<Real>       VectorView;
01498 typedef MatrixView_<Real>       MatrixView;
01500 typedef RowVectorView_<Real>    RowVectorView;
01501 
01503 typedef VectorView_<Complex>    ComplexVectorView;
01505 typedef MatrixView_<Complex>    ComplexMatrixView;
01507 typedef RowVectorView_<Complex> ComplexRowVectorView;
01508 
01510 typedef Vector_<float>          fVector;
01512 typedef Vector_<double>         dVector;
01514 typedef Vector_<fComplex>       fComplexVector;
01516 typedef Vector_<dComplex>       dComplexVector;
01517 
01519 typedef VectorView_<float>      fVectorView;
01521 typedef VectorView_<double>     dVectorView;
01523 typedef VectorView_<fComplex>   fComplexVectorView;
01525 typedef VectorView_<dComplex>   dComplexVectorView;
01526 
01528 typedef RowVector_<float>       fRowVector;
01530 typedef RowVector_<double>      dRowVector;
01532 typedef RowVector_<fComplex>    fComplexRowVector;
01534 typedef RowVector_<dComplex>    dComplexRowVector;
01535 
01537 typedef RowVectorView_<float>   fRowVectorView;
01539 typedef RowVectorView_<double>  dRowVectorView;
01541 typedef RowVectorView_<fComplex> fComplexRowVectorView;
01543 typedef RowVectorView_<dComplex> dComplexRowVectorView;
01544 
01546 typedef Matrix_<float>          fMatrix;
01548 typedef Matrix_<double>         dMatrix;
01550 typedef Matrix_<fComplex>       fComplexMatrix;
01552 typedef Matrix_<dComplex>       dComplexMatrix;
01553 
01555 typedef MatrixView_<float>      fMatrixView;
01557 typedef MatrixView_<double>     dMatrixView;
01559 typedef MatrixView_<fComplex>   fComplexMatrixView;
01561 typedef MatrixView_<dComplex>   dComplexMatrixView;
01564 } //namespace SimTK
01565 
01566 #endif //SimTK_SIMMATRIX_BIGMATRIX_H_
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines