Simbody
3.5
|
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_