Simbody  3.5
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  *                       Simbody(tm): SimTKcommon                             *
00006  * -------------------------------------------------------------------------- *
00007  * This is part of the SimTK biosimulation toolkit originating from           *
00008  * Simbios, the NIH National Center for Physics-Based Simulation of           *
00009  * Biological Structures at Stanford, funded under the NIH Roadmap for        *
00010  * Medical Research, grant U54 GM072970. See https://simtk.org/home/simbody.  *
00011  *                                                                            *
00012  * Portions copyright (c) 2005-12 Stanford University and the Authors.        *
00013  * Authors: Michael Sherman                                                   *
00014  * Contributors: Peter Eastman                                                *
00015  *                                                                            *
00016  * Licensed under the Apache License, Version 2.0 (the "License"); you may    *
00017  * not use this file except in compliance with the License. You may obtain a  *
00018  * copy of the License at http://www.apache.org/licenses/LICENSE-2.0.         *
00019  *                                                                            *
00020  * Unless required by applicable law or agreed to in writing, software        *
00021  * distributed under the License is distributed on an "AS IS" BASIS,          *
00022  * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.   *
00023  * See the License for the specific language governing permissions and        *
00024  * limitations under the License.                                             *
00025  * -------------------------------------------------------------------------- */
00026 
00031 #include "SimTKcommon/internal/common.h"
00032 
00033 
00034 namespace SimTK {
00035 
00036 // The following functions are used internally by Row.
00037 
00038 namespace Impl {
00039 
00040 // For those wimpy compilers that don't unroll short, constant-limit loops, 
00041 // Peter Eastman added these recursive template implementations of 
00042 // elementwise add, subtract, and copy. Sherm added multiply and divide.
00043 
00044 template <class E1, int S1, class E2, int S2> void
00045 conformingAdd(const Row<1,E1,S1>& r1, const Row<1,E2,S2>& r2, 
00046               Row<1,typename CNT<E1>::template Result<E2>::Add>& result) {
00047     result[0] = r1[0] + r2[0];
00048 }
00049 template <int N, class E1, int S1, class E2, int S2> void
00050 conformingAdd(const Row<N,E1,S1>& r1, const Row<N,E2,S2>& r2, 
00051               Row<N,typename CNT<E1>::template Result<E2>::Add>& result) {
00052     conformingAdd(reinterpret_cast<const Row<N-1,E1,S1>&>(r1), 
00053                   reinterpret_cast<const Row<N-1,E2,S2>&>(r2), 
00054                   reinterpret_cast<Row<N-1,typename CNT<E1>::
00055                               template Result<E2>::Add>&>(result));
00056     result[N-1] = r1[N-1] + r2[N-1];
00057 }
00058 
00059 template <class E1, int S1, class E2, int S2> void
00060 conformingSubtract(const Row<1,E1,S1>& r1, const Row<1,E2,S2>& r2, 
00061                    Row<1,typename CNT<E1>::template Result<E2>::Sub>& 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,
00066                    Row<N,typename CNT<E1>::template Result<E2>::Sub>& result) {
00067     conformingSubtract(reinterpret_cast<const Row<N-1,E1,S1>&>(r1), 
00068                        reinterpret_cast<const Row<N-1,E2,S2>&>(r2), 
00069                        reinterpret_cast<Row<N-1,typename CNT<E1>::
00070                                    template Result<E2>::Sub>&>(result));
00071     result[N-1] = r1[N-1] - r2[N-1];
00072 }
00073 
00074 template <class E1, int S1, class E2, int S2> void
00075 elementwiseMultiply(const Row<1,E1,S1>& r1, const Row<1,E2,S2>& r2, 
00076               Row<1,typename CNT<E1>::template Result<E2>::Mul>& result) {
00077     result[0] = r1[0] * r2[0];
00078 }
00079 template <int N, class E1, int S1, class E2, int S2> void
00080 elementwiseMultiply(const Row<N,E1,S1>& r1, const Row<N,E2,S2>& r2, 
00081               Row<N,typename CNT<E1>::template Result<E2>::Mul>& result) {
00082     elementwiseMultiply(reinterpret_cast<const Row<N-1,E1,S1>&>(r1), 
00083                         reinterpret_cast<const Row<N-1,E2,S2>&>(r2), 
00084                         reinterpret_cast<Row<N-1,typename CNT<E1>::
00085                                     template Result<E2>::Mul>&>(result));
00086     result[N-1] = r1[N-1] * r2[N-1];
00087 }
00088 
00089 template <class E1, int S1, class E2, int S2> void
00090 elementwiseDivide(const Row<1,E1,S1>& r1, const Row<1,E2,S2>& r2, 
00091               Row<1,typename CNT<E1>::template Result<E2>::Dvd>& result) {
00092     result[0] = r1[0] / r2[0];
00093 }
00094 template <int N, class E1, int S1, class E2, int S2> void
00095 elementwiseDivide(const Row<N,E1,S1>& r1, const Row<N,E2,S2>& r2, 
00096               Row<N,typename CNT<E1>::template Result<E2>::Dvd>& result) {
00097     elementwiseDivide(reinterpret_cast<const Row<N-1,E1,S1>&>(r1), 
00098                         reinterpret_cast<const Row<N-1,E2,S2>&>(r2), 
00099                         reinterpret_cast<Row<N-1,typename CNT<E1>::
00100                                     template Result<E2>::Dvd>&>(result));
00101     result[N-1] = r1[N-1] / r2[N-1];
00102 }
00103 
00104 template <class E1, int S1, class E2, int S2> void
00105 copy(Row<1,E1,S1>& r1, const Row<1,E2,S2>& r2) {
00106     r1[0] = r2[0];
00107 }
00108 template <int N, class E1, int S1, class E2, int S2> void
00109 copy(Row<N,E1,S1>& r1, const Row<N,E2,S2>& r2) {
00110     copy(reinterpret_cast<Row<N-1,E1,S1>&>(r1), 
00111          reinterpret_cast<const Row<N-1,E2,S2>&>(r2));
00112     r1[N-1] = r2[N-1];
00113 }
00114 
00115 }
00116 
00132 template <int N, class ELT, int STRIDE> class Row {
00133     typedef ELT                                 E;
00134     typedef typename CNT<E>::TNeg               ENeg;
00135     typedef typename CNT<E>::TWithoutNegator    EWithoutNegator;
00136     typedef typename CNT<E>::TReal              EReal;
00137     typedef typename CNT<E>::TImag              EImag;
00138     typedef typename CNT<E>::TComplex           EComplex;
00139     typedef typename CNT<E>::THerm              EHerm;
00140     typedef typename CNT<E>::TPosTrans          EPosTrans;
00141     typedef typename CNT<E>::TSqHermT           ESqHermT;
00142     typedef typename CNT<E>::TSqTHerm           ESqTHerm;
00143 
00144     typedef typename CNT<E>::TSqrt              ESqrt;
00145     typedef typename CNT<E>::TAbs               EAbs;
00146     typedef typename CNT<E>::TStandard          EStandard;
00147     typedef typename CNT<E>::TInvert            EInvert;
00148     typedef typename CNT<E>::TNormalize         ENormalize;
00149 
00150     typedef typename CNT<E>::Scalar             EScalar;
00151     typedef typename CNT<E>::ULessScalar        EULessScalar;
00152     typedef typename CNT<E>::Number             ENumber;
00153     typedef typename CNT<E>::StdNumber          EStdNumber;
00154     typedef typename CNT<E>::Precision          EPrecision;
00155     typedef typename CNT<E>::ScalarNormSq       EScalarNormSq;
00156 
00157 public:
00158 
00159     enum {
00160         NRows               = 1,
00161         NCols               = N,
00162         NPackedElements     = N,
00163         NActualElements     = N * STRIDE,   // includes trailing gap
00164         NActualScalars      = CNT<E>::NActualScalars * NActualElements,
00165         RowSpacing          = NActualElements,
00166         ColSpacing          = STRIDE,
00167         ImagOffset          = NTraits<ENumber>::ImagOffset,
00168         RealStrideFactor    = 1, // composite types don't change size when
00169                                  // cast from complex to real or imaginary
00170         ArgDepth            = ((int)CNT<E>::ArgDepth < (int)MAX_RESOLVED_DEPTH 
00171                                 ? CNT<E>::ArgDepth + 1 
00172                                 : MAX_RESOLVED_DEPTH),
00173         IsScalar            = 0,
00174         IsULessScalar       = 0,
00175         IsNumber            = 0,
00176         IsStdNumber         = 0,
00177         IsPrecision         = 0,
00178         SignInterpretation  = CNT<E>::SignInterpretation
00179     };
00180 
00181     typedef Row<N,E,STRIDE>                 T;
00182     typedef Row<N,ENeg,STRIDE>              TNeg;
00183     typedef Row<N,EWithoutNegator,STRIDE>   TWithoutNegator;
00184 
00185     typedef Row<N,EReal,STRIDE*CNT<E>::RealStrideFactor>         
00186                                             TReal;
00187     typedef Row<N,EImag,STRIDE*CNT<E>::RealStrideFactor>         
00188                                             TImag;
00189     typedef Row<N,EComplex,STRIDE>          TComplex;
00190     typedef Vec<N,EHerm,STRIDE>             THerm;
00191     typedef Vec<N,E,STRIDE>                 TPosTrans;
00192     typedef E                               TElement;
00193     typedef Row                             TRow;
00194     typedef E                               TCol;
00195 
00196     // These are the results of calculations, so are returned in new, packed
00197     // memory. Be sure to refer to element types here which are also packed.
00198     typedef Vec<N,ESqrt,1>                  TSqrt;      // Note stride
00199     typedef Row<N,EAbs,1>                   TAbs;       // Note stride
00200     typedef Row<N,EStandard,1>              TStandard;
00201     typedef Vec<N,EInvert,1>                TInvert;    // packed
00202     typedef Row<N,ENormalize,1>             TNormalize;
00203 
00204     typedef SymMat<N,ESqHermT>              TSqHermT;   // result of self outer product
00205     typedef EScalarNormSq                   TSqTHerm;   // result of self dot product
00206 
00207     // These recurse right down to the underlying scalar type no matter how
00208     // deep the elements are.
00209     typedef EScalar                         Scalar;
00210     typedef EULessScalar                    ULessScalar;
00211     typedef ENumber                         Number;
00212     typedef EStdNumber                      StdNumber;
00213     typedef EPrecision                      Precision;
00214     typedef EScalarNormSq                   ScalarNormSq;
00215 
00216     static int size() { return N; }
00217     static int nrow() { return 1; }
00218     static int ncol() { return N; }
00219 
00220 
00221     // Scalar norm square is sum( conjugate squares of all scalars )
00222     ScalarNormSq scalarNormSqr() const { 
00223         ScalarNormSq sum(0);
00224         for(int i=0;i<N;++i) sum += CNT<E>::scalarNormSqr(d[i*STRIDE]);
00225         return sum;
00226     }
00227 
00228     // sqrt() is elementwise square root; that is, the return value has the same
00229     // dimension as this Vec but with each element replaced by whatever it thinks
00230     // its square root is.
00231     TSqrt sqrt() const {
00232         TSqrt rsqrt;
00233         for(int i=0;i<N;++i) rsqrt[i] = CNT<E>::sqrt(d[i*STRIDE]);
00234         return rsqrt;
00235     }
00236 
00237     // abs() is elementwise absolute value; that is, the return value has the same
00238     // dimension as this Row but with each element replaced by whatever it thinks
00239     // its absolute value is.
00240     TAbs abs() const {
00241         TAbs rabs;
00242         for(int i=0;i<N;++i) rabs[i] = CNT<E>::abs(d[i*STRIDE]);
00243         return rabs;
00244     }
00245 
00246     TStandard standardize() const {
00247         TStandard rstd;
00248         for(int i=0;i<N;++i) rstd[i] = CNT<E>::standardize(d[i*STRIDE]);
00249         return rstd;
00250     }
00251 
00252     // Sum just adds up all the elements, getting rid of negators and
00253     // conjugates in the process.
00254     EStandard sum() const {
00255         E sum(0);
00256         for (int i=0;i<N;++i) sum += d[i*STRIDE];
00257         return CNT<E>::standardize(sum);
00258     }
00259 
00260     // This gives the resulting rowvector type when (v[i] op P) is applied to each element of v.
00261     // It is a row of length N, stride 1, and element types which are the regular
00262     // composite result of E op P. Typically P is a scalar type but it doesn't have to be.
00263     template <class P> struct EltResult { 
00264         typedef Row<N, typename CNT<E>::template Result<P>::Mul, 1> Mul;
00265         typedef Row<N, typename CNT<E>::template Result<P>::Dvd, 1> Dvd;
00266         typedef Row<N, typename CNT<E>::template Result<P>::Add, 1> Add;
00267         typedef Row<N, typename CNT<E>::template Result<P>::Sub, 1> Sub;
00268     };
00269 
00270     // This is the composite result for v op P where P is some kind of appropriately shaped
00271     // non-scalar type.
00272     template <class P> struct Result { 
00273         typedef MulCNTs<1,N,ArgDepth,Row,ColSpacing,RowSpacing,
00274             CNT<P>::NRows, CNT<P>::NCols, CNT<P>::ArgDepth,
00275             P, CNT<P>::ColSpacing, CNT<P>::RowSpacing> MulOp;
00276         typedef typename MulOp::Type Mul;
00277 
00278         typedef MulCNTsNonConforming<1,N,ArgDepth,Row,ColSpacing,RowSpacing,
00279             CNT<P>::NRows, CNT<P>::NCols, CNT<P>::ArgDepth,
00280             P, CNT<P>::ColSpacing, CNT<P>::RowSpacing> MulOpNonConforming;
00281         typedef typename MulOpNonConforming::Type MulNon;
00282 
00283 
00284         typedef DvdCNTs<1,N,ArgDepth,Row,ColSpacing,RowSpacing,
00285             CNT<P>::NRows, CNT<P>::NCols, CNT<P>::ArgDepth,
00286             P, CNT<P>::ColSpacing, CNT<P>::RowSpacing> DvdOp;
00287         typedef typename DvdOp::Type Dvd;
00288 
00289         typedef AddCNTs<1,N,ArgDepth,Row,ColSpacing,RowSpacing,
00290             CNT<P>::NRows, CNT<P>::NCols, CNT<P>::ArgDepth,
00291             P, CNT<P>::ColSpacing, CNT<P>::RowSpacing> AddOp;
00292         typedef typename AddOp::Type Add;
00293 
00294         typedef SubCNTs<1,N,ArgDepth,Row,ColSpacing,RowSpacing,
00295             CNT<P>::NRows, CNT<P>::NCols, CNT<P>::ArgDepth,
00296             P, CNT<P>::ColSpacing, CNT<P>::RowSpacing> SubOp;
00297         typedef typename SubOp::Type Sub;
00298     };
00299 
00300     // Shape-preserving element substitution (always packed)
00301     template <class P> struct Substitute {
00302         typedef Row<N,P> Type;
00303     };
00304 
00305     // Default construction initializes to NaN when debugging but
00306     // is left uninitialized otherwise.
00307     Row(){ 
00308     #ifndef NDEBUG
00309         setToNaN();
00310     #endif
00311     }
00312 
00313     // It's important not to use the default copy constructor or copy
00314     // assignment because the compiler doesn't understand that we may
00315     // have noncontiguous storage and will try to copy the whole array.
00316     Row(const Row& src) {
00317         Impl::copy(*this, src);
00318     }
00319     Row& operator=(const Row& src) {    // no harm if src and 'this' are the same
00320         Impl::copy(*this, src);
00321         return *this;
00322     }
00323 
00324     // We want an implicit conversion from a Row of the same length
00325     // and element type but with a different stride.
00326     template <int SS> Row(const Row<N,E,SS>& src) {
00327         Impl::copy(*this, src);
00328     }
00329 
00330     // We want an implicit conversion from a Row of the same length
00331     // and *negated* element type, possibly with a different stride.
00332     template <int SS> Row(const Row<N,ENeg,SS>& src) {
00333         Impl::copy(*this, src);
00334     }
00335 
00336     // Construct a Row from a Row of the same length, with any
00337     // stride. Works as long as the element types are compatible.
00338     template <class EE, int SS> explicit Row(const Row<N,EE,SS>& vv) {
00339         Impl::copy(*this, vv);
00340     }
00341 
00342     // Construction using an element assigns to each element.
00343     explicit Row(const E& e)
00344       { for (int i=0;i<N;++i) d[i*STRIDE]=e; }
00345 
00346     // Construction using a negated element assigns to each element.
00347     explicit Row(const ENeg& ne)
00348       { for (int i=0;i<N;++i) d[i*STRIDE]=ne; }
00349 
00350     // Given an int, turn it into a suitable floating point number
00351     // and then feed that to the above single-element constructor.
00352     explicit Row(int i) 
00353       { new (this) Row(E(Precision(i))); }
00354 
00355     // A bevy of constructors for Rows up to length 6.
00356     Row(const E& e0,const E& e1)
00357       { assert(N==2);(*this)[0]=e0;(*this)[1]=e1; }
00358     Row(const E& e0,const E& e1,const E& e2)
00359       { assert(N==3);(*this)[0]=e0;(*this)[1]=e1;(*this)[2]=e2; }
00360     Row(const E& e0,const E& e1,const E& e2,const E& e3)
00361       { assert(N==4);(*this)[0]=e0;(*this)[1]=e1;(*this)[2]=e2;(*this)[3]=e3; }
00362     Row(const E& e0,const E& e1,const E& e2,const E& e3,const E& e4)
00363       { assert(N==5);(*this)[0]=e0;(*this)[1]=e1;(*this)[2]=e2;
00364         (*this)[3]=e3;(*this)[4]=e4; }
00365     Row(const E& e0,const E& e1,const E& e2,const E& e3,const E& e4,const E& e5)
00366       { assert(N==6);(*this)[0]=e0;(*this)[1]=e1;(*this)[2]=e2;
00367         (*this)[3]=e3;(*this)[4]=e4;(*this)[5]=e5; }
00368     Row(const E& e0,const E& e1,const E& e2,const E& e3,const E& e4,const E& e5,const E& e6)
00369       { assert(N==7);(*this)[0]=e0;(*this)[1]=e1;(*this)[2]=e2;
00370         (*this)[3]=e3;(*this)[4]=e4;(*this)[5]=e5;(*this)[6]=e6; }
00371     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)
00372       { assert(N==8);(*this)[0]=e0;(*this)[1]=e1;(*this)[2]=e2;
00373         (*this)[3]=e3;(*this)[4]=e4;(*this)[5]=e5;(*this)[6]=e6;(*this)[7]=e7; }
00374     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)
00375       { assert(N==9);(*this)[0]=e0;(*this)[1]=e1;(*this)[2]=e2;
00376         (*this)[3]=e3;(*this)[4]=e4;(*this)[5]=e5;(*this)[6]=e6;(*this)[7]=e7;(*this)[8]=e8; }
00377 
00378     // Construction from a pointer to anything assumes we're pointing
00379     // at an element list of the right length.
00380     template <class EE> explicit Row(const EE* p)
00381       { assert(p); for(int i=0;i<N;++i) d[i*STRIDE]=p[i]; }
00382     template <class EE> Row& operator=(const EE* p)
00383       { assert(p); for(int i=0;i<N;++i) d[i*STRIDE]=p[i]; return *this; }
00384 
00385     // Conforming assignment ops.
00386     template <class EE, int SS> Row& operator=(const Row<N,EE,SS>& vv) {
00387         Impl::copy(*this, vv);
00388         return *this;
00389     }
00390     template <class EE, int SS> Row& operator+=(const Row<N,EE,SS>& r)
00391       { for(int i=0;i<N;++i) d[i*STRIDE] += r[i]; return *this; }
00392     template <class EE, int SS> Row& operator+=(const Row<N,negator<EE>,SS>& r)
00393       { for(int i=0;i<N;++i) d[i*STRIDE] -= -(r[i]); return *this; }
00394     template <class EE, int SS> Row& operator-=(const Row<N,EE,SS>& r)
00395       { for(int i=0;i<N;++i) d[i*STRIDE] -= r[i]; return *this; }
00396     template <class EE, int SS> Row& operator-=(const Row<N,negator<EE>,SS>& r)
00397       { for(int i=0;i<N;++i) d[i*STRIDE] += -(r[i]); return *this; }
00398 
00399     // Conforming binary ops with 'this' on left, producing new packed result.
00400     // Cases: r=r+r, r=r-r, s=r*v r=r*m
00401 
00403     template <class EE, int SS> Row<N,typename CNT<E>::template Result<EE>::Add>
00404     conformingAdd(const Row<N,EE,SS>& r) const {
00405         Row<N,typename CNT<E>::template Result<EE>::Add> result;
00406         Impl::conformingAdd(*this, r, result);
00407         return result;
00408     }
00409 
00411     template <class EE, int SS> Row<N,typename CNT<E>::template Result<EE>::Sub>
00412     conformingSubtract(const Row<N,EE,SS>& r) const {
00413         Row<N,typename CNT<E>::template Result<EE>::Sub> result;
00414         Impl::conformingSubtract(*this, r, result);
00415         return result;
00416     }
00417 
00419     template <class EE, int SS> typename CNT<E>::template Result<EE>::Mul
00420     conformingMultiply(const Vec<N,EE,SS>& r) const {
00421         return (*this)*r;
00422     }
00423 
00425     template <int MatNCol, class EE, int CS, int RS> 
00426     Row<MatNCol,typename CNT<E>::template Result<EE>::Mul>
00427     conformingMultiply(const Mat<N,MatNCol,EE,CS,RS>& m) const {
00428         Row<MatNCol,typename CNT<E>::template Result<EE>::Mul> result;
00429         for (int j=0;j<N;++j) result[j] = conformingMultiply(m(j));
00430         return result;
00431     }
00432 
00434     template <class EE, int SS> Row<N,typename CNT<E>::template Result<EE>::Mul>
00435     elementwiseMultiply(const Row<N,EE,SS>& r) const {
00436         Row<N,typename CNT<E>::template Result<EE>::Mul> result;
00437         Impl::elementwiseMultiply(*this, r, result);
00438         return result;
00439     }
00440 
00442     template <class EE, int SS> Row<N,typename CNT<E>::template Result<EE>::Dvd>
00443     elementwiseDivide(const Row<N,EE,SS>& r) const {
00444         Row<N,typename CNT<E>::template Result<EE>::Dvd> result;
00445         Impl::elementwiseDivide(*this, r, result);
00446         return result;
00447     }
00448 
00449     const E& operator[](int i) const { assert(0 <= i && i < N); return d[i*STRIDE]; }
00450     E&       operator[](int i)         { assert(0 <= i && i < N); return d[i*STRIDE]; }
00451     const E& operator()(int i) const { return (*this)[i]; }
00452     E&       operator()(int i)         { return (*this)[i]; }
00453 
00454     ScalarNormSq normSqr() const { return scalarNormSqr(); }
00455     typename CNT<ScalarNormSq>::TSqrt 
00456         norm() const { return CNT<ScalarNormSq>::sqrt(scalarNormSqr()); }
00457 
00458     // If the elements of this Row are scalars, the result is what you get by
00459     // dividing each element by the norm() calculated above. If the elements are
00460     // *not* scalars, then the elements are *separately* normalized. That means
00461     // you will get a different answer from Row<2,Row3>::normalize() than you
00462     // would from a Row<6>::normalize() containing the same scalars.
00463     //
00464     // Normalize returns a row of the same dimension but in new, packed storage
00465     // and with a return type that does not include negator<> even if the original
00466     // Row<> does, because we can eliminate the negation here almost for free.
00467     // But we can't standardize (change conjugate to complex) for free, so we'll retain
00468     // conjugates if there are any.
00469     TNormalize normalize() const {
00470         if (CNT<E>::IsScalar) {
00471             return castAwayNegatorIfAny() / (SignInterpretation*norm());
00472         } else {
00473             TNormalize elementwiseNormalized;
00474             for (int j=0; j<N; ++j) 
00475                 elementwiseNormalized[j] = CNT<E>::normalize((*this)[j]);
00476             return elementwiseNormalized;
00477         }
00478     }
00479 
00480     TInvert invert() const {assert(false); return TInvert();} // TODO default inversion
00481 
00482     const Row&   operator+() const { return *this; }
00483     const TNeg&  operator-() const { return negate(); }
00484     TNeg&        operator-()       { return updNegate(); }
00485     const THerm& operator~() const { return transpose(); }
00486     THerm&       operator~()       { return updTranspose(); }
00487 
00488     const TNeg&  negate() const { return *reinterpret_cast<const TNeg*>(this); }
00489     TNeg&        updNegate()    { return *reinterpret_cast<TNeg*>(this); }
00490 
00491     const THerm& transpose()    const { return *reinterpret_cast<const THerm*>(this); }
00492     THerm&       updTranspose()       { return *reinterpret_cast<THerm*>(this); }
00493 
00494     const TPosTrans& positionalTranspose() const
00495         { return *reinterpret_cast<const TPosTrans*>(this); }
00496     TPosTrans&       updPositionalTranspose()
00497         { return *reinterpret_cast<TPosTrans*>(this); }
00498 
00499     const TReal& real() const { return *reinterpret_cast<const TReal*>(this); }
00500     TReal&       real()       { return *reinterpret_cast<      TReal*>(this); }
00501 
00502     // Had to contort these routines to get them through VC++ 7.net
00503     const TImag& imag()    const { 
00504         const int offs = ImagOffset;
00505         const EImag* p = reinterpret_cast<const EImag*>(this);
00506         return *reinterpret_cast<const TImag*>(p+offs);
00507     }
00508     TImag& imag() { 
00509         const int offs = ImagOffset;
00510         EImag* p = reinterpret_cast<EImag*>(this);
00511         return *reinterpret_cast<TImag*>(p+offs);
00512     }
00513 
00514     const TWithoutNegator& castAwayNegatorIfAny() const {return *reinterpret_cast<const TWithoutNegator*>(this);}
00515     TWithoutNegator&       updCastAwayNegatorIfAny()    {return *reinterpret_cast<TWithoutNegator*>(this);}
00516 
00517 
00518     // These are elementwise binary operators, (this op ee) by default but 
00519     // (ee op this) if 'FromLeft' appears in the name. The result is a packed 
00520     // Row<N> but the element type may change. These are mostly used to 
00521     // implement global operators. We call these "scalar" operators but 
00522     // actually the "scalar" can be a composite type.
00523 
00524     //TODO: consider converting 'e' to Standard Numbers as precalculation and 
00525     // changing return type appropriately.
00526     template <class EE> Row<N, typename CNT<E>::template Result<EE>::Mul>
00527     scalarMultiply(const EE& e) const {
00528         Row<N, typename CNT<E>::template Result<EE>::Mul> result;
00529         for (int j=0; j<N; ++j) result[j] = (*this)[j] * e;
00530         return result;
00531     }
00532     template <class EE> Row<N, typename CNT<EE>::template Result<E>::Mul>
00533     scalarMultiplyFromLeft(const EE& e) const {
00534         Row<N, typename CNT<EE>::template Result<E>::Mul> result;
00535         for (int j=0; j<N; ++j) result[j] = e * (*this)[j];
00536         return result;
00537     }
00538 
00539     // TODO: should precalculate and store 1/e, while converting to Standard 
00540     // Numbers. Note that return type should change appropriately.
00541     template <class EE> Row<N, typename CNT<E>::template Result<EE>::Dvd>
00542     scalarDivide(const EE& e) const {
00543         Row<N, typename CNT<E>::template Result<EE>::Dvd> result;
00544         for (int j=0; j<N; ++j) result[j] = (*this)[j] / e;
00545         return result;
00546     }
00547     template <class EE> Row<N, typename CNT<EE>::template Result<E>::Dvd>
00548     scalarDivideFromLeft(const EE& e) const {
00549         Row<N, typename CNT<EE>::template Result<E>::Dvd> result;
00550         for (int j=0; j<N; ++j) result[j] = e / (*this)[j];
00551         return result;
00552     }
00553 
00554     template <class EE> Row<N, typename CNT<E>::template Result<EE>::Add>
00555     scalarAdd(const EE& e) const {
00556         Row<N, typename CNT<E>::template Result<EE>::Add> result;
00557         for (int j=0; j<N; ++j) result[j] = (*this)[j] + e;
00558         return result;
00559     }
00560     // Add is commutative, so no 'FromLeft'.
00561 
00562     template <class EE> Row<N, typename CNT<E>::template Result<EE>::Sub>
00563     scalarSubtract(const EE& e) const {
00564         Row<N, typename CNT<E>::template Result<EE>::Sub> result;
00565         for (int j=0; j<N; ++j) result[j] = (*this)[j] - e;
00566         return result;
00567     }
00568     template <class EE> Row<N, typename CNT<EE>::template Result<E>::Sub>
00569     scalarSubtractFromLeft(const EE& e) const {
00570         Row<N, typename CNT<EE>::template Result<E>::Sub> result;
00571         for (int j=0; j<N; ++j) result[j] = e - (*this)[j];
00572         return result;
00573     }
00574 
00575     // Generic assignments for any element type not listed explicitly, including scalars.
00576     // These are done repeatedly for each element and only work if the operation can
00577     // be performed leaving the original element type.
00578     template <class EE> Row& operator =(const EE& e) {return scalarEq(e);}
00579     template <class EE> Row& operator+=(const EE& e) {return scalarPlusEq(e);}
00580     template <class EE> Row& operator-=(const EE& e) {return scalarMinusEq(e);}
00581     template <class EE> Row& operator*=(const EE& e) {return scalarTimesEq(e);}
00582     template <class EE> Row& operator/=(const EE& e) {return scalarDivideEq(e);}
00583 
00584     // Generalized scalar assignment & computed assignment methods. These will work
00585     // for any assignment-compatible element, not just scalars.
00586     template <class EE> Row& scalarEq(const EE& ee)
00587       { for(int i=0;i<N;++i) d[i*STRIDE] = ee; return *this; }
00588     template <class EE> Row& scalarPlusEq(const EE& ee)
00589       { for(int i=0;i<N;++i) d[i*STRIDE] += ee; return *this; }
00590     template <class EE> Row& scalarMinusEq(const EE& ee)
00591       { for(int i=0;i<N;++i) d[i*STRIDE] -= ee; return *this; }
00592     template <class EE> Row& scalarMinusEqFromLeft(const EE& ee)
00593       { for(int i=0;i<N;++i) d[i*STRIDE] = ee - d[i*STRIDE]; return *this; }
00594     template <class EE> Row& scalarTimesEq(const EE& ee)
00595       { for(int i=0;i<N;++i) d[i*STRIDE] *= ee; return *this; }
00596     template <class EE> Row& scalarTimesEqFromLeft(const EE& ee)
00597       { for(int i=0;i<N;++i) d[i*STRIDE] = ee * d[i*STRIDE]; return *this; }
00598     template <class EE> Row& scalarDivideEq(const EE& ee)
00599       { for(int i=0;i<N;++i) d[i*STRIDE] /= ee; return *this; }
00600     template <class EE> Row& scalarDivideEqFromLeft(const EE& ee)
00601       { for(int i=0;i<N;++i) d[i*STRIDE] = ee / d[i*STRIDE]; return *this; }
00602 
00603 
00604     // Specialize for int to avoid warnings and ambiguities.
00605     Row& scalarEq(int ee)       {return scalarEq(Precision(ee));}
00606     Row& scalarPlusEq(int ee)   {return scalarPlusEq(Precision(ee));}
00607     Row& scalarMinusEq(int ee)  {return scalarMinusEq(Precision(ee));}
00608     Row& scalarTimesEq(int ee)  {return scalarTimesEq(Precision(ee));}
00609     Row& scalarDivideEq(int ee) {return scalarDivideEq(Precision(ee));}
00610     Row& scalarMinusEqFromLeft(int ee)  {return scalarMinusEqFromLeft(Precision(ee));}
00611     Row& scalarTimesEqFromLeft(int ee)  {return scalarTimesEqFromLeft(Precision(ee));}
00612     Row& scalarDivideEqFromLeft(int ee) {return scalarDivideEqFromLeft(Precision(ee));}
00613 
00616     void setToNaN() {
00617         (*this) = CNT<ELT>::getNaN();
00618     }
00619 
00621     void setToZero() {
00622         (*this) = ELT(0);
00623     }
00624 
00630     template <int NN>
00631     const Row<NN,ELT,STRIDE>& getSubRow(int j) const {
00632         assert(0 <= j && j + NN <= N);
00633         return Row<NN,ELT,STRIDE>::getAs(&(*this)[j]);
00634     }
00640     template <int NN>
00641     Row<NN,ELT,STRIDE>& updSubRow(int j) {
00642         assert(0 <= j && j + NN <= N);
00643         return Row<NN,ELT,STRIDE>::updAs(&(*this)[j]);
00644     }
00645 
00649     template <int NN>
00650     static const Row& getSubRow(const Row<NN,ELT,STRIDE>& r, int j) {
00651         assert(0 <= j && j + N <= NN);
00652         return getAs(&r[j]);
00653     }
00657     template <int NN>
00658     static Row& updSubRow(Row<NN,ELT,STRIDE>& r, int j) {
00659         assert(0 <= j && j + N <= NN);
00660         return updAs(&r[j]);
00661     }
00662 
00666     Row<N-1,ELT,1> drop1(int p) const {
00667         assert(0 <= p && p < N);
00668         Row<N-1,ELT,1> out;
00669         int nxt=0;
00670         for (int i=0; i<N-1; ++i, ++nxt) {
00671             if (nxt==p) ++nxt;  // skip the loser
00672             out[i] = (*this)[nxt];
00673         }
00674         return out;
00675     }
00676 
00680     template <class EE> Row<N+1,ELT,1> append1(const EE& v) const {
00681         Row<N+1,ELT,1> out;
00682         Row<N,ELT,1>::updAs(&out[0]) = (*this);
00683         out[N] = v;
00684         return out;
00685     }
00686 
00687 
00693     template <class EE> Row<N+1,ELT,1> insert1(int p, const EE& v) const {
00694         assert(0 <= p && p <= N);
00695         if (p==N) return append1(v);
00696         Row<N+1,ELT,1> out;
00697         int nxt=0;
00698         for (int i=0; i<N; ++i, ++nxt) {
00699             if (i==p) out[nxt++] = v;
00700             out[nxt] = (*this)[i];
00701         }
00702         return out;
00703     }
00704 
00707     static const Row& getAs(const ELT* p)  {return *reinterpret_cast<const Row*>(p);}
00710     static Row&       updAs(ELT* p)        {return *reinterpret_cast<Row*>(p);}
00711 
00715     static Row<N,ELT,1> getNaN() { return Row<N,ELT,1>(CNT<ELT>::getNaN()); }
00716 
00718     bool isNaN() const {
00719         for (int j=0; j<N; ++j)
00720             if (CNT<ELT>::isNaN((*this)[j]))
00721                 return true;
00722         return false;
00723     }
00724 
00727     bool isInf() const {
00728         bool seenInf = false;
00729         for (int j=0; j<N; ++j) {
00730             const ELT& e = (*this)[j];
00731             if (!CNT<ELT>::isFinite(e)) {
00732                 if (!CNT<ELT>::isInf(e)) 
00733                     return false; // something bad was found
00734                 seenInf = true; 
00735             }
00736         }
00737         return seenInf;
00738     }
00739 
00742     bool isFinite() const {
00743         for (int j=0; j<N; ++j)
00744             if (!CNT<ELT>::isFinite((*this)[j]))
00745                 return false;
00746         return true;
00747     }
00748 
00751     static double getDefaultTolerance() {return CNT<ELT>::getDefaultTolerance();}
00752 
00755     template <class E2, int CS2>
00756     bool isNumericallyEqual(const Row<N,E2,CS2>& r, double tol) const {
00757         for (int j=0; j<N; ++j)
00758             if (!CNT<ELT>::isNumericallyEqual((*this)(j), r(j), tol))
00759                 return false;
00760         return true;
00761     }
00762 
00766     template <class E2, int CS2>
00767     bool isNumericallyEqual(const Row<N,E2,CS2>& r) const {
00768         const double tol = std::max(getDefaultTolerance(),r.getDefaultTolerance());
00769         return isNumericallyEqual(r, tol);
00770     }
00771 
00776     bool isNumericallyEqual
00777        (const ELT& e,
00778         double     tol = getDefaultTolerance()) const 
00779     {
00780         for (int j=0; j<N; ++j)
00781             if (!CNT<ELT>::isNumericallyEqual((*this)(j), e, tol))
00782                 return false;
00783         return true;
00784     }
00785 private:
00786     ELT d[NActualElements];    // data
00787 };
00788 
00790 // Global operators involving two rows.    //
00791 //   v+v, v-v, v==v, v!=v                  //
00793 
00794 // v3 = v1 + v2 where all v's have the same length N. 
00795 template <int N, class E1, int S1, class E2, int S2> inline
00796 typename Row<N,E1,S1>::template Result< Row<N,E2,S2> >::Add
00797 operator+(const Row<N,E1,S1>& l, const Row<N,E2,S2>& r) { 
00798     return Row<N,E1,S1>::template Result< Row<N,E2,S2> >
00799         ::AddOp::perform(l,r);
00800 }
00801 
00802 // v3 = v1 - v2, similar to +
00803 template <int N, class E1, int S1, class E2, int S2> inline
00804 typename Row<N,E1,S1>::template Result< Row<N,E2,S2> >::Sub
00805 operator-(const Row<N,E1,S1>& l, const Row<N,E2,S2>& r) { 
00806     return Row<N,E1,S1>::template Result< Row<N,E2,S2> >
00807         ::SubOp::perform(l,r);
00808 }
00809 
00811 template <int N, class E1, int S1, class E2, int S2> inline bool
00812 operator==(const Row<N,E1,S1>& l, const Row<N,E2,S2>& r) { 
00813     for (int i=0; i < N; ++i) if (l[i] != r[i]) return false;
00814     return true;
00815 }
00817 template <int N, class E1, int S1, class E2, int S2> inline bool
00818 operator!=(const Row<N,E1,S1>& l, const Row<N,E2,S2>& r) {return !(l==r);} 
00819 
00821 template <int N, class E1, int S1, class E2, int S2> inline bool
00822 operator<(const Row<N,E1,S1>& l, const Row<N,E2,S2>& r) 
00823 {   for (int i=0; i < N; ++i) if (l[i] >= r[i]) return false;
00824     return true; }
00826 template <int N, class E1, int S1, class E2> inline bool
00827 operator<(const Row<N,E1,S1>& v, const E2& e) 
00828 {   for (int i=0; i < N; ++i) if (v[i] >= e) return false;
00829     return true; }
00830 
00832 template <int N, class E1, int S1, class E2, int S2> inline bool
00833 operator>(const Row<N,E1,S1>& l, const Row<N,E2,S2>& r) 
00834 {   for (int i=0; i < N; ++i) if (l[i] <= r[i]) return false;
00835     return true; }
00837 template <int N, class E1, int S1, class E2> inline bool
00838 operator>(const Row<N,E1,S1>& v, const E2& e) 
00839 {   for (int i=0; i < N; ++i) if (v[i] <= e) return false;
00840     return true; }
00841 
00844 template <int N, class E1, int S1, class E2, int S2> inline bool
00845 operator<=(const Row<N,E1,S1>& l, const Row<N,E2,S2>& r) 
00846 {   for (int i=0; i < N; ++i) if (l[i] > r[i]) return false;
00847     return true; }
00850 template <int N, class E1, int S1, class E2> inline bool
00851 operator<=(const Row<N,E1,S1>& v, const E2& e) 
00852 {   for (int i=0; i < N; ++i) if (v[i] > e) return false;
00853     return true; }
00854 
00857 template <int N, class E1, int S1, class E2, int S2> inline bool
00858 operator>=(const Row<N,E1,S1>& l, const Row<N,E2,S2>& r) 
00859 {   for (int i=0; i < N; ++i) if (l[i] < r[i]) return false;
00860     return true; }
00863 template <int N, class E1, int S1, class E2> inline bool
00864 operator>=(const Row<N,E1,S1>& v, const E2& e) 
00865 {   for (int i=0; i < N; ++i) if (v[i] < e) return false;
00866     return true; }
00867 
00869 // Global operators involving a row and a scalar. //
00871 
00872 // I haven't been able to figure out a nice way to templatize for the
00873 // built-in reals without introducing a lot of unwanted type matches
00874 // as well. So we'll just grind them out explicitly here.
00875 
00876 // SCALAR MULTIPLY
00877 
00878 // v = v*real, real*v 
00879 template <int N, class E, int S> inline
00880 typename Row<N,E,S>::template Result<float>::Mul
00881 operator*(const Row<N,E,S>& l, const float& r)
00882   { return Row<N,E,S>::template Result<float>::MulOp::perform(l,r); }
00883 template <int N, class E, int S> inline
00884 typename Row<N,E,S>::template Result<float>::Mul
00885 operator*(const float& l, const Row<N,E,S>& r) {return r*l;}
00886 
00887 template <int N, class E, int S> inline
00888 typename Row<N,E,S>::template Result<double>::Mul
00889 operator*(const Row<N,E,S>& l, const double& r)
00890   { return Row<N,E,S>::template Result<double>::MulOp::perform(l,r); }
00891 template <int N, class E, int S> inline
00892 typename Row<N,E,S>::template Result<double>::Mul
00893 operator*(const double& l, const Row<N,E,S>& r) {return r*l;}
00894 
00895 template <int N, class E, int S> inline
00896 typename Row<N,E,S>::template Result<long double>::Mul
00897 operator*(const Row<N,E,S>& l, const long double& r)
00898   { return Row<N,E,S>::template Result<long double>::MulOp::perform(l,r); }
00899 template <int N, class E, int S> inline
00900 typename Row<N,E,S>::template Result<long double>::Mul
00901 operator*(const long double& l, const Row<N,E,S>& r) {return r*l;}
00902 
00903 // v = v*int, int*v -- just convert int to v's precision float
00904 template <int N, class E, int S> inline
00905 typename Row<N,E,S>::template Result<typename CNT<E>::Precision>::Mul
00906 operator*(const Row<N,E,S>& l, int r) {return l * (typename CNT<E>::Precision)r;}
00907 template <int N, class E, int S> inline
00908 typename Row<N,E,S>::template Result<typename CNT<E>::Precision>::Mul
00909 operator*(int l, const Row<N,E,S>& r) {return r * (typename CNT<E>::Precision)l;}
00910 
00911 // Complex, conjugate, and negator are all easy to templatize.
00912 
00913 // v = v*complex, complex*v
00914 template <int N, class E, int S, class R> inline
00915 typename Row<N,E,S>::template Result<std::complex<R> >::Mul
00916 operator*(const Row<N,E,S>& l, const std::complex<R>& r)
00917   { return Row<N,E,S>::template Result<std::complex<R> >::MulOp::perform(l,r); }
00918 template <int N, class E, int S, class R> inline
00919 typename Row<N,E,S>::template Result<std::complex<R> >::Mul
00920 operator*(const std::complex<R>& l, const Row<N,E,S>& r) {return r*l;}
00921 
00922 // v = v*conjugate, conjugate*v (convert conjugate->complex)
00923 template <int N, class E, int S, class R> inline
00924 typename Row<N,E,S>::template Result<std::complex<R> >::Mul
00925 operator*(const Row<N,E,S>& l, const conjugate<R>& r) {return l*(std::complex<R>)r;}
00926 template <int N, class E, int S, class R> inline
00927 typename Row<N,E,S>::template Result<std::complex<R> >::Mul
00928 operator*(const conjugate<R>& l, const Row<N,E,S>& r) {return r*(std::complex<R>)l;}
00929 
00930 // v = v*negator, negator*v: convert negator to standard number
00931 template <int N, class E, int S, class R> inline
00932 typename Row<N,E,S>::template Result<typename negator<R>::StdNumber>::Mul
00933 operator*(const Row<N,E,S>& l, const negator<R>& r) {return l * (typename negator<R>::StdNumber)(R)r;}
00934 template <int N, class E, int S, class R> inline
00935 typename Row<N,E,S>::template Result<typename negator<R>::StdNumber>::Mul
00936 operator*(const negator<R>& l, const Row<N,E,S>& r) {return r * (typename negator<R>::StdNumber)(R)l;}
00937 
00938 
00939 // SCALAR DIVIDE. This is a scalar operation when the scalar is on the right,
00940 // but when it is on the left it means scalar * pseudoInverse(row), which is 
00941 // a vec.
00942 
00943 // v = v/real, real/v 
00944 template <int N, class E, int S> inline
00945 typename Row<N,E,S>::template Result<float>::Dvd
00946 operator/(const Row<N,E,S>& l, const float& r)
00947   { return Row<N,E,S>::template Result<float>::DvdOp::perform(l,r); }
00948 template <int N, class E, int S> inline
00949 typename CNT<float>::template Result<Row<N,E,S> >::Dvd
00950 operator/(const float& l, const Row<N,E,S>& r)
00951   { return CNT<float>::template Result<Row<N,E,S> >::DvdOp::perform(l,r); }
00952 
00953 template <int N, class E, int S> inline
00954 typename Row<N,E,S>::template Result<double>::Dvd
00955 operator/(const Row<N,E,S>& l, const double& r)
00956   { return Row<N,E,S>::template Result<double>::DvdOp::perform(l,r); }
00957 template <int N, class E, int S> inline
00958 typename CNT<double>::template Result<Row<N,E,S> >::Dvd
00959 operator/(const double& l, const Row<N,E,S>& r)
00960   { return CNT<double>::template Result<Row<N,E,S> >::DvdOp::perform(l,r); }
00961 
00962 template <int N, class E, int S> inline
00963 typename Row<N,E,S>::template Result<long double>::Dvd
00964 operator/(const Row<N,E,S>& l, const long double& r)
00965   { return Row<N,E,S>::template Result<long double>::DvdOp::perform(l,r); }
00966 template <int N, class E, int S> inline
00967 typename CNT<long double>::template Result<Row<N,E,S> >::Dvd
00968 operator/(const long double& l, const Row<N,E,S>& r)
00969   { return CNT<long double>::template Result<Row<N,E,S> >::DvdOp::perform(l,r); }
00970 
00971 // v = v/int, int/v -- just convert int to v's precision float
00972 template <int N, class E, int S> inline
00973 typename Row<N,E,S>::template Result<typename CNT<E>::Precision>::Dvd
00974 operator/(const Row<N,E,S>& l, int r) {return l / (typename CNT<E>::Precision)r;}
00975 template <int N, class E, int S> inline
00976 typename CNT<typename CNT<E>::Precision>::template Result<Row<N,E,S> >::Dvd
00977 operator/(int l, const Row<N,E,S>& r) {return (typename CNT<E>::Precision)l / r;}
00978 
00979 
00980 // Complex, conjugate, and negator are all easy to templatize.
00981 
00982 // v = v/complex, complex/v
00983 template <int N, class E, int S, class R> inline
00984 typename Row<N,E,S>::template Result<std::complex<R> >::Dvd
00985 operator/(const Row<N,E,S>& l, const std::complex<R>& r)
00986   { return Row<N,E,S>::template Result<std::complex<R> >::DvdOp::perform(l,r); }
00987 template <int N, class E, int S, class R> inline
00988 typename CNT<std::complex<R> >::template Result<Row<N,E,S> >::Dvd
00989 operator/(const std::complex<R>& l, const Row<N,E,S>& r)
00990   { return CNT<std::complex<R> >::template Result<Row<N,E,S> >::DvdOp::perform(l,r); }
00991 
00992 // v = v/conjugate, conjugate/v (convert conjugate->complex)
00993 template <int N, class E, int S, class R> inline
00994 typename Row<N,E,S>::template Result<std::complex<R> >::Dvd
00995 operator/(const Row<N,E,S>& l, const conjugate<R>& r) {return l/(std::complex<R>)r;}
00996 template <int N, class E, int S, class R> inline
00997 typename CNT<std::complex<R> >::template Result<Row<N,E,S> >::Dvd
00998 operator/(const conjugate<R>& l, const Row<N,E,S>& r) {return (std::complex<R>)l/r;}
00999 
01000 // v = v/negator, negator/v: convert negator to number
01001 template <int N, class E, int S, class R> inline
01002 typename Row<N,E,S>::template Result<typename negator<R>::StdNumber>::Dvd
01003 operator/(const Row<N,E,S>& l, const negator<R>& r) {return l/(typename negator<R>::StdNumber)(R)r;}
01004 template <int N, class E, int S, class R> inline
01005 typename CNT<R>::template Result<Row<N,E,S> >::Dvd
01006 operator/(const negator<R>& l, const Row<N,E,S>& r) {return (typename negator<R>::StdNumber)(R)l/r;}
01007 
01008 
01009 // Add and subtract are odd as scalar ops. They behave as though the
01010 // scalar stands for a vector each of whose elements is that scalar,
01011 // and then a normal vector add or subtract is done.
01012 
01013 // SCALAR ADD
01014 
01015 // v = v+real, real+v 
01016 template <int N, class E, int S> inline
01017 typename Row<N,E,S>::template Result<float>::Add
01018 operator+(const Row<N,E,S>& l, const float& r)
01019   { return Row<N,E,S>::template Result<float>::AddOp::perform(l,r); }
01020 template <int N, class E, int S> inline
01021 typename Row<N,E,S>::template Result<float>::Add
01022 operator+(const float& l, const Row<N,E,S>& r) {return r+l;}
01023 
01024 template <int N, class E, int S> inline
01025 typename Row<N,E,S>::template Result<double>::Add
01026 operator+(const Row<N,E,S>& l, const double& r)
01027   { return Row<N,E,S>::template Result<double>::AddOp::perform(l,r); }
01028 template <int N, class E, int S> inline
01029 typename Row<N,E,S>::template Result<double>::Add
01030 operator+(const double& l, const Row<N,E,S>& r) {return r+l;}
01031 
01032 template <int N, class E, int S> inline
01033 typename Row<N,E,S>::template Result<long double>::Add
01034 operator+(const Row<N,E,S>& l, const long double& r)
01035   { return Row<N,E,S>::template Result<long double>::AddOp::perform(l,r); }
01036 template <int N, class E, int S> inline
01037 typename Row<N,E,S>::template Result<long double>::Add
01038 operator+(const long double& l, const Row<N,E,S>& r) {return r+l;}
01039 
01040 // v = v+int, int+v -- just convert int to v's precision float
01041 template <int N, class E, int S> inline
01042 typename Row<N,E,S>::template Result<typename CNT<E>::Precision>::Add
01043 operator+(const Row<N,E,S>& l, int r) {return l + (typename CNT<E>::Precision)r;}
01044 template <int N, class E, int S> inline
01045 typename Row<N,E,S>::template Result<typename CNT<E>::Precision>::Add
01046 operator+(int l, const Row<N,E,S>& r) {return r + (typename CNT<E>::Precision)l;}
01047 
01048 // Complex, conjugate, and negator are all easy to templatize.
01049 
01050 // v = v+complex, complex+v
01051 template <int N, class E, int S, class R> inline
01052 typename Row<N,E,S>::template Result<std::complex<R> >::Add
01053 operator+(const Row<N,E,S>& l, const std::complex<R>& r)
01054   { return Row<N,E,S>::template Result<std::complex<R> >::AddOp::perform(l,r); }
01055 template <int N, class E, int S, class R> inline
01056 typename Row<N,E,S>::template Result<std::complex<R> >::Add
01057 operator+(const std::complex<R>& l, const Row<N,E,S>& r) {return r+l;}
01058 
01059 // v = v+conjugate, conjugate+v (convert conjugate->complex)
01060 template <int N, class E, int S, class R> inline
01061 typename Row<N,E,S>::template Result<std::complex<R> >::Add
01062 operator+(const Row<N,E,S>& l, const conjugate<R>& r) {return l+(std::complex<R>)r;}
01063 template <int N, class E, int S, class R> inline
01064 typename Row<N,E,S>::template Result<std::complex<R> >::Add
01065 operator+(const conjugate<R>& l, const Row<N,E,S>& r) {return r+(std::complex<R>)l;}
01066 
01067 // v = v+negator, negator+v: convert negator to standard number
01068 template <int N, class E, int S, class R> inline
01069 typename Row<N,E,S>::template Result<typename negator<R>::StdNumber>::Add
01070 operator+(const Row<N,E,S>& l, const negator<R>& r) {return l + (typename negator<R>::StdNumber)(R)r;}
01071 template <int N, class E, int S, class R> inline
01072 typename Row<N,E,S>::template Result<typename negator<R>::StdNumber>::Add
01073 operator+(const negator<R>& l, const Row<N,E,S>& r) {return r + (typename negator<R>::StdNumber)(R)l;}
01074 
01075 // SCALAR SUBTRACT -- careful, not commutative.
01076 
01077 // v = v-real, real-v 
01078 template <int N, class E, int S> inline
01079 typename Row<N,E,S>::template Result<float>::Sub
01080 operator-(const Row<N,E,S>& l, const float& r)
01081   { return Row<N,E,S>::template Result<float>::SubOp::perform(l,r); }
01082 template <int N, class E, int S> inline
01083 typename CNT<float>::template Result<Row<N,E,S> >::Sub
01084 operator-(const float& l, const Row<N,E,S>& r)
01085   { return CNT<float>::template Result<Row<N,E,S> >::SubOp::perform(l,r); }
01086 
01087 template <int N, class E, int S> inline
01088 typename Row<N,E,S>::template Result<double>::Sub
01089 operator-(const Row<N,E,S>& l, const double& r)
01090   { return Row<N,E,S>::template Result<double>::SubOp::perform(l,r); }
01091 template <int N, class E, int S> inline
01092 typename CNT<double>::template Result<Row<N,E,S> >::Sub
01093 operator-(const double& l, const Row<N,E,S>& r)
01094   { return CNT<double>::template Result<Row<N,E,S> >::SubOp::perform(l,r); }
01095 
01096 template <int N, class E, int S> inline
01097 typename Row<N,E,S>::template Result<long double>::Sub
01098 operator-(const Row<N,E,S>& l, const long double& r)
01099   { return Row<N,E,S>::template Result<long double>::SubOp::perform(l,r); }
01100 template <int N, class E, int S> inline
01101 typename CNT<long double>::template Result<Row<N,E,S> >::Sub
01102 operator-(const long double& l, const Row<N,E,S>& r)
01103   { return CNT<long double>::template Result<Row<N,E,S> >::SubOp::perform(l,r); }
01104 
01105 // v = v-int, int-v // just convert int to v's precision float
01106 template <int N, class E, int S> inline
01107 typename Row<N,E,S>::template Result<typename CNT<E>::Precision>::Sub
01108 operator-(const Row<N,E,S>& l, int r) {return l - (typename CNT<E>::Precision)r;}
01109 template <int N, class E, int S> inline
01110 typename CNT<typename CNT<E>::Precision>::template Result<Row<N,E,S> >::Sub
01111 operator-(int l, const Row<N,E,S>& r) {return (typename CNT<E>::Precision)l - r;}
01112 
01113 
01114 // Complex, conjugate, and negator are all easy to templatize.
01115 
01116 // v = v-complex, complex-v
01117 template <int N, class E, int S, class R> inline
01118 typename Row<N,E,S>::template Result<std::complex<R> >::Sub
01119 operator-(const Row<N,E,S>& l, const std::complex<R>& r)
01120   { return Row<N,E,S>::template Result<std::complex<R> >::SubOp::perform(l,r); }
01121 template <int N, class E, int S, class R> inline
01122 typename CNT<std::complex<R> >::template Result<Row<N,E,S> >::Sub
01123 operator-(const std::complex<R>& l, const Row<N,E,S>& r)
01124   { return CNT<std::complex<R> >::template Result<Row<N,E,S> >::SubOp::perform(l,r); }
01125 
01126 // v = v-conjugate, conjugate-v (convert conjugate->complex)
01127 template <int N, class E, int S, class R> inline
01128 typename Row<N,E,S>::template Result<std::complex<R> >::Sub
01129 operator-(const Row<N,E,S>& l, const conjugate<R>& r) {return l-(std::complex<R>)r;}
01130 template <int N, class E, int S, class R> inline
01131 typename CNT<std::complex<R> >::template Result<Row<N,E,S> >::Sub
01132 operator-(const conjugate<R>& l, const Row<N,E,S>& r) {return (std::complex<R>)l-r;}
01133 
01134 // v = v-negator, negator-v: convert negator to standard number
01135 template <int N, class E, int S, class R> inline
01136 typename Row<N,E,S>::template Result<typename negator<R>::StdNumber>::Sub
01137 operator-(const Row<N,E,S>& l, const negator<R>& r) {return l-(typename negator<R>::StdNumber)(R)r;}
01138 template <int N, class E, int S, class R> inline
01139 typename CNT<R>::template Result<Row<N,E,S> >::Sub
01140 operator-(const negator<R>& l, const Row<N,E,S>& r) {return (typename negator<R>::StdNumber)(R)l-r;}
01141 
01142 
01143 // Row I/O
01144 template <int N, class E, int S, class CHAR, class TRAITS> inline
01145 std::basic_ostream<CHAR,TRAITS>&
01146 operator<<(std::basic_ostream<CHAR,TRAITS>& o, const Row<N,E,S>& v) {
01147     o << "[" << v[0]; for(int i=1;i<N;++i) o<<','<<v[i]; o<<']'; return o;
01148 }
01149 
01152 template <int N, class E, int S, class CHAR, class TRAITS> inline
01153 std::basic_istream<CHAR,TRAITS>&
01154 operator>>(std::basic_istream<CHAR,TRAITS>& is, Row<N,E,S>& v) {
01155     CHAR openBracket, closeBracket;
01156     is >> openBracket; if (is.fail()) return is;
01157     if (openBracket==CHAR('('))
01158         closeBracket = CHAR(')');
01159     else if (openBracket==CHAR('['))
01160         closeBracket = CHAR(']');
01161     else {
01162         closeBracket = CHAR(0);
01163         is.unget(); if (is.fail()) return is;
01164     }
01165 
01166     for (int i=0; i < N; ++i) {
01167         is >> v[i];
01168         if (is.fail()) return is;
01169         if (i != N-1) {
01170             CHAR c; is >> c; if (is.fail()) return is;
01171             if (c != ',') is.unget();
01172             if (is.fail()) return is;
01173         }
01174     }
01175 
01176     // Get the closing bracket if there was an opening one. If we don't
01177     // see the expected character we'll set the fail bit in the istream.
01178     if (closeBracket != CHAR(0)) {
01179         CHAR closer; is >> closer; if (is.fail()) return is;
01180         if (closer != closeBracket) {
01181             is.unget(); if (is.fail()) return is;
01182             is.setstate( std::ios::failbit );
01183         }
01184     }
01185 
01186     return is;
01187 }
01188 
01189 } //namespace SimTK
01190 
01191 
01192 #endif //SimTK_SIMMATRIX_SMALLMATRIX_ROW_H_
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines