00001 #ifndef SimTK_SIMMATRIX_SMALLMATRIX_ROW_H_
00002 #define SimTK_SIMMATRIX_SMALLMATRIX_ROW_H_
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00039 #include "SimTKcommon/internal/common.h"
00040
00041
00042 namespace SimTK {
00043
00044
00045
00046 namespace Impl {
00047
00048
00049
00050
00051 template <class E1, int S1, class E2, int S2> void
00052 conformingAdd(const Row<1,E1,S1>& r1, const Row<1,E2,S2>& r2, Row<1,typename CNT<E1>::template Result<E2>::Add>& result) {
00053 result[0] = r1[0] + r2[0];
00054 }
00055 template <int N, class E1, int S1, class E2, int S2> void
00056 conformingAdd(const Row<N,E1,S1>& r1, const Row<N,E2,S2>& r2, Row<N,typename CNT<E1>::template Result<E2>::Add>& result) {
00057 conformingAdd(reinterpret_cast<const Row<N-1,E1,S1>&>(r1), reinterpret_cast<const Row<N-1,E2,S2>&>(r2), reinterpret_cast<Row<N-1,typename CNT<E1>::template Result<E2>::Add>&>(result));
00058 result[N-1] = r1[N-1] + r2[N-1];
00059 }
00060 template <class E1, int S1, class E2, int S2> void
00061 conformingSubtract(const Row<1,E1,S1>& r1, const Row<1,E2,S2>& r2, Row<1,typename CNT<E1>::template Result<E2>::Add>& result) {
00062 result[0] = r1[0] - r2[0];
00063 }
00064 template <int N, class E1, int S1, class E2, int S2> void
00065 conformingSubtract(const Row<N,E1,S1>& r1, const Row<N,E2,S2>& r2, Row<N,typename CNT<E1>::template Result<E2>::Add>& result) {
00066 conformingSubtract(reinterpret_cast<const Row<N-1,E1,S1>&>(r1), reinterpret_cast<const Row<N-1,E2,S2>&>(r2), reinterpret_cast<Row<N-1,typename CNT<E1>::template Result<E2>::Add>&>(result));
00067 result[N-1] = r1[N-1] - r2[N-1];
00068 }
00069 template <class E1, int S1, class E2, int S2> void
00070 copy(Row<1,E1,S1>& r1, const Row<1,E2,S2>& r2) {
00071 r1[0] = r2[0];
00072 }
00073 template <int N, class E1, int S1, class E2, int S2> void
00074 copy(Row<N,E1,S1>& r1, const Row<N,E2,S2>& r2) {
00075 copy(reinterpret_cast<Row<N-1,E1,S1>&>(r1), reinterpret_cast<const Row<N-1,E2,S2>&>(r2));
00076 r1[N-1] = r2[N-1];
00077 }
00078
00079 }
00080
00082 template <int N, class ELT, int STRIDE> class Row {
00083 typedef ELT E;
00084 typedef typename CNT<E>::TNeg ENeg;
00085 typedef typename CNT<E>::TWithoutNegator EWithoutNegator;
00086 typedef typename CNT<E>::TReal EReal;
00087 typedef typename CNT<E>::TImag EImag;
00088 typedef typename CNT<E>::TComplex EComplex;
00089 typedef typename CNT<E>::THerm EHerm;
00090 typedef typename CNT<E>::TPosTrans EPosTrans;
00091 typedef typename CNT<E>::TSqHermT ESqHermT;
00092 typedef typename CNT<E>::TSqTHerm ESqTHerm;
00093
00094 typedef typename CNT<E>::TSqrt ESqrt;
00095 typedef typename CNT<E>::TAbs EAbs;
00096 typedef typename CNT<E>::TStandard EStandard;
00097 typedef typename CNT<E>::TInvert EInvert;
00098 typedef typename CNT<E>::TNormalize ENormalize;
00099
00100 typedef typename CNT<E>::Scalar EScalar;
00101 typedef typename CNT<E>::ULessScalar EULessScalar;
00102 typedef typename CNT<E>::Number ENumber;
00103 typedef typename CNT<E>::StdNumber EStdNumber;
00104 typedef typename CNT<E>::Precision EPrecision;
00105 typedef typename CNT<E>::ScalarNormSq EScalarNormSq;
00106
00107 public:
00108
00109 enum {
00110 NRows = 1,
00111 NCols = N,
00112 NPackedElements = N,
00113 NActualElements = N * STRIDE,
00114 NActualScalars = CNT<E>::NActualScalars * NActualElements,
00115 RowSpacing = NActualElements,
00116 ColSpacing = STRIDE,
00117 ImagOffset = NTraits<ENumber>::ImagOffset,
00118 RealStrideFactor = 1,
00119
00120 ArgDepth = ((int)CNT<E>::ArgDepth < (int)MAX_RESOLVED_DEPTH
00121 ? CNT<E>::ArgDepth + 1
00122 : MAX_RESOLVED_DEPTH),
00123 IsScalar = 0,
00124 IsULessScalar = 0,
00125 IsNumber = 0,
00126 IsStdNumber = 0,
00127 IsPrecision = 0,
00128 SignInterpretation = CNT<E>::SignInterpretation
00129 };
00130
00131 typedef Row<N,E,STRIDE> T;
00132 typedef Row<N,ENeg,STRIDE> TNeg;
00133 typedef Row<N,EWithoutNegator,STRIDE> TWithoutNegator;
00134
00135 typedef Row<N,EReal,STRIDE*CNT<E>::RealStrideFactor>
00136 TReal;
00137 typedef Row<N,EImag,STRIDE*CNT<E>::RealStrideFactor>
00138 TImag;
00139 typedef Row<N,EComplex,STRIDE> TComplex;
00140 typedef Vec<N,EHerm,STRIDE> THerm;
00141 typedef Vec<N,E,STRIDE> TPosTrans;
00142 typedef E TElement;
00143 typedef Row TRow;
00144 typedef E TCol;
00145
00146
00147
00148 typedef Vec<N,ESqrt,1> TSqrt;
00149 typedef Row<N,EAbs,1> TAbs;
00150 typedef Row<N,EStandard,1> TStandard;
00151 typedef Vec<N,EInvert,1> TInvert;
00152 typedef Row<N,ENormalize,1> TNormalize;
00153
00154 typedef SymMat<N,ESqHermT> TSqHermT;
00155 typedef EScalarNormSq TSqTHerm;
00156
00157
00158
00159 typedef EScalar Scalar;
00160 typedef EULessScalar ULessScalar;
00161 typedef ENumber Number;
00162 typedef EStdNumber StdNumber;
00163 typedef EPrecision Precision;
00164 typedef EScalarNormSq ScalarNormSq;
00165
00166 int size() const { return N; }
00167 int nrow() const { return 1; }
00168 int ncol() const { return N; }
00169
00170
00171
00172 ScalarNormSq scalarNormSqr() const {
00173 ScalarNormSq sum(0);
00174 for(int i=0;i<N;++i) sum += CNT<E>::scalarNormSqr(d[i*STRIDE]);
00175 return sum;
00176 }
00177
00178
00179
00180
00181 TSqrt sqrt() const {
00182 TSqrt rsqrt;
00183 for(int i=0;i<N;++i) rsqrt[i] = CNT<E>::sqrt(d[i*STRIDE]);
00184 return rsqrt;
00185 }
00186
00187
00188
00189
00190 TAbs abs() const {
00191 TAbs rabs;
00192 for(int i=0;i<N;++i) rabs[i] = CNT<E>::abs(d[i*STRIDE]);
00193 return rabs;
00194 }
00195
00196 TStandard standardize() const {
00197 TStandard rstd;
00198 for(int i=0;i<N;++i) rstd[i] = CNT<E>::standardize(d[i*STRIDE]);
00199 return rstd;
00200 }
00201
00202
00203
00204 EStandard sum() const {
00205 E sum(0);
00206 for (int i=0;i<N;++i) sum += d[i*STRIDE];
00207 return CNT<E>::standardize(sum);
00208 }
00209
00210
00211
00212
00213 template <class P> struct EltResult {
00214 typedef Row<N, typename CNT<E>::template Result<P>::Mul, 1> Mul;
00215 typedef Row<N, typename CNT<E>::template Result<P>::Dvd, 1> Dvd;
00216 typedef Row<N, typename CNT<E>::template Result<P>::Add, 1> Add;
00217 typedef Row<N, typename CNT<E>::template Result<P>::Sub, 1> Sub;
00218 };
00219
00220
00221
00222 template <class P> struct Result {
00223 typedef MulCNTs<1,N,ArgDepth,Row,ColSpacing,RowSpacing,
00224 CNT<P>::NRows, CNT<P>::NCols, CNT<P>::ArgDepth,
00225 P, CNT<P>::ColSpacing, CNT<P>::RowSpacing> MulOp;
00226 typedef typename MulOp::Type Mul;
00227
00228 typedef MulCNTsNonConforming<1,N,ArgDepth,Row,ColSpacing,RowSpacing,
00229 CNT<P>::NRows, CNT<P>::NCols, CNT<P>::ArgDepth,
00230 P, CNT<P>::ColSpacing, CNT<P>::RowSpacing> MulOpNonConforming;
00231 typedef typename MulOpNonConforming::Type MulNon;
00232
00233
00234 typedef DvdCNTs<1,N,ArgDepth,Row,ColSpacing,RowSpacing,
00235 CNT<P>::NRows, CNT<P>::NCols, CNT<P>::ArgDepth,
00236 P, CNT<P>::ColSpacing, CNT<P>::RowSpacing> DvdOp;
00237 typedef typename DvdOp::Type Dvd;
00238
00239 typedef AddCNTs<1,N,ArgDepth,Row,ColSpacing,RowSpacing,
00240 CNT<P>::NRows, CNT<P>::NCols, CNT<P>::ArgDepth,
00241 P, CNT<P>::ColSpacing, CNT<P>::RowSpacing> AddOp;
00242 typedef typename AddOp::Type Add;
00243
00244 typedef SubCNTs<1,N,ArgDepth,Row,ColSpacing,RowSpacing,
00245 CNT<P>::NRows, CNT<P>::NCols, CNT<P>::ArgDepth,
00246 P, CNT<P>::ColSpacing, CNT<P>::RowSpacing> SubOp;
00247 typedef typename SubOp::Type Sub;
00248 };
00249
00250
00251 template <class P> struct Substitute {
00252 typedef Row<N,P> Type;
00253 };
00254
00255
00256
00257 Row(){
00258 #ifndef NDEBUG
00259 setToNaN();
00260 #endif
00261 }
00262
00263
00264
00265
00266 Row(const Row& src) {
00267 Impl::copy(*this, src);
00268 }
00269 Row& operator=(const Row& src) {
00270 Impl::copy(*this, src);
00271 return *this;
00272 }
00273
00274
00275
00276 template <int SS> Row(const Row<N,E,SS>& src) {
00277 Impl::copy(*this, src);
00278 }
00279
00280
00281
00282 template <int SS> Row(const Row<N,ENeg,SS>& src) {
00283 Impl::copy(*this, src);
00284 }
00285
00286
00287
00288 template <class EE, int SS> explicit Row(const Row<N,EE,SS>& vv) {
00289 Impl::copy(*this, vv);
00290 }
00291
00292
00293 explicit Row(const ELT& e)
00294 { for (int i=0;i<N;++i) d[i*STRIDE]=e; }
00295
00296
00297 Row(const E& e0,const E& e1)
00298 { assert(N==2);(*this)[0]=e0;(*this)[1]=e1; }
00299 Row(const E& e0,const E& e1,const E& e2)
00300 { assert(N==3);(*this)[0]=e0;(*this)[1]=e1;(*this)[2]=e2; }
00301 Row(const E& e0,const E& e1,const E& e2,const E& e3)
00302 { assert(N==4);(*this)[0]=e0;(*this)[1]=e1;(*this)[2]=e2;(*this)[3]=e3; }
00303 Row(const E& e0,const E& e1,const E& e2,const E& e3,const E& e4)
00304 { assert(N==5);(*this)[0]=e0;(*this)[1]=e1;(*this)[2]=e2;
00305 (*this)[3]=e3;(*this)[4]=e4; }
00306 Row(const E& e0,const E& e1,const E& e2,const E& e3,const E& e4,const E& e5)
00307 { assert(N==6);(*this)[0]=e0;(*this)[1]=e1;(*this)[2]=e2;
00308 (*this)[3]=e3;(*this)[4]=e4;(*this)[5]=e5; }
00309 Row(const E& e0,const E& e1,const E& e2,const E& e3,const E& e4,const E& e5,const E& e6)
00310 { assert(N==7);(*this)[0]=e0;(*this)[1]=e1;(*this)[2]=e2;
00311 (*this)[3]=e3;(*this)[4]=e4;(*this)[5]=e5;(*this)[6]=e6; }
00312 Row(const E& e0,const E& e1,const E& e2,const E& e3,const E& e4,const E& e5,const E& e6,const E& e7)
00313 { assert(N==8);(*this)[0]=e0;(*this)[1]=e1;(*this)[2]=e2;
00314 (*this)[3]=e3;(*this)[4]=e4;(*this)[5]=e5;(*this)[6]=e6;(*this)[7]=e7; }
00315 Row(const E& e0,const E& e1,const E& e2,const E& e3,const E& e4,const E& e5,const E& e6,const E& e7,const E& e8)
00316 { assert(N==9);(*this)[0]=e0;(*this)[1]=e1;(*this)[2]=e2;
00317 (*this)[3]=e3;(*this)[4]=e4;(*this)[5]=e5;(*this)[6]=e6;(*this)[7]=e7;(*this)[8]=e8; }
00318
00319
00320
00321 template <class EE> explicit Row(const EE* p)
00322 { assert(p); for(int i=0;i<N;++i) d[i*STRIDE]=p[i]; }
00323 template <class EE> Row& operator=(const EE* p)
00324 { assert(p); for(int i=0;i<N;++i) d[i*STRIDE]=p[i]; return *this; }
00325
00326
00327 template <class EE, int SS> Row& operator=(const Row<N,EE,SS>& vv) {
00328 Impl::copy(*this, vv);
00329 return *this;
00330 }
00331 template <class EE, int SS> Row& operator+=(const Row<N,EE,SS>& r)
00332 { for(int i=0;i<N;++i) d[i*STRIDE] += r[i]; return *this; }
00333 template <class EE, int SS> Row& operator+=(const Row<N,negator<EE>,SS>& r)
00334 { for(int i=0;i<N;++i) d[i*STRIDE] -= -(r[i]); return *this; }
00335 template <class EE, int SS> Row& operator-=(const Row<N,EE,SS>& r)
00336 { for(int i=0;i<N;++i) d[i*STRIDE] -= r[i]; return *this; }
00337 template <class EE, int SS> Row& operator-=(const Row<N,negator<EE>,SS>& r)
00338 { for(int i=0;i<N;++i) d[i*STRIDE] += -(r[i]); return *this; }
00339
00340
00341
00342 template <class EE, int SS> Row<N,typename CNT<E>::template Result<EE>::Add>
00343 conformingAdd(const Row<N,EE,SS>& r) const {
00344 Row<N,typename CNT<E>::template Result<EE>::Add> result;
00345 Impl::conformingAdd(*this, r, result);
00346 return result;
00347 }
00348 template <class EE, int SS> Row<N,typename CNT<E>::template Result<EE>::Sub>
00349 conformingSubtract(const Row<N,EE,SS>& r) const {
00350 Row<N,typename CNT<E>::template Result<EE>::Sub> result;
00351 Impl::conformingSubtract(*this, r, result);
00352 return result;
00353 }
00354
00355
00356 template <class EE, int SS> typename CNT<E>::template Result<EE>::Mul
00357 conformingMultiply(const Vec<N,EE,SS>& r) const {
00358 return (*this)*r;
00359 }
00360
00361
00362 template <int MatNCol, class EE, int CS, int RS>
00363 Row<MatNCol,typename CNT<E>::template Result<EE>::Mul>
00364 conformingMultiply(const Mat<N,MatNCol,EE,CS,RS>& m) const {
00365 Row<MatNCol,typename CNT<E>::template Result<EE>::Mul> result;
00366 for (int j=0;j<N;++j) result[j] = conformingMultiply(m(j));
00367 return result;
00368 }
00369
00370 const E& operator[](int i) const { assert(0 <= i && i < N); return d[i*STRIDE]; }
00371 E& operator[](int i) { assert(0 <= i && i < N); return d[i*STRIDE]; }
00372 const E& operator()(int i) const { return (*this)[i]; }
00373 E& operator()(int i) { return (*this)[i]; }
00374
00375 ScalarNormSq normSqr() const { return scalarNormSqr(); }
00376 typename CNT<ScalarNormSq>::TSqrt
00377 norm() const { return CNT<ScalarNormSq>::sqrt(scalarNormSqr()); }
00378
00379
00380
00381
00382
00383
00384
00385
00386
00387
00388
00389
00390 TNormalize normalize() const {
00391 if (CNT<E>::IsScalar) {
00392 return castAwayNegatorIfAny() / (SignInterpretation*norm());
00393 } else {
00394 TNormalize elementwiseNormalized;
00395 for (int j=0; j<N; ++j)
00396 elementwiseNormalized[j] = CNT<E>::normalize((*this)[j]);
00397 return elementwiseNormalized;
00398 }
00399 }
00400
00401 TInvert invert() const {assert(false); return TInvert();}
00402
00403 const Row& operator+() const { return *this; }
00404 const TNeg& operator-() const { return negate(); }
00405 TNeg& operator-() { return updNegate(); }
00406 const THerm& operator~() const { return transpose(); }
00407 THerm& operator~() { return updTranspose(); }
00408
00409 const TNeg& negate() const { return *reinterpret_cast<const TNeg*>(this); }
00410 TNeg& updNegate() { return *reinterpret_cast<TNeg*>(this); }
00411
00412 const THerm& transpose() const { return *reinterpret_cast<const THerm*>(this); }
00413 THerm& updTranspose() { return *reinterpret_cast<THerm*>(this); }
00414
00415 const TPosTrans& positionalTranspose() const
00416 { return *reinterpret_cast<const TPosTrans*>(this); }
00417 TPosTrans& updPositionalTranspose()
00418 { return *reinterpret_cast<TPosTrans*>(this); }
00419
00420 const TReal& real() const { return *reinterpret_cast<const TReal*>(this); }
00421 TReal& real() { return *reinterpret_cast< TReal*>(this); }
00422
00423
00424 const TImag& imag() const {
00425 const int offs = ImagOffset;
00426 const EImag* p = reinterpret_cast<const EImag*>(this);
00427 return *reinterpret_cast<const TImag*>(p+offs);
00428 }
00429 TImag& imag() {
00430 const int offs = ImagOffset;
00431 EImag* p = reinterpret_cast<EImag*>(this);
00432 return *reinterpret_cast<TImag*>(p+offs);
00433 }
00434
00435 const TWithoutNegator& castAwayNegatorIfAny() const {return *reinterpret_cast<const TWithoutNegator*>(this);}
00436 TWithoutNegator& updCastAwayNegatorIfAny() {return *reinterpret_cast<TWithoutNegator*>(this);}
00437
00438
00439
00440
00441
00442
00443
00444
00445 template <class EE> Row<N, typename CNT<E>::template Result<EE>::Mul>
00446 scalarMultiply(const EE& e) const {
00447 Row<N, typename CNT<E>::template Result<EE>::Mul> result;
00448 for (int j=0; j<N; ++j) result[j] = (*this)[j] * e;
00449 return result;
00450 }
00451 template <class EE> Row<N, typename CNT<EE>::template Result<E>::Mul>
00452 scalarMultiplyFromLeft(const EE& e) const {
00453 Row<N, typename CNT<EE>::template Result<E>::Mul> result;
00454 for (int j=0; j<N; ++j) result[j] = e * (*this)[j];
00455 return result;
00456 }
00457
00458
00459
00460 template <class EE> Row<N, typename CNT<E>::template Result<EE>::Dvd>
00461 scalarDivide(const EE& e) const {
00462 Row<N, typename CNT<E>::template Result<EE>::Dvd> result;
00463 for (int j=0; j<N; ++j) result[j] = (*this)[j] / e;
00464 return result;
00465 }
00466 template <class EE> Row<N, typename CNT<EE>::template Result<E>::Dvd>
00467 scalarDivideFromLeft(const EE& e) const {
00468 Row<N, typename CNT<EE>::template Result<E>::Dvd> result;
00469 for (int j=0; j<N; ++j) result[j] = e / (*this)[j];
00470 return result;
00471 }
00472
00473 template <class EE> Row<N, typename CNT<E>::template Result<EE>::Add>
00474 scalarAdd(const EE& e) const {
00475 Row<N, typename CNT<E>::template Result<EE>::Add> result;
00476 for (int j=0; j<N; ++j) result[j] = (*this)[j] + e;
00477 return result;
00478 }
00479
00480
00481 template <class EE> Row<N, typename CNT<E>::template Result<EE>::Sub>
00482 scalarSubtract(const EE& e) const {
00483 Row<N, typename CNT<E>::template Result<EE>::Sub> result;
00484 for (int j=0; j<N; ++j) result[j] = (*this)[j] - e;
00485 return result;
00486 }
00487 template <class EE> Row<N, typename CNT<EE>::template Result<E>::Sub>
00488 scalarSubtractFromLeft(const EE& e) const {
00489 Row<N, typename CNT<EE>::template Result<E>::Sub> result;
00490 for (int j=0; j<N; ++j) result[j] = e - (*this)[j];
00491 return result;
00492 }
00493
00494
00495
00496
00497 template <class EE> Row& operator =(const EE& e) {return scalarEq(e);}
00498 template <class EE> Row& operator+=(const EE& e) {return scalarPlusEq(e);}
00499 template <class EE> Row& operator-=(const EE& e) {return scalarMinusEq(e);}
00500 template <class EE> Row& operator*=(const EE& e) {return scalarTimesEq(e);}
00501 template <class EE> Row& operator/=(const EE& e) {return scalarDivideEq(e);}
00502
00503
00504
00505 template <class EE> Row& scalarEq(const EE& ee)
00506 { for(int i=0;i<N;++i) d[i*STRIDE] = ee; return *this; }
00507 template <class EE> Row& scalarPlusEq(const EE& ee)
00508 { for(int i=0;i<N;++i) d[i*STRIDE] += ee; return *this; }
00509 template <class EE> Row& scalarMinusEq(const EE& ee)
00510 { for(int i=0;i<N;++i) d[i*STRIDE] -= ee; return *this; }
00511 template <class EE> Row& scalarInverseMinusEq(const EE& ee)
00512 { for(int i=0;i<N;++i) d[i*STRIDE] = ee - d[i*STRIDE]; return *this; }
00513 template <class EE> Row& scalarTimesEq(const EE& ee)
00514 { for(int i=0;i<N;++i) d[i*STRIDE] *= ee; return *this; }
00515 template <class EE> Row& scalarDivideEq(const EE& ee)
00516 { for(int i=0;i<N;++i) d[i*STRIDE] /= ee; return *this; }
00517 template <class EE> Row& scalarInverseDivideEq(const EE& ee)
00518 { for(int i=0;i<N;++i) d[i*STRIDE] = ee / d[i*STRIDE]; return *this; }
00519
00520 void setToNaN() {
00521 (*this) = CNT<ELT>::getNaN();
00522 }
00523
00524 void setToZero() {
00525 (*this) = ELT(0);
00526 }
00527
00528
00529
00530 template <int NN>
00531 const Row<NN,ELT,STRIDE>& getSubRow(int j) const {
00532 assert(0 <= j && j + NN <= N);
00533 return Row<NN,ELT,STRIDE>::getAs(&(*this)[j]);
00534 }
00535 template <int NN>
00536 Row<NN,ELT,STRIDE>& updSubRow(int j) {
00537 assert(0 <= j && j + NN <= N);
00538 return Row<NN,ELT,STRIDE>::updAs(&(*this)[j]);
00539 }
00540
00541
00542
00543
00544 Row<N-1,ELT,1> drop1(int p) const {
00545 assert(0 <= p && p < N);
00546 Row<N-1,ELT,1> out;
00547 int nxt=0;
00548 for (int i=0; i<N-1; ++i, ++nxt) {
00549 if (nxt==p) ++nxt;
00550 out[i] = (*this)[nxt];
00551 }
00552 return out;
00553 }
00554
00555
00556
00557
00558 template <class EE> Row<N+1,ELT,1> append1(const EE& v) const {
00559 Row<N+1,ELT,1> out;
00560 Row<N,ELT,1>::updAs(&out[0]) = (*this);
00561 out[N] = v;
00562 return out;
00563 }
00564
00565
00566
00567
00568
00569
00570
00571 template <class EE> Row<N+1,ELT,1> insert1(int p, const EE& v) const {
00572 assert(0 <= p && p <= N);
00573 if (p==N) return append1(v);
00574 Row<N+1,ELT,1> out;
00575 int nxt=0;
00576 for (int i=0; i<N; ++i, ++nxt) {
00577 if (i==p) out[nxt++] = v;
00578 out[nxt] = (*this)[i];
00579 }
00580 return out;
00581 }
00582
00583
00584 static const Row& getAs(const ELT* p) {return *reinterpret_cast<const Row*>(p);}
00585 static Row& updAs(ELT* p) {return *reinterpret_cast<Row*>(p);}
00586
00587
00588 template <int NN>
00589 static const Row& getSubRow(const Row<NN,ELT,STRIDE>& r, int j) {
00590 assert(0 <= j && j + N <= NN);
00591 return getAs(&r[j]);
00592 }
00593 template <int NN>
00594 static Row& updSubRow(Row<NN,ELT,STRIDE>& r, int j) {
00595 assert(0 <= j && j + N <= NN);
00596 return updAs(&r[j]);
00597 }
00598
00599 static Row<N,ELT,1> getNaN() { return Row<N,ELT,1>(CNT<ELT>::getNaN()); }
00600
00602 bool isNaN() const {
00603 for (int j=0; j<N; ++j)
00604 if (CNT<ELT>::isNaN((*this)[j]))
00605 return true;
00606 return false;
00607 }
00608
00611 bool isInf() const {
00612 bool seenInf = false;
00613 for (int j=0; j<N; ++j) {
00614 const ELT& e = (*this)[j];
00615 if (!CNT<ELT>::isFinite(e)) {
00616 if (!CNT<ELT>::isInf(e))
00617 return false;
00618 seenInf = true;
00619 }
00620 }
00621 return seenInf;
00622 }
00623
00625 bool isFinite() const {
00626 for (int j=0; j<N; ++j)
00627 if (!CNT<ELT>::isFinite((*this)[j]))
00628 return false;
00629 return true;
00630 }
00631
00634 static double getDefaultTolerance() {return CNT<ELT>::getDefaultTolerance();}
00635
00638 template <class E2, int CS2>
00639 bool isNumericallyEqual(const Row<N,E2,CS2>& r, double tol) const {
00640 for (int j=0; j<N; ++j)
00641 if (!CNT<ELT>::isNumericallyEqual((*this)(j), r(j), tol))
00642 return false;
00643 return true;
00644 }
00645
00649 template <class E2, int CS2>
00650 bool isNumericallyEqual(const Row<N,E2,CS2>& r) const {
00651 const double tol = std::max(getDefaultTolerance(),r.getDefaultTolerance());
00652 return isNumericallyEqual(r, tol);
00653 }
00654
00659 bool isNumericallyEqual
00660 (const ELT& e,
00661 double tol = getDefaultTolerance()) const
00662 {
00663 for (int j=0; j<N; ++j)
00664 if (!CNT<ELT>::isNumericallyEqual((*this)(j), e, tol))
00665 return false;
00666 return true;
00667 }
00668 private:
00669 ELT d[NActualElements];
00670 };
00671
00673
00674
00676
00677
00678 template <int N, class E1, int S1, class E2, int S2> inline
00679 typename Row<N,E1,S1>::template Result< Row<N,E2,S2> >::Add
00680 operator+(const Row<N,E1,S1>& l, const Row<N,E2,S2>& r) {
00681 return Row<N,E1,S1>::template Result< Row<N,E2,S2> >
00682 ::AddOp::perform(l,r);
00683 }
00684
00685
00686 template <int N, class E1, int S1, class E2, int S2> inline
00687 typename Row<N,E1,S1>::template Result< Row<N,E2,S2> >::Sub
00688 operator-(const Row<N,E1,S1>& l, const Row<N,E2,S2>& r) {
00689 return Row<N,E1,S1>::template Result< Row<N,E2,S2> >
00690 ::SubOp::perform(l,r);
00691 }
00692
00693
00694 template <int N, class E1, int S1, class E2, int S2> inline bool
00695 operator==(const Row<N,E1,S1>& l, const Row<N,E2,S2>& r) {
00696 for (int i=0; i < N; ++i)
00697 if (l[i] != r[i]) return false;
00698 return true;
00699 }
00700
00701
00702 template <int N, class E1, int S1, class E2, int S2> inline bool
00703 operator!=(const Row<N,E1,S1>& l, const Row<N,E2,S2>& r) {return !(l==r);}
00704
00705
00707
00709
00710
00711
00712
00713
00714
00715
00716
00717 template <int N, class E, int S> inline
00718 typename Row<N,E,S>::template Result<float>::Mul
00719 operator*(const Row<N,E,S>& l, const float& r)
00720 { return Row<N,E,S>::template Result<float>::MulOp::perform(l,r); }
00721 template <int N, class E, int S> inline
00722 typename Row<N,E,S>::template Result<float>::Mul
00723 operator*(const float& l, const Row<N,E,S>& r) {return r*l;}
00724
00725 template <int N, class E, int S> inline
00726 typename Row<N,E,S>::template Result<double>::Mul
00727 operator*(const Row<N,E,S>& l, const double& r)
00728 { return Row<N,E,S>::template Result<double>::MulOp::perform(l,r); }
00729 template <int N, class E, int S> inline
00730 typename Row<N,E,S>::template Result<double>::Mul
00731 operator*(const double& l, const Row<N,E,S>& r) {return r*l;}
00732
00733 template <int N, class E, int S> inline
00734 typename Row<N,E,S>::template Result<long double>::Mul
00735 operator*(const Row<N,E,S>& l, const long double& r)
00736 { return Row<N,E,S>::template Result<long double>::MulOp::perform(l,r); }
00737 template <int N, class E, int S> inline
00738 typename Row<N,E,S>::template Result<long double>::Mul
00739 operator*(const long double& l, const Row<N,E,S>& r) {return r*l;}
00740
00741
00742 template <int N, class E, int S> inline
00743 typename Row<N,E,S>::template Result<typename CNT<E>::Precision>::Mul
00744 operator*(const Row<N,E,S>& l, int r) {return l * (typename CNT<E>::Precision)r;}
00745 template <int N, class E, int S> inline
00746 typename Row<N,E,S>::template Result<typename CNT<E>::Precision>::Mul
00747 operator*(int l, const Row<N,E,S>& r) {return r * (typename CNT<E>::Precision)l;}
00748
00749
00750
00751
00752 template <int N, class E, int S, class R> inline
00753 typename Row<N,E,S>::template Result<std::complex<R> >::Mul
00754 operator*(const Row<N,E,S>& l, const std::complex<R>& r)
00755 { return Row<N,E,S>::template Result<std::complex<R> >::MulOp::perform(l,r); }
00756 template <int N, class E, int S, class R> inline
00757 typename Row<N,E,S>::template Result<std::complex<R> >::Mul
00758 operator*(const std::complex<R>& l, const Row<N,E,S>& r) {return r*l;}
00759
00760
00761 template <int N, class E, int S, class R> inline
00762 typename Row<N,E,S>::template Result<std::complex<R> >::Mul
00763 operator*(const Row<N,E,S>& l, const conjugate<R>& r) {return l*(std::complex<R>)r;}
00764 template <int N, class E, int S, class R> inline
00765 typename Row<N,E,S>::template Result<std::complex<R> >::Mul
00766 operator*(const conjugate<R>& l, const Row<N,E,S>& r) {return r*(std::complex<R>)l;}
00767
00768
00769 template <int N, class E, int S, class R> inline
00770 typename Row<N,E,S>::template Result<typename negator<R>::StdNumber>::Mul
00771 operator*(const Row<N,E,S>& l, const negator<R>& r) {return l * (typename negator<R>::StdNumber)(R)r;}
00772 template <int N, class E, int S, class R> inline
00773 typename Row<N,E,S>::template Result<typename negator<R>::StdNumber>::Mul
00774 operator*(const negator<R>& l, const Row<N,E,S>& r) {return r * (typename negator<R>::StdNumber)(R)l;}
00775
00776
00777
00778
00779
00780
00781 template <int N, class E, int S> inline
00782 typename Row<N,E,S>::template Result<float>::Dvd
00783 operator/(const Row<N,E,S>& l, const float& r)
00784 { return Row<N,E,S>::template Result<float>::DvdOp::perform(l,r); }
00785 template <int N, class E, int S> inline
00786 typename CNT<float>::template Result<Row<N,E,S> >::Dvd
00787 operator/(const float& l, const Row<N,E,S>& r)
00788 { return CNT<float>::template Result<Row<N,E,S> >::DvdOp::perform(l,r); }
00789
00790 template <int N, class E, int S> inline
00791 typename Row<N,E,S>::template Result<double>::Dvd
00792 operator/(const Row<N,E,S>& l, const double& r)
00793 { return Row<N,E,S>::template Result<double>::DvdOp::perform(l,r); }
00794 template <int N, class E, int S> inline
00795 typename CNT<double>::template Result<Row<N,E,S> >::Dvd
00796 operator/(const double& l, const Row<N,E,S>& r)
00797 { return CNT<double>::template Result<Row<N,E,S> >::DvdOp::perform(l,r); }
00798
00799 template <int N, class E, int S> inline
00800 typename Row<N,E,S>::template Result<long double>::Dvd
00801 operator/(const Row<N,E,S>& l, const long double& r)
00802 { return Row<N,E,S>::template Result<long double>::DvdOp::perform(l,r); }
00803 template <int N, class E, int S> inline
00804 typename CNT<long double>::template Result<Row<N,E,S> >::Dvd
00805 operator/(const long double& l, const Row<N,E,S>& r)
00806 { return CNT<long double>::template Result<Row<N,E,S> >::DvdOp::perform(l,r); }
00807
00808
00809 template <int N, class E, int S> inline
00810 typename Row<N,E,S>::template Result<typename CNT<E>::Precision>::Dvd
00811 operator/(const Row<N,E,S>& l, int r) {return l / (typename CNT<E>::Precision)r;}
00812 template <int N, class E, int S> inline
00813 typename CNT<typename CNT<E>::Precision>::template Result<Row<N,E,S> >::Dvd
00814 operator/(int l, const Row<N,E,S>& r) {return (typename CNT<E>::Precision)l / r;}
00815
00816
00817
00818
00819
00820 template <int N, class E, int S, class R> inline
00821 typename Row<N,E,S>::template Result<std::complex<R> >::Dvd
00822 operator/(const Row<N,E,S>& l, const std::complex<R>& r)
00823 { return Row<N,E,S>::template Result<std::complex<R> >::DvdOp::perform(l,r); }
00824 template <int N, class E, int S, class R> inline
00825 typename CNT<std::complex<R> >::template Result<Row<N,E,S> >::Dvd
00826 operator/(const std::complex<R>& l, const Row<N,E,S>& r)
00827 { return CNT<std::complex<R> >::template Result<Row<N,E,S> >::DvdOp::perform(l,r); }
00828
00829
00830 template <int N, class E, int S, class R> inline
00831 typename Row<N,E,S>::template Result<std::complex<R> >::Dvd
00832 operator/(const Row<N,E,S>& l, const conjugate<R>& r) {return l/(std::complex<R>)r;}
00833 template <int N, class E, int S, class R> inline
00834 typename CNT<std::complex<R> >::template Result<Row<N,E,S> >::Dvd
00835 operator/(const conjugate<R>& l, const Row<N,E,S>& r) {return (std::complex<R>)l/r;}
00836
00837
00838 template <int N, class E, int S, class R> inline
00839 typename Row<N,E,S>::template Result<typename negator<R>::StdNumber>::Dvd
00840 operator/(const Row<N,E,S>& l, const negator<R>& r) {return l/(typename negator<R>::StdNumber)(R)r;}
00841 template <int N, class E, int S, class R> inline
00842 typename CNT<R>::template Result<Row<N,E,S> >::Dvd
00843 operator/(const negator<R>& l, const Row<N,E,S>& r) {return (typename negator<R>::StdNumber)(R)l/r;}
00844
00845
00846
00847
00848
00849
00850
00851
00852
00853 template <int N, class E, int S> inline
00854 typename Row<N,E,S>::template Result<float>::Add
00855 operator+(const Row<N,E,S>& l, const float& r)
00856 { return Row<N,E,S>::template Result<float>::AddOp::perform(l,r); }
00857 template <int N, class E, int S> inline
00858 typename Row<N,E,S>::template Result<float>::Add
00859 operator+(const float& l, const Row<N,E,S>& r) {return r+l;}
00860
00861 template <int N, class E, int S> inline
00862 typename Row<N,E,S>::template Result<double>::Add
00863 operator+(const Row<N,E,S>& l, const double& r)
00864 { return Row<N,E,S>::template Result<double>::AddOp::perform(l,r); }
00865 template <int N, class E, int S> inline
00866 typename Row<N,E,S>::template Result<double>::Add
00867 operator+(const double& l, const Row<N,E,S>& r) {return r+l;}
00868
00869 template <int N, class E, int S> inline
00870 typename Row<N,E,S>::template Result<long double>::Add
00871 operator+(const Row<N,E,S>& l, const long double& r)
00872 { return Row<N,E,S>::template Result<long double>::AddOp::perform(l,r); }
00873 template <int N, class E, int S> inline
00874 typename Row<N,E,S>::template Result<long double>::Add
00875 operator+(const long double& l, const Row<N,E,S>& r) {return r+l;}
00876
00877
00878 template <int N, class E, int S> inline
00879 typename Row<N,E,S>::template Result<typename CNT<E>::Precision>::Add
00880 operator+(const Row<N,E,S>& l, int r) {return l + (typename CNT<E>::Precision)r;}
00881 template <int N, class E, int S> inline
00882 typename Row<N,E,S>::template Result<typename CNT<E>::Precision>::Add
00883 operator+(int l, const Row<N,E,S>& r) {return r + (typename CNT<E>::Precision)l;}
00884
00885
00886
00887
00888 template <int N, class E, int S, class R> inline
00889 typename Row<N,E,S>::template Result<std::complex<R> >::Add
00890 operator+(const Row<N,E,S>& l, const std::complex<R>& r)
00891 { return Row<N,E,S>::template Result<std::complex<R> >::AddOp::perform(l,r); }
00892 template <int N, class E, int S, class R> inline
00893 typename Row<N,E,S>::template Result<std::complex<R> >::Add
00894 operator+(const std::complex<R>& l, const Row<N,E,S>& r) {return r+l;}
00895
00896
00897 template <int N, class E, int S, class R> inline
00898 typename Row<N,E,S>::template Result<std::complex<R> >::Add
00899 operator+(const Row<N,E,S>& l, const conjugate<R>& r) {return l+(std::complex<R>)r;}
00900 template <int N, class E, int S, class R> inline
00901 typename Row<N,E,S>::template Result<std::complex<R> >::Add
00902 operator+(const conjugate<R>& l, const Row<N,E,S>& r) {return r+(std::complex<R>)l;}
00903
00904
00905 template <int N, class E, int S, class R> inline
00906 typename Row<N,E,S>::template Result<typename negator<R>::StdNumber>::Add
00907 operator+(const Row<N,E,S>& l, const negator<R>& r) {return l + (typename negator<R>::StdNumber)(R)r;}
00908 template <int N, class E, int S, class R> inline
00909 typename Row<N,E,S>::template Result<typename negator<R>::StdNumber>::Add
00910 operator+(const negator<R>& l, const Row<N,E,S>& r) {return r + (typename negator<R>::StdNumber)(R)l;}
00911
00912
00913
00914
00915 template <int N, class E, int S> inline
00916 typename Row<N,E,S>::template Result<float>::Sub
00917 operator-(const Row<N,E,S>& l, const float& r)
00918 { return Row<N,E,S>::template Result<float>::SubOp::perform(l,r); }
00919 template <int N, class E, int S> inline
00920 typename CNT<float>::template Result<Row<N,E,S> >::Sub
00921 operator-(const float& l, const Row<N,E,S>& r)
00922 { return CNT<float>::template Result<Row<N,E,S> >::SubOp::perform(l,r); }
00923
00924 template <int N, class E, int S> inline
00925 typename Row<N,E,S>::template Result<double>::Sub
00926 operator-(const Row<N,E,S>& l, const double& r)
00927 { return Row<N,E,S>::template Result<double>::SubOp::perform(l,r); }
00928 template <int N, class E, int S> inline
00929 typename CNT<double>::template Result<Row<N,E,S> >::Sub
00930 operator-(const double& l, const Row<N,E,S>& r)
00931 { return CNT<double>::template Result<Row<N,E,S> >::SubOp::perform(l,r); }
00932
00933 template <int N, class E, int S> inline
00934 typename Row<N,E,S>::template Result<long double>::Sub
00935 operator-(const Row<N,E,S>& l, const long double& r)
00936 { return Row<N,E,S>::template Result<long double>::SubOp::perform(l,r); }
00937 template <int N, class E, int S> inline
00938 typename CNT<long double>::template Result<Row<N,E,S> >::Sub
00939 operator-(const long double& l, const Row<N,E,S>& r)
00940 { return CNT<long double>::template Result<Row<N,E,S> >::SubOp::perform(l,r); }
00941
00942
00943 template <int N, class E, int S> inline
00944 typename Row<N,E,S>::template Result<typename CNT<E>::Precision>::Sub
00945 operator-(const Row<N,E,S>& l, int r) {return l - (typename CNT<E>::Precision)r;}
00946 template <int N, class E, int S> inline
00947 typename CNT<typename CNT<E>::Precision>::template Result<Row<N,E,S> >::Sub
00948 operator-(int l, const Row<N,E,S>& r) {return (typename CNT<E>::Precision)l - r;}
00949
00950
00951
00952
00953
00954 template <int N, class E, int S, class R> inline
00955 typename Row<N,E,S>::template Result<std::complex<R> >::Sub
00956 operator-(const Row<N,E,S>& l, const std::complex<R>& r)
00957 { return Row<N,E,S>::template Result<std::complex<R> >::SubOp::perform(l,r); }
00958 template <int N, class E, int S, class R> inline
00959 typename CNT<std::complex<R> >::template Result<Row<N,E,S> >::Sub
00960 operator-(const std::complex<R>& l, const Row<N,E,S>& r)
00961 { return CNT<std::complex<R> >::template Result<Row<N,E,S> >::SubOp::perform(l,r); }
00962
00963
00964 template <int N, class E, int S, class R> inline
00965 typename Row<N,E,S>::template Result<std::complex<R> >::Sub
00966 operator-(const Row<N,E,S>& l, const conjugate<R>& r) {return l-(std::complex<R>)r;}
00967 template <int N, class E, int S, class R> inline
00968 typename CNT<std::complex<R> >::template Result<Row<N,E,S> >::Sub
00969 operator-(const conjugate<R>& l, const Row<N,E,S>& r) {return (std::complex<R>)l-r;}
00970
00971
00972 template <int N, class E, int S, class R> inline
00973 typename Row<N,E,S>::template Result<typename negator<R>::StdNumber>::Sub
00974 operator-(const Row<N,E,S>& l, const negator<R>& r) {return l-(typename negator<R>::StdNumber)(R)r;}
00975 template <int N, class E, int S, class R> inline
00976 typename CNT<R>::template Result<Row<N,E,S> >::Sub
00977 operator-(const negator<R>& l, const Row<N,E,S>& r) {return (typename negator<R>::StdNumber)(R)l-r;}
00978
00979
00980
00981 template <int N, class E, int S, class CHAR, class TRAITS> inline
00982 std::basic_ostream<CHAR,TRAITS>&
00983 operator<<(std::basic_ostream<CHAR,TRAITS>& o, const Row<N,E,S>& v) {
00984 o << "[" << v[0]; for(int i=1;i<N;++i) o<<','<<v[i]; o<<']'; return o;
00985 }
00986
00987 template <int N, class E, int S, class CHAR, class TRAITS> inline
00988 std::basic_istream<CHAR,TRAITS>&
00989 operator>>(std::basic_istream<CHAR,TRAITS>& is, Row<N,E,S>& v) {
00990
00991 assert(false);
00992 return is;
00993 }
00994
00995 }
00996
00997
00998 #endif //SimTK_SIMMATRIX_SMALLMATRIX_ROW_H_