Simbody  3.5
Row.h
Go to the documentation of this file.
1 #ifndef SimTK_SIMMATRIX_SMALLMATRIX_ROW_H_
2 #define SimTK_SIMMATRIX_SMALLMATRIX_ROW_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 
34 namespace SimTK {
35 
36 // The following functions are used internally by Row.
37 
38 namespace Impl {
39 
40 // For those wimpy compilers that don't unroll short, constant-limit loops,
41 // Peter Eastman added these recursive template implementations of
42 // elementwise add, subtract, and copy. Sherm added multiply and divide.
43 
44 template <class E1, int S1, class E2, int S2> void
45 conformingAdd(const Row<1,E1,S1>& r1, const Row<1,E2,S2>& r2,
46  Row<1,typename CNT<E1>::template Result<E2>::Add>& result) {
47  result[0] = r1[0] + r2[0];
48 }
49 template <int N, class E1, int S1, class E2, int S2> void
50 conformingAdd(const Row<N,E1,S1>& r1, const Row<N,E2,S2>& r2,
51  Row<N,typename CNT<E1>::template Result<E2>::Add>& result) {
52  conformingAdd(reinterpret_cast<const Row<N-1,E1,S1>&>(r1),
53  reinterpret_cast<const Row<N-1,E2,S2>&>(r2),
54  reinterpret_cast<Row<N-1,typename CNT<E1>::
55  template Result<E2>::Add>&>(result));
56  result[N-1] = r1[N-1] + r2[N-1];
57 }
58 
59 template <class E1, int S1, class E2, int S2> void
61  Row<1,typename CNT<E1>::template Result<E2>::Sub>& result) {
62  result[0] = r1[0] - r2[0];
63 }
64 template <int N, class E1, int S1, class E2, int S2> void
66  Row<N,typename CNT<E1>::template Result<E2>::Sub>& result) {
67  conformingSubtract(reinterpret_cast<const Row<N-1,E1,S1>&>(r1),
68  reinterpret_cast<const Row<N-1,E2,S2>&>(r2),
69  reinterpret_cast<Row<N-1,typename CNT<E1>::
70  template Result<E2>::Sub>&>(result));
71  result[N-1] = r1[N-1] - r2[N-1];
72 }
73 
74 template <class E1, int S1, class E2, int S2> void
76  Row<1,typename CNT<E1>::template Result<E2>::Mul>& result) {
77  result[0] = r1[0] * r2[0];
78 }
79 template <int N, class E1, int S1, class E2, int S2> void
81  Row<N,typename CNT<E1>::template Result<E2>::Mul>& result) {
82  elementwiseMultiply(reinterpret_cast<const Row<N-1,E1,S1>&>(r1),
83  reinterpret_cast<const Row<N-1,E2,S2>&>(r2),
84  reinterpret_cast<Row<N-1,typename CNT<E1>::
85  template Result<E2>::Mul>&>(result));
86  result[N-1] = r1[N-1] * r2[N-1];
87 }
88 
89 template <class E1, int S1, class E2, int S2> void
91  Row<1,typename CNT<E1>::template Result<E2>::Dvd>& result) {
92  result[0] = r1[0] / r2[0];
93 }
94 template <int N, class E1, int S1, class E2, int S2> void
96  Row<N,typename CNT<E1>::template Result<E2>::Dvd>& result) {
97  elementwiseDivide(reinterpret_cast<const Row<N-1,E1,S1>&>(r1),
98  reinterpret_cast<const Row<N-1,E2,S2>&>(r2),
99  reinterpret_cast<Row<N-1,typename CNT<E1>::
100  template Result<E2>::Dvd>&>(result));
101  result[N-1] = r1[N-1] / r2[N-1];
102 }
103 
104 template <class E1, int S1, class E2, int S2> void
105 copy(Row<1,E1,S1>& r1, const Row<1,E2,S2>& r2) {
106  r1[0] = r2[0];
107 }
108 template <int N, class E1, int S1, class E2, int S2> void
109 copy(Row<N,E1,S1>& r1, const Row<N,E2,S2>& r2) {
110  copy(reinterpret_cast<Row<N-1,E1,S1>&>(r1),
111  reinterpret_cast<const Row<N-1,E2,S2>&>(r2));
112  r1[N-1] = r2[N-1];
113 }
114 
115 }
116 
132 template <int N, class ELT, int STRIDE> class Row {
133  typedef ELT E;
134  typedef typename CNT<E>::TNeg ENeg;
135  typedef typename CNT<E>::TWithoutNegator EWithoutNegator;
136  typedef typename CNT<E>::TReal EReal;
137  typedef typename CNT<E>::TImag EImag;
138  typedef typename CNT<E>::TComplex EComplex;
139  typedef typename CNT<E>::THerm EHerm;
140  typedef typename CNT<E>::TPosTrans EPosTrans;
141  typedef typename CNT<E>::TSqHermT ESqHermT;
142  typedef typename CNT<E>::TSqTHerm ESqTHerm;
143 
144  typedef typename CNT<E>::TSqrt ESqrt;
145  typedef typename CNT<E>::TAbs EAbs;
146  typedef typename CNT<E>::TStandard EStandard;
147  typedef typename CNT<E>::TInvert EInvert;
148  typedef typename CNT<E>::TNormalize ENormalize;
149 
150  typedef typename CNT<E>::Scalar EScalar;
151  typedef typename CNT<E>::ULessScalar EULessScalar;
152  typedef typename CNT<E>::Number ENumber;
153  typedef typename CNT<E>::StdNumber EStdNumber;
154  typedef typename CNT<E>::Precision EPrecision;
155  typedef typename CNT<E>::ScalarNormSq EScalarNormSq;
156 
157 public:
158 
159  enum {
160  NRows = 1,
161  NCols = N,
162  NPackedElements = N,
163  NActualElements = N * STRIDE, // includes trailing gap
164  NActualScalars = CNT<E>::NActualScalars * NActualElements,
165  RowSpacing = NActualElements,
166  ColSpacing = STRIDE,
168  RealStrideFactor = 1, // composite types don't change size when
169  // cast from complex to real or imaginary
170  ArgDepth = ((int)CNT<E>::ArgDepth < (int)MAX_RESOLVED_DEPTH
171  ? CNT<E>::ArgDepth + 1
173  IsScalar = 0,
174  IsULessScalar = 0,
175  IsNumber = 0,
176  IsStdNumber = 0,
177  IsPrecision = 0,
178  SignInterpretation = CNT<E>::SignInterpretation
179  };
180 
184 
192  typedef E TElement;
193  typedef Row TRow;
194  typedef E TCol;
195 
196  // These are the results of calculations, so are returned in new, packed
197  // memory. Be sure to refer to element types here which are also packed.
198  typedef Vec<N,ESqrt,1> TSqrt; // Note stride
199  typedef Row<N,EAbs,1> TAbs; // Note stride
201  typedef Vec<N,EInvert,1> TInvert; // packed
203 
204  typedef SymMat<N,ESqHermT> TSqHermT; // result of self outer product
205  typedef EScalarNormSq TSqTHerm; // result of self dot product
206 
207  // These recurse right down to the underlying scalar type no matter how
208  // deep the elements are.
209  typedef EScalar Scalar;
210  typedef EULessScalar ULessScalar;
211  typedef ENumber Number;
212  typedef EStdNumber StdNumber;
213  typedef EPrecision Precision;
214  typedef EScalarNormSq ScalarNormSq;
215 
216  static int size() { return N; }
217  static int nrow() { return 1; }
218  static int ncol() { return N; }
219 
220 
221  // Scalar norm square is sum( conjugate squares of all scalars )
222  ScalarNormSq scalarNormSqr() const {
223  ScalarNormSq sum(0);
224  for(int i=0;i<N;++i) sum += CNT<E>::scalarNormSqr(d[i*STRIDE]);
225  return sum;
226  }
227 
228  // sqrt() is elementwise square root; that is, the return value has the same
229  // dimension as this Vec but with each element replaced by whatever it thinks
230  // its square root is.
231  TSqrt sqrt() const {
232  TSqrt rsqrt;
233  for(int i=0;i<N;++i) rsqrt[i] = CNT<E>::sqrt(d[i*STRIDE]);
234  return rsqrt;
235  }
236 
237  // abs() is elementwise absolute value; that is, the return value has the same
238  // dimension as this Row but with each element replaced by whatever it thinks
239  // its absolute value is.
240  TAbs abs() const {
241  TAbs rabs;
242  for(int i=0;i<N;++i) rabs[i] = CNT<E>::abs(d[i*STRIDE]);
243  return rabs;
244  }
245 
246  TStandard standardize() const {
247  TStandard rstd;
248  for(int i=0;i<N;++i) rstd[i] = CNT<E>::standardize(d[i*STRIDE]);
249  return rstd;
250  }
251 
252  // Sum just adds up all the elements, getting rid of negators and
253  // conjugates in the process.
254  EStandard sum() const {
255  E sum(0);
256  for (int i=0;i<N;++i) sum += d[i*STRIDE];
257  return CNT<E>::standardize(sum);
258  }
259 
260  // This gives the resulting rowvector type when (v[i] op P) is applied to each element of v.
261  // It is a row of length N, stride 1, and element types which are the regular
262  // composite result of E op P. Typically P is a scalar type but it doesn't have to be.
263  template <class P> struct EltResult {
268  };
269 
270  // This is the composite result for v op P where P is some kind of appropriately shaped
271  // non-scalar type.
272  template <class P> struct Result {
273  typedef MulCNTs<1,N,ArgDepth,Row,ColSpacing,RowSpacing,
276  typedef typename MulOp::Type Mul;
277 
278  typedef MulCNTsNonConforming<1,N,ArgDepth,Row,ColSpacing,RowSpacing,
279  CNT<P>::NRows, CNT<P>::NCols, CNT<P>::ArgDepth,
280  P, CNT<P>::ColSpacing, CNT<P>::RowSpacing> MulOpNonConforming;
281  typedef typename MulOpNonConforming::Type MulNon;
282 
283 
284  typedef DvdCNTs<1,N,ArgDepth,Row,ColSpacing,RowSpacing,
285  CNT<P>::NRows, CNT<P>::NCols, CNT<P>::ArgDepth,
286  P, CNT<P>::ColSpacing, CNT<P>::RowSpacing> DvdOp;
287  typedef typename DvdOp::Type Dvd;
288 
289  typedef AddCNTs<1,N,ArgDepth,Row,ColSpacing,RowSpacing,
290  CNT<P>::NRows, CNT<P>::NCols, CNT<P>::ArgDepth,
291  P, CNT<P>::ColSpacing, CNT<P>::RowSpacing> AddOp;
292  typedef typename AddOp::Type Add;
293 
294  typedef SubCNTs<1,N,ArgDepth,Row,ColSpacing,RowSpacing,
295  CNT<P>::NRows, CNT<P>::NCols, CNT<P>::ArgDepth,
296  P, CNT<P>::ColSpacing, CNT<P>::RowSpacing> SubOp;
297  typedef typename SubOp::Type Sub;
298  };
299 
300  // Shape-preserving element substitution (always packed)
301  template <class P> struct Substitute {
302  typedef Row<N,P> Type;
303  };
304 
305  // Default construction initializes to NaN when debugging but
306  // is left uninitialized otherwise.
307  Row(){
308  #ifndef NDEBUG
309  setToNaN();
310  #endif
311  }
312 
313  // It's important not to use the default copy constructor or copy
314  // assignment because the compiler doesn't understand that we may
315  // have noncontiguous storage and will try to copy the whole array.
316  Row(const Row& src) {
317  Impl::copy(*this, src);
318  }
319  Row& operator=(const Row& src) { // no harm if src and 'this' are the same
320  Impl::copy(*this, src);
321  return *this;
322  }
323 
324  // We want an implicit conversion from a Row of the same length
325  // and element type but with a different stride.
326  template <int SS> Row(const Row<N,E,SS>& src) {
327  Impl::copy(*this, src);
328  }
329 
330  // We want an implicit conversion from a Row of the same length
331  // and *negated* element type, possibly with a different stride.
332  template <int SS> Row(const Row<N,ENeg,SS>& src) {
333  Impl::copy(*this, src);
334  }
335 
336  // Construct a Row from a Row of the same length, with any
337  // stride. Works as long as the element types are compatible.
338  template <class EE, int SS> explicit Row(const Row<N,EE,SS>& vv) {
339  Impl::copy(*this, vv);
340  }
341 
342  // Construction using an element assigns to each element.
343  explicit Row(const E& e)
344  { for (int i=0;i<N;++i) d[i*STRIDE]=e; }
345 
346  // Construction using a negated element assigns to each element.
347  explicit Row(const ENeg& ne)
348  { for (int i=0;i<N;++i) d[i*STRIDE]=ne; }
349 
350  // Given an int, turn it into a suitable floating point number
351  // and then feed that to the above single-element constructor.
352  explicit Row(int i)
353  { new (this) Row(E(Precision(i))); }
354 
355  // A bevy of constructors for Rows up to length 6.
356  Row(const E& e0,const E& e1)
357  { assert(N==2);(*this)[0]=e0;(*this)[1]=e1; }
358  Row(const E& e0,const E& e1,const E& e2)
359  { assert(N==3);(*this)[0]=e0;(*this)[1]=e1;(*this)[2]=e2; }
360  Row(const E& e0,const E& e1,const E& e2,const E& e3)
361  { assert(N==4);(*this)[0]=e0;(*this)[1]=e1;(*this)[2]=e2;(*this)[3]=e3; }
362  Row(const E& e0,const E& e1,const E& e2,const E& e3,const E& e4)
363  { assert(N==5);(*this)[0]=e0;(*this)[1]=e1;(*this)[2]=e2;
364  (*this)[3]=e3;(*this)[4]=e4; }
365  Row(const E& e0,const E& e1,const E& e2,const E& e3,const E& e4,const E& e5)
366  { assert(N==6);(*this)[0]=e0;(*this)[1]=e1;(*this)[2]=e2;
367  (*this)[3]=e3;(*this)[4]=e4;(*this)[5]=e5; }
368  Row(const E& e0,const E& e1,const E& e2,const E& e3,const E& e4,const E& e5,const E& e6)
369  { assert(N==7);(*this)[0]=e0;(*this)[1]=e1;(*this)[2]=e2;
370  (*this)[3]=e3;(*this)[4]=e4;(*this)[5]=e5;(*this)[6]=e6; }
371  Row(const E& e0,const E& e1,const E& e2,const E& e3,const E& e4,const E& e5,const E& e6,const E& e7)
372  { assert(N==8);(*this)[0]=e0;(*this)[1]=e1;(*this)[2]=e2;
373  (*this)[3]=e3;(*this)[4]=e4;(*this)[5]=e5;(*this)[6]=e6;(*this)[7]=e7; }
374  Row(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)
375  { assert(N==9);(*this)[0]=e0;(*this)[1]=e1;(*this)[2]=e2;
376  (*this)[3]=e3;(*this)[4]=e4;(*this)[5]=e5;(*this)[6]=e6;(*this)[7]=e7;(*this)[8]=e8; }
377 
378  // Construction from a pointer to anything assumes we're pointing
379  // at an element list of the right length.
380  template <class EE> explicit Row(const EE* p)
381  { assert(p); for(int i=0;i<N;++i) d[i*STRIDE]=p[i]; }
382  template <class EE> Row& operator=(const EE* p)
383  { assert(p); for(int i=0;i<N;++i) d[i*STRIDE]=p[i]; return *this; }
384 
385  // Conforming assignment ops.
386  template <class EE, int SS> Row& operator=(const Row<N,EE,SS>& vv) {
387  Impl::copy(*this, vv);
388  return *this;
389  }
390  template <class EE, int SS> Row& operator+=(const Row<N,EE,SS>& r)
391  { for(int i=0;i<N;++i) d[i*STRIDE] += r[i]; return *this; }
392  template <class EE, int SS> Row& operator+=(const Row<N,negator<EE>,SS>& r)
393  { for(int i=0;i<N;++i) d[i*STRIDE] -= -(r[i]); return *this; }
394  template <class EE, int SS> Row& operator-=(const Row<N,EE,SS>& r)
395  { for(int i=0;i<N;++i) d[i*STRIDE] -= r[i]; return *this; }
396  template <class EE, int SS> Row& operator-=(const Row<N,negator<EE>,SS>& r)
397  { for(int i=0;i<N;++i) d[i*STRIDE] += -(r[i]); return *this; }
398 
399  // Conforming binary ops with 'this' on left, producing new packed result.
400  // Cases: r=r+r, r=r-r, s=r*v r=r*m
401 
403  template <class EE, int SS> Row<N,typename CNT<E>::template Result<EE>::Add>
404  conformingAdd(const Row<N,EE,SS>& r) const {
405  Row<N,typename CNT<E>::template Result<EE>::Add> result;
406  Impl::conformingAdd(*this, r, result);
407  return result;
408  }
409 
411  template <class EE, int SS> Row<N,typename CNT<E>::template Result<EE>::Sub>
413  Row<N,typename CNT<E>::template Result<EE>::Sub> result;
414  Impl::conformingSubtract(*this, r, result);
415  return result;
416  }
417 
419  template <class EE, int SS> typename CNT<E>::template Result<EE>::Mul
421  return (*this)*r;
422  }
423 
425  template <int MatNCol, class EE, int CS, int RS>
428  Row<MatNCol,typename CNT<E>::template Result<EE>::Mul> result;
429  for (int j=0;j<N;++j) result[j] = conformingMultiply(m(j));
430  return result;
431  }
432 
434  template <class EE, int SS> Row<N,typename CNT<E>::template Result<EE>::Mul>
436  Row<N,typename CNT<E>::template Result<EE>::Mul> result;
437  Impl::elementwiseMultiply(*this, r, result);
438  return result;
439  }
440 
442  template <class EE, int SS> Row<N,typename CNT<E>::template Result<EE>::Dvd>
443  elementwiseDivide(const Row<N,EE,SS>& r) const {
444  Row<N,typename CNT<E>::template Result<EE>::Dvd> result;
445  Impl::elementwiseDivide(*this, r, result);
446  return result;
447  }
448 
449  const E& operator[](int i) const { assert(0 <= i && i < N); return d[i*STRIDE]; }
450  E& operator[](int i) { assert(0 <= i && i < N); return d[i*STRIDE]; }
451  const E& operator()(int i) const { return (*this)[i]; }
452  E& operator()(int i) { return (*this)[i]; }
453 
454  ScalarNormSq normSqr() const { return scalarNormSqr(); }
455  typename CNT<ScalarNormSq>::TSqrt
456  norm() const { return CNT<ScalarNormSq>::sqrt(scalarNormSqr()); }
457 
458  // If the elements of this Row are scalars, the result is what you get by
459  // dividing each element by the norm() calculated above. If the elements are
460  // *not* scalars, then the elements are *separately* normalized. That means
461  // you will get a different answer from Row<2,Row3>::normalize() than you
462  // would from a Row<6>::normalize() containing the same scalars.
463  //
464  // Normalize returns a row of the same dimension but in new, packed storage
465  // and with a return type that does not include negator<> even if the original
466  // Row<> does, because we can eliminate the negation here almost for free.
467  // But we can't standardize (change conjugate to complex) for free, so we'll retain
468  // conjugates if there are any.
469  TNormalize normalize() const {
470  if (CNT<E>::IsScalar) {
471  return castAwayNegatorIfAny() / (SignInterpretation*norm());
472  } else {
473  TNormalize elementwiseNormalized;
474  for (int j=0; j<N; ++j)
475  elementwiseNormalized[j] = CNT<E>::normalize((*this)[j]);
476  return elementwiseNormalized;
477  }
478  }
479 
480  TInvert invert() const {assert(false); return TInvert();} // TODO default inversion
481 
482  const Row& operator+() const { return *this; }
483  const TNeg& operator-() const { return negate(); }
484  TNeg& operator-() { return updNegate(); }
485  const THerm& operator~() const { return transpose(); }
486  THerm& operator~() { return updTranspose(); }
487 
488  const TNeg& negate() const { return *reinterpret_cast<const TNeg*>(this); }
489  TNeg& updNegate() { return *reinterpret_cast<TNeg*>(this); }
490 
491  const THerm& transpose() const { return *reinterpret_cast<const THerm*>(this); }
492  THerm& updTranspose() { return *reinterpret_cast<THerm*>(this); }
493 
494  const TPosTrans& positionalTranspose() const
495  { return *reinterpret_cast<const TPosTrans*>(this); }
497  { return *reinterpret_cast<TPosTrans*>(this); }
498 
499  const TReal& real() const { return *reinterpret_cast<const TReal*>(this); }
500  TReal& real() { return *reinterpret_cast< TReal*>(this); }
501 
502  // Had to contort these routines to get them through VC++ 7.net
503  const TImag& imag() const {
504  const int offs = ImagOffset;
505  const EImag* p = reinterpret_cast<const EImag*>(this);
506  return *reinterpret_cast<const TImag*>(p+offs);
507  }
508  TImag& imag() {
509  const int offs = ImagOffset;
510  EImag* p = reinterpret_cast<EImag*>(this);
511  return *reinterpret_cast<TImag*>(p+offs);
512  }
513 
514  const TWithoutNegator& castAwayNegatorIfAny() const {return *reinterpret_cast<const TWithoutNegator*>(this);}
515  TWithoutNegator& updCastAwayNegatorIfAny() {return *reinterpret_cast<TWithoutNegator*>(this);}
516 
517 
518  // These are elementwise binary operators, (this op ee) by default but
519  // (ee op this) if 'FromLeft' appears in the name. The result is a packed
520  // Row<N> but the element type may change. These are mostly used to
521  // implement global operators. We call these "scalar" operators but
522  // actually the "scalar" can be a composite type.
523 
524  //TODO: consider converting 'e' to Standard Numbers as precalculation and
525  // changing return type appropriately.
526  template <class EE> Row<N, typename CNT<E>::template Result<EE>::Mul>
527  scalarMultiply(const EE& e) const {
528  Row<N, typename CNT<E>::template Result<EE>::Mul> result;
529  for (int j=0; j<N; ++j) result[j] = (*this)[j] * e;
530  return result;
531  }
532  template <class EE> Row<N, typename CNT<EE>::template Result<E>::Mul>
533  scalarMultiplyFromLeft(const EE& e) const {
534  Row<N, typename CNT<EE>::template Result<E>::Mul> result;
535  for (int j=0; j<N; ++j) result[j] = e * (*this)[j];
536  return result;
537  }
538 
539  // TODO: should precalculate and store 1/e, while converting to Standard
540  // Numbers. Note that return type should change appropriately.
541  template <class EE> Row<N, typename CNT<E>::template Result<EE>::Dvd>
542  scalarDivide(const EE& e) const {
543  Row<N, typename CNT<E>::template Result<EE>::Dvd> result;
544  for (int j=0; j<N; ++j) result[j] = (*this)[j] / e;
545  return result;
546  }
547  template <class EE> Row<N, typename CNT<EE>::template Result<E>::Dvd>
548  scalarDivideFromLeft(const EE& e) const {
549  Row<N, typename CNT<EE>::template Result<E>::Dvd> result;
550  for (int j=0; j<N; ++j) result[j] = e / (*this)[j];
551  return result;
552  }
553 
554  template <class EE> Row<N, typename CNT<E>::template Result<EE>::Add>
555  scalarAdd(const EE& e) const {
556  Row<N, typename CNT<E>::template Result<EE>::Add> result;
557  for (int j=0; j<N; ++j) result[j] = (*this)[j] + e;
558  return result;
559  }
560  // Add is commutative, so no 'FromLeft'.
561 
562  template <class EE> Row<N, typename CNT<E>::template Result<EE>::Sub>
563  scalarSubtract(const EE& e) const {
564  Row<N, typename CNT<E>::template Result<EE>::Sub> result;
565  for (int j=0; j<N; ++j) result[j] = (*this)[j] - e;
566  return result;
567  }
568  template <class EE> Row<N, typename CNT<EE>::template Result<E>::Sub>
569  scalarSubtractFromLeft(const EE& e) const {
570  Row<N, typename CNT<EE>::template Result<E>::Sub> result;
571  for (int j=0; j<N; ++j) result[j] = e - (*this)[j];
572  return result;
573  }
574 
575  // Generic assignments for any element type not listed explicitly, including scalars.
576  // These are done repeatedly for each element and only work if the operation can
577  // be performed leaving the original element type.
578  template <class EE> Row& operator =(const EE& e) {return scalarEq(e);}
579  template <class EE> Row& operator+=(const EE& e) {return scalarPlusEq(e);}
580  template <class EE> Row& operator-=(const EE& e) {return scalarMinusEq(e);}
581  template <class EE> Row& operator*=(const EE& e) {return scalarTimesEq(e);}
582  template <class EE> Row& operator/=(const EE& e) {return scalarDivideEq(e);}
583 
584  // Generalized scalar assignment & computed assignment methods. These will work
585  // for any assignment-compatible element, not just scalars.
586  template <class EE> Row& scalarEq(const EE& ee)
587  { for(int i=0;i<N;++i) d[i*STRIDE] = ee; return *this; }
588  template <class EE> Row& scalarPlusEq(const EE& ee)
589  { for(int i=0;i<N;++i) d[i*STRIDE] += ee; return *this; }
590  template <class EE> Row& scalarMinusEq(const EE& ee)
591  { for(int i=0;i<N;++i) d[i*STRIDE] -= ee; return *this; }
592  template <class EE> Row& scalarMinusEqFromLeft(const EE& ee)
593  { for(int i=0;i<N;++i) d[i*STRIDE] = ee - d[i*STRIDE]; return *this; }
594  template <class EE> Row& scalarTimesEq(const EE& ee)
595  { for(int i=0;i<N;++i) d[i*STRIDE] *= ee; return *this; }
596  template <class EE> Row& scalarTimesEqFromLeft(const EE& ee)
597  { for(int i=0;i<N;++i) d[i*STRIDE] = ee * d[i*STRIDE]; return *this; }
598  template <class EE> Row& scalarDivideEq(const EE& ee)
599  { for(int i=0;i<N;++i) d[i*STRIDE] /= ee; return *this; }
600  template <class EE> Row& scalarDivideEqFromLeft(const EE& ee)
601  { for(int i=0;i<N;++i) d[i*STRIDE] = ee / d[i*STRIDE]; return *this; }
602 
603 
604  // Specialize for int to avoid warnings and ambiguities.
605  Row& scalarEq(int ee) {return scalarEq(Precision(ee));}
606  Row& scalarPlusEq(int ee) {return scalarPlusEq(Precision(ee));}
607  Row& scalarMinusEq(int ee) {return scalarMinusEq(Precision(ee));}
608  Row& scalarTimesEq(int ee) {return scalarTimesEq(Precision(ee));}
609  Row& scalarDivideEq(int ee) {return scalarDivideEq(Precision(ee));}
610  Row& scalarMinusEqFromLeft(int ee) {return scalarMinusEqFromLeft(Precision(ee));}
611  Row& scalarTimesEqFromLeft(int ee) {return scalarTimesEqFromLeft(Precision(ee));}
612  Row& scalarDivideEqFromLeft(int ee) {return scalarDivideEqFromLeft(Precision(ee));}
613 
616  void setToNaN() {
617  (*this) = CNT<ELT>::getNaN();
618  }
619 
621  void setToZero() {
622  (*this) = ELT(0);
623  }
624 
630  template <int NN>
631  const Row<NN,ELT,STRIDE>& getSubRow(int j) const {
632  assert(0 <= j && j + NN <= N);
633  return Row<NN,ELT,STRIDE>::getAs(&(*this)[j]);
634  }
640  template <int NN>
642  assert(0 <= j && j + NN <= N);
643  return Row<NN,ELT,STRIDE>::updAs(&(*this)[j]);
644  }
645 
649  template <int NN>
650  static const Row& getSubRow(const Row<NN,ELT,STRIDE>& r, int j) {
651  assert(0 <= j && j + N <= NN);
652  return getAs(&r[j]);
653  }
657  template <int NN>
658  static Row& updSubRow(Row<NN,ELT,STRIDE>& r, int j) {
659  assert(0 <= j && j + N <= NN);
660  return updAs(&r[j]);
661  }
662 
666  Row<N-1,ELT,1> drop1(int p) const {
667  assert(0 <= p && p < N);
668  Row<N-1,ELT,1> out;
669  int nxt=0;
670  for (int i=0; i<N-1; ++i, ++nxt) {
671  if (nxt==p) ++nxt; // skip the loser
672  out[i] = (*this)[nxt];
673  }
674  return out;
675  }
676 
680  template <class EE> Row<N+1,ELT,1> append1(const EE& v) const {
681  Row<N+1,ELT,1> out;
682  Row<N,ELT,1>::updAs(&out[0]) = (*this);
683  out[N] = v;
684  return out;
685  }
686 
687 
693  template <class EE> Row<N+1,ELT,1> insert1(int p, const EE& v) const {
694  assert(0 <= p && p <= N);
695  if (p==N) return append1(v);
696  Row<N+1,ELT,1> out;
697  int nxt=0;
698  for (int i=0; i<N; ++i, ++nxt) {
699  if (i==p) out[nxt++] = v;
700  out[nxt] = (*this)[i];
701  }
702  return out;
703  }
704 
707  static const Row& getAs(const ELT* p) {return *reinterpret_cast<const Row*>(p);}
710  static Row& updAs(ELT* p) {return *reinterpret_cast<Row*>(p);}
711 
716 
718  bool isNaN() const {
719  for (int j=0; j<N; ++j)
720  if (CNT<ELT>::isNaN((*this)[j]))
721  return true;
722  return false;
723  }
724 
727  bool isInf() const {
728  bool seenInf = false;
729  for (int j=0; j<N; ++j) {
730  const ELT& e = (*this)[j];
731  if (!CNT<ELT>::isFinite(e)) {
732  if (!CNT<ELT>::isInf(e))
733  return false; // something bad was found
734  seenInf = true;
735  }
736  }
737  return seenInf;
738  }
739 
742  bool isFinite() const {
743  for (int j=0; j<N; ++j)
744  if (!CNT<ELT>::isFinite((*this)[j]))
745  return false;
746  return true;
747  }
748 
752 
755  template <class E2, int CS2>
756  bool isNumericallyEqual(const Row<N,E2,CS2>& r, double tol) const {
757  for (int j=0; j<N; ++j)
758  if (!CNT<ELT>::isNumericallyEqual((*this)(j), r(j), tol))
759  return false;
760  return true;
761  }
762 
766  template <class E2, int CS2>
767  bool isNumericallyEqual(const Row<N,E2,CS2>& r) const {
768  const double tol = std::max(getDefaultTolerance(),r.getDefaultTolerance());
769  return isNumericallyEqual(r, tol);
770  }
771 
776  bool isNumericallyEqual
777  (const ELT& e,
778  double tol = getDefaultTolerance()) const
779  {
780  for (int j=0; j<N; ++j)
781  if (!CNT<ELT>::isNumericallyEqual((*this)(j), e, tol))
782  return false;
783  return true;
784  }
785 private:
786  ELT d[NActualElements]; // data
787 };
788 
790 // Global operators involving two rows. //
791 // v+v, v-v, v==v, v!=v //
793 
794 // v3 = v1 + v2 where all v's have the same length N.
795 template <int N, class E1, int S1, class E2, int S2> inline
796 typename Row<N,E1,S1>::template Result< Row<N,E2,S2> >::Add
797 operator+(const Row<N,E1,S1>& l, const Row<N,E2,S2>& r) {
798  return Row<N,E1,S1>::template Result< Row<N,E2,S2> >
799  ::AddOp::perform(l,r);
800 }
801 
802 // v3 = v1 - v2, similar to +
803 template <int N, class E1, int S1, class E2, int S2> inline
804 typename Row<N,E1,S1>::template Result< Row<N,E2,S2> >::Sub
805 operator-(const Row<N,E1,S1>& l, const Row<N,E2,S2>& r) {
806  return Row<N,E1,S1>::template Result< Row<N,E2,S2> >
807  ::SubOp::perform(l,r);
808 }
809 
811 template <int N, class E1, int S1, class E2, int S2> inline bool
812 operator==(const Row<N,E1,S1>& l, const Row<N,E2,S2>& r) {
813  for (int i=0; i < N; ++i) if (l[i] != r[i]) return false;
814  return true;
815 }
817 template <int N, class E1, int S1, class E2, int S2> inline bool
818 operator!=(const Row<N,E1,S1>& l, const Row<N,E2,S2>& r) {return !(l==r);}
819 
821 template <int N, class E1, int S1, class E2, int S2> inline bool
822 operator<(const Row<N,E1,S1>& l, const Row<N,E2,S2>& r)
823 { for (int i=0; i < N; ++i) if (l[i] >= r[i]) return false;
824  return true; }
826 template <int N, class E1, int S1, class E2> inline bool
827 operator<(const Row<N,E1,S1>& v, const E2& e)
828 { for (int i=0; i < N; ++i) if (v[i] >= e) return false;
829  return true; }
830 
832 template <int N, class E1, int S1, class E2, int S2> inline bool
833 operator>(const Row<N,E1,S1>& l, const Row<N,E2,S2>& r)
834 { for (int i=0; i < N; ++i) if (l[i] <= r[i]) return false;
835  return true; }
837 template <int N, class E1, int S1, class E2> inline bool
838 operator>(const Row<N,E1,S1>& v, const E2& e)
839 { for (int i=0; i < N; ++i) if (v[i] <= e) return false;
840  return true; }
841 
844 template <int N, class E1, int S1, class E2, int S2> inline bool
845 operator<=(const Row<N,E1,S1>& l, const Row<N,E2,S2>& r)
846 { for (int i=0; i < N; ++i) if (l[i] > r[i]) return false;
847  return true; }
850 template <int N, class E1, int S1, class E2> inline bool
851 operator<=(const Row<N,E1,S1>& v, const E2& e)
852 { for (int i=0; i < N; ++i) if (v[i] > e) return false;
853  return true; }
854 
857 template <int N, class E1, int S1, class E2, int S2> inline bool
859 { for (int i=0; i < N; ++i) if (l[i] < r[i]) return false;
860  return true; }
863 template <int N, class E1, int S1, class E2> inline bool
864 operator>=(const Row<N,E1,S1>& v, const E2& e)
865 { for (int i=0; i < N; ++i) if (v[i] < e) return false;
866  return true; }
867 
869 // Global operators involving a row and a scalar. //
871 
872 // I haven't been able to figure out a nice way to templatize for the
873 // built-in reals without introducing a lot of unwanted type matches
874 // as well. So we'll just grind them out explicitly here.
875 
876 // SCALAR MULTIPLY
877 
878 // v = v*real, real*v
879 template <int N, class E, int S> inline
880 typename Row<N,E,S>::template Result<float>::Mul
881 operator*(const Row<N,E,S>& l, const float& r)
882  { return Row<N,E,S>::template Result<float>::MulOp::perform(l,r); }
883 template <int N, class E, int S> inline
884 typename Row<N,E,S>::template Result<float>::Mul
885 operator*(const float& l, const Row<N,E,S>& r) {return r*l;}
886 
887 template <int N, class E, int S> inline
888 typename Row<N,E,S>::template Result<double>::Mul
889 operator*(const Row<N,E,S>& l, const double& r)
890  { return Row<N,E,S>::template Result<double>::MulOp::perform(l,r); }
891 template <int N, class E, int S> inline
892 typename Row<N,E,S>::template Result<double>::Mul
893 operator*(const double& l, const Row<N,E,S>& r) {return r*l;}
894 
895 template <int N, class E, int S> inline
896 typename Row<N,E,S>::template Result<long double>::Mul
897 operator*(const Row<N,E,S>& l, const long double& r)
898  { return Row<N,E,S>::template Result<long double>::MulOp::perform(l,r); }
899 template <int N, class E, int S> inline
900 typename Row<N,E,S>::template Result<long double>::Mul
901 operator*(const long double& l, const Row<N,E,S>& r) {return r*l;}
902 
903 // v = v*int, int*v -- just convert int to v's precision float
904 template <int N, class E, int S> inline
905 typename Row<N,E,S>::template Result<typename CNT<E>::Precision>::Mul
906 operator*(const Row<N,E,S>& l, int r) {return l * (typename CNT<E>::Precision)r;}
907 template <int N, class E, int S> inline
908 typename Row<N,E,S>::template Result<typename CNT<E>::Precision>::Mul
909 operator*(int l, const Row<N,E,S>& r) {return r * (typename CNT<E>::Precision)l;}
910 
911 // Complex, conjugate, and negator are all easy to templatize.
912 
913 // v = v*complex, complex*v
914 template <int N, class E, int S, class R> inline
915 typename Row<N,E,S>::template Result<std::complex<R> >::Mul
916 operator*(const Row<N,E,S>& l, const std::complex<R>& r)
917  { return Row<N,E,S>::template Result<std::complex<R> >::MulOp::perform(l,r); }
918 template <int N, class E, int S, class R> inline
919 typename Row<N,E,S>::template Result<std::complex<R> >::Mul
920 operator*(const std::complex<R>& l, const Row<N,E,S>& r) {return r*l;}
921 
922 // v = v*conjugate, conjugate*v (convert conjugate->complex)
923 template <int N, class E, int S, class R> inline
924 typename Row<N,E,S>::template Result<std::complex<R> >::Mul
925 operator*(const Row<N,E,S>& l, const conjugate<R>& r) {return l*(std::complex<R>)r;}
926 template <int N, class E, int S, class R> inline
927 typename Row<N,E,S>::template Result<std::complex<R> >::Mul
928 operator*(const conjugate<R>& l, const Row<N,E,S>& r) {return r*(std::complex<R>)l;}
929 
930 // v = v*negator, negator*v: convert negator to standard number
931 template <int N, class E, int S, class R> inline
932 typename Row<N,E,S>::template Result<typename negator<R>::StdNumber>::Mul
933 operator*(const Row<N,E,S>& l, const negator<R>& r) {return l * (typename negator<R>::StdNumber)(R)r;}
934 template <int N, class E, int S, class R> inline
935 typename Row<N,E,S>::template Result<typename negator<R>::StdNumber>::Mul
936 operator*(const negator<R>& l, const Row<N,E,S>& r) {return r * (typename negator<R>::StdNumber)(R)l;}
937 
938 
939 // SCALAR DIVIDE. This is a scalar operation when the scalar is on the right,
940 // but when it is on the left it means scalar * pseudoInverse(row), which is
941 // a vec.
942 
943 // v = v/real, real/v
944 template <int N, class E, int S> inline
945 typename Row<N,E,S>::template Result<float>::Dvd
946 operator/(const Row<N,E,S>& l, const float& r)
947  { return Row<N,E,S>::template Result<float>::DvdOp::perform(l,r); }
948 template <int N, class E, int S> inline
949 typename CNT<float>::template Result<Row<N,E,S> >::Dvd
950 operator/(const float& l, const Row<N,E,S>& r)
951  { return CNT<float>::template Result<Row<N,E,S> >::DvdOp::perform(l,r); }
952 
953 template <int N, class E, int S> inline
954 typename Row<N,E,S>::template Result<double>::Dvd
955 operator/(const Row<N,E,S>& l, const double& r)
956  { return Row<N,E,S>::template Result<double>::DvdOp::perform(l,r); }
957 template <int N, class E, int S> inline
958 typename CNT<double>::template Result<Row<N,E,S> >::Dvd
959 operator/(const double& l, const Row<N,E,S>& r)
960  { return CNT<double>::template Result<Row<N,E,S> >::DvdOp::perform(l,r); }
961 
962 template <int N, class E, int S> inline
963 typename Row<N,E,S>::template Result<long double>::Dvd
964 operator/(const Row<N,E,S>& l, const long double& r)
965  { return Row<N,E,S>::template Result<long double>::DvdOp::perform(l,r); }
966 template <int N, class E, int S> inline
967 typename CNT<long double>::template Result<Row<N,E,S> >::Dvd
968 operator/(const long double& l, const Row<N,E,S>& r)
969  { return CNT<long double>::template Result<Row<N,E,S> >::DvdOp::perform(l,r); }
970 
971 // v = v/int, int/v -- just convert int to v's precision float
972 template <int N, class E, int S> inline
973 typename Row<N,E,S>::template Result<typename CNT<E>::Precision>::Dvd
974 operator/(const Row<N,E,S>& l, int r) {return l / (typename CNT<E>::Precision)r;}
975 template <int N, class E, int S> inline
976 typename CNT<typename CNT<E>::Precision>::template Result<Row<N,E,S> >::Dvd
977 operator/(int l, const Row<N,E,S>& r) {return (typename CNT<E>::Precision)l / r;}
978 
979 
980 // Complex, conjugate, and negator are all easy to templatize.
981 
982 // v = v/complex, complex/v
983 template <int N, class E, int S, class R> inline
984 typename Row<N,E,S>::template Result<std::complex<R> >::Dvd
985 operator/(const Row<N,E,S>& l, const std::complex<R>& r)
986  { return Row<N,E,S>::template Result<std::complex<R> >::DvdOp::perform(l,r); }
987 template <int N, class E, int S, class R> inline
988 typename CNT<std::complex<R> >::template Result<Row<N,E,S> >::Dvd
989 operator/(const std::complex<R>& l, const Row<N,E,S>& r)
990  { return CNT<std::complex<R> >::template Result<Row<N,E,S> >::DvdOp::perform(l,r); }
991 
992 // v = v/conjugate, conjugate/v (convert conjugate->complex)
993 template <int N, class E, int S, class R> inline
994 typename Row<N,E,S>::template Result<std::complex<R> >::Dvd
995 operator/(const Row<N,E,S>& l, const conjugate<R>& r) {return l/(std::complex<R>)r;}
996 template <int N, class E, int S, class R> inline
997 typename CNT<std::complex<R> >::template Result<Row<N,E,S> >::Dvd
998 operator/(const conjugate<R>& l, const Row<N,E,S>& r) {return (std::complex<R>)l/r;}
999 
1000 // v = v/negator, negator/v: convert negator to number
1001 template <int N, class E, int S, class R> inline
1002 typename Row<N,E,S>::template Result<typename negator<R>::StdNumber>::Dvd
1003 operator/(const Row<N,E,S>& l, const negator<R>& r) {return l/(typename negator<R>::StdNumber)(R)r;}
1004 template <int N, class E, int S, class R> inline
1005 typename CNT<R>::template Result<Row<N,E,S> >::Dvd
1006 operator/(const negator<R>& l, const Row<N,E,S>& r) {return (typename negator<R>::StdNumber)(R)l/r;}
1007 
1008 
1009 // Add and subtract are odd as scalar ops. They behave as though the
1010 // scalar stands for a vector each of whose elements is that scalar,
1011 // and then a normal vector add or subtract is done.
1012 
1013 // SCALAR ADD
1014 
1015 // v = v+real, real+v
1016 template <int N, class E, int S> inline
1017 typename Row<N,E,S>::template Result<float>::Add
1018 operator+(const Row<N,E,S>& l, const float& r)
1019  { return Row<N,E,S>::template Result<float>::AddOp::perform(l,r); }
1020 template <int N, class E, int S> inline
1021 typename Row<N,E,S>::template Result<float>::Add
1022 operator+(const float& l, const Row<N,E,S>& r) {return r+l;}
1023 
1024 template <int N, class E, int S> inline
1025 typename Row<N,E,S>::template Result<double>::Add
1026 operator+(const Row<N,E,S>& l, const double& r)
1027  { return Row<N,E,S>::template Result<double>::AddOp::perform(l,r); }
1028 template <int N, class E, int S> inline
1029 typename Row<N,E,S>::template Result<double>::Add
1030 operator+(const double& l, const Row<N,E,S>& r) {return r+l;}
1031 
1032 template <int N, class E, int S> inline
1033 typename Row<N,E,S>::template Result<long double>::Add
1034 operator+(const Row<N,E,S>& l, const long double& r)
1035  { return Row<N,E,S>::template Result<long double>::AddOp::perform(l,r); }
1036 template <int N, class E, int S> inline
1037 typename Row<N,E,S>::template Result<long double>::Add
1038 operator+(const long double& l, const Row<N,E,S>& r) {return r+l;}
1039 
1040 // v = v+int, int+v -- just convert int to v's precision float
1041 template <int N, class E, int S> inline
1042 typename Row<N,E,S>::template Result<typename CNT<E>::Precision>::Add
1043 operator+(const Row<N,E,S>& l, int r) {return l + (typename CNT<E>::Precision)r;}
1044 template <int N, class E, int S> inline
1045 typename Row<N,E,S>::template Result<typename CNT<E>::Precision>::Add
1046 operator+(int l, const Row<N,E,S>& r) {return r + (typename CNT<E>::Precision)l;}
1047 
1048 // Complex, conjugate, and negator are all easy to templatize.
1049 
1050 // v = v+complex, complex+v
1051 template <int N, class E, int S, class R> inline
1052 typename Row<N,E,S>::template Result<std::complex<R> >::Add
1053 operator+(const Row<N,E,S>& l, const std::complex<R>& r)
1054  { return Row<N,E,S>::template Result<std::complex<R> >::AddOp::perform(l,r); }
1055 template <int N, class E, int S, class R> inline
1056 typename Row<N,E,S>::template Result<std::complex<R> >::Add
1057 operator+(const std::complex<R>& l, const Row<N,E,S>& r) {return r+l;}
1058 
1059 // v = v+conjugate, conjugate+v (convert conjugate->complex)
1060 template <int N, class E, int S, class R> inline
1061 typename Row<N,E,S>::template Result<std::complex<R> >::Add
1062 operator+(const Row<N,E,S>& l, const conjugate<R>& r) {return l+(std::complex<R>)r;}
1063 template <int N, class E, int S, class R> inline
1064 typename Row<N,E,S>::template Result<std::complex<R> >::Add
1065 operator+(const conjugate<R>& l, const Row<N,E,S>& r) {return r+(std::complex<R>)l;}
1066 
1067 // v = v+negator, negator+v: convert negator to standard number
1068 template <int N, class E, int S, class R> inline
1069 typename Row<N,E,S>::template Result<typename negator<R>::StdNumber>::Add
1070 operator+(const Row<N,E,S>& l, const negator<R>& r) {return l + (typename negator<R>::StdNumber)(R)r;}
1071 template <int N, class E, int S, class R> inline
1072 typename Row<N,E,S>::template Result<typename negator<R>::StdNumber>::Add
1073 operator+(const negator<R>& l, const Row<N,E,S>& r) {return r + (typename negator<R>::StdNumber)(R)l;}
1074 
1075 // SCALAR SUBTRACT -- careful, not commutative.
1076 
1077 // v = v-real, real-v
1078 template <int N, class E, int S> inline
1079 typename Row<N,E,S>::template Result<float>::Sub
1080 operator-(const Row<N,E,S>& l, const float& r)
1081  { return Row<N,E,S>::template Result<float>::SubOp::perform(l,r); }
1082 template <int N, class E, int S> inline
1083 typename CNT<float>::template Result<Row<N,E,S> >::Sub
1084 operator-(const float& l, const Row<N,E,S>& r)
1085  { return CNT<float>::template Result<Row<N,E,S> >::SubOp::perform(l,r); }
1086 
1087 template <int N, class E, int S> inline
1088 typename Row<N,E,S>::template Result<double>::Sub
1089 operator-(const Row<N,E,S>& l, const double& r)
1090  { return Row<N,E,S>::template Result<double>::SubOp::perform(l,r); }
1091 template <int N, class E, int S> inline
1092 typename CNT<double>::template Result<Row<N,E,S> >::Sub
1093 operator-(const double& l, const Row<N,E,S>& r)
1094  { return CNT<double>::template Result<Row<N,E,S> >::SubOp::perform(l,r); }
1095 
1096 template <int N, class E, int S> inline
1097 typename Row<N,E,S>::template Result<long double>::Sub
1098 operator-(const Row<N,E,S>& l, const long double& r)
1099  { return Row<N,E,S>::template Result<long double>::SubOp::perform(l,r); }
1100 template <int N, class E, int S> inline
1101 typename CNT<long double>::template Result<Row<N,E,S> >::Sub
1102 operator-(const long double& l, const Row<N,E,S>& r)
1103  { return CNT<long double>::template Result<Row<N,E,S> >::SubOp::perform(l,r); }
1104 
1105 // v = v-int, int-v // just convert int to v's precision float
1106 template <int N, class E, int S> inline
1107 typename Row<N,E,S>::template Result<typename CNT<E>::Precision>::Sub
1108 operator-(const Row<N,E,S>& l, int r) {return l - (typename CNT<E>::Precision)r;}
1109 template <int N, class E, int S> inline
1110 typename CNT<typename CNT<E>::Precision>::template Result<Row<N,E,S> >::Sub
1111 operator-(int l, const Row<N,E,S>& r) {return (typename CNT<E>::Precision)l - r;}
1112 
1113 
1114 // Complex, conjugate, and negator are all easy to templatize.
1115 
1116 // v = v-complex, complex-v
1117 template <int N, class E, int S, class R> inline
1118 typename Row<N,E,S>::template Result<std::complex<R> >::Sub
1119 operator-(const Row<N,E,S>& l, const std::complex<R>& r)
1120  { return Row<N,E,S>::template Result<std::complex<R> >::SubOp::perform(l,r); }
1121 template <int N, class E, int S, class R> inline
1122 typename CNT<std::complex<R> >::template Result<Row<N,E,S> >::Sub
1123 operator-(const std::complex<R>& l, const Row<N,E,S>& r)
1124  { return CNT<std::complex<R> >::template Result<Row<N,E,S> >::SubOp::perform(l,r); }
1125 
1126 // v = v-conjugate, conjugate-v (convert conjugate->complex)
1127 template <int N, class E, int S, class R> inline
1128 typename Row<N,E,S>::template Result<std::complex<R> >::Sub
1129 operator-(const Row<N,E,S>& l, const conjugate<R>& r) {return l-(std::complex<R>)r;}
1130 template <int N, class E, int S, class R> inline
1131 typename CNT<std::complex<R> >::template Result<Row<N,E,S> >::Sub
1132 operator-(const conjugate<R>& l, const Row<N,E,S>& r) {return (std::complex<R>)l-r;}
1133 
1134 // v = v-negator, negator-v: convert negator to standard number
1135 template <int N, class E, int S, class R> inline
1136 typename Row<N,E,S>::template Result<typename negator<R>::StdNumber>::Sub
1137 operator-(const Row<N,E,S>& l, const negator<R>& r) {return l-(typename negator<R>::StdNumber)(R)r;}
1138 template <int N, class E, int S, class R> inline
1139 typename CNT<R>::template Result<Row<N,E,S> >::Sub
1140 operator-(const negator<R>& l, const Row<N,E,S>& r) {return (typename negator<R>::StdNumber)(R)l-r;}
1141 
1142 
1143 // Row I/O
1144 template <int N, class E, int S, class CHAR, class TRAITS> inline
1145 std::basic_ostream<CHAR,TRAITS>&
1146 operator<<(std::basic_ostream<CHAR,TRAITS>& o, const Row<N,E,S>& v) {
1147  o << "[" << v[0]; for(int i=1;i<N;++i) o<<','<<v[i]; o<<']'; return o;
1148 }
1149 
1152 template <int N, class E, int S, class CHAR, class TRAITS> inline
1153 std::basic_istream<CHAR,TRAITS>&
1154 operator>>(std::basic_istream<CHAR,TRAITS>& is, Row<N,E,S>& v) {
1155  CHAR openBracket, closeBracket;
1156  is >> openBracket; if (is.fail()) return is;
1157  if (openBracket==CHAR('('))
1158  closeBracket = CHAR(')');
1159  else if (openBracket==CHAR('['))
1160  closeBracket = CHAR(']');
1161  else {
1162  closeBracket = CHAR(0);
1163  is.unget(); if (is.fail()) return is;
1164  }
1165 
1166  for (int i=0; i < N; ++i) {
1167  is >> v[i];
1168  if (is.fail()) return is;
1169  if (i != N-1) {
1170  CHAR c; is >> c; if (is.fail()) return is;
1171  if (c != ',') is.unget();
1172  if (is.fail()) return is;
1173  }
1174  }
1175 
1176  // Get the closing bracket if there was an opening one. If we don't
1177  // see the expected character we'll set the fail bit in the istream.
1178  if (closeBracket != CHAR(0)) {
1179  CHAR closer; is >> closer; if (is.fail()) return is;
1180  if (closer != closeBracket) {
1181  is.unget(); if (is.fail()) return is;
1182  is.setstate( std::ios::failbit );
1183  }
1184  }
1185 
1186  return is;
1187 }
1188 
1189 } //namespace SimTK
1190 
1191 
1192 #endif //SimTK_SIMMATRIX_SMALLMATRIX_ROW_H_
AddCNTs< 1, N, ArgDepth, Row, ColSpacing, RowSpacing, CNT< P >::NRows, CNT< P >::NCols, CNT< P >::ArgDepth, P, CNT< P >::ColSpacing, CNT< P >::RowSpacing > AddOp
Definition: Row.h:291
AddOp::Type Add
Definition: Row.h:292
Vec< N, EHerm, STRIDE > THerm
Definition: Row.h:190
Matrix_< E > operator/(const MatrixBase< E > &l, const typename CNT< E >::StdNumber &r)
Definition: BigMatrix.h:613
TNeg & updNegate()
Definition: Row.h:489
Row< N, typename CNT< E >::template Result< EE >::Mul > scalarMultiply(const EE &e) const
Definition: Row.h:527
Row & scalarPlusEq(int ee)
Definition: Row.h:606
PhiMatrixTranspose transpose(const PhiMatrix &phi)
Definition: SpatialAlgebra.h:720
Row & scalarMinusEq(int ee)
Definition: Row.h:607
EPrecision Precision
Definition: Row.h:213
Row & scalarTimesEqFromLeft(int ee)
Definition: Row.h:611
K::ScalarNormSq ScalarNormSq
Definition: CompositeNumericalTypes.h:166
Row(const ENeg &ne)
Definition: Row.h:347
Row< N, typename CNT< E >::template Result< EE >::Add > conformingAdd(const Row< N, EE, SS > &r) const
Vector addition – use operator+ instead.
Definition: Row.h:404
K::ULessScalar ULessScalar
Definition: CompositeNumericalTypes.h:161
TImag & imag()
Definition: Row.h:508
DvdOp::Type Dvd
Definition: Row.h:287
EULessScalar ULessScalar
Definition: Row.h:210
Row TRow
Definition: Row.h:193
K::TReal TReal
Definition: CompositeNumericalTypes.h:141
Vec< N, EInvert, 1 > TInvert
Definition: Row.h:201
EStandard sum() const
Definition: Row.h:254
Row & scalarDivideEq(const EE &ee)
Definition: Row.h:598
Row & scalarEq(const EE &ee)
Definition: Row.h:586
Row< N, typename CNT< EE >::template Result< E >::Sub > scalarSubtractFromLeft(const EE &e) const
Definition: Row.h:569
Row< N, E, STRIDE > T
Definition: Row.h:181
Row(const E &e0, const E &e1, const E &e2, const E &e3, const E &e4, const E &e5)
Definition: Row.h:365
This is a small, fixed-size symmetric or Hermitian matrix designed for no-overhead inline computation...
Definition: SimTKcommon/include/SimTKcommon/internal/common.h:608
Row & scalarTimesEq(const EE &ee)
Definition: Row.h:594
Row & scalarDivideEqFromLeft(const EE &ee)
Definition: Row.h:600
Row & operator=(const Row< N, EE, SS > &vv)
Definition: Row.h:386
TNeg & operator-()
Definition: Row.h:484
EScalar Scalar
Definition: Row.h:209
Row & scalarPlusEq(const EE &ee)
Definition: Row.h:588
Row & operator+=(const Row< N, EE, SS > &r)
Definition: Row.h:390
Row< N, typename CNT< E >::template Result< EE >::Mul > elementwiseMultiply(const Row< N, EE, SS > &r) const
Elementwise multiply (Matlab .
Definition: Row.h:435
TAbs abs() const
Definition: Row.h:240
MulOpNonConforming::Type MulNon
Definition: Row.h:281
Row & operator-=(const Row< N, negator< EE >, SS > &r)
Definition: Row.h:396
THerm & updTranspose()
Definition: Row.h:492
This is the top-level SimTK namespace into which all SimTK names are placed to avoid collision with o...
Definition: Assembler.h:37
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
SymMat< N, ESqHermT > TSqHermT
Definition: Row.h:204
static TSqrt sqrt(const K &t)
Definition: CompositeNumericalTypes.h:239
Row(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: Row.h:374
K::Scalar Scalar
Definition: CompositeNumericalTypes.h:160
Row< N, EWithoutNegator, STRIDE > TWithoutNegator
Definition: Row.h:183
Row & operator=(const EE *p)
Definition: Row.h:382
K::TNormalize TNormalize
Definition: CompositeNumericalTypes.h:158
Matrix_< typename CNT< E1 >::template Result< E2 >::Sub > operator-(const MatrixBase< E1 > &l, const MatrixBase< E2 > &r)
Definition: BigMatrix.h:584
ENumber Number
Definition: Row.h:211
K::TImag TImag
Definition: CompositeNumericalTypes.h:142
MulCNTsNonConforming< 1, N, ArgDepth, Row, ColSpacing, RowSpacing, CNT< P >::NRows, CNT< P >::NCols, CNT< P >::ArgDepth, P, CNT< P >::ColSpacing, CNT< P >::RowSpacing > MulOpNonConforming
Definition: Row.h:280
Row & operator-=(const EE &e)
Definition: Row.h:580
static int size()
Definition: Row.h:216
Row< N, typename CNT< E >::template Result< P >::Add, 1 > Add
Definition: Row.h:266
CNT< ScalarNormSq >::TSqrt norm() const
Definition: Row.h:456
std::basic_istream< CHAR, TRAITS > & operator>>(std::basic_istream< CHAR, TRAITS > &is, conjugate< R > &c)
Definition: conjugate.h:800
Row< N, EReal, STRIDE *CNT< E >::RealStrideFactor > TReal
Definition: Row.h:186
const E & operator()(int i) const
Definition: Row.h:451
MulCNTs< 1, N, ArgDepth, Row, ColSpacing, RowSpacing, CNT< P >::NRows, CNT< P >::NCols, CNT< P >::ArgDepth, P, CNT< P >::ColSpacing, CNT< P >::RowSpacing > MulOp
Definition: Row.h:275
negator<N>, where N is a number type (real, complex, conjugate), is represented in memory identically...
Definition: String.h:44
Definition: CompositeNumericalTypes.h:120
Row & operator/=(const EE &e)
Definition: Row.h:582
static double getDefaultTolerance()
Definition: CompositeNumericalTypes.h:269
Row(const E &e0, const E &e1, const E &e2, const E &e3, const E &e4, const E &e5, const E &e6)
Definition: Row.h:368
bool isInf() const
Return true if any element of this Row contains a +Infinity or -Infinity somewhere but no element con...
Definition: Row.h:727
TPosTrans & updPositionalTranspose()
Definition: Row.h:496
bool operator==(const PhiMatrix &p1, const PhiMatrix &p2)
Definition: SpatialAlgebra.h:774
ScalarNormSq scalarNormSqr() const
Definition: Row.h:222
Row & operator=(const Row &src)
Definition: Row.h:319
static TStandard standardize(const K &t)
Definition: CompositeNumericalTypes.h:241
Row< N+1, ELT, 1 > insert1(int p, const EE &v) const
Return a row one larger than this one by inserting an element before the indicated one...
Definition: Row.h:693
const THerm & operator~() const
Definition: Row.h:485
Row(int i)
Definition: Row.h:352
static Row & updSubRow(Row< NN, ELT, STRIDE > &r, int j)
Extract a subvector of type Row from a longer one that has the same element type and stride...
Definition: Row.h:658
Row & scalarMinusEqFromLeft(int ee)
Definition: Row.h:610
Row< N, typename CNT< EE >::template Result< E >::Dvd > scalarDivideFromLeft(const EE &e) const
Definition: Row.h:548
Row(const Row< N, EE, SS > &vv)
Definition: Row.h:338
Row< N, P > Type
Definition: Row.h:302
ELEM sum(const VectorBase< ELEM > &v)
Definition: VectorMath.h:147
K::TSqTHerm TSqTHerm
Definition: CompositeNumericalTypes.h:147
Row< N, EAbs, 1 > TAbs
Definition: Row.h:199
Row(const E &e0, const E &e1, const E &e2)
Definition: Row.h:358
TNormalize normalize() const
Definition: Row.h:469
THerm & operator~()
Definition: Row.h:486
Row(const EE *p)
Definition: Row.h:380
DvdCNTs< 1, N, ArgDepth, Row, ColSpacing, RowSpacing, CNT< P >::NRows, CNT< P >::NCols, CNT< P >::ArgDepth, P, CNT< P >::ColSpacing, CNT< P >::RowSpacing > DvdOp
Definition: Row.h:286
Row< N, typename CNT< E >::template Result< EE >::Sub > conformingSubtract(const Row< N, EE, SS > &r) const
Vector subtraction – use operator- instead.
Definition: Row.h:412
This is a fixed-length column vector designed for no-overhead inline computation. ...
Definition: SimTKcommon/include/SimTKcommon/internal/common.h:605
Row< N, EComplex, STRIDE > TComplex
Definition: Row.h:189
Row & operator-=(const Row< N, EE, SS > &r)
Definition: Row.h:394
K::Precision Precision
Definition: CompositeNumericalTypes.h:164
static Row & updAs(ELT *p)
Recast a writable ordinary C++ array E[] to a writable Row<N,E,S>; assumes compatible length...
Definition: Row.h:710
Row< N, typename CNT< EE >::template Result< E >::Mul > scalarMultiplyFromLeft(const EE &e) const
Definition: Row.h:533
SubOp::Type Sub
Definition: Row.h:297
Matrix_< E > operator*(const MatrixBase< E > &l, const typename CNT< E >::StdNumber &r)
Definition: BigMatrix.h:605
Row & scalarTimesEq(int ee)
Definition: Row.h:608
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
TStandard standardize() const
Definition: Row.h:246
bool isNaN() const
Return true if any element of this Row contains a NaN anywhere.
Definition: Row.h:718
CNT< E >::template Result< EE >::Mul conformingMultiply(const Vec< N, EE, SS > &r) const
Same as dot product (s = row*col) – use operator* or dot() instead.
Definition: Row.h:420
ELEM max(const VectorBase< ELEM > &v)
Definition: VectorMath.h:251
Row< N, typename CNT< E >::template Result< P >::Sub, 1 > Sub
Definition: Row.h:267
Row & scalarMinusEq(const EE &ee)
Definition: Row.h:590
E & operator()(int i)
Definition: Row.h:452
static double getDefaultTolerance()
For approximate comparisions, the default tolerance to use for a vector is the same as its elements&#39; ...
Definition: Row.h:751
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
Row< MatNCol, typename CNT< E >::template Result< EE >::Mul > conformingMultiply(const Mat< N, MatNCol, EE, CS, RS > &m) const
Row times a conforming matrix, row=row*mat – use operator* instead.
Definition: Row.h:427
EScalarNormSq ScalarNormSq
Definition: Row.h:214
TSqrt sqrt() const
Definition: Row.h:231
MulOp::Type Mul
Definition: Row.h:276
K::TPosTrans TPosTrans
Definition: CompositeNumericalTypes.h:145
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
Row()
Definition: Row.h:307
Row(const Row &src)
Definition: Row.h:316
TWithoutNegator & updCastAwayNegatorIfAny()
Definition: Row.h:515
const Row & operator+() const
Definition: Row.h:482
Row< N, typename CNT< E >::template Result< EE >::Dvd > scalarDivide(const EE &e) const
Definition: Row.h:542
K::StdNumber StdNumber
Definition: CompositeNumericalTypes.h:163
Row & scalarEq(int ee)
Definition: Row.h:605
RowVectorBase< typename CNT< ELEM >::TAbs > abs(const RowVectorBase< ELEM > &v)
Definition: VectorMath.h:120
Row< N, typename CNT< E >::template Result< EE >::Dvd > elementwiseDivide(const Row< N, EE, SS > &r) const
Elementwise divide (Matlab .
Definition: Row.h:443
Definition: Row.h:301
bool operator!=(const conjugate< R > &a, const float &b)
Definition: conjugate.h:859
Row< N, typename CNT< E >::template Result< EE >::Add > scalarAdd(const EE &e) const
Definition: Row.h:555
Row< N, typename CNT< E >::template Result< P >::Mul, 1 > Mul
Definition: Row.h:264
Specialized information about Composite Numerical Types which allows us to define appropriate templat...
Definition: CompositeNumericalTypes.h:136
Row & scalarMinusEqFromLeft(const EE &ee)
Definition: Row.h:592
bool isFinite() const
Return true if no element of this Row contains an Infinity or a NaN anywhere.
Definition: Row.h:742
TInvert invert() const
Definition: Row.h:480
const TNeg & negate() const
Definition: Row.h:488
This is a fixed-length row vector designed for no-overhead inline computation.
Definition: SimTKcommon/include/SimTKcommon/internal/common.h:606
Row< N, EStandard, 1 > TStandard
Definition: Row.h:200
Row< N, EImag, STRIDE *CNT< E >::RealStrideFactor > TImag
Definition: Row.h:188
E & operator[](int i)
Definition: Row.h:450
const Real E
e = Real(exp(1))
const TWithoutNegator & castAwayNegatorIfAny() const
Definition: Row.h:514
Mandatory first inclusion for any Simbody source or header file.
bool isNumericallyEqual(const Row< N, E2, CS2 > &r) const
Test whether this row vector is numerically equal to some other row with the same shape...
Definition: Row.h:767
Row< N, typename CNT< E >::template Result< P >::Dvd, 1 > Dvd
Definition: Row.h:265
static const Row & getAs(const ELT *p)
Recast an ordinary C++ array E[] to a const Row<N,E,S>; assumes compatible length, stride, and packing.
Definition: Row.h:707
Row & operator+=(const Row< N, negator< EE >, SS > &r)
Definition: Row.h:392
K::TNeg TNeg
Definition: CompositeNumericalTypes.h:139
Row< N-1, ELT, 1 > drop1(int p) const
Return a row one smaller than this one by dropping the element at the indicated position p...
Definition: Row.h:666
static Row< N, ELT, 1 > getNaN()
Return a Row of the same length and element type as this one but with all elements set to NaN...
Definition: Row.h:715
K::TStandard TStandard
Definition: CompositeNumericalTypes.h:156
Row(const Row< N, ENeg, SS > &src)
Definition: Row.h:332
Row & operator*=(const EE &e)
Definition: Row.h:581
void copy(Row< 1, E1, S1 > &r1, const Row< 1, E2, S2 > &r2)
Definition: Row.h:105
K::TWithoutNegator TWithoutNegator
Definition: CompositeNumericalTypes.h:140
const Row< NN, ELT, STRIDE > & getSubRow(int j) const
Extract a const reference to a sub-Row with size known at compile time.
Definition: Row.h:631
Row & scalarTimesEqFromLeft(const EE &ee)
Definition: Row.h:596
Row< N, ENormalize, 1 > TNormalize
Definition: Row.h:202
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
Row< N, ENeg, STRIDE > TNeg
Definition: Row.h:182
Row & scalarDivideEq(int ee)
Definition: Row.h:609
Matrix_< typename CNT< E1 >::template Result< E2 >::Add > operator+(const MatrixBase< E1 > &l, const MatrixBase< E2 > &r)
Definition: BigMatrix.h:568
Row(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: Row.h:371
SubCNTs< 1, N, ArgDepth, Row, ColSpacing, RowSpacing, CNT< P >::NRows, CNT< P >::NCols, CNT< P >::ArgDepth, P, CNT< P >::ColSpacing, CNT< P >::RowSpacing > SubOp
Definition: Row.h:296
void setToZero()
Set every scalar in this Row to zero.
Definition: Row.h:621
Row & scalarDivideEqFromLeft(int ee)
Definition: Row.h:612
EStdNumber StdNumber
Definition: Row.h:212
static int nrow()
Definition: Row.h:217
const THerm & transpose() const
Definition: Row.h:491
E TElement
Definition: Row.h:192
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
E TCol
Definition: Row.h:194
K::TComplex TComplex
Definition: CompositeNumericalTypes.h:143
bool isNumericallyEqual(const Row< N, E2, CS2 > &r, double tol) const
Test whether this row is numerically equal to some other row with the same shape, using a specified t...
Definition: Row.h:756
TReal & real()
Definition: Row.h:500
K::Number Number
Definition: CompositeNumericalTypes.h:162
Definition: Row.h:263
const TReal & real() const
Definition: Row.h:499
static K getNaN()
Definition: CompositeNumericalTypes.h:246
Definition: Row.h:272
Row(const E &e0, const E &e1)
Definition: Row.h:356
Vec< N, ESqrt, 1 > TSqrt
Definition: Row.h:198
Row(const E &e)
Definition: Row.h:343
ScalarNormSq normSqr() const
Definition: Row.h:454
K::TSqHermT TSqHermT
Definition: CompositeNumericalTypes.h:146
const E & operator[](int i) const
Definition: Row.h:449
const TImag & imag() const
Definition: Row.h:503
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
const TPosTrans & positionalTranspose() const
Definition: Row.h:494
K::THerm THerm
Definition: CompositeNumericalTypes.h:144
const TNeg & operator-() const
Definition: Row.h:483
Row(const E &e0, const E &e1, const E &e2, const E &e3)
Definition: Row.h:360
Definition: negator.h:64
Row & operator+=(const EE &e)
Definition: Row.h:579
static int ncol()
Definition: Row.h:218
Row< NN, ELT, STRIDE > & updSubRow(int j)
Extract a writable reference to a sub-Row with size known at compile time.
Definition: Row.h:641
Row(const E &e0, const E &e1, const E &e2, const E &e3, const E &e4)
Definition: Row.h:362
void setToNaN()
Set every scalar in this Row to NaN; this is the default initial value in Debug builds, but not in Release.
Definition: Row.h:616
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
Row(const Row< N, E, SS > &src)
Definition: Row.h:326
Vec< N, E, STRIDE > TPosTrans
Definition: Row.h:191
Row< N, typename CNT< E >::template Result< EE >::Sub > scalarSubtract(const EE &e) const
Definition: Row.h:563
Row< N+1, ELT, 1 > append1(const EE &v) const
Return a row one larger than this one by adding an element to the end.
Definition: Row.h:680
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
EScalarNormSq TSqTHerm
Definition: Row.h:205
static const Row & getSubRow(const Row< NN, ELT, STRIDE > &r, int j)
Extract a subvector of type Row from a longer one that has the same element type and stride...
Definition: Row.h:650