Row.h

Go to the documentation of this file.
00001 #ifndef SimTK_SIMMATRIX_SMALLMATRIX_ROW_H_
00002 #define SimTK_SIMMATRIX_SMALLMATRIX_ROW_H_
00003 
00004 /* -------------------------------------------------------------------------- *
00005  *                      SimTK Core: SimTK Simmatrix(tm)                       *
00006  * -------------------------------------------------------------------------- *
00007  * This is part of the SimTK Core biosimulation toolkit originating from      *
00008  * Simbios, the NIH National Center for Physics-Based Simulation of           *
00009  * Biological Structures at Stanford, funded under the NIH Roadmap for        *
00010  * Medical Research, grant U54 GM072970. See https://simtk.org.               *
00011  *                                                                            *
00012  * Portions copyright (c) 2005-8 Stanford University and the Authors.         *
00013  * Authors: Michael Sherman                                                   *
00014  * Contributors: Peter Eastman                                                *
00015  *                                                                            *
00016  * Permission is hereby granted, free of charge, to any person obtaining a    *
00017  * copy of this software and associated documentation files (the "Software"), *
00018  * to deal in the Software without restriction, including without limitation  *
00019  * the rights to use, copy, modify, merge, publish, distribute, sublicense,   *
00020  * and/or sell copies of the Software, and to permit persons to whom the      *
00021  * Software is furnished to do so, subject to the following conditions:       *
00022  *                                                                            *
00023  * The above copyright notice and this permission notice shall be included in *
00024  * all copies or substantial portions of the Software.                        *
00025  *                                                                            *
00026  * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR *
00027  * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,   *
00028  * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL    *
00029  * THE AUTHORS, CONTRIBUTORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM,    *
00030  * DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR      *
00031  * OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE  *
00032  * USE OR OTHER DEALINGS IN THE SOFTWARE.                                     *
00033  * -------------------------------------------------------------------------- */
00034 
00039 #include "SimTKcommon/internal/common.h"
00040 
00041 
00042 namespace SimTK {
00043 
00044 // The following functions are used internally by Row.
00045 
00046 namespace Impl {
00047 
00048 // For those wimpy compilers that don't unroll short, constant-limit loops, Peter Eastman added these
00049 // recursive template implementations of add and subtract.
00050 
00051 template <class E1, int S1, class E2, int S2> void
00052 conformingAdd(const Row<1,E1,S1>& r1, const Row<1,E2,S2>& r2, Row<1,typename CNT<E1>::template Result<E2>::Add>& result) {
00053     result[0] = r1[0] + r2[0];
00054 }
00055 template <int N, class E1, int S1, class E2, int S2> void
00056 conformingAdd(const Row<N,E1,S1>& r1, const Row<N,E2,S2>& r2, Row<N,typename CNT<E1>::template Result<E2>::Add>& result) {
00057     conformingAdd(reinterpret_cast<const Row<N-1,E1,S1>&>(r1), reinterpret_cast<const Row<N-1,E2,S2>&>(r2), reinterpret_cast<Row<N-1,typename CNT<E1>::template Result<E2>::Add>&>(result));
00058     result[N-1] = r1[N-1] + r2[N-1];
00059 }
00060 template <class E1, int S1, class E2, int S2> void
00061 conformingSubtract(const Row<1,E1,S1>& r1, const Row<1,E2,S2>& r2, Row<1,typename CNT<E1>::template Result<E2>::Add>& result) {
00062     result[0] = r1[0] - r2[0];
00063 }
00064 template <int N, class E1, int S1, class E2, int S2> void
00065 conformingSubtract(const Row<N,E1,S1>& r1, const Row<N,E2,S2>& r2, Row<N,typename CNT<E1>::template Result<E2>::Add>& result) {
00066     conformingSubtract(reinterpret_cast<const Row<N-1,E1,S1>&>(r1), reinterpret_cast<const Row<N-1,E2,S2>&>(r2), reinterpret_cast<Row<N-1,typename CNT<E1>::template Result<E2>::Add>&>(result));
00067     result[N-1] = r1[N-1] - r2[N-1];
00068 }
00069 template <class E1, int S1, class E2, int S2> void
00070 copy(Row<1,E1,S1>& r1, const Row<1,E2,S2>& r2) {
00071     r1[0] = r2[0];
00072 }
00073 template <int N, class E1, int S1, class E2, int S2> void
00074 copy(Row<N,E1,S1>& r1, const Row<N,E2,S2>& r2) {
00075     copy(reinterpret_cast<Row<N-1,E1,S1>&>(r1), reinterpret_cast<const Row<N-1,E2,S2>&>(r2));
00076     r1[N-1] = r2[N-1];
00077 }
00078 
00079 }
00080 
00082 template <int N, class ELT, int STRIDE> class Row {
00083     typedef ELT                                 E;
00084     typedef typename CNT<E>::TNeg               ENeg;
00085     typedef typename CNT<E>::TWithoutNegator    EWithoutNegator;
00086     typedef typename CNT<E>::TReal              EReal;
00087     typedef typename CNT<E>::TImag              EImag;
00088     typedef typename CNT<E>::TComplex           EComplex;
00089     typedef typename CNT<E>::THerm              EHerm;
00090     typedef typename CNT<E>::TPosTrans          EPosTrans;
00091     typedef typename CNT<E>::TSqHermT           ESqHermT;
00092     typedef typename CNT<E>::TSqTHerm           ESqTHerm;
00093 
00094     typedef typename CNT<E>::TSqrt              ESqrt;
00095     typedef typename CNT<E>::TAbs               EAbs;
00096     typedef typename CNT<E>::TStandard          EStandard;
00097     typedef typename CNT<E>::TInvert            EInvert;
00098     typedef typename CNT<E>::TNormalize         ENormalize;
00099 
00100     typedef typename CNT<E>::Scalar             EScalar;
00101     typedef typename CNT<E>::ULessScalar        EULessScalar;
00102     typedef typename CNT<E>::Number             ENumber;
00103     typedef typename CNT<E>::StdNumber          EStdNumber;
00104     typedef typename CNT<E>::Precision          EPrecision;
00105     typedef typename CNT<E>::ScalarNormSq       EScalarNormSq;
00106 
00107 public:
00108 
00109     enum {
00110         NRows               = 1,
00111         NCols               = N,
00112         NPackedElements     = N,
00113         NActualElements     = N * STRIDE,   // includes trailing gap
00114         NActualScalars      = CNT<E>::NActualScalars * NActualElements,
00115         RowSpacing          = NActualElements,
00116         ColSpacing          = STRIDE,
00117         ImagOffset          = NTraits<ENumber>::ImagOffset,
00118         RealStrideFactor    = 1, // composite types don't change size when
00119                                  // cast from complex to real or imaginary
00120         ArgDepth            = ((int)CNT<E>::ArgDepth < (int)MAX_RESOLVED_DEPTH 
00121                                 ? CNT<E>::ArgDepth + 1 
00122                                 : MAX_RESOLVED_DEPTH),
00123         IsScalar            = 0,
00124         IsULessScalar       = 0,
00125         IsNumber            = 0,
00126         IsStdNumber         = 0,
00127         IsPrecision         = 0,
00128         SignInterpretation  = CNT<E>::SignInterpretation
00129     };
00130 
00131     typedef Row<N,E,STRIDE>                 T;
00132     typedef Row<N,ENeg,STRIDE>              TNeg;
00133     typedef Row<N,EWithoutNegator,STRIDE>   TWithoutNegator;
00134 
00135     typedef Row<N,EReal,STRIDE*CNT<E>::RealStrideFactor>         
00136                                             TReal;
00137     typedef Row<N,EImag,STRIDE*CNT<E>::RealStrideFactor>         
00138                                             TImag;
00139     typedef Row<N,EComplex,STRIDE>          TComplex;
00140     typedef Vec<N,EHerm,STRIDE>             THerm;
00141     typedef Vec<N,E,STRIDE>                 TPosTrans;
00142     typedef E                               TElement;
00143     typedef Row                             TRow;
00144     typedef E                               TCol;
00145 
00146     // These are the results of calculations, so are returned in new, packed
00147     // memory. Be sure to refer to element types here which are also packed.
00148     typedef Vec<N,ESqrt,1>                  TSqrt;      // Note stride
00149     typedef Row<N,EAbs,1>                   TAbs;       // Note stride
00150     typedef Row<N,EStandard,1>              TStandard;
00151     typedef Vec<N,EInvert,1>                TInvert;    // packed
00152     typedef Row<N,ENormalize,1>             TNormalize;
00153 
00154     typedef SymMat<N,ESqHermT>              TSqHermT;   // result of self outer product
00155     typedef EScalarNormSq                   TSqTHerm;   // result of self dot product
00156 
00157     // These recurse right down to the underlying scalar type no matter how
00158     // deep the elements are.
00159     typedef EScalar                         Scalar;
00160     typedef EULessScalar                    ULessScalar;
00161     typedef ENumber                         Number;
00162     typedef EStdNumber                      StdNumber;
00163     typedef EPrecision                      Precision;
00164     typedef EScalarNormSq                   ScalarNormSq;
00165 
00166     int size()   const  { return N; }
00167     int nrow()   const  { return 1; }
00168     int ncol()   const  { return N; }
00169 
00170 
00171     // Scalar norm square is sum( conjugate squares of all scalars )
00172     ScalarNormSq scalarNormSqr() const { 
00173         ScalarNormSq sum(0);
00174         for(int i=0;i<N;++i) sum += CNT<E>::scalarNormSqr(d[i*STRIDE]);
00175         return sum;
00176     }
00177 
00178     // sqrt() is elementwise square root; that is, the return value has the same
00179     // dimension as this Vec but with each element replaced by whatever it thinks
00180     // its square root is.
00181     TSqrt sqrt() const {
00182         TSqrt rsqrt;
00183         for(int i=0;i<N;++i) rsqrt[i] = CNT<E>::sqrt(d[i*STRIDE]);
00184         return rsqrt;
00185     }
00186 
00187     // abs() is elementwise absolute value; that is, the return value has the same
00188     // dimension as this Row but with each element replaced by whatever it thinks
00189     // its absolute value is.
00190     TAbs abs() const {
00191         TAbs rabs;
00192         for(int i=0;i<N;++i) rabs[i] = CNT<E>::abs(d[i*STRIDE]);
00193         return rabs;
00194     }
00195 
00196     TStandard standardize() const {
00197         TStandard rstd;
00198         for(int i=0;i<N;++i) rstd[i] = CNT<E>::standardize(d[i*STRIDE]);
00199         return rstd;
00200     }
00201 
00202     // Sum just adds up all the elements, getting rid of negators and
00203     // conjugates in the process.
00204     EStandard sum() const {
00205         E sum(0);
00206         for (int i=0;i<N;++i) sum += d[i*STRIDE];
00207         return CNT<E>::standardize(sum);
00208     }
00209 
00210     // This gives the resulting rowvector type when (v[i] op P) is applied to each element of v.
00211     // It is a row of length N, stride 1, and element types which are the regular
00212     // composite result of E op P. Typically P is a scalar type but it doesn't have to be.
00213     template <class P> struct EltResult { 
00214         typedef Row<N, typename CNT<E>::template Result<P>::Mul, 1> Mul;
00215         typedef Row<N, typename CNT<E>::template Result<P>::Dvd, 1> Dvd;
00216         typedef Row<N, typename CNT<E>::template Result<P>::Add, 1> Add;
00217         typedef Row<N, typename CNT<E>::template Result<P>::Sub, 1> Sub;
00218     };
00219 
00220     // This is the composite result for v op P where P is some kind of appropriately shaped
00221     // non-scalar type.
00222     template <class P> struct Result { 
00223         typedef MulCNTs<1,N,ArgDepth,Row,ColSpacing,RowSpacing,
00224             CNT<P>::NRows, CNT<P>::NCols, CNT<P>::ArgDepth,
00225             P, CNT<P>::ColSpacing, CNT<P>::RowSpacing> MulOp;
00226         typedef typename MulOp::Type Mul;
00227 
00228         typedef MulCNTsNonConforming<1,N,ArgDepth,Row,ColSpacing,RowSpacing,
00229             CNT<P>::NRows, CNT<P>::NCols, CNT<P>::ArgDepth,
00230             P, CNT<P>::ColSpacing, CNT<P>::RowSpacing> MulOpNonConforming;
00231         typedef typename MulOpNonConforming::Type MulNon;
00232 
00233 
00234         typedef DvdCNTs<1,N,ArgDepth,Row,ColSpacing,RowSpacing,
00235             CNT<P>::NRows, CNT<P>::NCols, CNT<P>::ArgDepth,
00236             P, CNT<P>::ColSpacing, CNT<P>::RowSpacing> DvdOp;
00237         typedef typename DvdOp::Type Dvd;
00238 
00239         typedef AddCNTs<1,N,ArgDepth,Row,ColSpacing,RowSpacing,
00240             CNT<P>::NRows, CNT<P>::NCols, CNT<P>::ArgDepth,
00241             P, CNT<P>::ColSpacing, CNT<P>::RowSpacing> AddOp;
00242         typedef typename AddOp::Type Add;
00243 
00244         typedef SubCNTs<1,N,ArgDepth,Row,ColSpacing,RowSpacing,
00245             CNT<P>::NRows, CNT<P>::NCols, CNT<P>::ArgDepth,
00246             P, CNT<P>::ColSpacing, CNT<P>::RowSpacing> SubOp;
00247         typedef typename SubOp::Type Sub;
00248     };
00249 
00250     // Shape-preserving element substitution (always packed)
00251     template <class P> struct Substitute {
00252         typedef Row<N,P> Type;
00253     };
00254 
00255     // Default construction initializes to NaN when debugging but
00256     // is left uninitialized otherwise.
00257     Row(){ 
00258     #ifndef NDEBUG
00259         setToNaN();
00260     #endif
00261     }
00262 
00263     // It's important not to use the default copy constructor or copy
00264     // assignment because the compiler doesn't understand that we may
00265     // have noncontiguous storage and will try to copy the whole array.
00266     Row(const Row& src) {
00267         Impl::copy(*this, src);
00268     }
00269     Row& operator=(const Row& src) {    // no harm if src and 'this' are the same
00270         Impl::copy(*this, src);
00271         return *this;
00272     }
00273 
00274     // We want an implicit conversion from a Row of the same length
00275     // and element type but with a different stride.
00276     template <int SS> Row(const Row<N,E,SS>& src) {
00277         Impl::copy(*this, src);
00278     }
00279 
00280     // We want an implicit conversion from a Row of the same length
00281     // and *negated* element type, possibly with a different stride.
00282     template <int SS> Row(const Row<N,ENeg,SS>& src) {
00283         Impl::copy(*this, src);
00284     }
00285 
00286     // Construct a Row from a Row of the same length, with any
00287     // stride. Works as long as the element types are compatible.
00288     template <class EE, int SS> explicit Row(const Row<N,EE,SS>& vv) {
00289         Impl::copy(*this, vv);
00290     }
00291 
00292     // Construction using an element assigns to each element.
00293     explicit Row(const ELT& e)
00294       { for (int i=0;i<N;++i) d[i*STRIDE]=e; }
00295 
00296     // A bevy of constructors for Rows up to length 6.
00297     Row(const E& e0,const E& e1)
00298       { assert(N==2);(*this)[0]=e0;(*this)[1]=e1; }
00299     Row(const E& e0,const E& e1,const E& e2)
00300       { assert(N==3);(*this)[0]=e0;(*this)[1]=e1;(*this)[2]=e2; }
00301     Row(const E& e0,const E& e1,const E& e2,const E& e3)
00302       { assert(N==4);(*this)[0]=e0;(*this)[1]=e1;(*this)[2]=e2;(*this)[3]=e3; }
00303     Row(const E& e0,const E& e1,const E& e2,const E& e3,const E& e4)
00304       { assert(N==5);(*this)[0]=e0;(*this)[1]=e1;(*this)[2]=e2;
00305         (*this)[3]=e3;(*this)[4]=e4; }
00306     Row(const E& e0,const E& e1,const E& e2,const E& e3,const E& e4,const E& e5)
00307       { assert(N==6);(*this)[0]=e0;(*this)[1]=e1;(*this)[2]=e2;
00308         (*this)[3]=e3;(*this)[4]=e4;(*this)[5]=e5; }
00309     Row(const E& e0,const E& e1,const E& e2,const E& e3,const E& e4,const E& e5,const E& e6)
00310       { assert(N==7);(*this)[0]=e0;(*this)[1]=e1;(*this)[2]=e2;
00311         (*this)[3]=e3;(*this)[4]=e4;(*this)[5]=e5;(*this)[6]=e6; }
00312     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)
00313       { assert(N==8);(*this)[0]=e0;(*this)[1]=e1;(*this)[2]=e2;
00314         (*this)[3]=e3;(*this)[4]=e4;(*this)[5]=e5;(*this)[6]=e6;(*this)[7]=e7; }
00315     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)
00316       { assert(N==9);(*this)[0]=e0;(*this)[1]=e1;(*this)[2]=e2;
00317         (*this)[3]=e3;(*this)[4]=e4;(*this)[5]=e5;(*this)[6]=e6;(*this)[7]=e7;(*this)[8]=e8; }
00318 
00319     // Construction from a pointer to anything assumes we're pointing
00320     // at an element list of the right length.
00321     template <class EE> explicit Row(const EE* p)
00322       { assert(p); for(int i=0;i<N;++i) d[i*STRIDE]=p[i]; }
00323     template <class EE> Row& operator=(const EE* p)
00324       { assert(p); for(int i=0;i<N;++i) d[i*STRIDE]=p[i]; return *this; }
00325 
00326     // Conforming assignment ops.
00327     template <class EE, int SS> Row& operator=(const Row<N,EE,SS>& vv) {
00328         Impl::copy(*this, vv);
00329         return *this;
00330     }
00331     template <class EE, int SS> Row& operator+=(const Row<N,EE,SS>& r)
00332       { for(int i=0;i<N;++i) d[i*STRIDE] += r[i]; return *this; }
00333     template <class EE, int SS> Row& operator+=(const Row<N,negator<EE>,SS>& r)
00334       { for(int i=0;i<N;++i) d[i*STRIDE] -= -(r[i]); return *this; }
00335     template <class EE, int SS> Row& operator-=(const Row<N,EE,SS>& r)
00336       { for(int i=0;i<N;++i) d[i*STRIDE] -= r[i]; return *this; }
00337     template <class EE, int SS> Row& operator-=(const Row<N,negator<EE>,SS>& r)
00338       { for(int i=0;i<N;++i) d[i*STRIDE] += -(r[i]); return *this; }
00339 
00340     // Conforming binary ops with 'this' on left, producing new packed result.
00341     // Cases: r=r+r, r=r-r, s=r*v r=r*m
00342     template <class EE, int SS> Row<N,typename CNT<E>::template Result<EE>::Add>
00343     conformingAdd(const Row<N,EE,SS>& r) const {
00344         Row<N,typename CNT<E>::template Result<EE>::Add> result;
00345         Impl::conformingAdd(*this, r, result);
00346         return result;
00347     }
00348     template <class EE, int SS> Row<N,typename CNT<E>::template Result<EE>::Sub>
00349     conformingSubtract(const Row<N,EE,SS>& r) const {
00350         Row<N,typename CNT<E>::template Result<EE>::Sub> result;
00351         Impl::conformingSubtract(*this, r, result);
00352         return result;
00353     }
00354 
00355     // dot product
00356     template <class EE, int SS> typename CNT<E>::template Result<EE>::Mul
00357     conformingMultiply(const Vec<N,EE,SS>& r) const {
00358         return (*this)*r;
00359     }
00360 
00361     // row=row*mat
00362     template <int MatNCol, class EE, int CS, int RS> 
00363     Row<MatNCol,typename CNT<E>::template Result<EE>::Mul>
00364     conformingMultiply(const Mat<N,MatNCol,EE,CS,RS>& m) const {
00365         Row<MatNCol,typename CNT<E>::template Result<EE>::Mul> result;
00366         for (int j=0;j<N;++j) result[j] = conformingMultiply(m(j));
00367         return result;
00368     }
00369 
00370     const E& operator[](int i) const { assert(0 <= i && i < N); return d[i*STRIDE]; }
00371     E&       operator[](int i)       { assert(0 <= i && i < N); return d[i*STRIDE]; }
00372     const E& operator()(int i) const { return (*this)[i]; }
00373     E&       operator()(int i)       { return (*this)[i]; }
00374 
00375     ScalarNormSq normSqr() const { return scalarNormSqr(); }
00376     typename CNT<ScalarNormSq>::TSqrt 
00377         norm() const { return CNT<ScalarNormSq>::sqrt(scalarNormSqr()); }
00378 
00379     // If the elements of this Row are scalars, the result is what you get by
00380     // dividing each element by the norm() calculated above. If the elements are
00381     // *not* scalars, then the elements are *separately* normalized. That means
00382     // you will get a different answer from Row<2,Row3>::normalize() than you
00383     // would from a Row<6>::normalize() containing the same scalars.
00384     //
00385     // Normalize returns a row of the same dimension but in new, packed storage
00386     // and with a return type that does not include negator<> even if the original
00387     // Row<> does, because we can eliminate the negation here almost for free.
00388     // But we can't standardize (change conjugate to complex) for free, so we'll retain
00389     // conjugates if there are any.
00390     TNormalize normalize() const {
00391         if (CNT<E>::IsScalar) {
00392             return castAwayNegatorIfAny() / (SignInterpretation*norm());
00393         } else {
00394             TNormalize elementwiseNormalized;
00395             for (int j=0; j<N; ++j) 
00396                 elementwiseNormalized[j] = CNT<E>::normalize((*this)[j]);
00397             return elementwiseNormalized;
00398         }
00399     }
00400 
00401     TInvert invert() const {assert(false); return TInvert();} // TODO default inversion
00402 
00403     const Row&   operator+() const { return *this; }
00404     const TNeg&  operator-() const { return negate(); }
00405     TNeg&        operator-()       { return updNegate(); }
00406     const THerm& operator~() const { return transpose(); }
00407     THerm&       operator~()       { return updTranspose(); }
00408 
00409     const TNeg&  negate() const { return *reinterpret_cast<const TNeg*>(this); }
00410     TNeg&        updNegate()    { return *reinterpret_cast<TNeg*>(this); }
00411 
00412     const THerm& transpose()    const { return *reinterpret_cast<const THerm*>(this); }
00413     THerm&       updTranspose()       { return *reinterpret_cast<THerm*>(this); }
00414 
00415     const TPosTrans& positionalTranspose() const
00416         { return *reinterpret_cast<const TPosTrans*>(this); }
00417     TPosTrans&       updPositionalTranspose()
00418         { return *reinterpret_cast<TPosTrans*>(this); }
00419 
00420     const TReal& real() const { return *reinterpret_cast<const TReal*>(this); }
00421     TReal&       real()       { return *reinterpret_cast<      TReal*>(this); }
00422 
00423     // Had to contort these routines to get them through VC++ 7.net
00424     const TImag& imag()    const { 
00425         const int offs = ImagOffset;
00426         const EImag* p = reinterpret_cast<const EImag*>(this);
00427         return *reinterpret_cast<const TImag*>(p+offs);
00428     }
00429     TImag& imag() { 
00430         const int offs = ImagOffset;
00431         EImag* p = reinterpret_cast<EImag*>(this);
00432         return *reinterpret_cast<TImag*>(p+offs);
00433     }
00434 
00435     const TWithoutNegator& castAwayNegatorIfAny() const {return *reinterpret_cast<const TWithoutNegator*>(this);}
00436     TWithoutNegator&       updCastAwayNegatorIfAny()    {return *reinterpret_cast<TWithoutNegator*>(this);}
00437 
00438 
00439     // These are elementwise binary operators, (this op ee) by default but (ee op this) if
00440     // 'FromLeft' appears in the name. The result is a packed Row<N> but the element type
00441     // may change. These are mostly used to implement global operators.
00442 
00443     //TODO: consider converting 'e' to Standard Numbers as precalculation and changing
00444     // return type appropriately.
00445     template <class EE> Row<N, typename CNT<E>::template Result<EE>::Mul>
00446     scalarMultiply(const EE& e) const {
00447         Row<N, typename CNT<E>::template Result<EE>::Mul> result;
00448         for (int j=0; j<N; ++j) result[j] = (*this)[j] * e;
00449         return result;
00450     }
00451     template <class EE> Row<N, typename CNT<EE>::template Result<E>::Mul>
00452     scalarMultiplyFromLeft(const EE& e) const {
00453         Row<N, typename CNT<EE>::template Result<E>::Mul> result;
00454         for (int j=0; j<N; ++j) result[j] = e * (*this)[j];
00455         return result;
00456     }
00457 
00458     // TODO: should precalculate and store 1/e, while converting to Standard Numbers. Note
00459     // that return type should change appropriately.
00460     template <class EE> Row<N, typename CNT<E>::template Result<EE>::Dvd>
00461     scalarDivide(const EE& e) const {
00462         Row<N, typename CNT<E>::template Result<EE>::Dvd> result;
00463         for (int j=0; j<N; ++j) result[j] = (*this)[j] / e;
00464         return result;
00465     }
00466     template <class EE> Row<N, typename CNT<EE>::template Result<E>::Dvd>
00467     scalarDivideFromLeft(const EE& e) const {
00468         Row<N, typename CNT<EE>::template Result<E>::Dvd> result;
00469         for (int j=0; j<N; ++j) result[j] = e / (*this)[j];
00470         return result;
00471     }
00472 
00473     template <class EE> Row<N, typename CNT<E>::template Result<EE>::Add>
00474     scalarAdd(const EE& e) const {
00475         Row<N, typename CNT<E>::template Result<EE>::Add> result;
00476         for (int j=0; j<N; ++j) result[j] = (*this)[j] + e;
00477         return result;
00478     }
00479     // Add is commutative, so no 'FromLeft'.
00480 
00481     template <class EE> Row<N, typename CNT<E>::template Result<EE>::Sub>
00482     scalarSubtract(const EE& e) const {
00483         Row<N, typename CNT<E>::template Result<EE>::Sub> result;
00484         for (int j=0; j<N; ++j) result[j] = (*this)[j] - e;
00485         return result;
00486     }
00487     template <class EE> Row<N, typename CNT<EE>::template Result<E>::Sub>
00488     scalarSubtractFromLeft(const EE& e) const {
00489         Row<N, typename CNT<EE>::template Result<E>::Sub> result;
00490         for (int j=0; j<N; ++j) result[j] = e - (*this)[j];
00491         return result;
00492     }
00493 
00494     // Generic assignments for any element type not listed explicitly, including scalars.
00495     // These are done repeatedly for each element and only work if the operation can
00496     // be performed leaving the original element type.
00497     template <class EE> Row& operator =(const EE& e) {return scalarEq(e);}
00498     template <class EE> Row& operator+=(const EE& e) {return scalarPlusEq(e);}
00499     template <class EE> Row& operator-=(const EE& e) {return scalarMinusEq(e);}
00500     template <class EE> Row& operator*=(const EE& e) {return scalarTimesEq(e);}
00501     template <class EE> Row& operator/=(const EE& e) {return scalarDivideEq(e);}
00502 
00503     // Generalized scalar assignment & computed assignment methods. These will work
00504     // for any assignment-compatible element, not just scalars.
00505     template <class EE> Row& scalarEq(const EE& ee)
00506       { for(int i=0;i<N;++i) d[i*STRIDE] = ee; return *this; }
00507     template <class EE> Row& scalarPlusEq(const EE& ee)
00508       { for(int i=0;i<N;++i) d[i*STRIDE] += ee; return *this; }
00509     template <class EE> Row& scalarMinusEq(const EE& ee)
00510       { for(int i=0;i<N;++i) d[i*STRIDE] -= ee; return *this; }
00511     template <class EE> Row& scalarInverseMinusEq(const EE& ee)
00512       { for(int i=0;i<N;++i) d[i*STRIDE] = ee - d[i*STRIDE]; return *this; }
00513     template <class EE> Row& scalarTimesEq(const EE& ee)
00514       { for(int i=0;i<N;++i) d[i*STRIDE] *= ee; return *this; }
00515     template <class EE> Row& scalarDivideEq(const EE& ee)
00516       { for(int i=0;i<N;++i) d[i*STRIDE] /= ee; return *this; }
00517     template <class EE> Row& scalarInverseDivideEq(const EE& ee)
00518       { for(int i=0;i<N;++i) d[i*STRIDE] = ee / d[i*STRIDE]; return *this; }
00519 
00520     void setToNaN() {
00521         (*this) = CNT<ELT>::getNaN();
00522     }
00523 
00524     void setToZero() {
00525         (*this) = ELT(0);
00526     }
00527 
00528     // Extract a sub-Row with size known at compile time. These have to be
00529     // called with explicit template arguments, e.g. getSubRow<3>(j).
00530     template <int NN>
00531     const Row<NN,ELT,STRIDE>& getSubRow(int j) const {
00532         assert(0 <= j && j + NN <= N);
00533         return Row<NN,ELT,STRIDE>::getAs(&(*this)[j]);
00534     }
00535     template <int NN>
00536     Row<NN,ELT,STRIDE>& updSubRow(int j) {
00537         assert(0 <= j && j + NN <= N);
00538         return Row<NN,ELT,STRIDE>::updAs(&(*this)[j]);
00539     }
00540 
00541     // Return a row one smaller than this one by dropping the element
00542     // at the indicated position p. The result is packed but has same
00543     // element type as this one.
00544     Row<N-1,ELT,1> drop1(int p) const {
00545         assert(0 <= p && p < N);
00546         Row<N-1,ELT,1> out;
00547         int nxt=0;
00548         for (int i=0; i<N-1; ++i, ++nxt) {
00549             if (nxt==p) ++nxt;  // skip the loser
00550             out[i] = (*this)[nxt];
00551         }
00552         return out;
00553     }
00554 
00555     // Return a vector one larger than this one by adding an element
00556     // to the end. The result is packed but has same element type as
00557     // this one. Works for any assignment compatible element.
00558     template <class EE> Row<N+1,ELT,1> append1(const EE& v) const {
00559         Row<N+1,ELT,1> out;
00560         Row<N,ELT,1>::updAs(&out[0]) = (*this);
00561         out[N] = v;
00562         return out;
00563     }
00564 
00565 
00566     // Return a vector one larger than this one by inserting an element
00567     // *before* the indicated one. The result is packed but has same element type as
00568     // this one. Works for any assignment compatible element. The index
00569     // can be one greater than normally allowed in which case the element
00570     // is appended.
00571     template <class EE> Row<N+1,ELT,1> insert1(int p, const EE& v) const {
00572         assert(0 <= p && p <= N);
00573         if (p==N) return append1(v);
00574         Row<N+1,ELT,1> out;
00575         int nxt=0;
00576         for (int i=0; i<N; ++i, ++nxt) {
00577             if (i==p) out[nxt++] = v;
00578             out[nxt] = (*this)[i];
00579         }
00580         return out;
00581     }
00582 
00583     // These assume we are given a pointer to d[0] of a Row<N,E,S> like this one.
00584     static const Row& getAs(const ELT* p)  {return *reinterpret_cast<const Row*>(p);}
00585     static Row&       updAs(ELT* p)        {return *reinterpret_cast<Row*>(p);}
00586 
00587     // Extract a subrow from a longer one. Element type and stride must match.
00588     template <int NN>
00589     static const Row& getSubRow(const Row<NN,ELT,STRIDE>& r, int j) {
00590         assert(0 <= j && j + N <= NN);
00591         return getAs(&r[j]);
00592     }
00593     template <int NN>
00594     static Row& updSubRow(Row<NN,ELT,STRIDE>& r, int j) {
00595         assert(0 <= j && j + N <= NN);
00596         return updAs(&r[j]);
00597     }
00598 
00599     static Row<N,ELT,1> getNaN() { return Row<N,ELT,1>(CNT<ELT>::getNaN()); }
00600 
00602     bool isNaN() const {
00603         for (int j=0; j<N; ++j)
00604             if (CNT<ELT>::isNaN((*this)[j]))
00605                 return true;
00606         return false;
00607     }
00608 
00611     bool isInf() const {
00612         bool seenInf = false;
00613         for (int j=0; j<N; ++j) {
00614             const ELT& e = (*this)[j];
00615             if (!CNT<ELT>::isFinite(e)) {
00616                 if (!CNT<ELT>::isInf(e)) 
00617                     return false; // something bad was found
00618                 seenInf = true; 
00619             }
00620         }
00621         return seenInf;
00622     }
00623 
00625     bool isFinite() const {
00626         for (int j=0; j<N; ++j)
00627             if (!CNT<ELT>::isFinite((*this)[j]))
00628                 return false;
00629         return true;
00630     }
00631 
00634     static double getDefaultTolerance() {return CNT<ELT>::getDefaultTolerance();}
00635 
00638     template <class E2, int CS2>
00639     bool isNumericallyEqual(const Row<N,E2,CS2>& r, double tol) const {
00640         for (int j=0; j<N; ++j)
00641             if (!CNT<ELT>::isNumericallyEqual((*this)(j), r(j), tol))
00642                 return false;
00643         return true;
00644     }
00645 
00649     template <class E2, int CS2>
00650     bool isNumericallyEqual(const Row<N,E2,CS2>& r) const {
00651         const double tol = std::max(getDefaultTolerance(),r.getDefaultTolerance());
00652         return isNumericallyEqual(r, tol);
00653     }
00654 
00659     bool isNumericallyEqual
00660        (const ELT& e,
00661         double     tol = getDefaultTolerance()) const 
00662     {
00663         for (int j=0; j<N; ++j)
00664             if (!CNT<ELT>::isNumericallyEqual((*this)(j), e, tol))
00665                 return false;
00666         return true;
00667     }
00668 private:
00669     ELT d[NActualElements];    // data
00670 };
00671 
00673 // Global operators involving two rows.    //
00674 //   v+v, v-v, v==v, v!=v                  //
00676 
00677 // v3 = v1 + v2 where all v's have the same length N. 
00678 template <int N, class E1, int S1, class E2, int S2> inline
00679 typename Row<N,E1,S1>::template Result< Row<N,E2,S2> >::Add
00680 operator+(const Row<N,E1,S1>& l, const Row<N,E2,S2>& r) { 
00681     return Row<N,E1,S1>::template Result< Row<N,E2,S2> >
00682         ::AddOp::perform(l,r);
00683 }
00684 
00685 // v3 = v1 - v2, similar to +
00686 template <int N, class E1, int S1, class E2, int S2> inline
00687 typename Row<N,E1,S1>::template Result< Row<N,E2,S2> >::Sub
00688 operator-(const Row<N,E1,S1>& l, const Row<N,E2,S2>& r) { 
00689     return Row<N,E1,S1>::template Result< Row<N,E2,S2> >
00690         ::SubOp::perform(l,r);
00691 }
00692 
00693 // bool = v1 == v2, v1 and v2 have the same length M
00694 template <int N, class E1, int S1, class E2, int S2> inline bool
00695 operator==(const Row<N,E1,S1>& l, const Row<N,E2,S2>& r) { 
00696     for (int i=0; i < N; ++i)
00697         if (l[i] != r[i]) return false;
00698     return true;
00699 }
00700 
00701 // bool = v1 != v2, v1 and v2 have the same length M
00702 template <int N, class E1, int S1, class E2, int S2> inline bool
00703 operator!=(const Row<N,E1,S1>& l, const Row<N,E2,S2>& r) {return !(l==r);} 
00704 
00705 
00707 // Global operators involving a row and a scalar. //
00709 
00710 // I haven't been able to figure out a nice way to templatize for the
00711 // built-in reals without introducing a lot of unwanted type matches
00712 // as well. So we'll just grind them out explicitly here.
00713 
00714 // SCALAR MULTIPLY
00715 
00716 // v = v*real, real*v 
00717 template <int N, class E, int S> inline
00718 typename Row<N,E,S>::template Result<float>::Mul
00719 operator*(const Row<N,E,S>& l, const float& r)
00720   { return Row<N,E,S>::template Result<float>::MulOp::perform(l,r); }
00721 template <int N, class E, int S> inline
00722 typename Row<N,E,S>::template Result<float>::Mul
00723 operator*(const float& l, const Row<N,E,S>& r) {return r*l;}
00724 
00725 template <int N, class E, int S> inline
00726 typename Row<N,E,S>::template Result<double>::Mul
00727 operator*(const Row<N,E,S>& l, const double& r)
00728   { return Row<N,E,S>::template Result<double>::MulOp::perform(l,r); }
00729 template <int N, class E, int S> inline
00730 typename Row<N,E,S>::template Result<double>::Mul
00731 operator*(const double& l, const Row<N,E,S>& r) {return r*l;}
00732 
00733 template <int N, class E, int S> inline
00734 typename Row<N,E,S>::template Result<long double>::Mul
00735 operator*(const Row<N,E,S>& l, const long double& r)
00736   { return Row<N,E,S>::template Result<long double>::MulOp::perform(l,r); }
00737 template <int N, class E, int S> inline
00738 typename Row<N,E,S>::template Result<long double>::Mul
00739 operator*(const long double& l, const Row<N,E,S>& r) {return r*l;}
00740 
00741 // v = v*int, int*v -- just convert int to v's precision float
00742 template <int N, class E, int S> inline
00743 typename Row<N,E,S>::template Result<typename CNT<E>::Precision>::Mul
00744 operator*(const Row<N,E,S>& l, int r) {return l * (typename CNT<E>::Precision)r;}
00745 template <int N, class E, int S> inline
00746 typename Row<N,E,S>::template Result<typename CNT<E>::Precision>::Mul
00747 operator*(int l, const Row<N,E,S>& r) {return r * (typename CNT<E>::Precision)l;}
00748 
00749 // Complex, conjugate, and negator are all easy to templatize.
00750 
00751 // v = v*complex, complex*v
00752 template <int N, class E, int S, class R> inline
00753 typename Row<N,E,S>::template Result<std::complex<R> >::Mul
00754 operator*(const Row<N,E,S>& l, const std::complex<R>& r)
00755   { return Row<N,E,S>::template Result<std::complex<R> >::MulOp::perform(l,r); }
00756 template <int N, class E, int S, class R> inline
00757 typename Row<N,E,S>::template Result<std::complex<R> >::Mul
00758 operator*(const std::complex<R>& l, const Row<N,E,S>& r) {return r*l;}
00759 
00760 // v = v*conjugate, conjugate*v (convert conjugate->complex)
00761 template <int N, class E, int S, class R> inline
00762 typename Row<N,E,S>::template Result<std::complex<R> >::Mul
00763 operator*(const Row<N,E,S>& l, const conjugate<R>& r) {return l*(std::complex<R>)r;}
00764 template <int N, class E, int S, class R> inline
00765 typename Row<N,E,S>::template Result<std::complex<R> >::Mul
00766 operator*(const conjugate<R>& l, const Row<N,E,S>& r) {return r*(std::complex<R>)l;}
00767 
00768 // v = v*negator, negator*v: convert negator to standard number
00769 template <int N, class E, int S, class R> inline
00770 typename Row<N,E,S>::template Result<typename negator<R>::StdNumber>::Mul
00771 operator*(const Row<N,E,S>& l, const negator<R>& r) {return l * (typename negator<R>::StdNumber)(R)r;}
00772 template <int N, class E, int S, class R> inline
00773 typename Row<N,E,S>::template Result<typename negator<R>::StdNumber>::Mul
00774 operator*(const negator<R>& l, const Row<N,E,S>& r) {return r * (typename negator<R>::StdNumber)(R)l;}
00775 
00776 
00777 // SCALAR DIVIDE. This is a scalar operation when the scalar is on the right,
00778 // but when it is on the left it means scalar * pseudoInverse(row), which is a vec.
00779 
00780 // v = v/real, real/v 
00781 template <int N, class E, int S> inline
00782 typename Row<N,E,S>::template Result<float>::Dvd
00783 operator/(const Row<N,E,S>& l, const float& r)
00784   { return Row<N,E,S>::template Result<float>::DvdOp::perform(l,r); }
00785 template <int N, class E, int S> inline
00786 typename CNT<float>::template Result<Row<N,E,S> >::Dvd
00787 operator/(const float& l, const Row<N,E,S>& r)
00788   { return CNT<float>::template Result<Row<N,E,S> >::DvdOp::perform(l,r); }
00789 
00790 template <int N, class E, int S> inline
00791 typename Row<N,E,S>::template Result<double>::Dvd
00792 operator/(const Row<N,E,S>& l, const double& r)
00793   { return Row<N,E,S>::template Result<double>::DvdOp::perform(l,r); }
00794 template <int N, class E, int S> inline
00795 typename CNT<double>::template Result<Row<N,E,S> >::Dvd
00796 operator/(const double& l, const Row<N,E,S>& r)
00797   { return CNT<double>::template Result<Row<N,E,S> >::DvdOp::perform(l,r); }
00798 
00799 template <int N, class E, int S> inline
00800 typename Row<N,E,S>::template Result<long double>::Dvd
00801 operator/(const Row<N,E,S>& l, const long double& r)
00802   { return Row<N,E,S>::template Result<long double>::DvdOp::perform(l,r); }
00803 template <int N, class E, int S> inline
00804 typename CNT<long double>::template Result<Row<N,E,S> >::Dvd
00805 operator/(const long double& l, const Row<N,E,S>& r)
00806   { return CNT<long double>::template Result<Row<N,E,S> >::DvdOp::perform(l,r); }
00807 
00808 // v = v/int, int/v -- just convert int to v's precision float
00809 template <int N, class E, int S> inline
00810 typename Row<N,E,S>::template Result<typename CNT<E>::Precision>::Dvd
00811 operator/(const Row<N,E,S>& l, int r) {return l / (typename CNT<E>::Precision)r;}
00812 template <int N, class E, int S> inline
00813 typename CNT<typename CNT<E>::Precision>::template Result<Row<N,E,S> >::Dvd
00814 operator/(int l, const Row<N,E,S>& r) {return (typename CNT<E>::Precision)l / r;}
00815 
00816 
00817 // Complex, conjugate, and negator are all easy to templatize.
00818 
00819 // v = v/complex, complex/v
00820 template <int N, class E, int S, class R> inline
00821 typename Row<N,E,S>::template Result<std::complex<R> >::Dvd
00822 operator/(const Row<N,E,S>& l, const std::complex<R>& r)
00823   { return Row<N,E,S>::template Result<std::complex<R> >::DvdOp::perform(l,r); }
00824 template <int N, class E, int S, class R> inline
00825 typename CNT<std::complex<R> >::template Result<Row<N,E,S> >::Dvd
00826 operator/(const std::complex<R>& l, const Row<N,E,S>& r)
00827   { return CNT<std::complex<R> >::template Result<Row<N,E,S> >::DvdOp::perform(l,r); }
00828 
00829 // v = v/conjugate, conjugate/v (convert conjugate->complex)
00830 template <int N, class E, int S, class R> inline
00831 typename Row<N,E,S>::template Result<std::complex<R> >::Dvd
00832 operator/(const Row<N,E,S>& l, const conjugate<R>& r) {return l/(std::complex<R>)r;}
00833 template <int N, class E, int S, class R> inline
00834 typename CNT<std::complex<R> >::template Result<Row<N,E,S> >::Dvd
00835 operator/(const conjugate<R>& l, const Row<N,E,S>& r) {return (std::complex<R>)l/r;}
00836 
00837 // v = v/negator, negator/v: convert negator to number
00838 template <int N, class E, int S, class R> inline
00839 typename Row<N,E,S>::template Result<typename negator<R>::StdNumber>::Dvd
00840 operator/(const Row<N,E,S>& l, const negator<R>& r) {return l/(typename negator<R>::StdNumber)(R)r;}
00841 template <int N, class E, int S, class R> inline
00842 typename CNT<R>::template Result<Row<N,E,S> >::Dvd
00843 operator/(const negator<R>& l, const Row<N,E,S>& r) {return (typename negator<R>::StdNumber)(R)l/r;}
00844 
00845 
00846 // Add and subtract are odd as scalar ops. They behave as though the
00847 // scalar stands for a vector each of whose elements is that scalar,
00848 // and then a normal vector add or subtract is done.
00849 
00850 // SCALAR ADD
00851 
00852 // v = v+real, real+v 
00853 template <int N, class E, int S> inline
00854 typename Row<N,E,S>::template Result<float>::Add
00855 operator+(const Row<N,E,S>& l, const float& r)
00856   { return Row<N,E,S>::template Result<float>::AddOp::perform(l,r); }
00857 template <int N, class E, int S> inline
00858 typename Row<N,E,S>::template Result<float>::Add
00859 operator+(const float& l, const Row<N,E,S>& r) {return r+l;}
00860 
00861 template <int N, class E, int S> inline
00862 typename Row<N,E,S>::template Result<double>::Add
00863 operator+(const Row<N,E,S>& l, const double& r)
00864   { return Row<N,E,S>::template Result<double>::AddOp::perform(l,r); }
00865 template <int N, class E, int S> inline
00866 typename Row<N,E,S>::template Result<double>::Add
00867 operator+(const double& l, const Row<N,E,S>& r) {return r+l;}
00868 
00869 template <int N, class E, int S> inline
00870 typename Row<N,E,S>::template Result<long double>::Add
00871 operator+(const Row<N,E,S>& l, const long double& r)
00872   { return Row<N,E,S>::template Result<long double>::AddOp::perform(l,r); }
00873 template <int N, class E, int S> inline
00874 typename Row<N,E,S>::template Result<long double>::Add
00875 operator+(const long double& l, const Row<N,E,S>& r) {return r+l;}
00876 
00877 // v = v+int, int+v -- just convert int to v's precision float
00878 template <int N, class E, int S> inline
00879 typename Row<N,E,S>::template Result<typename CNT<E>::Precision>::Add
00880 operator+(const Row<N,E,S>& l, int r) {return l + (typename CNT<E>::Precision)r;}
00881 template <int N, class E, int S> inline
00882 typename Row<N,E,S>::template Result<typename CNT<E>::Precision>::Add
00883 operator+(int l, const Row<N,E,S>& r) {return r + (typename CNT<E>::Precision)l;}
00884 
00885 // Complex, conjugate, and negator are all easy to templatize.
00886 
00887 // v = v+complex, complex+v
00888 template <int N, class E, int S, class R> inline
00889 typename Row<N,E,S>::template Result<std::complex<R> >::Add
00890 operator+(const Row<N,E,S>& l, const std::complex<R>& r)
00891   { return Row<N,E,S>::template Result<std::complex<R> >::AddOp::perform(l,r); }
00892 template <int N, class E, int S, class R> inline
00893 typename Row<N,E,S>::template Result<std::complex<R> >::Add
00894 operator+(const std::complex<R>& l, const Row<N,E,S>& r) {return r+l;}
00895 
00896 // v = v+conjugate, conjugate+v (convert conjugate->complex)
00897 template <int N, class E, int S, class R> inline
00898 typename Row<N,E,S>::template Result<std::complex<R> >::Add
00899 operator+(const Row<N,E,S>& l, const conjugate<R>& r) {return l+(std::complex<R>)r;}
00900 template <int N, class E, int S, class R> inline
00901 typename Row<N,E,S>::template Result<std::complex<R> >::Add
00902 operator+(const conjugate<R>& l, const Row<N,E,S>& r) {return r+(std::complex<R>)l;}
00903 
00904 // v = v+negator, negator+v: convert negator to standard number
00905 template <int N, class E, int S, class R> inline
00906 typename Row<N,E,S>::template Result<typename negator<R>::StdNumber>::Add
00907 operator+(const Row<N,E,S>& l, const negator<R>& r) {return l + (typename negator<R>::StdNumber)(R)r;}
00908 template <int N, class E, int S, class R> inline
00909 typename Row<N,E,S>::template Result<typename negator<R>::StdNumber>::Add
00910 operator+(const negator<R>& l, const Row<N,E,S>& r) {return r + (typename negator<R>::StdNumber)(R)l;}
00911 
00912 // SCALAR SUBTRACT -- careful, not commutative.
00913 
00914 // v = v-real, real-v 
00915 template <int N, class E, int S> inline
00916 typename Row<N,E,S>::template Result<float>::Sub
00917 operator-(const Row<N,E,S>& l, const float& r)
00918   { return Row<N,E,S>::template Result<float>::SubOp::perform(l,r); }
00919 template <int N, class E, int S> inline
00920 typename CNT<float>::template Result<Row<N,E,S> >::Sub
00921 operator-(const float& l, const Row<N,E,S>& r)
00922   { return CNT<float>::template Result<Row<N,E,S> >::SubOp::perform(l,r); }
00923 
00924 template <int N, class E, int S> inline
00925 typename Row<N,E,S>::template Result<double>::Sub
00926 operator-(const Row<N,E,S>& l, const double& r)
00927   { return Row<N,E,S>::template Result<double>::SubOp::perform(l,r); }
00928 template <int N, class E, int S> inline
00929 typename CNT<double>::template Result<Row<N,E,S> >::Sub
00930 operator-(const double& l, const Row<N,E,S>& r)
00931   { return CNT<double>::template Result<Row<N,E,S> >::SubOp::perform(l,r); }
00932 
00933 template <int N, class E, int S> inline
00934 typename Row<N,E,S>::template Result<long double>::Sub
00935 operator-(const Row<N,E,S>& l, const long double& r)
00936   { return Row<N,E,S>::template Result<long double>::SubOp::perform(l,r); }
00937 template <int N, class E, int S> inline
00938 typename CNT<long double>::template Result<Row<N,E,S> >::Sub
00939 operator-(const long double& l, const Row<N,E,S>& r)
00940   { return CNT<long double>::template Result<Row<N,E,S> >::SubOp::perform(l,r); }
00941 
00942 // v = v-int, int-v // just convert int to v's precision float
00943 template <int N, class E, int S> inline
00944 typename Row<N,E,S>::template Result<typename CNT<E>::Precision>::Sub
00945 operator-(const Row<N,E,S>& l, int r) {return l - (typename CNT<E>::Precision)r;}
00946 template <int N, class E, int S> inline
00947 typename CNT<typename CNT<E>::Precision>::template Result<Row<N,E,S> >::Sub
00948 operator-(int l, const Row<N,E,S>& r) {return (typename CNT<E>::Precision)l - r;}
00949 
00950 
00951 // Complex, conjugate, and negator are all easy to templatize.
00952 
00953 // v = v-complex, complex-v
00954 template <int N, class E, int S, class R> inline
00955 typename Row<N,E,S>::template Result<std::complex<R> >::Sub
00956 operator-(const Row<N,E,S>& l, const std::complex<R>& r)
00957   { return Row<N,E,S>::template Result<std::complex<R> >::SubOp::perform(l,r); }
00958 template <int N, class E, int S, class R> inline
00959 typename CNT<std::complex<R> >::template Result<Row<N,E,S> >::Sub
00960 operator-(const std::complex<R>& l, const Row<N,E,S>& r)
00961   { return CNT<std::complex<R> >::template Result<Row<N,E,S> >::SubOp::perform(l,r); }
00962 
00963 // v = v-conjugate, conjugate-v (convert conjugate->complex)
00964 template <int N, class E, int S, class R> inline
00965 typename Row<N,E,S>::template Result<std::complex<R> >::Sub
00966 operator-(const Row<N,E,S>& l, const conjugate<R>& r) {return l-(std::complex<R>)r;}
00967 template <int N, class E, int S, class R> inline
00968 typename CNT<std::complex<R> >::template Result<Row<N,E,S> >::Sub
00969 operator-(const conjugate<R>& l, const Row<N,E,S>& r) {return (std::complex<R>)l-r;}
00970 
00971 // v = v-negator, negator-v: convert negator to standard number
00972 template <int N, class E, int S, class R> inline
00973 typename Row<N,E,S>::template Result<typename negator<R>::StdNumber>::Sub
00974 operator-(const Row<N,E,S>& l, const negator<R>& r) {return l-(typename negator<R>::StdNumber)(R)r;}
00975 template <int N, class E, int S, class R> inline
00976 typename CNT<R>::template Result<Row<N,E,S> >::Sub
00977 operator-(const negator<R>& l, const Row<N,E,S>& r) {return (typename negator<R>::StdNumber)(R)l-r;}
00978 
00979 
00980 // Row I/O
00981 template <int N, class E, int S, class CHAR, class TRAITS> inline
00982 std::basic_ostream<CHAR,TRAITS>&
00983 operator<<(std::basic_ostream<CHAR,TRAITS>& o, const Row<N,E,S>& v) {
00984     o << "[" << v[0]; for(int i=1;i<N;++i) o<<','<<v[i]; o<<']'; return o;
00985 }
00986 
00987 template <int N, class E, int S, class CHAR, class TRAITS> inline
00988 std::basic_istream<CHAR,TRAITS>&
00989 operator>>(std::basic_istream<CHAR,TRAITS>& is, Row<N,E,S>& v) {
00990     // TODO: not sure how to do Row input yet
00991     assert(false);
00992     return is;
00993 }
00994 
00995 } //namespace SimTK
00996 
00997 
00998 #endif //SimTK_SIMMATRIX_SMALLMATRIX_ROW_H_

Generated by  doxygen 1.6.2