Simbody  3.4
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
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 
31 #include "SimTKcommon/internal/common.h"
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
60 conformingSubtract(const Row<1,E1,S1>& r1, const Row<1,E2,S2>& r2,
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
65 conformingSubtract(const Row<N,E1,S1>& r1, const Row<N,E2,S2>& r2,
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
75 elementwiseMultiply(const Row<1,E1,S1>& r1, const Row<1,E2,S2>& r2,
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
80 elementwiseMultiply(const Row<N,E1,S1>& r1, const Row<N,E2,S2>& r2,
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
90 elementwiseDivide(const Row<1,E1,S1>& r1, const Row<1,E2,S2>& r2,
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
95 elementwiseDivide(const Row<N,E1,S1>& r1, const Row<N,E2,S2>& r2,
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,
167  ImagOffset = NTraits<ENumber>::ImagOffset,
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
172  : MAX_RESOLVED_DEPTH),
173  IsScalar = 0,
174  IsULessScalar = 0,
175  IsNumber = 0,
176  IsStdNumber = 0,
177  IsPrecision = 0,
178  SignInterpretation = CNT<E>::SignInterpretation
179  };
180 
181  typedef Row<N,E,STRIDE> T;
182  typedef Row<N,ENeg,STRIDE> TNeg;
183  typedef Row<N,EWithoutNegator,STRIDE> TWithoutNegator;
184 
185  typedef Row<N,EReal,STRIDE*CNT<E>::RealStrideFactor>
186  TReal;
187  typedef Row<N,EImag,STRIDE*CNT<E>::RealStrideFactor>
188  TImag;
189  typedef Row<N,EComplex,STRIDE> TComplex;
190  typedef Vec<N,EHerm,STRIDE> THerm;
191  typedef Vec<N,E,STRIDE> TPosTrans;
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
200  typedef Row<N,EStandard,1> TStandard;
201  typedef Vec<N,EInvert,1> TInvert; // packed
202  typedef Row<N,ENormalize,1> TNormalize;
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 {
264  typedef Row<N, typename CNT<E>::template Result<P>::Mul, 1> Mul;
265  typedef Row<N, typename CNT<E>::template Result<P>::Dvd, 1> Dvd;
266  typedef Row<N, typename CNT<E>::template Result<P>::Add, 1> Add;
267  typedef Row<N, typename CNT<E>::template Result<P>::Sub, 1> Sub;
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,
274  CNT<P>::NRows, CNT<P>::NCols, CNT<P>::ArgDepth,
275  P, CNT<P>::ColSpacing, CNT<P>::RowSpacing> MulOp;
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>
412  conformingSubtract(const Row<N,EE,SS>& r) const {
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
420  conformingMultiply(const Vec<N,EE,SS>& r) const {
421  return (*this)*r;
422  }
423 
425  template <int MatNCol, class EE, int CS, int RS>
426  Row<MatNCol,typename CNT<E>::template Result<EE>::Mul>
427  conformingMultiply(const Mat<N,MatNCol,EE,CS,RS>& m) const {
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>
435  elementwiseMultiply(const Row<N,EE,SS>& r) const {
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); }
496  TPosTrans& updPositionalTranspose()
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>
641  Row<NN,ELT,STRIDE>& updSubRow(int j) {
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 
715  static Row<N,ELT,1> getNaN() { return Row<N,ELT,1>(CNT<ELT>::getNaN()); }
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 
751  static double getDefaultTolerance() {return CNT<ELT>::getDefaultTolerance();}
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
858 operator>=(const Row<N,E1,S1>& l, const Row<N,E2,S2>& r)
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_