Simbody
3.5
|
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_