Simbody  3.5
Vec.h
Go to the documentation of this file.
1 #ifndef SimTK_SIMMATRIX_SMALLMATRIX_VEC_H_
2 #define SimTK_SIMMATRIX_SMALLMATRIX_VEC_H_
3 
4 /* -------------------------------------------------------------------------- *
5  * Simbody(tm): SimTKcommon *
6  * -------------------------------------------------------------------------- *
7  * This is part of the SimTK biosimulation toolkit originating from *
8  * Simbios, the NIH National Center for Physics-Based Simulation of *
9  * Biological Structures at Stanford, funded under the NIH Roadmap for *
10  * Medical Research, grant U54 GM072970. See https://simtk.org/home/simbody. *
11  * *
12  * Portions copyright (c) 2005-12 Stanford University and the Authors. *
13  * Authors: Michael Sherman *
14  * Contributors: Peter Eastman *
15  * *
16  * Licensed under the Apache License, Version 2.0 (the "License"); you may *
17  * not use this file except in compliance with the License. You may obtain a *
18  * copy of the License at http://www.apache.org/licenses/LICENSE-2.0. *
19  * *
20  * Unless required by applicable law or agreed to in writing, software *
21  * distributed under the License is distributed on an "AS IS" BASIS, *
22  * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. *
23  * See the License for the specific language governing permissions and *
24  * limitations under the License. *
25  * -------------------------------------------------------------------------- */
26 
32 
33 namespace SimTK {
34 
35 
36 // The following functions are used internally by Vec.
37 
38 // Hide from Doxygen.
40 namespace Impl {
41 
42 // For those wimpy compilers that don't unroll short, constant-limit loops,
43 // Peter Eastman added these recursive template implementations of
44 // elementwise add, subtract, and copy. Sherm added multiply and divide.
45 
46 template <class E1, int S1, class E2, int S2> void
47 conformingAdd(const Vec<1,E1,S1>& r1, const Vec<1,E2,S2>& r2,
48  Vec<1,typename CNT<E1>::template Result<E2>::Add>& result) {
49  result[0] = r1[0] + r2[0];
50 }
51 template <int N, class E1, int S1, class E2, int S2> void
52 conformingAdd(const Vec<N,E1,S1>& r1, const Vec<N,E2,S2>& r2,
53  Vec<N,typename CNT<E1>::template Result<E2>::Add>& result) {
54  conformingAdd(reinterpret_cast<const Vec<N-1,E1,S1>&>(r1),
55  reinterpret_cast<const Vec<N-1,E2,S2>&>(r2),
56  reinterpret_cast<Vec<N-1,typename CNT<E1>::
57  template Result<E2>::Add>&>(result));
58  result[N-1] = r1[N-1] + r2[N-1];
59 }
60 
61 template <class E1, int S1, class E2, int S2> void
62 conformingSubtract(const Vec<1,E1,S1>& r1, const Vec<1,E2,S2>& r2,
63  Vec<1,typename CNT<E1>::template Result<E2>::Sub>& result) {
64  result[0] = r1[0] - r2[0];
65 }
66 template <int N, class E1, int S1, class E2, int S2> void
67 conformingSubtract(const Vec<N,E1,S1>& r1, const Vec<N,E2,S2>& r2,
68  Vec<N,typename CNT<E1>::template Result<E2>::Sub>& result) {
69  conformingSubtract(reinterpret_cast<const Vec<N-1,E1,S1>&>(r1),
70  reinterpret_cast<const Vec<N-1,E2,S2>&>(r2),
71  reinterpret_cast<Vec<N-1,typename CNT<E1>::
72  template Result<E2>::Sub>&>(result));
73  result[N-1] = r1[N-1] - r2[N-1];
74 }
75 
76 template <class E1, int S1, class E2, int S2> void
77 elementwiseMultiply(const Vec<1,E1,S1>& r1, const Vec<1,E2,S2>& r2,
78  Vec<1,typename CNT<E1>::template Result<E2>::Mul>& result) {
79  result[0] = r1[0] * r2[0];
80 }
81 template <int N, class E1, int S1, class E2, int S2> void
82 elementwiseMultiply(const Vec<N,E1,S1>& r1, const Vec<N,E2,S2>& r2,
83  Vec<N,typename CNT<E1>::template Result<E2>::Mul>& result) {
84  elementwiseMultiply(reinterpret_cast<const Vec<N-1,E1,S1>&>(r1),
85  reinterpret_cast<const Vec<N-1,E2,S2>&>(r2),
86  reinterpret_cast<Vec<N-1,typename CNT<E1>::
87  template Result<E2>::Mul>&>(result));
88  result[N-1] = r1[N-1] * r2[N-1];
89 }
90 
91 template <class E1, int S1, class E2, int S2> void
92 elementwiseDivide(const Vec<1,E1,S1>& r1, const Vec<1,E2,S2>& r2,
93  Vec<1,typename CNT<E1>::template Result<E2>::Dvd>& result) {
94  result[0] = r1[0] / r2[0];
95 }
96 template <int N, class E1, int S1, class E2, int S2> void
97 elementwiseDivide(const Vec<N,E1,S1>& r1, const Vec<N,E2,S2>& r2,
98  Vec<N,typename CNT<E1>::template Result<E2>::Dvd>& result) {
99  elementwiseDivide(reinterpret_cast<const Vec<N-1,E1,S1>&>(r1),
100  reinterpret_cast<const Vec<N-1,E2,S2>&>(r2),
101  reinterpret_cast<Vec<N-1,typename CNT<E1>::
102  template Result<E2>::Dvd>&>(result));
103  result[N-1] = r1[N-1] / r2[N-1];
104 }
105 
106 template <class E1, int S1, class E2, int S2> void
107 copy(Vec<1,E1,S1>& r1, const Vec<1,E2,S2>& r2) {
108  r1[0] = r2[0];
109 }
110 template <int N, class E1, int S1, class E2, int S2> void
111 copy(Vec<N,E1,S1>& r1, const Vec<N,E2,S2>& r2) {
112  copy(reinterpret_cast<Vec<N-1,E1,S1>&>(r1),
113  reinterpret_cast<const Vec<N-1,E2,S2>&>(r2));
114  r1[N-1] = r2[N-1];
115 }
116 
117 }
183 template <int M, class ELT, int STRIDE>
184 class Vec {
185 public:
191  typedef ELT E;
193  typedef typename CNT<E>::TNeg ENeg;
198  typedef typename CNT<E>::TReal EReal;
202  typedef typename CNT<E>::TImag EImag;
205  typedef typename CNT<E>::TComplex EComplex;
207  typedef typename CNT<E>::THerm EHerm;
209  typedef typename CNT<E>::TPosTrans EPosTrans;
212  typedef typename CNT<E>::TSqHermT ESqHermT;
214  typedef typename CNT<E>::TSqTHerm ESqTHerm;
216  typedef typename CNT<E>::TSqrt ESqrt;
218  typedef typename CNT<E>::TAbs EAbs;
221  typedef typename CNT<E>::TStandard EStandard;
224  typedef typename CNT<E>::TInvert EInvert;
226  typedef typename CNT<E>::TNormalize ENormalize;
227 
228  typedef typename CNT<E>::Scalar EScalar;
230  typedef typename CNT<E>::Number ENumber;
231  typedef typename CNT<E>::StdNumber EStdNumber;
232  typedef typename CNT<E>::Precision EPrecision;
234 
237  enum {
238  NRows = M,
239  NCols = 1,
240  NPackedElements = M,
241  NActualElements = M * STRIDE, // includes trailing gap
242  NActualScalars = CNT<E>::NActualScalars * NActualElements,
243  RowSpacing = STRIDE,
244  ColSpacing = NActualElements,
246  RealStrideFactor = 1, // composite types don't change size when
247  // cast from complex to real or imaginary
248  ArgDepth = ((int)CNT<E>::ArgDepth < (int)MAX_RESOLVED_DEPTH
249  ? CNT<E>::ArgDepth + 1
251  IsScalar = 0,
252  IsULessScalar = 0,
253  IsNumber = 0,
254  IsStdNumber = 0,
255  IsPrecision = 0,
256  SignInterpretation = CNT<E>::SignInterpretation
257  };
258 
259  // These are reinterpretations of the current data, so have the
260  // same packing (stride).
261 
287  typedef E TElement;
289  typedef E TRow;
291  typedef Vec TCol;
292 
293  // These are the results of calculations, so are returned in new, packed
294  // memory. Be sure to refer to element types here which are also packed.
295  typedef Vec<M,ESqrt,1> TSqrt; // Note stride
296  typedef Vec<M,EAbs,1> TAbs; // Note stride
300 
301  typedef ESqHermT TSqHermT; // result of self dot product
302  typedef SymMat<M,ESqTHerm> TSqTHerm; // result of self outer product
303 
304  // These recurse right down to the underlying scalar type no matter how
305  // deep the elements are.
306  typedef EScalar Scalar;
307  typedef EULessScalar ULessScalar;
308  typedef ENumber Number;
309  typedef EStdNumber StdNumber;
310  typedef EPrecision Precision;
311  typedef EScalarNormSq ScalarNormSq;
316  static int size() { return M; }
318  static int nrow() { return M; }
320  static int ncol() { return 1; }
321 
322 
325  ScalarNormSq scalarNormSqr() const {
326  ScalarNormSq sum(0);
327  for(int i=0;i<M;++i) sum += CNT<E>::scalarNormSqr(d[i*STRIDE]);
328  return sum;
329  }
330 
335  TSqrt sqrt() const {
336  TSqrt vsqrt;
337  for(int i=0;i<M;++i) vsqrt[i] = CNT<E>::sqrt(d[i*STRIDE]);
338  return vsqrt;
339  }
340 
345  TAbs abs() const {
346  TAbs vabs;
347  for(int i=0;i<M;++i) vabs[i] = CNT<E>::abs(d[i*STRIDE]);
348  return vabs;
349  }
350 
355  TStandard standardize() const {
356  TStandard vstd;
357  for(int i=0;i<M;++i) vstd[i] = CNT<E>::standardize(d[i*STRIDE]);
358  return vstd;
359  }
360 
364  EStandard sum() const {
365  E sum(0);
366  for (int i=0;i<M;++i) sum += d[i*STRIDE];
367  return CNT<E>::standardize(sum);
368  }
369 
370 
371  // This gives the resulting vector type when (v[i] op P) is applied to
372  // each element of v. It is a vector of length M, stride 1, and element
373  // types which are the regular composite result of E op P. Typically P is
374  // a scalar type but it doesn't have to be.
375  template <class P> struct EltResult {
380  };
381 
382  // This is the composite result for v op P where P is some kind of
383  // appropriately shaped non-scalar type.
384  template <class P> struct Result {
385  typedef MulCNTs<M,1,ArgDepth,Vec,ColSpacing,RowSpacing,
388  typedef typename MulOp::Type Mul;
389 
390  typedef MulCNTsNonConforming<M,1,ArgDepth,Vec,ColSpacing,RowSpacing,
391  CNT<P>::NRows, CNT<P>::NCols, CNT<P>::ArgDepth,
392  P, CNT<P>::ColSpacing, CNT<P>::RowSpacing> MulOpNonConforming;
393  typedef typename MulOpNonConforming::Type MulNon;
394 
395  typedef DvdCNTs<M,1,ArgDepth,Vec,ColSpacing,RowSpacing,
396  CNT<P>::NRows, CNT<P>::NCols, CNT<P>::ArgDepth,
397  P, CNT<P>::ColSpacing, CNT<P>::RowSpacing> DvdOp;
398  typedef typename DvdOp::Type Dvd;
399 
400  typedef AddCNTs<M,1,ArgDepth,Vec,ColSpacing,RowSpacing,
401  CNT<P>::NRows, CNT<P>::NCols, CNT<P>::ArgDepth,
402  P, CNT<P>::ColSpacing, CNT<P>::RowSpacing> AddOp;
403  typedef typename AddOp::Type Add;
404 
405  typedef SubCNTs<M,1,ArgDepth,Vec,ColSpacing,RowSpacing,
406  CNT<P>::NRows, CNT<P>::NCols, CNT<P>::ArgDepth,
407  P, CNT<P>::ColSpacing, CNT<P>::RowSpacing> SubOp;
408  typedef typename SubOp::Type Sub;
409  };
410 
416  template <class P> struct Substitute {
417  typedef Vec<M,P> Type;
418  };
419 
423  Vec(){
424  #ifndef NDEBUG
425  setToNaN();
426  #endif
427  }
428 
429  // It's important not to use the default copy constructor or copy
430  // assignment because the compiler doesn't understand that we may
431  // have noncontiguous storage and will try to copy the whole array.
432 
436  Vec(const Vec& src) {
437  Impl::copy(*this, src);
438  }
443  Vec& operator=(const Vec& src) {
444  Impl::copy(*this, src);
445  return *this;
446  }
447 
450  template <int SS> Vec(const Vec<M,E,SS>& src) {
451  Impl::copy(*this, src);
452  }
453 
456  template <int SS> Vec(const Vec<M,ENeg,SS>& src) {
457  Impl::copy(*this, src);
458  }
459 
462  template <class EE, int SS> explicit Vec(const Vec<M,EE,SS>& src) {
463  Impl::copy(*this, src);
464  }
465 
468  explicit Vec(const E& e) {for (int i=0;i<M;++i) d[i*STRIDE]=e;}
469 
474  explicit Vec(const ENeg& ne) {
475  const E e = ne; // requires floating point negation
476  for (int i=0;i<M;++i) d[i*STRIDE]=e;
477  }
478 
483  explicit Vec(int i) {new (this) Vec(E(Precision(i)));}
484 
485  // A bevy of constructors for Vecs up to length 9.
486 
488  Vec(const E& e0,const E& e1)
489  { assert(M==2);(*this)[0]=e0;(*this)[1]=e1; }
490  Vec(const E& e0,const E& e1,const E& e2)
491  { assert(M==3);(*this)[0]=e0;(*this)[1]=e1;(*this)[2]=e2; }
492  Vec(const E& e0,const E& e1,const E& e2,const E& e3)
493  { assert(M==4);(*this)[0]=e0;(*this)[1]=e1;(*this)[2]=e2;(*this)[3]=e3; }
494  Vec(const E& e0,const E& e1,const E& e2,const E& e3,const E& e4)
495  { assert(M==5);(*this)[0]=e0;(*this)[1]=e1;(*this)[2]=e2;
496  (*this)[3]=e3;(*this)[4]=e4; }
497  Vec(const E& e0,const E& e1,const E& e2,const E& e3,const E& e4,const E& e5)
498  { assert(M==6);(*this)[0]=e0;(*this)[1]=e1;(*this)[2]=e2;
499  (*this)[3]=e3;(*this)[4]=e4;(*this)[5]=e5; }
500  Vec(const E& e0,const E& e1,const E& e2,const E& e3,const E& e4,const E& e5, const E& e6)
501  { assert(M==7);(*this)[0]=e0;(*this)[1]=e1;(*this)[2]=e2;
502  (*this)[3]=e3;(*this)[4]=e4;(*this)[5]=e5;(*this)[6]=e6; }
503  Vec(const E& e0,const E& e1,const E& e2,const E& e3,const E& e4,const E& e5, const E& e6, const E& e7)
504  { assert(M==8);(*this)[0]=e0;(*this)[1]=e1;(*this)[2]=e2;
505  (*this)[3]=e3;(*this)[4]=e4;(*this)[5]=e5;(*this)[6]=e6;(*this)[7]=e7; }
506  Vec(const E& e0,const E& e1,const E& e2,const E& e3,const E& e4,const E& e5, const E& e6, const E& e7, const E& e8)
507  { assert(M==9);(*this)[0]=e0;(*this)[1]=e1;(*this)[2]=e2;
508  (*this)[3]=e3;(*this)[4]=e4;(*this)[5]=e5;(*this)[6]=e6;(*this)[7]=e7;(*this)[8]=e8; }
509 
514  template <class EE> explicit Vec(const EE* p)
515  { assert(p); for(int i=0;i<M;++i) d[i*STRIDE]=p[i]; }
516 
521  template <class EE> Vec& operator=(const EE* p)
522  { assert(p); for(int i=0;i<M;++i) d[i*STRIDE]=p[i]; return *this; }
523 
526  template <class EE, int SS> Vec& operator=(const Vec<M,EE,SS>& vv)
527  { Impl::copy(*this, vv); return *this; }
528 
531  template <class EE, int SS> Vec& operator+=(const Vec<M,EE,SS>& r)
532  { for(int i=0;i<M;++i) d[i*STRIDE] += r[i]; return *this; }
536  template <class EE, int SS> Vec& operator+=(const Vec<M,negator<EE>,SS>& r)
537  { for(int i=0;i<M;++i) d[i*STRIDE] -= -(r[i]); return *this; }
538 
541  template <class EE, int SS> Vec& operator-=(const Vec<M,EE,SS>& r)
542  { for(int i=0;i<M;++i) d[i*STRIDE] -= r[i]; return *this; }
546  template <class EE, int SS> Vec& operator-=(const Vec<M,negator<EE>,SS>& r)
547  { for(int i=0;i<M;++i) d[i*STRIDE] += -(r[i]); return *this; }
548 
549  // Conforming binary ops with 'this' on left, producing new packed result.
550  // Cases: v=v+v, v=v-v, m=v*r
551 
553  template <class EE, int SS> Vec<M,typename CNT<E>::template Result<EE>::Add>
554  conformingAdd(const Vec<M,EE,SS>& r) const {
555  Vec<M,typename CNT<E>::template Result<EE>::Add> result;
556  Impl::conformingAdd(*this, r, result);
557  return result;
558  }
560  template <class EE, int SS> Vec<M,typename CNT<E>::template Result<EE>::Sub>
562  Vec<M,typename CNT<E>::template Result<EE>::Sub> result;
563  Impl::conformingSubtract(*this, r, result);
564  return result;
565  }
566 
569  template <class EE, int SS> Mat<M,M,typename CNT<E>::template Result<EE>::Mul>
571  Mat<M,M,typename CNT<E>::template Result<EE>::Mul> result;
572  for (int j=0;j<M;++j) result(j) = scalarMultiply(r(j));
573  return result;
574  }
575 
577  template <class EE, int SS> Vec<M,typename CNT<E>::template Result<EE>::Mul>
579  Vec<M,typename CNT<E>::template Result<EE>::Mul> result;
580  Impl::elementwiseMultiply(*this, r, result);
581  return result;
582  }
584  template <class EE, int SS> Vec<M,typename CNT<E>::template Result<EE>::Dvd>
585  elementwiseDivide(const Vec<M,EE,SS>& r) const {
586  Vec<M,typename CNT<E>::template Result<EE>::Dvd> result;
587  Impl::elementwiseDivide(*this, r, result);
588  return result;
589  }
590 
594  const E& operator[](int i) const
595  { assert(0 <= i && i < M); return d[i*STRIDE]; }
597  const E& operator()(int i) const {return (*this)[i];}
598 
602  E& operator[](int i) {assert(0 <= i && i < M); return d[i*STRIDE];}
604  E& operator()(int i) {return (*this)[i];}
605 
606  ScalarNormSq normSqr() const { return scalarNormSqr(); }
607  typename CNT<ScalarNormSq>::TSqrt
608  norm() const { return CNT<ScalarNormSq>::sqrt(scalarNormSqr()); }
609 
621  TNormalize normalize() const {
622  if (CNT<E>::IsScalar) {
623  return castAwayNegatorIfAny() / (SignInterpretation*norm());
624  } else {
625  TNormalize elementwiseNormalized;
626  for (int i=0; i<M; ++i)
627  elementwiseNormalized[i] = CNT<E>::normalize((*this)[i]);
628  return elementwiseNormalized;
629  }
630  }
631 
633  TInvert invert() const {assert(false); return TInvert();} // TODO default inversion
634 
636  const Vec& operator+() const { return *this; }
640  const TNeg& operator-() const { return negate(); }
643  TNeg& operator-() { return updNegate(); }
647  const THerm& operator~() const { return transpose(); }
651  THerm& operator~() { return updTranspose(); }
652 
654  const TNeg& negate() const { return *reinterpret_cast<const TNeg*>(this); }
657  TNeg& updNegate() { return *reinterpret_cast< TNeg*>(this); }
658 
660  const THerm& transpose() const { return *reinterpret_cast<const THerm*>(this); }
663  THerm& updTranspose() { return *reinterpret_cast< THerm*>(this); }
664 
669  const TPosTrans& positionalTranspose() const
670  { return *reinterpret_cast<const TPosTrans*>(this); }
673  { return *reinterpret_cast<TPosTrans*>(this); }
674 
679  const TReal& real() const { return *reinterpret_cast<const TReal*>(this); }
682  TReal& real() { return *reinterpret_cast< TReal*>(this); }
683 
684  // Had to contort these next two routines to get them through VC++ 7.net
685 
690  const TImag& imag() const {
691  const int offs = ImagOffset;
692  const EImag* p = reinterpret_cast<const EImag*>(this);
693  return *reinterpret_cast<const TImag*>(p+offs);
694  }
697  TImag& imag() {
698  const int offs = ImagOffset;
699  EImag* p = reinterpret_cast<EImag*>(this);
700  return *reinterpret_cast<TImag*>(p+offs);
701  }
702 
706  const TWithoutNegator& castAwayNegatorIfAny() const
707  { return *reinterpret_cast<const TWithoutNegator*>(this); }
710  TWithoutNegator& updCastAwayNegatorIfAny()
711  { return *reinterpret_cast<TWithoutNegator*>(this); }
712 
713  // These are elementwise binary operators, (this op ee) by default but
714  // (ee op this) if 'FromLeft' appears in the name. The result is a packed
715  // Vec<M> but the element type may change. These are mostly used to
716  // implement global operators. We call these "scalar" operators but
717  // actually the "scalar" can be a composite type.
718 
719  //TODO: consider converting 'e' to Standard Numbers as precalculation and
720  // changing return type appropriately.
721  template <class EE> Vec<M, typename CNT<E>::template Result<EE>::Mul>
722  scalarMultiply(const EE& e) const {
723  Vec<M, typename CNT<E>::template Result<EE>::Mul> result;
724  for (int i=0; i<M; ++i) result[i] = (*this)[i] * e;
725  return result;
726  }
727  template <class EE> Vec<M, typename CNT<EE>::template Result<E>::Mul>
728  scalarMultiplyFromLeft(const EE& e) const {
729  Vec<M, typename CNT<EE>::template Result<E>::Mul> result;
730  for (int i=0; i<M; ++i) result[i] = e * (*this)[i];
731  return result;
732  }
733 
734  // TODO: should precalculate and store 1/e, while converting to Standard
735  // Numbers. Note that return type should change appropriately.
736  template <class EE> Vec<M, typename CNT<E>::template Result<EE>::Dvd>
737  scalarDivide(const EE& e) const {
738  Vec<M, typename CNT<E>::template Result<EE>::Dvd> result;
739  for (int i=0; i<M; ++i) result[i] = (*this)[i] / e;
740  return result;
741  }
742  template <class EE> Vec<M, typename CNT<EE>::template Result<E>::Dvd>
743  scalarDivideFromLeft(const EE& e) const {
744  Vec<M, typename CNT<EE>::template Result<E>::Dvd> result;
745  for (int i=0; i<M; ++i) result[i] = e / (*this)[i];
746  return result;
747  }
748 
749  template <class EE> Vec<M, typename CNT<E>::template Result<EE>::Add>
750  scalarAdd(const EE& e) const {
751  Vec<M, typename CNT<E>::template Result<EE>::Add> result;
752  for (int i=0; i<M; ++i) result[i] = (*this)[i] + e;
753  return result;
754  }
755  // Add is commutative, so no 'FromLeft'.
756 
757  template <class EE> Vec<M, typename CNT<E>::template Result<EE>::Sub>
758  scalarSubtract(const EE& e) const {
759  Vec<M, typename CNT<E>::template Result<EE>::Sub> result;
760  for (int i=0; i<M; ++i) result[i] = (*this)[i] - e;
761  return result;
762  }
763  template <class EE> Vec<M, typename CNT<EE>::template Result<E>::Sub>
764  scalarSubtractFromLeft(const EE& e) const {
765  Vec<M, typename CNT<EE>::template Result<E>::Sub> result;
766  for (int i=0; i<M; ++i) result[i] = e - (*this)[i];
767  return result;
768  }
769 
770  // Generic assignments for any element type not listed explicitly, including scalars.
771  // These are done repeatedly for each element and only work if the operation can
772  // be performed leaving the original element type.
773  template <class EE> Vec& operator =(const EE& e) {return scalarEq(e);}
774  template <class EE> Vec& operator+=(const EE& e) {return scalarPlusEq(e);}
775  template <class EE> Vec& operator-=(const EE& e) {return scalarMinusEq(e);}
776  template <class EE> Vec& operator*=(const EE& e) {return scalarTimesEq(e);}
777  template <class EE> Vec& operator/=(const EE& e) {return scalarDivideEq(e);}
778 
779  // Generalized element assignment & computed assignment methods. These will work
780  // for any assignment-compatible element, not just scalars.
781  template <class EE> Vec& scalarEq(const EE& ee)
782  { for(int i=0;i<M;++i) d[i*STRIDE] = ee; return *this; }
783  template <class EE> Vec& scalarPlusEq(const EE& ee)
784  { for(int i=0;i<M;++i) d[i*STRIDE] += ee; return *this; }
785  template <class EE> Vec& scalarMinusEq(const EE& ee)
786  { for(int i=0;i<M;++i) d[i*STRIDE] -= ee; return *this; }
787  template <class EE> Vec& scalarMinusEqFromLeft(const EE& ee)
788  { for(int i=0;i<M;++i) d[i*STRIDE] = ee - d[i*STRIDE]; return *this; }
789  template <class EE> Vec& scalarTimesEq(const EE& ee)
790  { for(int i=0;i<M;++i) d[i*STRIDE] *= ee; return *this; }
791  template <class EE> Vec& scalarTimesEqFromLeft(const EE& ee)
792  { for(int i=0;i<M;++i) d[i*STRIDE] = ee * d[i*STRIDE]; return *this; }
793  template <class EE> Vec& scalarDivideEq(const EE& ee)
794  { for(int i=0;i<M;++i) d[i*STRIDE] /= ee; return *this; }
795  template <class EE> Vec& scalarDivideEqFromLeft(const EE& ee)
796  { for(int i=0;i<M;++i) d[i*STRIDE] = ee / d[i*STRIDE]; return *this; }
797 
798  // Specialize for int to avoid warnings and ambiguities.
799  Vec& scalarEq(int ee) {return scalarEq(Precision(ee));}
800  Vec& scalarPlusEq(int ee) {return scalarPlusEq(Precision(ee));}
801  Vec& scalarMinusEq(int ee) {return scalarMinusEq(Precision(ee));}
802  Vec& scalarTimesEq(int ee) {return scalarTimesEq(Precision(ee));}
803  Vec& scalarDivideEq(int ee) {return scalarDivideEq(Precision(ee));}
804  Vec& scalarMinusEqFromLeft(int ee) {return scalarMinusEqFromLeft(Precision(ee));}
805  Vec& scalarTimesEqFromLeft(int ee) {return scalarTimesEqFromLeft(Precision(ee));}
806  Vec& scalarDivideEqFromLeft(int ee) {return scalarDivideEqFromLeft(Precision(ee));}
807 
810  void setToNaN() {
811  (*this) = CNT<ELT>::getNaN();
812  }
813 
815  void setToZero() {
816  (*this) = ELT(0);
817  }
818 
824  template <int MM>
825  const Vec<MM,ELT,STRIDE>& getSubVec(int i) const {
826  assert(0 <= i && i + MM <= M);
827  return Vec<MM,ELT,STRIDE>::getAs(&(*this)[i]);
828  }
834  template <int MM>
836  assert(0 <= i && i + MM <= M);
837  return Vec<MM,ELT,STRIDE>::updAs(&(*this)[i]);
838  }
839 
840 
844  template <int MM>
845  static const Vec& getSubVec(const Vec<MM,ELT,STRIDE>& v, int i) {
846  assert(0 <= i && i + M <= MM);
847  return getAs(&v[i]);
848  }
852  template <int MM>
853  static Vec& updSubVec(Vec<MM,ELT,STRIDE>& v, int i) {
854  assert(0 <= i && i + M <= MM);
855  return updAs(&v[i]);
856  }
857 
861  Vec<M-1,ELT,1> drop1(int p) const {
862  assert(0 <= p && p < M);
863  Vec<M-1,ELT,1> out;
864  int nxt=0;
865  for (int i=0; i<M-1; ++i, ++nxt) {
866  if (nxt==p) ++nxt; // skip the loser
867  out[i] = (*this)[nxt];
868  }
869  return out;
870  }
871 
875  template <class EE> Vec<M+1,ELT,1> append1(const EE& v) const {
876  Vec<M+1,ELT,1> out;
877  Vec<M,ELT,1>::updAs(&out[0]) = (*this);
878  out[M] = v;
879  return out;
880  }
881 
882 
888  template <class EE> Vec<M+1,ELT,1> insert1(int p, const EE& v) const {
889  assert(0 <= p && p <= M);
890  if (p==M) return append1(v);
891  Vec<M+1,ELT,1> out;
892  int nxt=0;
893  for (int i=0; i<M; ++i, ++nxt) {
894  if (i==p) out[nxt++] = v;
895  out[nxt] = (*this)[i];
896  }
897  return out;
898  }
899 
902  static const Vec& getAs(const ELT* p)
903  { return *reinterpret_cast<const Vec*>(p); }
906  static Vec& updAs(ELT* p)
907  { return *reinterpret_cast<Vec*>(p); }
908 
909 
914 
916  bool isNaN() const {
917  for (int i=0; i<M; ++i)
918  if (CNT<ELT>::isNaN((*this)[i]))
919  return true;
920  return false;
921  }
922 
925  bool isInf() const {
926  bool seenInf = false;
927  for (int i=0; i<M; ++i) {
928  const ELT& e = (*this)[i];
929  if (!CNT<ELT>::isFinite(e)) {
930  if (!CNT<ELT>::isInf(e))
931  return false; // something bad was found
932  seenInf = true;
933  }
934  }
935  return seenInf;
936  }
937 
940  bool isFinite() const {
941  for (int i=0; i<M; ++i)
942  if (!CNT<ELT>::isFinite((*this)[i]))
943  return false;
944  return true;
945  }
946 
950 
953  template <class E2, int RS2>
954  bool isNumericallyEqual(const Vec<M,E2,RS2>& v, double tol) const {
955  for (int i=0; i<M; ++i)
956  if (!CNT<ELT>::isNumericallyEqual((*this)[i], v[i], tol))
957  return false;
958  return true;
959  }
960 
964  template <class E2, int RS2>
965  bool isNumericallyEqual(const Vec<M,E2,RS2>& v) const {
966  const double tol = std::max(getDefaultTolerance(),v.getDefaultTolerance());
967  return isNumericallyEqual(v, tol);
968  }
969 
974  bool isNumericallyEqual
975  (const ELT& e,
976  double tol = getDefaultTolerance()) const
977  {
978  for (int i=0; i<M; ++i)
979  if (!CNT<ELT>::isNumericallyEqual((*this)[i], e, tol))
980  return false;
981  return true;
982  }
983 
984  // Functions to be used for Scripting in MATLAB and languages that do not support operator overloading
986  std::string toString() const {
987  std::stringstream stream;
988  stream << (*this);
989  return stream.str();
990  }
991 
993  void set(int i, const E& value)
994  { (*this)[i] = value; }
995 
997  const E& get(int i) const
998  { return operator[](i); }
999 
1000 private:
1001  // TODO: should be an array of scalars rather than elements to control
1002  // packing more carefully.
1003  ELT d[NActualElements]; // data
1004 };
1005 
1007 // Global operators involving two vectors. //
1008 // v+v, v-v, v==v, v!=v //
1010 
1011 // v3 = v1 + v2 where all v's have the same length M.
1012 template <int M, class E1, int S1, class E2, int S2> inline
1013 typename Vec<M,E1,S1>::template Result< Vec<M,E2,S2> >::Add
1014 operator+(const Vec<M,E1,S1>& l, const Vec<M,E2,S2>& r) {
1015  return Vec<M,E1,S1>::template Result< Vec<M,E2,S2> >
1016  ::AddOp::perform(l,r);
1017 }
1018 
1019 // v3 = v1 - v2, similar to +
1020 template <int M, class E1, int S1, class E2, int S2> inline
1021 typename Vec<M,E1,S1>::template Result< Vec<M,E2,S2> >::Sub
1022 operator-(const Vec<M,E1,S1>& l, const Vec<M,E2,S2>& r) {
1023  return Vec<M,E1,S1>::template Result< Vec<M,E2,S2> >
1024  ::SubOp::perform(l,r);
1025 }
1026 
1028 template <int M, class E1, int S1, class E2, int S2> inline bool
1030 { for (int i=0; i < M; ++i) if (l[i] != r[i]) return false;
1031  return true; }
1033 template <int M, class E1, int S1, class E2, int S2> inline bool
1034 operator!=(const Vec<M,E1,S1>& l, const Vec<M,E2,S2>& r) {return !(l==r);}
1035 
1037 template <int M, class E1, int S1, class E2> inline bool
1038 operator==(const Vec<M,E1,S1>& v, const E2& e)
1039 { for (int i=0; i < M; ++i) if (v[i] != e) return false;
1040  return true; }
1042 template <int M, class E1, int S1, class E2> inline bool
1043 operator!=(const Vec<M,E1,S1>& v, const E2& e) {return !(v==e);}
1044 
1046 template <int M, class E1, int S1, class E2, int S2> inline bool
1047 operator<(const Vec<M,E1,S1>& l, const Vec<M,E2,S2>& r)
1048 { for (int i=0; i < M; ++i) if (l[i] >= r[i]) return false;
1049  return true; }
1051 template <int M, class E1, int S1, class E2> inline bool
1052 operator<(const Vec<M,E1,S1>& v, const E2& e)
1053 { for (int i=0; i < M; ++i) if (v[i] >= e) return false;
1054  return true; }
1055 
1057 template <int M, class E1, int S1, class E2, int S2> inline bool
1059 { for (int i=0; i < M; ++i) if (l[i] <= r[i]) return false;
1060  return true; }
1062 template <int M, class E1, int S1, class E2> inline bool
1063 operator>(const Vec<M,E1,S1>& v, const E2& e)
1064 { for (int i=0; i < M; ++i) if (v[i] <= e) return false;
1065  return true; }
1066 
1069 template <int M, class E1, int S1, class E2, int S2> inline bool
1070 operator<=(const Vec<M,E1,S1>& l, const Vec<M,E2,S2>& r)
1071 { for (int i=0; i < M; ++i) if (l[i] > r[i]) return false;
1072  return true; }
1075 template <int M, class E1, int S1, class E2> inline bool
1076 operator<=(const Vec<M,E1,S1>& v, const E2& e)
1077 { for (int i=0; i < M; ++i) if (v[i] > e) return false;
1078  return true; }
1079 
1082 template <int M, class E1, int S1, class E2, int S2> inline bool
1084 { for (int i=0; i < M; ++i) if (l[i] < r[i]) return false;
1085  return true; }
1088 template <int M, class E1, int S1, class E2> inline bool
1089 operator>=(const Vec<M,E1,S1>& v, const E2& e)
1090 { for (int i=0; i < M; ++i) if (v[i] < e) return false;
1091  return true; }
1092 
1094 // Global operators involving a vector and a scalar. //
1096 
1097 // I haven't been able to figure out a nice way to templatize for the
1098 // built-in reals without introducing a lot of unwanted type matches
1099 // as well. So we'll just grind them out explicitly here.
1100 
1101 // SCALAR MULTIPLY
1102 
1103 // v = v*real, real*v
1104 template <int M, class E, int S> inline
1105 typename Vec<M,E,S>::template Result<float>::Mul
1106 operator*(const Vec<M,E,S>& l, const float& r)
1107  { return Vec<M,E,S>::template Result<float>::MulOp::perform(l,r); }
1108 template <int M, class E, int S> inline
1109 typename Vec<M,E,S>::template Result<float>::Mul
1110 operator*(const float& l, const Vec<M,E,S>& r) {return r*l;}
1111 
1112 template <int M, class E, int S> inline
1113 typename Vec<M,E,S>::template Result<double>::Mul
1114 operator*(const Vec<M,E,S>& l, const double& r)
1115  { return Vec<M,E,S>::template Result<double>::MulOp::perform(l,r); }
1116 template <int M, class E, int S> inline
1117 typename Vec<M,E,S>::template Result<double>::Mul
1118 operator*(const double& l, const Vec<M,E,S>& r) {return r*l;}
1119 
1120 template <int M, class E, int S> inline
1121 typename Vec<M,E,S>::template Result<long double>::Mul
1122 operator*(const Vec<M,E,S>& l, const long double& r)
1123  { return Vec<M,E,S>::template Result<long double>::MulOp::perform(l,r); }
1124 template <int M, class E, int S> inline
1125 typename Vec<M,E,S>::template Result<long double>::Mul
1126 operator*(const long double& l, const Vec<M,E,S>& r) {return r*l;}
1127 
1128 // v = v*int, int*v -- just convert int to v's precision float
1129 template <int M, class E, int S> inline
1130 typename Vec<M,E,S>::template Result<typename CNT<E>::Precision>::Mul
1131 operator*(const Vec<M,E,S>& l, int r) {return l * (typename CNT<E>::Precision)r;}
1132 template <int M, class E, int S> inline
1133 typename Vec<M,E,S>::template Result<typename CNT<E>::Precision>::Mul
1134 operator*(int l, const Vec<M,E,S>& r) {return r * (typename CNT<E>::Precision)l;}
1135 
1136 // Complex, conjugate, and negator are all easy to templatize.
1137 
1138 // v = v*complex, complex*v
1139 template <int M, class E, int S, class R> inline
1140 typename Vec<M,E,S>::template Result<std::complex<R> >::Mul
1141 operator*(const Vec<M,E,S>& l, const std::complex<R>& r)
1142  { return Vec<M,E,S>::template Result<std::complex<R> >::MulOp::perform(l,r); }
1143 template <int M, class E, int S, class R> inline
1144 typename Vec<M,E,S>::template Result<std::complex<R> >::Mul
1145 operator*(const std::complex<R>& l, const Vec<M,E,S>& r) {return r*l;}
1146 
1147 // v = v*conjugate, conjugate*v (convert conjugate->complex)
1148 template <int M, class E, int S, class R> inline
1149 typename Vec<M,E,S>::template Result<std::complex<R> >::Mul
1150 operator*(const Vec<M,E,S>& l, const conjugate<R>& r) {return l*(std::complex<R>)r;}
1151 template <int M, class E, int S, class R> inline
1152 typename Vec<M,E,S>::template Result<std::complex<R> >::Mul
1153 operator*(const conjugate<R>& l, const Vec<M,E,S>& r) {return r*(std::complex<R>)l;}
1154 
1155 // v = v*negator, negator*v: convert negator to standard number
1156 template <int M, class E, int S, class R> inline
1157 typename Vec<M,E,S>::template Result<typename negator<R>::StdNumber>::Mul
1158 operator*(const Vec<M,E,S>& l, const negator<R>& r) {return l * (typename negator<R>::StdNumber)(R)r;}
1159 template <int M, class E, int S, class R> inline
1160 typename Vec<M,E,S>::template Result<typename negator<R>::StdNumber>::Mul
1161 operator*(const negator<R>& l, const Vec<M,E,S>& r) {return r * (typename negator<R>::StdNumber)(R)l;}
1162 
1163 
1164 // SCALAR DIVIDE. This is a scalar operation when the scalar is on the right,
1165 // but when it is on the left it means scalar * pseudoInverse(vec), which is
1166 // a row.
1167 
1168 // v = v/real, real/v
1169 template <int M, class E, int S> inline
1170 typename Vec<M,E,S>::template Result<float>::Dvd
1171 operator/(const Vec<M,E,S>& l, const float& r)
1172  { return Vec<M,E,S>::template Result<float>::DvdOp::perform(l,r); }
1173 template <int M, class E, int S> inline
1174 typename CNT<float>::template Result<Vec<M,E,S> >::Dvd
1175 operator/(const float& l, const Vec<M,E,S>& r)
1176  { return CNT<float>::template Result<Vec<M,E,S> >::DvdOp::perform(l,r); }
1177 
1178 template <int M, class E, int S> inline
1179 typename Vec<M,E,S>::template Result<double>::Dvd
1180 operator/(const Vec<M,E,S>& l, const double& r)
1181  { return Vec<M,E,S>::template Result<double>::DvdOp::perform(l,r); }
1182 template <int M, class E, int S> inline
1183 typename CNT<double>::template Result<Vec<M,E,S> >::Dvd
1184 operator/(const double& l, const Vec<M,E,S>& r)
1185  { return CNT<double>::template Result<Vec<M,E,S> >::DvdOp::perform(l,r); }
1186 
1187 template <int M, class E, int S> inline
1188 typename Vec<M,E,S>::template Result<long double>::Dvd
1189 operator/(const Vec<M,E,S>& l, const long double& r)
1190  { return Vec<M,E,S>::template Result<long double>::DvdOp::perform(l,r); }
1191 template <int M, class E, int S> inline
1192 typename CNT<long double>::template Result<Vec<M,E,S> >::Dvd
1193 operator/(const long double& l, const Vec<M,E,S>& r)
1194  { return CNT<long double>::template Result<Vec<M,E,S> >::DvdOp::perform(l,r); }
1195 
1196 // v = v/int, int/v -- just convert int to v's precision float
1197 template <int M, class E, int S> inline
1198 typename Vec<M,E,S>::template Result<typename CNT<E>::Precision>::Dvd
1199 operator/(const Vec<M,E,S>& l, int r) {return l / (typename CNT<E>::Precision)r;}
1200 template <int M, class E, int S> inline
1201 typename CNT<typename CNT<E>::Precision>::template Result<Vec<M,E,S> >::Dvd
1202 operator/(int l, const Vec<M,E,S>& r) {return (typename CNT<E>::Precision)l / r;}
1203 
1204 
1205 // Complex, conjugate, and negator are all easy to templatize.
1206 
1207 // v = v/complex, complex/v
1208 template <int M, class E, int S, class R> inline
1209 typename Vec<M,E,S>::template Result<std::complex<R> >::Dvd
1210 operator/(const Vec<M,E,S>& l, const std::complex<R>& r)
1211  { return Vec<M,E,S>::template Result<std::complex<R> >::DvdOp::perform(l,r); }
1212 template <int M, class E, int S, class R> inline
1213 typename CNT<std::complex<R> >::template Result<Vec<M,E,S> >::Dvd
1214 operator/(const std::complex<R>& l, const Vec<M,E,S>& r)
1215  { return CNT<std::complex<R> >::template Result<Vec<M,E,S> >::DvdOp::perform(l,r); }
1216 
1217 // v = v/conjugate, conjugate/v (convert conjugate->complex)
1218 template <int M, class E, int S, class R> inline
1219 typename Vec<M,E,S>::template Result<std::complex<R> >::Dvd
1220 operator/(const Vec<M,E,S>& l, const conjugate<R>& r) {return l/(std::complex<R>)r;}
1221 template <int M, class E, int S, class R> inline
1222 typename CNT<std::complex<R> >::template Result<Vec<M,E,S> >::Dvd
1223 operator/(const conjugate<R>& l, const Vec<M,E,S>& r) {return (std::complex<R>)l/r;}
1224 
1225 // v = v/negator, negator/v: convert negator to number
1226 template <int M, class E, int S, class R> inline
1227 typename Vec<M,E,S>::template Result<typename negator<R>::StdNumber>::Dvd
1228 operator/(const Vec<M,E,S>& l, const negator<R>& r) {return l/(typename negator<R>::StdNumber)(R)r;}
1229 template <int M, class E, int S, class R> inline
1230 typename CNT<R>::template Result<Vec<M,E,S> >::Dvd
1231 operator/(const negator<R>& l, const Vec<M,E,S>& r) {return (typename negator<R>::StdNumber)(R)l/r;}
1232 
1233 
1234 // Add and subtract are odd as scalar ops. They behave as though the
1235 // scalar stands for a vector each of whose elements is that scalar,
1236 // and then a normal vector add or subtract is done.
1237 
1238 // SCALAR ADD
1239 
1240 // v = v+real, real+v
1241 template <int M, class E, int S> inline
1242 typename Vec<M,E,S>::template Result<float>::Add
1243 operator+(const Vec<M,E,S>& l, const float& r)
1244  { return Vec<M,E,S>::template Result<float>::AddOp::perform(l,r); }
1245 template <int M, class E, int S> inline
1246 typename Vec<M,E,S>::template Result<float>::Add
1247 operator+(const float& l, const Vec<M,E,S>& r) {return r+l;}
1248 
1249 template <int M, class E, int S> inline
1250 typename Vec<M,E,S>::template Result<double>::Add
1251 operator+(const Vec<M,E,S>& l, const double& r)
1252  { return Vec<M,E,S>::template Result<double>::AddOp::perform(l,r); }
1253 template <int M, class E, int S> inline
1254 typename Vec<M,E,S>::template Result<double>::Add
1255 operator+(const double& l, const Vec<M,E,S>& r) {return r+l;}
1256 
1257 template <int M, class E, int S> inline
1258 typename Vec<M,E,S>::template Result<long double>::Add
1259 operator+(const Vec<M,E,S>& l, const long double& r)
1260  { return Vec<M,E,S>::template Result<long double>::AddOp::perform(l,r); }
1261 template <int M, class E, int S> inline
1262 typename Vec<M,E,S>::template Result<long double>::Add
1263 operator+(const long double& l, const Vec<M,E,S>& r) {return r+l;}
1264 
1265 // v = v+int, int+v -- just convert int to v's precision float
1266 template <int M, class E, int S> inline
1267 typename Vec<M,E,S>::template Result<typename CNT<E>::Precision>::Add
1268 operator+(const Vec<M,E,S>& l, int r) {return l + (typename CNT<E>::Precision)r;}
1269 template <int M, class E, int S> inline
1270 typename Vec<M,E,S>::template Result<typename CNT<E>::Precision>::Add
1271 operator+(int l, const Vec<M,E,S>& r) {return r + (typename CNT<E>::Precision)l;}
1272 
1273 // Complex, conjugate, and negator are all easy to templatize.
1274 
1275 // v = v+complex, complex+v
1276 template <int M, class E, int S, class R> inline
1277 typename Vec<M,E,S>::template Result<std::complex<R> >::Add
1278 operator+(const Vec<M,E,S>& l, const std::complex<R>& r)
1279  { return Vec<M,E,S>::template Result<std::complex<R> >::AddOp::perform(l,r); }
1280 template <int M, class E, int S, class R> inline
1281 typename Vec<M,E,S>::template Result<std::complex<R> >::Add
1282 operator+(const std::complex<R>& l, const Vec<M,E,S>& r) {return r+l;}
1283 
1284 // v = v+conjugate, conjugate+v (convert conjugate->complex)
1285 template <int M, class E, int S, class R> inline
1286 typename Vec<M,E,S>::template Result<std::complex<R> >::Add
1287 operator+(const Vec<M,E,S>& l, const conjugate<R>& r) {return l+(std::complex<R>)r;}
1288 template <int M, class E, int S, class R> inline
1289 typename Vec<M,E,S>::template Result<std::complex<R> >::Add
1290 operator+(const conjugate<R>& l, const Vec<M,E,S>& r) {return r+(std::complex<R>)l;}
1291 
1292 // v = v+negator, negator+v: convert negator to standard number
1293 template <int M, class E, int S, class R> inline
1294 typename Vec<M,E,S>::template Result<typename negator<R>::StdNumber>::Add
1295 operator+(const Vec<M,E,S>& l, const negator<R>& r) {return l + (typename negator<R>::StdNumber)(R)r;}
1296 template <int M, class E, int S, class R> inline
1297 typename Vec<M,E,S>::template Result<typename negator<R>::StdNumber>::Add
1298 operator+(const negator<R>& l, const Vec<M,E,S>& r) {return r + (typename negator<R>::StdNumber)(R)l;}
1299 
1300 // SCALAR SUBTRACT -- careful, not commutative.
1301 
1302 // v = v-real, real-v
1303 template <int M, class E, int S> inline
1304 typename Vec<M,E,S>::template Result<float>::Sub
1305 operator-(const Vec<M,E,S>& l, const float& r)
1306  { return Vec<M,E,S>::template Result<float>::SubOp::perform(l,r); }
1307 template <int M, class E, int S> inline
1308 typename CNT<float>::template Result<Vec<M,E,S> >::Sub
1309 operator-(const float& l, const Vec<M,E,S>& r)
1310  { return CNT<float>::template Result<Vec<M,E,S> >::SubOp::perform(l,r); }
1311 
1312 template <int M, class E, int S> inline
1313 typename Vec<M,E,S>::template Result<double>::Sub
1314 operator-(const Vec<M,E,S>& l, const double& r)
1315  { return Vec<M,E,S>::template Result<double>::SubOp::perform(l,r); }
1316 template <int M, class E, int S> inline
1317 typename CNT<double>::template Result<Vec<M,E,S> >::Sub
1318 operator-(const double& l, const Vec<M,E,S>& r)
1319  { return CNT<double>::template Result<Vec<M,E,S> >::SubOp::perform(l,r); }
1320 
1321 template <int M, class E, int S> inline
1322 typename Vec<M,E,S>::template Result<long double>::Sub
1323 operator-(const Vec<M,E,S>& l, const long double& r)
1324  { return Vec<M,E,S>::template Result<long double>::SubOp::perform(l,r); }
1325 template <int M, class E, int S> inline
1326 typename CNT<long double>::template Result<Vec<M,E,S> >::Sub
1327 operator-(const long double& l, const Vec<M,E,S>& r)
1328  { return CNT<long double>::template Result<Vec<M,E,S> >::SubOp::perform(l,r); }
1329 
1330 // v = v-int, int-v // just convert int to v's precision float
1331 template <int M, class E, int S> inline
1332 typename Vec<M,E,S>::template Result<typename CNT<E>::Precision>::Sub
1333 operator-(const Vec<M,E,S>& l, int r) {return l - (typename CNT<E>::Precision)r;}
1334 template <int M, class E, int S> inline
1335 typename CNT<typename CNT<E>::Precision>::template Result<Vec<M,E,S> >::Sub
1336 operator-(int l, const Vec<M,E,S>& r) {return (typename CNT<E>::Precision)l - r;}
1337 
1338 
1339 // Complex, conjugate, and negator are all easy to templatize.
1340 
1341 // v = v-complex, complex-v
1342 template <int M, class E, int S, class R> inline
1343 typename Vec<M,E,S>::template Result<std::complex<R> >::Sub
1344 operator-(const Vec<M,E,S>& l, const std::complex<R>& r)
1345  { return Vec<M,E,S>::template Result<std::complex<R> >::SubOp::perform(l,r); }
1346 template <int M, class E, int S, class R> inline
1347 typename CNT<std::complex<R> >::template Result<Vec<M,E,S> >::Sub
1348 operator-(const std::complex<R>& l, const Vec<M,E,S>& r)
1349  { return CNT<std::complex<R> >::template Result<Vec<M,E,S> >::SubOp::perform(l,r); }
1350 
1351 // v = v-conjugate, conjugate-v (convert conjugate->complex)
1352 template <int M, class E, int S, class R> inline
1353 typename Vec<M,E,S>::template Result<std::complex<R> >::Sub
1354 operator-(const Vec<M,E,S>& l, const conjugate<R>& r) {return l-(std::complex<R>)r;}
1355 template <int M, class E, int S, class R> inline
1356 typename CNT<std::complex<R> >::template Result<Vec<M,E,S> >::Sub
1357 operator-(const conjugate<R>& l, const Vec<M,E,S>& r) {return (std::complex<R>)l-r;}
1358 
1359 // v = v-negator, negator-v: convert negator to standard number
1360 template <int M, class E, int S, class R> inline
1361 typename Vec<M,E,S>::template Result<typename negator<R>::StdNumber>::Sub
1362 operator-(const Vec<M,E,S>& l, const negator<R>& r) {return l-(typename negator<R>::StdNumber)(R)r;}
1363 template <int M, class E, int S, class R> inline
1364 typename CNT<R>::template Result<Vec<M,E,S> >::Sub
1365 operator-(const negator<R>& l, const Vec<M,E,S>& r) {return (typename negator<R>::StdNumber)(R)l-r;}
1366 
1367 // Vec I/O
1368 template <int M, class E, int S, class CHAR, class TRAITS> inline
1369 std::basic_ostream<CHAR,TRAITS>&
1370 operator<<(std::basic_ostream<CHAR,TRAITS>& o, const Vec<M,E,S>& v) {
1371  o << "~[" << v[0]; for(int i=1;i<M;++i) o<<','<<v[i]; o<<']'; return o;
1372 }
1373 
1376 template <int M, class E, int S, class CHAR, class TRAITS> inline
1377 std::basic_istream<CHAR,TRAITS>&
1378 operator>>(std::basic_istream<CHAR,TRAITS>& is, Vec<M,E,S>& v) {
1379  CHAR tilde;
1380  is >> tilde; if (is.fail()) return is;
1381  if (tilde != CHAR('~')) {
1382  tilde = CHAR(0);
1383  is.unget(); if (is.fail()) return is;
1384  }
1385 
1386  CHAR openBracket, closeBracket;
1387  is >> openBracket; if (is.fail()) return is;
1388  if (openBracket==CHAR('('))
1389  closeBracket = CHAR(')');
1390  else if (openBracket==CHAR('['))
1391  closeBracket = CHAR(']');
1392  else {
1393  closeBracket = CHAR(0);
1394  is.unget(); if (is.fail()) return is;
1395  }
1396 
1397  // If we saw a "~" but then we didn't see any brackets, that's an
1398  // error. Set the fail bit and return.
1399  if (tilde != CHAR(0) && closeBracket == CHAR(0)) {
1400  is.setstate( std::ios::failbit );
1401  return is;
1402  }
1403 
1404  for (int i=0; i < M; ++i) {
1405  is >> v[i];
1406  if (is.fail()) return is;
1407  if (i != M-1) {
1408  CHAR c; is >> c; if (is.fail()) return is;
1409  if (c != ',') is.unget();
1410  if (is.fail()) return is;
1411  }
1412  }
1413 
1414  // Get the closing bracket if there was an opening one. If we don't
1415  // see the expected character we'll set the fail bit in the istream.
1416  if (closeBracket != CHAR(0)) {
1417  CHAR closer; is >> closer; if (is.fail()) return is;
1418  if (closer != closeBracket) {
1419  is.unget(); if (is.fail()) return is;
1420  is.setstate( std::ios::failbit );
1421  }
1422  }
1423 
1424  return is;
1425 }
1426 
1427 } //namespace SimTK
1428 
1429 
1430 #endif //SimTK_SIMMATRIX_SMALLMATRIX_VEC_H_
Matrix_< E > operator/(const MatrixBase< E > &l, const typename CNT< E >::StdNumber &r)
Definition: BigMatrix.h:613
TImag & imag()
Recast to show only the imaginary portion of this Vec and return a writable reference.
Definition: Vec.h:697
PhiMatrixTranspose transpose(const PhiMatrix &phi)
Definition: SpatialAlgebra.h:720
std::string toString() const
Print Vec into a string and return it.
Definition: Vec.h:986
bool isFinite() const
Return true if no element of this Vec contains an Infinity or a NaN anywhere.
Definition: Vec.h:940
TAbs abs() const
Elementwise absolute value; that is, the return value has the same dimension as this Vec but with eac...
Definition: Vec.h:345
CNT< E >::TSqHermT ESqHermT
Type of the expression ~E*E (default vector and matrix square; symmetric).
Definition: Vec.h:212
K::ScalarNormSq ScalarNormSq
Definition: CompositeNumericalTypes.h:166
SubOp::Type Sub
Definition: Vec.h:408
static Vec< M, ELT, 1 > getNaN()
Return a Vec of the same length and element type as this one but with all elements set to NaN...
Definition: Vec.h:913
K::ULessScalar ULessScalar
Definition: CompositeNumericalTypes.h:161
TNeg & operator-()
Recast to negated type and return a writable reference; writing to this will cause the negated result...
Definition: Vec.h:643
static int size()
The number of elements in this Vec (note that stride does not affect this number.) ...
Definition: Vec.h:316
Vec(const E &e0, const E &e1, const E &e2, const E &e3, const E &e4, const E &e5, const E &e6, const E &e7, const E &e8)
Definition: Vec.h:506
Vec< M, P > Type
Definition: Vec.h:417
K::TReal TReal
Definition: CompositeNumericalTypes.h:141
CNT< E >::TSqrt ESqrt
Type required to hold the result of sqrt(E).
Definition: Vec.h:216
const TImag & imag() const
Return a reference to the imaginary portion of this Vec if it has complex elements; otherwise the typ...
Definition: Vec.h:690
ScalarNormSq scalarNormSqr() const
Scalar norm square is sum( conjugate squares of all underlying scalars ), where conjugate square of s...
Definition: Vec.h:325
Vec< M, typename CNT< E >::template Result< EE >::Dvd > elementwiseDivide(const Vec< M, EE, SS > &r) const
Elementwise divide (Matlab " ./ " operator).
Definition: Vec.h:585
This is a small, fixed-size symmetric or Hermitian matrix designed for no-overhead inline computation...
Definition: SimTKcommon/include/SimTKcommon/internal/common.h:608
const E & operator[](int i) const
Select an element of this Vec and return a const reference to it.
Definition: Vec.h:594
SymMat< M, ESqTHerm > TSqTHerm
These compile-time constants are required of every Composite Numerical Type (CNT).
Definition: Vec.h:302
Vec< M, typename CNT< EE >::template Result< E >::Sub > scalarSubtractFromLeft(const EE &e) const
Definition: Vec.h:764
Vec & scalarEq(const EE &ee)
Definition: Vec.h:781
CNT< E >::TImag EImag
Type showing the imaginary part of an element of this Vec as real, if elements are complex; otherwise...
Definition: Vec.h:202
MulOp::Type Mul
Definition: Vec.h:388
CNT< E >::TReal EReal
Type showing just the real part of an element of this Vec if elements are complex; otherwise just the...
Definition: Vec.h:198
TSqrt sqrt() const
Elementwise square root; that is, the return value has the same length as this Vec but with each elem...
Definition: Vec.h:335
Vec & scalarTimesEq(const EE &ee)
Definition: Vec.h:789
Vec & scalarTimesEqFromLeft(int ee)
Definition: Vec.h:805
This is the top-level SimTK namespace into which all SimTK names are placed to avoid collision with o...
Definition: Assembler.h:37
bool isNaN() const
Return true if any element of this Vec contains a NaN anywhere.
Definition: Vec.h:916
NTraits< N >::StdNumber StdNumber
Definition: negator.h:107
SimTK::conjugate<R> should be instantiated only for float, double, long double.
Definition: String.h:45
K::TSqrt TSqrt
Definition: CompositeNumericalTypes.h:154
Vec< M-1, ELT, 1 > drop1(int p) const
Return a vector one smaller than this one by dropping the element at the indicated position p...
Definition: Vec.h:861
TWithoutNegator & updCastAwayNegatorIfAny()
Recast to remove negators from this Vec&#39;s type if present and return a writable reference.
Definition: Vec.h:710
Vec(const E &e)
Construction from a single value of this Vec&#39;s element type assigns that value to each element...
Definition: Vec.h:468
TInvert invert() const
This method is not supported for Vec objects.
Definition: Vec.h:633
static TSqrt sqrt(const K &t)
Definition: CompositeNumericalTypes.h:239
static Vec & updAs(ELT *p)
Recast a writable ordinary C++ array E[] to a writable Vec<M,E,S>; assumes compatible length...
Definition: Vec.h:906
Vec< M, ESqrt, 1 > TSqrt
These compile-time constants are required of every Composite Numerical Type (CNT).
Definition: Vec.h:295
Vec & operator=(const EE *p)
Assignment to a pointer to elements of any type EE assumes we&#39;re pointing at a C++ array of EE&#39;s of t...
Definition: Vec.h:521
Definition: Vec.h:375
TNeg & updNegate()
Non-operator version of unary negation; recasts and returns a writable reference. ...
Definition: Vec.h:657
K::Scalar Scalar
Definition: CompositeNumericalTypes.h:160
ESqHermT TSqHermT
These compile-time constants are required of every Composite Numerical Type (CNT).
Definition: Vec.h:301
Vec & operator+=(const EE &e)
Definition: Vec.h:774
CNT< E >::TNormalize ENormalize
Packed type that can hold the value returned from normalize(E).
Definition: Vec.h:226
K::TNormalize TNormalize
Definition: CompositeNumericalTypes.h:158
Vec< M, typename CNT< EE >::template Result< E >::Mul > scalarMultiplyFromLeft(const EE &e) const
Definition: Vec.h:728
Matrix_< typename CNT< E1 >::template Result< E2 >::Sub > operator-(const MatrixBase< E1 > &l, const MatrixBase< E2 > &r)
Definition: BigMatrix.h:584
E TElement
Element type of this Vec.
Definition: Vec.h:287
Vec< M, typename CNT< E >::template Result< EE >::Mul > elementwiseMultiply(const Vec< M, EE, SS > &r) const
Elementwise multiply (Matlab " .* " operator).
Definition: Vec.h:578
Vec(const E &e0, const E &e1, const E &e2, const E &e3, const E &e4, const E &e5, const E &e6)
Definition: Vec.h:500
static int nrow()
The number of rows in a Vec is the number of elements.
Definition: Vec.h:318
E & operator[](int i)
Select an element of this Vec and return a writable reference to it.
Definition: Vec.h:602
Vec(const E &e0, const E &e1, const E &e2, const E &e3, const E &e4, const E &e5, const E &e6, const E &e7)
Definition: Vec.h:503
K::TImag TImag
Definition: CompositeNumericalTypes.h:142
CNT< E >::StdNumber EStdNumber
These compile-time constants are required of every Composite Numerical Type (CNT).
Definition: Vec.h:231
Vec(const E &e0, const E &e1)
Construct a Vec<2,E> from two elements of type E, etc.
Definition: Vec.h:488
TStandard standardize() const
Return a copy of this Vec but with the underlying scalar type converted (if necessary) to one of the ...
Definition: Vec.h:355
Vec< M, typename CNT< E >::template Result< P >::Mul, 1 > Mul
Definition: Vec.h:376
Vec & operator-=(const EE &e)
Definition: Vec.h:775
EStdNumber StdNumber
These compile-time constants are required of every Composite Numerical Type (CNT).
Definition: Vec.h:309
AddCNTs< M, 1, ArgDepth, Vec, ColSpacing, RowSpacing, CNT< P >::NRows, CNT< P >::NCols, CNT< P >::ArgDepth, P, CNT< P >::ColSpacing, CNT< P >::RowSpacing > AddOp
Definition: Vec.h:402
std::basic_istream< CHAR, TRAITS > & operator>>(std::basic_istream< CHAR, TRAITS > &is, conjugate< R > &c)
Definition: conjugate.h:800
Vec(const EE *p)
Construction from a pointer to elements of any type EE assumes we&#39;re pointing at a C++ array of EE&#39;s ...
Definition: Vec.h:514
Vec< M, EAbs, 1 > TAbs
These compile-time constants are required of every Composite Numerical Type (CNT).
Definition: Vec.h:296
negator<N>, where N is a number type (real, complex, conjugate), is represented in memory identically...
Definition: String.h:44
Definition: CompositeNumericalTypes.h:120
static double getDefaultTolerance()
Definition: CompositeNumericalTypes.h:269
CNT< E >::THerm EHerm
Type of the Hermitian transpose of an element of this Vec.
Definition: Vec.h:207
Vec< M+1, ELT, 1 > insert1(int p, const EE &v) const
Return a vector one larger than this one by inserting an element before the indicated one...
Definition: Vec.h:888
const TNeg & operator-() const
Unary minus recasts this Vec to a type that has the opposite interpretation of the sign but is otherw...
Definition: Vec.h:640
SubCNTs< M, 1, ArgDepth, Vec, ColSpacing, RowSpacing, CNT< P >::NRows, CNT< P >::NCols, CNT< P >::ArgDepth, P, CNT< P >::ColSpacing, CNT< P >::RowSpacing > SubOp
Definition: Vec.h:407
TPosTrans & updPositionalTranspose()
Positional transpose returning a writable reference.
Definition: Vec.h:672
Vec(const Vec &src)
Copy constructor copies the logically-included elements from the source Vec; gaps due to stride are n...
Definition: Vec.h:436
bool operator==(const PhiMatrix &p1, const PhiMatrix &p2)
Definition: SpatialAlgebra.h:774
CNT< E >::TStandard EStandard
Return type of standardize(E) method; a packed type that can hold the value of an element after elimi...
Definition: Vec.h:221
static TStandard standardize(const K &t)
Definition: CompositeNumericalTypes.h:241
Vec(const E &e0, const E &e1, const E &e2, const E &e3)
Definition: Vec.h:492
Vec< M, EReal, STRIDE *CNT< E >::RealStrideFactor > TReal
Type of this Vec cast to show only the real part of its element; this might affect the stride...
Definition: Vec.h:273
CNT< E >::Scalar EScalar
These compile-time constants are required of every Composite Numerical Type (CNT).
Definition: Vec.h:228
EPrecision Precision
These compile-time constants are required of every Composite Numerical Type (CNT).
Definition: Vec.h:310
CNT< E >::TPosTrans EPosTrans
Type of a positional transpose of an element of this Vec.
Definition: Vec.h:209
CNT< E >::TSqTHerm ESqTHerm
Type of the expression E*~E ("row square"; symmetric).
Definition: Vec.h:214
Vec< M, typename CNT< E >::template Result< P >::Add, 1 > Add
Definition: Vec.h:378
Vec< M, EStandard, 1 > TStandard
These compile-time constants are required of every Composite Numerical Type (CNT).
Definition: Vec.h:297
Vec & scalarMinusEqFromLeft(int ee)
Definition: Vec.h:804
Vec & scalarPlusEq(const EE &ee)
Definition: Vec.h:783
const E & operator()(int i) const
Same as const operator[] above.
Definition: Vec.h:597
Vec & operator+=(const Vec< M, negator< EE >, SS > &r)
Add in a conforming Vec, of any negated element type and stride, provided that the element types are ...
Definition: Vec.h:536
Vec< M, E, STRIDE > T
The type of this Vec.
Definition: Vec.h:263
CNT< E >::ScalarNormSq EScalarNormSq
These compile-time constants are required of every Composite Numerical Type (CNT).
Definition: Vec.h:233
Mat< M, M, typename CNT< E >::template Result< EE >::Mul > conformingMultiply(const Row< M, EE, SS > &r) const
Same as outer product (m = col*row) – use operator* or outer() instead.
Definition: Vec.h:570
ELEM sum(const VectorBase< ELEM > &v)
Definition: VectorMath.h:147
K::TSqTHerm TSqTHerm
Definition: CompositeNumericalTypes.h:147
Vec(const Vec< M, ENeg, SS > &src)
This is an implicit conversion from a Vec of the same length and negated element type (possibly with ...
Definition: Vec.h:456
EStandard sum() const
Sum just adds up all the elements into a single return element that is the same type as this Vec&#39;s el...
Definition: Vec.h:364
This is a fixed-length column vector designed for no-overhead inline computation. ...
Definition: SimTKcommon/include/SimTKcommon/internal/common.h:605
const Vec & operator+() const
Unary plus does nothing.
Definition: Vec.h:636
Row< M, EInvert, 1 > TInvert
These compile-time constants are required of every Composite Numerical Type (CNT).
Definition: Vec.h:298
const TNeg & negate() const
Non-operator version of unary negation; just a recast.
Definition: Vec.h:654
Vec(const Vec< M, EE, SS > &src)
Construct a Vec from a Vec of the same length, with any stride.
Definition: Vec.h:462
const TReal & real() const
Return a reference to the real portion of this Vec if it has complex elements; otherwise the type doe...
Definition: Vec.h:679
const THerm & operator~() const
The Hermitian transpose operator recasts this Vec to a type that specifies the opposite storage order...
Definition: Vec.h:647
Vec< M, typename CNT< E >::template Result< EE >::Mul > scalarMultiply(const EE &e) const
Definition: Vec.h:722
static double getDefaultTolerance()
For approximate comparisions, the default tolerance to use for a vector is the same as its elements&#39; ...
Definition: Vec.h:949
K::Precision Precision
Definition: CompositeNumericalTypes.h:164
Vec & scalarEq(int ee)
Definition: Vec.h:799
Row< M, EHerm, STRIDE > THerm
Type of this Vec after casting to its Hermitian transpose; that is, the Vec turns into a Row and each...
Definition: Vec.h:282
ENumber Number
These compile-time constants are required of every Composite Numerical Type (CNT).
Definition: Vec.h:308
void setToNaN()
Set every scalar in this Vec to NaN; this is the default initial value in Debug builds, but not in Release.
Definition: Vec.h:810
Matrix_< E > operator*(const MatrixBase< E > &l, const typename CNT< E >::StdNumber &r)
Definition: BigMatrix.h:605
Vec< M, typename CNT< E >::template Result< EE >::Add > scalarAdd(const EE &e) const
Definition: Vec.h:750
void setToZero()
Set every scalar in this Vec to zero.
Definition: Vec.h:815
void elementwiseDivide(const Row< 1, E1, S1 > &r1, const Row< 1, E2, S2 > &r2, Row< 1, typename CNT< E1 >::template Result< E2 >::Dvd > &result)
Definition: Row.h:90
K::TInvert TInvert
Definition: CompositeNumericalTypes.h:157
THerm & operator~()
Recast to Hermitian transposed type and return a writable reference; the effect is that writing to el...
Definition: Vec.h:651
bool isNumericallyEqual(const Vec< M, E2, RS2 > &v, double tol) const
Test whether this vector is numerically equal to some other vector with the same shape, using a specified tolerance.
Definition: Vec.h:954
Vec & operator=(const Vec &src)
Copy assignment operator copies the logically-included elements from the source Vec; gaps due to stri...
Definition: Vec.h:443
Vec(const E &e0, const E &e1, const E &e2)
Definition: Vec.h:490
MulOpNonConforming::Type MulNon
Definition: Vec.h:393
ELEM max(const VectorBase< ELEM > &v)
Definition: VectorMath.h:251
Vec & scalarMinusEq(int ee)
Definition: Vec.h:801
Shape-preserving element substitution (always packed).
Definition: Vec.h:416
Vec< M, typename CNT< E >::template Result< EE >::Add > conformingAdd(const Vec< M, EE, SS > &r) const
Vector addition – use operator+ instead.
Definition: Vec.h:554
MulCNTsNonConforming< M, 1, ArgDepth, Vec, ColSpacing, RowSpacing, CNT< P >::NRows, CNT< P >::NCols, CNT< P >::ArgDepth, P, CNT< P >::ColSpacing, CNT< P >::RowSpacing > MulOpNonConforming
Definition: Vec.h:392
TNormalize normalize() const
If the elements of this Vec are scalars, the result is what you get by dividing each element by the n...
Definition: Vec.h:621
Vec & operator/=(const EE &e)
Definition: Vec.h:777
static const Vec & getSubVec(const Vec< MM, ELT, STRIDE > &v, int i)
Extract a subvector of type Vec from a longer one that has the same element type and stride...
Definition: Vec.h:845
bool operator>=(const Row< N, E1, S1 > &l, const Row< N, E2, S2 > &r)
bool = v1[i] >= v2[i], for all elements i This is not the same as !(v1<v2).
Definition: Row.h:858
CNT< E >::TWithoutNegator EWithoutNegator
Element type, stripped of negator<> if it has one.
Definition: Vec.h:195
Definition: Vec.h:384
CNT< E >::TComplex EComplex
Type that elements would have if complex, if E is currently real; otherwise just the element type E...
Definition: Vec.h:205
Vec(const ENeg &ne)
Construction from a single value of this Vec&#39;s negated element type assigns that value to each elemen...
Definition: Vec.h:474
K::TPosTrans TPosTrans
Definition: CompositeNumericalTypes.h:145
Vec & operator*=(const EE &e)
Definition: Vec.h:776
void elementwiseMultiply(const Row< 1, E1, S1 > &r1, const Row< 1, E2, S2 > &r2, Row< 1, typename CNT< E1 >::template Result< E2 >::Mul > &result)
Definition: Row.h:75
float norm(const conjugate< float > &c)
Definition: conjugate.h:775
Vec< M, EImag, STRIDE *CNT< E >::RealStrideFactor > TImag
Type of this Vec cast to show only the imaginary part of its element; this might affect the stride...
Definition: Vec.h:277
TReal & real()
Recast to show only the real portion of this Vec and return a writable reference. ...
Definition: Vec.h:682
const TPosTrans & positionalTranspose() const
Positional transpose turns this Vec into a Row but does not transpose the individual elements...
Definition: Vec.h:669
CNT< E >::TAbs EAbs
Type required to hold the result of abs(E).
Definition: Vec.h:218
K::StdNumber StdNumber
Definition: CompositeNumericalTypes.h:163
RowVectorBase< typename CNT< ELEM >::TAbs > abs(const RowVectorBase< ELEM > &v)
Definition: VectorMath.h:120
Vec & scalarTimesEqFromLeft(const EE &ee)
Definition: Vec.h:791
bool operator!=(const conjugate< R > &a, const float &b)
Definition: conjugate.h:859
E TRow
Type of a row of this CNT object (for a Vec, just its element type).
Definition: Vec.h:289
Specialized information about Composite Numerical Types which allows us to define appropriate templat...
Definition: CompositeNumericalTypes.h:136
Vec & scalarPlusEq(int ee)
Definition: Vec.h:800
THerm & updTranspose()
Non-operator version of Hermitian transpose; recasts and returns a writable reference.
Definition: Vec.h:663
DvdCNTs< M, 1, ArgDepth, Vec, ColSpacing, RowSpacing, CNT< P >::NRows, CNT< P >::NCols, CNT< P >::ArgDepth, P, CNT< P >::ColSpacing, CNT< P >::RowSpacing > DvdOp
Definition: Vec.h:397
Vec & scalarDivideEq(int ee)
Definition: Vec.h:803
CNT< E >::Precision EPrecision
These compile-time constants are required of every Composite Numerical Type (CNT).
Definition: Vec.h:232
This is a fixed-length row vector designed for no-overhead inline computation.
Definition: SimTKcommon/include/SimTKcommon/internal/common.h:606
Vec & scalarMinusEqFromLeft(const EE &ee)
Definition: Vec.h:787
const Real E
e = Real(exp(1))
Mandatory first inclusion for any Simbody source or header file.
Vec< M, ENormalize, 1 > TNormalize
These compile-time constants are required of every Composite Numerical Type (CNT).
Definition: Vec.h:299
E & operator()(int i)
Same as non-const operator[] above.
Definition: Vec.h:604
Vec< M, typename CNT< E >::template Result< EE >::Sub > scalarSubtract(const EE &e) const
Definition: Vec.h:758
K::TNeg TNeg
Definition: CompositeNumericalTypes.h:139
K::TStandard TStandard
Definition: CompositeNumericalTypes.h:156
void copy(Row< 1, E1, S1 > &r1, const Row< 1, E2, S2 > &r2)
Definition: Row.h:105
ScalarNormSq normSqr() const
Definition: Vec.h:606
K::TWithoutNegator TWithoutNegator
Definition: CompositeNumericalTypes.h:140
Vec & operator=(const Vec< M, EE, SS > &vv)
Assignment to a conforming Vec, of any element type and stride, provided that the element types are a...
Definition: Vec.h:526
void conformingSubtract(const Row< 1, E1, S1 > &r1, const Row< 1, E2, S2 > &r2, Row< 1, typename CNT< E1 >::template Result< E2 >::Sub > &result)
Definition: Row.h:60
const THerm & transpose() const
Non-operator version of Hermitian transpose; just a recast.
Definition: Vec.h:660
CNT< ScalarNormSq >::TSqrt norm() const
Definition: Vec.h:608
Vec< MM, ELT, STRIDE > & updSubVec(int i)
Extract a writable reference to a sub-Vec with size known at compile time.
Definition: Vec.h:835
Vec< M, typename CNT< E >::template Result< P >::Sub, 1 > Sub
Definition: Vec.h:379
CNT< E >::TInvert EInvert
Packed type that can hold the value returned from invert(E), the inverse type of an element...
Definition: Vec.h:224
Matrix_< typename CNT< E1 >::template Result< E2 >::Add > operator+(const MatrixBase< E1 > &l, const MatrixBase< E2 > &r)
Definition: BigMatrix.h:568
Vec< M, EComplex, STRIDE > TComplex
These compile-time constants are required of every Composite Numerical Type (CNT).
Definition: Vec.h:278
Vec< M+1, ELT, 1 > append1(const EE &v) const
Return a vector one larger than this one by adding an element to the end.
Definition: Vec.h:875
const Vec< MM, ELT, STRIDE > & getSubVec(int i) const
Extract a const reference to a sub-Vec with size known at compile time.
Definition: Vec.h:825
This class represents a small matrix whose size is known at compile time, containing elements of any ...
Definition: SimTKcommon/include/SimTKcommon/internal/common.h:607
Vec & scalarDivideEqFromLeft(int ee)
Definition: Vec.h:806
K::TComplex TComplex
Definition: CompositeNumericalTypes.h:143
Vec & scalarTimesEq(int ee)
Definition: Vec.h:802
K::Number Number
Definition: CompositeNumericalTypes.h:162
Vec & operator+=(const Vec< M, EE, SS > &r)
Add in a conforming Vec, of any element type and stride, provided that the element types are addition...
Definition: Vec.h:531
Row< M, E, STRIDE > TPosTrans
Type of this Vec after casting to its positional transpose; that is, the Vec turns into a Row but the...
Definition: Vec.h:285
bool isInf() const
Return true if any element of this Vec contains a +Infinity or -Infinity somewhere but no element con...
Definition: Vec.h:925
Vec & scalarDivideEq(const EE &ee)
Definition: Vec.h:793
static K getNaN()
Definition: CompositeNumericalTypes.h:246
Vec & operator-=(const Vec< M, negator< EE >, SS > &r)
Subtract off a conforming Vec, of any negated element type and stride, provided that the element type...
Definition: Vec.h:546
static Vec & updSubVec(Vec< MM, ELT, STRIDE > &v, int i)
Extract a subvector of type Vec from a longer one that has the same element type and stride...
Definition: Vec.h:853
Vec(const E &e0, const E &e1, const E &e2, const E &e3, const E &e4)
Definition: Vec.h:494
Vec< M, ENeg, STRIDE > TNeg
Type this Vec would have if its elements were interpreted as negated.
Definition: Vec.h:266
CNT< E >::Number ENumber
These compile-time constants are required of every Composite Numerical Type (CNT).
Definition: Vec.h:230
ELT E
Element type of this Vec.
Definition: Vec.h:191
EScalarNormSq ScalarNormSq
These compile-time constants are required of every Composite Numerical Type (CNT).
Definition: Vec.h:311
K::TSqHermT TSqHermT
Definition: CompositeNumericalTypes.h:146
bool isNumericallyEqual(const Vec< M, E2, RS2 > &v) const
Test whether this vector is numerically equal to some other vector with the same shape, using a default tolerance which is the looser of the default tolerances of the two objects being compared.
Definition: Vec.h:965
bool operator>(const Row< N, E1, S1 > &l, const Row< N, E2, S2 > &r)
bool = v1[i] > v2[i], for all elements i
Definition: Row.h:833
Vec(const E &e0, const E &e1, const E &e2, const E &e3, const E &e4, const E &e5)
Definition: Vec.h:497
CNT< E >::ULessScalar EULessScalar
These compile-time constants are required of every Composite Numerical Type (CNT).
Definition: Vec.h:229
AddOp::Type Add
Definition: Vec.h:403
Vec TCol
Type of a column of this CNT object (for a Vec, the whole thing).
Definition: Vec.h:291
EULessScalar ULessScalar
These compile-time constants are required of every Composite Numerical Type (CNT).
Definition: Vec.h:307
K::THerm THerm
Definition: CompositeNumericalTypes.h:144
Vec()
Default construction initializes Vec&#39;s elements to NaN when debugging but leaves them uninitialized g...
Definition: Vec.h:423
Vec< M, typename CNT< E >::template Result< EE >::Sub > conformingSubtract(const Vec< M, EE, SS > &r) const
Vector subtraction – use operator- instead.
Definition: Vec.h:561
Vec< M, typename CNT< E >::template Result< EE >::Dvd > scalarDivide(const EE &e) const
Definition: Vec.h:737
Vec & scalarMinusEq(const EE &ee)
Definition: Vec.h:785
static int ncol()
The number of columns in a Vec is always 1.
Definition: Vec.h:320
EScalar Scalar
These compile-time constants are required of every Composite Numerical Type (CNT).
Definition: Vec.h:306
Definition: negator.h:64
Vec(const Vec< M, E, SS > &src)
This is an implicit conversion from a Vec of the same length and element type but with a different st...
Definition: Vec.h:450
Vec & scalarDivideEqFromLeft(const EE &ee)
Definition: Vec.h:795
const TWithoutNegator & castAwayNegatorIfAny() const
Recast to remove negators from this Vec&#39;s type if present; this is handy for simplifying operations w...
Definition: Vec.h:706
Vec(int i)
Given an int value, turn it into a suitable floating point number, convert that to element type E and...
Definition: Vec.h:483
CNT< E >::TNeg ENeg
Negated version of this Vec&#39;s element type; ENeg==negator< E >.
Definition: Vec.h:193
void conformingAdd(const Row< 1, E1, S1 > &r1, const Row< 1, E2, S2 > &r2, Row< 1, typename CNT< E1 >::template Result< E2 >::Add > &result)
Definition: Row.h:45
Vec< M, typename CNT< EE >::template Result< E >::Dvd > scalarDivideFromLeft(const EE &e) const
Definition: Vec.h:743
Vec & operator-=(const Vec< M, EE, SS > &r)
Subtract off a conforming Vec, of any element type and stride, provided that the element types are ad...
Definition: Vec.h:541
DvdOp::Type Dvd
Definition: Vec.h:398
Vec< M, EWithoutNegator, STRIDE > TWithoutNegator
Type of this Vec with negator removed from its element type, if the element is negated.
Definition: Vec.h:269
Vec< M, typename CNT< E >::template Result< P >::Dvd, 1 > Dvd
Definition: Vec.h:377
static const Vec & getAs(const ELT *p)
Recast an ordinary C++ array E[] to a const Vec<M,E,S>; assumes compatible length, stride, and packing.
Definition: Vec.h:902
K::TAbs TAbs
Definition: CompositeNumericalTypes.h:155
bool isNumericallyEqual(const float &a, const float &b, double tol=RTraits< float >::getDefaultTolerance())
Compare two floats for approximate equality.
Definition: NTraits.h:313
MulCNTs< M, 1, ArgDepth, Vec, ColSpacing, RowSpacing, CNT< P >::NRows, CNT< P >::NCols, CNT< P >::ArgDepth, P, CNT< P >::ColSpacing, CNT< P >::RowSpacing > MulOp
Definition: Vec.h:387