Simbody  3.6
Mat.h
Go to the documentation of this file.
1 #ifndef SimTK_SIMMATRIX_SMALLMATRIX_MAT_H_
2 #define SimTK_SIMMATRIX_SMALLMATRIX_MAT_H_
3 
4 /* -------------------------------------------------------------------------- *
5  * Simbody(tm): SimTKcommon *
6  * -------------------------------------------------------------------------- *
7  * This is part of the SimTK biosimulation toolkit originating from *
8  * Simbios, the NIH National Center for Physics-Based Simulation of *
9  * Biological Structures at Stanford, funded under the NIH Roadmap for *
10  * Medical Research, grant U54 GM072970. See https://simtk.org/home/simbody. *
11  * *
12  * Portions copyright (c) 2005-12 Stanford University and the Authors. *
13  * Authors: Michael Sherman *
14  * Contributors: *
15  * *
16  * Licensed under the Apache License, Version 2.0 (the "License"); you may *
17  * not use this file except in compliance with the License. You may obtain a *
18  * copy of the License at http://www.apache.org/licenses/LICENSE-2.0. *
19  * *
20  * Unless required by applicable law or agreed to in writing, software *
21  * distributed under the License is distributed on an "AS IS" BASIS, *
22  * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. *
23  * See the License for the specific language governing permissions and *
24  * limitations under the License. *
25  * -------------------------------------------------------------------------- */
26 
32 
33 namespace SimTK {
34 
97 template <int M, int N, class ELT, int CS, int RS> class Mat {
98 public:
99  typedef ELT E;
100  typedef typename CNT<E>::TNeg ENeg;
102  typedef typename CNT<E>::TReal EReal;
103  typedef typename CNT<E>::TImag EImag;
104  typedef typename CNT<E>::TComplex EComplex;
105  typedef typename CNT<E>::THerm EHerm;
106  typedef typename CNT<E>::TPosTrans EPosTrans;
107  typedef typename CNT<E>::TSqHermT ESqHermT;
108  typedef typename CNT<E>::TSqTHerm ESqTHerm;
109 
110  typedef typename CNT<E>::TSqrt ESqrt;
111  typedef typename CNT<E>::TAbs EAbs;
112  typedef typename CNT<E>::TStandard EStandard;
113  typedef typename CNT<E>::TInvert EInvert;
114  typedef typename CNT<E>::TNormalize ENormalize;
115 
116  typedef typename CNT<E>::Scalar EScalar;
118  typedef typename CNT<E>::Number ENumber;
119  typedef typename CNT<E>::StdNumber EStdNumber;
120  typedef typename CNT<E>::Precision EPrecision;
122 
124  #ifndef SWIG
125  enum {
126  NRows = M,
127  NCols = N,
128  MinDim = N < M ? N : M,
129  MaxDim = N > M ? N : M,
133  NActualElements = (N-1)*CS + (M-1)*RS + 1,
136  RealStrideFactor = 1, // composite types don't change size when
137  // cast from complex to real or imaginary
139  ? CNT<E>::ArgDepth + 1
141  IsScalar = 0,
143  IsNumber = 0,
147  };
148  #endif
149 
153 
161  typedef E TElement;
162  typedef Row<N,E,CS> TRow;
163  typedef Vec<M,E,RS> TCol;
165 
166  // These are the results of calculations, so are returned in new, packed
167  // memory. Be sure to refer to element types here which are also packed.
168  typedef Mat<M,N,ESqrt,M,1> TSqrt; // Note strides are packed
169  typedef Mat<M,N,EAbs,M,1> TAbs; // Note strides are packed
171  typedef Mat<N,M,EInvert,N,1> TInvert; // like THerm but packed
173 
174  typedef SymMat<N,ESqHermT> TSqHermT; // ~Mat*Mat
175  typedef SymMat<M,ESqTHerm> TSqTHerm; // Mat*~Mat
176 
177  // Here the elements are copied unchanged but the result matrix
178  // is an ordinary packed, column order matrix.
180  typedef Mat<M-1,N,E,M-1,1> TDropRow;
181  typedef Mat<M,N-1,E,M,1> TDropCol;
182  typedef Mat<M-1,N-1,E,M-1,1> TDropRowCol;
186 
187  typedef EScalar Scalar;
188  typedef EULessScalar ULessScalar;
189  typedef ENumber Number;
190  typedef EStdNumber StdNumber;
191  typedef EPrecision Precision;
192  typedef EScalarNormSq ScalarNormSq;
193 
194  typedef THerm TransposeType; // TODO
195 
197  static int size() { return M*N; }
200  static int nrow() { return M; }
203  static int ncol() { return N; }
204 
210  ScalarNormSq scalarNormSqr() const {
211  ScalarNormSq sum(0);
212  for(int j=0;j<N;++j) sum += CNT<TCol>::scalarNormSqr((*this)(j));
213  return sum;
214  }
215 
219  TSqrt sqrt() const {
220  TSqrt msqrt;
221  for(int j=0;j<N;++j) msqrt(j) = (*this)(j).sqrt();
222  return msqrt;
223  }
224 
228  TAbs abs() const {
229  TAbs mabs;
230  for(int j=0;j<N;++j) mabs(j) = (*this)(j).abs();
231  return mabs;
232  }
233 
234  TStandard standardize() const {
235  TStandard mstd;
236  for(int j=0;j<N;++j) mstd(j) = (*this)(j).standardize();
237  return mstd;
238  }
239 
240  // This gives the resulting matrix type when (m(i,j) op P) is applied to each element.
241  // It is an MxN vector with default column and row spacing, and element types which
242  // are the regular composite result of E op P. Typically P is a scalar type but
243  // it doesn't have to be.
244  template <class P> struct EltResult {
249  };
250 
251  // This is the composite result for m op P where P is some kind of appropriately shaped
252  // non-scalar type.
253  template <class P> struct Result {
254  typedef MulCNTs<M,N,ArgDepth,Mat,ColSpacing,RowSpacing,
257  typedef typename MulOp::Type Mul;
258 
259  typedef MulCNTsNonConforming<M,N,ArgDepth,Mat,ColSpacing,RowSpacing,
260  CNT<P>::NRows, CNT<P>::NCols, CNT<P>::ArgDepth,
261  P, CNT<P>::ColSpacing, CNT<P>::RowSpacing> MulOpNonConforming;
262  typedef typename MulOpNonConforming::Type MulNon;
263 
264  typedef DvdCNTs<M,N,ArgDepth,Mat,ColSpacing,RowSpacing,
265  CNT<P>::NRows, CNT<P>::NCols, CNT<P>::ArgDepth,
266  P, CNT<P>::ColSpacing, CNT<P>::RowSpacing> DvdOp;
267  typedef typename DvdOp::Type Dvd;
268 
269  typedef AddCNTs<M,N,ArgDepth,Mat,ColSpacing,RowSpacing,
270  CNT<P>::NRows, CNT<P>::NCols, CNT<P>::ArgDepth,
271  P, CNT<P>::ColSpacing, CNT<P>::RowSpacing> AddOp;
272  typedef typename AddOp::Type Add;
273 
274  typedef SubCNTs<M,N,ArgDepth,Mat,ColSpacing,RowSpacing,
275  CNT<P>::NRows, CNT<P>::NCols, CNT<P>::ArgDepth,
276  P, CNT<P>::ColSpacing, CNT<P>::RowSpacing> SubOp;
277  typedef typename SubOp::Type Sub;
278  };
279 
280  // Shape-preserving element substitution (always packed)
281  template <class P> struct Substitute {
282  typedef Mat<M,N,P> Type;
283  };
284 
287  Mat(){
288  #ifndef NDEBUG
289  setToNaN();
290  #endif
291  }
292 
293  // It's important not to use the default copy constructor or copy
294  // assignment because the compiler doesn't understand that we may
295  // have noncontiguous storage and will try to copy the whole array.
296 
299  Mat(const Mat& src) {
300  for (int j=0; j<N; ++j)
301  (*this)(j) = src(j);
302  }
306  Mat& operator=(const Mat& src) {
307  for (int j=0; j<N; ++j)
308  (*this)(j) = src(j); // no harm if src and 'this' are the same
309  return *this;
310  }
311 
316  explicit Mat(const SymMat<M, ELT>& src) {
317  updDiag() = src.diag();
318  for (int j = 0; j < M; ++j)
319  for (int i = j+1; i < M; ++i) {
320  (*this)(i, j) = src.getEltLower(i, j);
321  (*this)(j, i) = src.getEltUpper(j, i);
322  }
323  }
324 
329  template <int CSS, int RSS>
330  Mat(const Mat<M,N,E,CSS,RSS>& src) {
331  for (int j=0; j<N; ++j)
332  (*this)(j) = src(j);
333  }
334 
340  template <int CSS, int RSS>
342  for (int j=0; j<N; ++j)
343  (*this)(j) = src(j);
344  }
345 
353  template <class EE, int CSS, int RSS>
354  explicit Mat(const Mat<M,N,EE,CSS,RSS>& mm)
355  { for (int j=0;j<N;++j) (*this)(j) = mm(j);}
356 
360  explicit Mat(const E& e)
361  { for (int j=0;j<N;++j) (*this)(j) = E(0); diag()=e; }
362 
367  explicit Mat(const ENeg& e)
368  { for (int j=0;j<N;++j) (*this)(j) = E(0); diag()=e; }
369 
376  explicit Mat(int i)
377  { new (this) Mat(E(Precision(i))); }
378 
379  // A bevy of constructors from individual exact-match elements IN ROW ORDER.
380  Mat(const E& e0,const E& e1)
381  {assert(M*N==2);d[rIx(0)]=e0;d[rIx(1)]=e1;}
382  Mat(const E& e0,const E& e1,const E& e2)
383  {assert(M*N==3);d[rIx(0)]=e0;d[rIx(1)]=e1;d[rIx(2)]=e2;}
384  Mat(const E& e0,const E& e1,const E& e2,const E& e3)
385  {assert(M*N==4);d[rIx(0)]=e0;d[rIx(1)]=e1;d[rIx(2)]=e2;d[rIx(3)]=e3;}
386  Mat(const E& e0,const E& e1,const E& e2,const E& e3,const E& e4)
387  {assert(M*N==5);d[rIx(0)]=e0;d[rIx(1)]=e1;d[rIx(2)]=e2;d[rIx(3)]=e3;d[rIx(4)]=e4;}
388  Mat(const E& e0,const E& e1,const E& e2,const E& e3,const E& e4,
389  const E& e5)
390  {assert(M*N==6);d[rIx(0)]=e0;d[rIx(1)]=e1;d[rIx(2)]=e2;d[rIx(3)]=e3;d[rIx(4)]=e4;
391  d[rIx(5)]=e5;}
392  Mat(const E& e0,const E& e1,const E& e2,const E& e3,const E& e4,
393  const E& e5,const E& e6)
394  {assert(M*N==7);d[rIx(0)]=e0;d[rIx(1)]=e1;d[rIx(2)]=e2;d[rIx(3)]=e3;d[rIx(4)]=e4;
395  d[rIx(5)]=e5;d[rIx(6)]=e6;}
396  Mat(const E& e0,const E& e1,const E& e2,const E& e3,const E& e4,
397  const E& e5,const E& e6,const E& e7)
398  {assert(M*N==8);d[rIx(0)]=e0;d[rIx(1)]=e1;d[rIx(2)]=e2;d[rIx(3)]=e3;d[rIx(4)]=e4;
399  d[rIx(5)]=e5;d[rIx(6)]=e6;d[rIx(7)]=e7;}
400  Mat(const E& e0,const E& e1,const E& e2,const E& e3,const E& e4,
401  const E& e5,const E& e6,const E& e7,const E& e8)
402  {assert(M*N==9);d[rIx(0)]=e0;d[rIx(1)]=e1;d[rIx(2)]=e2;d[rIx(3)]=e3;d[rIx(4)]=e4;
403  d[rIx(5)]=e5;d[rIx(6)]=e6;d[rIx(7)]=e7;d[rIx(8)]=e8;}
404  Mat(const E& e0,const E& e1,const E& e2,const E& e3,const E& e4,
405  const E& e5,const E& e6,const E& e7,const E& e8,const E& e9)
406  {assert(M*N==10);d[rIx(0)]=e0;d[rIx(1)]=e1;d[rIx(2)]=e2;d[rIx(3)]=e3;d[rIx(4)]=e4;
407  d[rIx(5)]=e5;d[rIx(6)]=e6;d[rIx(7)]=e7;d[rIx(8)]=e8;d[rIx(9)]=e9;}
408  Mat(const E& e0,const E& e1,const E& e2,const E& e3,const E& e4,
409  const E& e5,const E& e6,const E& e7,const E& e8,const E& e9,
410  const E& e10)
411  {assert(M*N==11);d[rIx(0)]=e0;d[rIx(1)]=e1;d[rIx(2)]=e2;d[rIx(3)]=e3;d[rIx(4)]=e4;
412  d[rIx(5)]=e5;d[rIx(6)]=e6;d[rIx(7)]=e7;d[rIx(8)]=e8;d[rIx(9)]=e9;d[rIx(10)]=e10;}
413  Mat(const E& e0,const E& e1,const E& e2,const E& e3,const E& e4,
414  const E& e5,const E& e6,const E& e7,const E& e8,const E& e9,
415  const E& e10, const E& e11)
416  {assert(M*N==12);d[rIx(0)]=e0;d[rIx(1)]=e1;d[rIx(2)]=e2;d[rIx(3)]=e3;d[rIx(4)]=e4;
417  d[rIx(5)]=e5;d[rIx(6)]=e6;d[rIx(7)]=e7;d[rIx(8)]=e8;d[rIx(9)]=e9;d[rIx(10)]=e10;
418  d[rIx(11)]=e11;}
419  Mat(const E& e0,const E& e1,const E& e2,const E& e3,const E& e4,
420  const E& e5,const E& e6,const E& e7,const E& e8,const E& e9,
421  const E& e10, const E& e11, const E& e12)
422  {assert(M*N==13);d[rIx(0)]=e0;d[rIx(1)]=e1;d[rIx(2)]=e2;d[rIx(3)]=e3;d[rIx(4)]=e4;
423  d[rIx(5)]=e5;d[rIx(6)]=e6;d[rIx(7)]=e7;d[rIx(8)]=e8;d[rIx(9)]=e9;d[rIx(10)]=e10;
424  d[rIx(11)]=e11;d[rIx(12)]=e12;}
425  Mat(const E& e0,const E& e1,const E& e2,const E& e3,const E& e4,
426  const E& e5,const E& e6,const E& e7,const E& e8,const E& e9,
427  const E& e10, const E& e11, const E& e12, const E& e13)
428  {assert(M*N==14);d[rIx(0)]=e0;d[rIx(1)]=e1;d[rIx(2)]=e2;d[rIx(3)]=e3;d[rIx(4)]=e4;
429  d[rIx(5)]=e5;d[rIx(6)]=e6;d[rIx(7)]=e7;d[rIx(8)]=e8;d[rIx(9)]=e9;d[rIx(10)]=e10;
430  d[rIx(11)]=e11;d[rIx(12)]=e12;d[rIx(13)]=e13;}
431  Mat(const E& e0,const E& e1,const E& e2,const E& e3,const E& e4,
432  const E& e5,const E& e6,const E& e7,const E& e8,const E& e9,
433  const E& e10, const E& e11, const E& e12, const E& e13, const E& e14)
434  {assert(M*N==15);d[rIx(0)]=e0;d[rIx(1)]=e1;d[rIx(2)]=e2;d[rIx(3)]=e3;d[rIx(4)]=e4;
435  d[rIx(5)]=e5;d[rIx(6)]=e6;d[rIx(7)]=e7;d[rIx(8)]=e8;d[rIx(9)]=e9;d[rIx(10)]=e10;
436  d[rIx(11)]=e11;d[rIx(12)]=e12;d[rIx(13)]=e13;d[rIx(14)]=e14;}
437  Mat(const E& e0,const E& e1,const E& e2,const E& e3,const E& e4,
438  const E& e5,const E& e6,const E& e7,const E& e8,const E& e9,
439  const E& e10, const E& e11, const E& e12, const E& e13, const E& e14,
440  const E& e15)
441  {assert(M*N==16);d[rIx(0)]=e0;d[rIx(1)]=e1;d[rIx(2)]=e2;d[rIx(3)]=e3;d[rIx(4)]=e4;
442  d[rIx(5)]=e5;d[rIx(6)]=e6;d[rIx(7)]=e7;d[rIx(8)]=e8;d[rIx(9)]=e9;d[rIx(10)]=e10;
443  d[rIx(11)]=e11;d[rIx(12)]=e12;d[rIx(13)]=e13;d[rIx(14)]=e14;d[rIx(15)]=e15;}
444 
445  // Construction from 1-6 *exact match* Rows
446  explicit Mat(const TRow& r0)
447  { assert(M==1); (*this)[0]=r0; }
448  Mat(const TRow& r0,const TRow& r1)
449  { assert(M==2);(*this)[0]=r0;(*this)[1]=r1; }
450  Mat(const TRow& r0,const TRow& r1,const TRow& r2)
451  { assert(M==3);(*this)[0]=r0;(*this)[1]=r1;(*this)[2]=r2; }
452  Mat(const TRow& r0,const TRow& r1,const TRow& r2,
453  const TRow& r3)
454  { assert(M==4);(*this)[0]=r0;(*this)[1]=r1;(*this)[2]=r2;(*this)[3]=r3; }
455  Mat(const TRow& r0,const TRow& r1,const TRow& r2,
456  const TRow& r3,const TRow& r4)
457  { assert(M==5);(*this)[0]=r0;(*this)[1]=r1;(*this)[2]=r2;
458  (*this)[3]=r3;(*this)[4]=r4; }
459  Mat(const TRow& r0,const TRow& r1,const TRow& r2,
460  const TRow& r3,const TRow& r4,const TRow& r5)
461  { assert(M==6);(*this)[0]=r0;(*this)[1]=r1;(*this)[2]=r2;
462  (*this)[3]=r3;(*this)[4]=r4;(*this)[5]=r5; }
463 
464  // Construction from 1-6 *compatible* Rows
465  template <class EE, int SS> explicit Mat(const Row<N,EE,SS>& r0)
466  { assert(M==1); (*this)[0]=r0; }
467  template <class EE, int SS> Mat(const Row<N,EE,SS>& r0,const Row<N,EE,SS>& r1)
468  { assert(M==2);(*this)[0]=r0;(*this)[1]=r1; }
469  template <class EE, int SS>
470  Mat(const Row<N,EE,SS>& r0,const Row<N,EE,SS>& r1,const Row<N,EE,SS>& r2)
471  { assert(M==3);(*this)[0]=r0;(*this)[1]=r1;(*this)[2]=r2; }
472  template <class EE, int SS>
473  Mat(const Row<N,EE,SS>& r0,const Row<N,EE,SS>& r1,const Row<N,EE,SS>& r2,
474  const Row<N,EE,SS>& r3)
475  { assert(M==4);(*this)[0]=r0;(*this)[1]=r1;(*this)[2]=r2;(*this)[3]=r3; }
476  template <class EE, int SS>
477  Mat(const Row<N,EE,SS>& r0,const Row<N,EE,SS>& r1,const Row<N,EE,SS>& r2,
478  const Row<N,EE,SS>& r3,const Row<N,EE,SS>& r4)
479  { assert(M==5);(*this)[0]=r0;(*this)[1]=r1;(*this)[2]=r2;
480  (*this)[3]=r3;(*this)[4]=r4; }
481  template <class EE, int SS>
482  Mat(const Row<N,EE,SS>& r0,const Row<N,EE,SS>& r1,const Row<N,EE,SS>& r2,
483  const Row<N,EE,SS>& r3,const Row<N,EE,SS>& r4,const Row<N,EE,SS>& r5)
484  { assert(M==6);(*this)[0]=r0;(*this)[1]=r1;(*this)[2]=r2;
485  (*this)[3]=r3;(*this)[4]=r4;(*this)[5]=r5; }
486 
487 
488  // Construction from 1-6 *exact match* Vecs
489  explicit Mat(const TCol& r0)
490  { assert(N==1); (*this)(0)=r0; }
491  Mat(const TCol& r0,const TCol& r1)
492  { assert(N==2);(*this)(0)=r0;(*this)(1)=r1; }
493  Mat(const TCol& r0,const TCol& r1,const TCol& r2)
494  { assert(N==3);(*this)(0)=r0;(*this)(1)=r1;(*this)(2)=r2; }
495  Mat(const TCol& r0,const TCol& r1,const TCol& r2,
496  const TCol& r3)
497  { assert(N==4);(*this)(0)=r0;(*this)(1)=r1;(*this)(2)=r2;(*this)(3)=r3; }
498  Mat(const TCol& r0,const TCol& r1,const TCol& r2,
499  const TCol& r3,const TCol& r4)
500  { assert(N==5);(*this)(0)=r0;(*this)(1)=r1;(*this)(2)=r2;
501  (*this)(3)=r3;(*this)(4)=r4; }
502  Mat(const TCol& r0,const TCol& r1,const TCol& r2,
503  const TCol& r3,const TCol& r4,const TCol& r5)
504  { assert(N==6);(*this)(0)=r0;(*this)(1)=r1;(*this)(2)=r2;
505  (*this)(3)=r3;(*this)(4)=r4;(*this)(5)=r5; }
506 
507  // Construction from 1-6 *compatible* Vecs
508  template <class EE, int SS> explicit Mat(const Vec<M,EE,SS>& r0)
509  { assert(N==1); (*this)(0)=r0; }
510  template <class EE, int SS> Mat(const Vec<M,EE,SS>& r0,const Vec<M,EE,SS>& r1)
511  { assert(N==2);(*this)(0)=r0;(*this)(1)=r1; }
512  template <class EE, int SS>
513  Mat(const Vec<M,EE,SS>& r0,const Vec<M,EE,SS>& r1,const Vec<M,EE,SS>& r2)
514  { assert(N==3);(*this)(0)=r0;(*this)(1)=r1;(*this)(2)=r2; }
515  template <class EE, int SS>
516  Mat(const Vec<M,EE,SS>& r0,const Vec<M,EE,SS>& r1,const Vec<M,EE,SS>& r2,
517  const Vec<M,EE,SS>& r3)
518  { assert(N==4);(*this)(0)=r0;(*this)(1)=r1;(*this)(2)=r2;(*this)(3)=r3; }
519  template <class EE, int SS>
520  Mat(const Vec<M,EE,SS>& r0,const Vec<M,EE,SS>& r1,const Vec<M,EE,SS>& r2,
521  const Vec<M,EE,SS>& r3,const Vec<M,EE,SS>& r4)
522  { assert(N==5);(*this)(0)=r0;(*this)(1)=r1;(*this)(2)=r2;
523  (*this)(3)=r3;(*this)(4)=r4; }
524  template <class EE, int SS>
525  Mat(const Vec<M,EE,SS>& r0,const Vec<M,EE,SS>& r1,const Vec<M,EE,SS>& r2,
526  const Vec<M,EE,SS>& r3,const Vec<M,EE,SS>& r4,const Vec<M,EE,SS>& r5)
527  { assert(N==6);(*this)(0)=r0;(*this)(1)=r1;(*this)(2)=r2;
528  (*this)(3)=r3;(*this)(4)=r4;(*this)(5)=r5; }
529 
530  // Construction from a pointer to anything assumes we're pointing
531  // at a packed element list of the right length, in row order.
532  template <class EE> explicit Mat(const EE* p)
533  { assert(p); for(int i=0;i<M;++i) (*this)[i]=&p[i*N]; }
534 
535  // Assignment works similarly to copy -- if the lengths match,
536  // go element-by-element. Otherwise, zero and then copy to each
537  // diagonal element.
538  template <class EE, int CSS, int RSS> Mat& operator=(const Mat<M,N,EE,CSS,RSS>& mm) {
539  for (int j=0; j<N; ++j) (*this)(j) = mm(j);
540  return *this;
541  }
542 
543  template <class EE> Mat& operator=(const EE* p) {
544  assert(p); for(int i=0;i<M;++i) (*this)[i]=&p[i*N];
545  return *this;
546  }
547 
548  // Assignment ops
549  template <class EE, int CSS, int RSS> Mat&
551  for (int j=0; j<N; ++j) (*this)(j) += mm(j);
552  return *this;
553  }
554  template <class EE, int CSS, int RSS> Mat&
555  operator+=(const Mat<M,N,negator<EE>,CSS,RSS>& mm) {
556  for (int j=0; j<N; ++j) (*this)(j) -= -(mm(j));
557  return *this;
558  }
559 
560  template <class EE, int CSS, int RSS> Mat&
562  for (int j=0; j<N; ++j) (*this)(j) -= mm(j);
563  return *this;
564  }
565  template <class EE, int CSS, int RSS> Mat&
566  operator-=(const Mat<M,N,negator<EE>,CSS,RSS>& mm) {
567  for (int j=0; j<N; ++j) (*this)(j) += -(mm(j));
568  return *this;
569  }
570 
571  // In place matrix multiply can only be done when the RHS matrix is square of dimension
572  // N (i.e., number of columns), and the elements are also *= compatible.
573  template <class EE, int CSS, int RSS> Mat&
575  const Mat t(*this);
576  for (int j=0; j<N; ++j)
577  for (int i=0; i<M; ++i)
578  (*this)(i,j) = t[i] * mm(j);
579  return *this;
580  }
581 
582  // Conforming binary ops with 'this' on left, producing new packed result.
583  // Cases: m=m+-m, m=m+-sy, m=m*m, m=m*sy, v=m*v
584 
585  // m= this + m
586  template <class E2, int CS2, int RS2>
587  typename Result<Mat<M,N,E2,CS2,RS2> >::Add
589  typename Result<Mat<M,N,E2,CS2,RS2> >::Add result;
590  for (int j=0;j<N;++j) result(j) = (*this)(j) + r(j);
591  return result;
592  }
593  // m= this - m
594  template <class E2, int CS2, int RS2>
595  typename Result<Mat<M,N,E2,CS2,RS2> >::Sub
597  typename Result<Mat<M,N,E2,CS2,RS2> >::Sub result;
598  for (int j=0;j<N;++j) result(j) = (*this)(j) - r(j);
599  return result;
600  }
601  // m= m - this
602  template <class E2, int CS2, int RS2>
605  return l.conformingSubtract(*this);
606  }
607 
608  // m= this .* m
609  template <class E2, int CS2, int RS2>
610  typename EltResult<E2>::Mul
612  typename EltResult<E2>::Mul result;
613  for (int j=0;j<N;++j)
614  result(j) = (*this)(j).elementwiseMultiply(r(j));
615  return result;
616  }
617 
618  // m= this ./ m
619  template <class E2, int CS2, int RS2>
620  typename EltResult<E2>::Dvd
622  typename EltResult<E2>::Dvd result;
623  for (int j=0;j<N;++j)
624  result(j) = (*this)(j).elementwiseDivide(r(j));
625  return result;
626  }
627 
628  // We always punt to the SymMat since it knows better what to do.
629  // m = this + sym
630  template <class E2, int RS2>
631  typename Result<SymMat<M,E2,RS2> >::Add
632  conformingAdd(const SymMat<M,E2,RS2>& sy) const {
633  assert(M==N);
634  return sy.conformingAdd(*this);
635  }
636  // m= this - sym
637  template <class E2, int RS2>
638  typename Result<SymMat<M,E2,RS2> >::Sub
640  assert(M==N);
641  return sy.conformingSubtractFromLeft(*this);
642  }
643  // m= sym - this
644  template <class E2, int RS2>
647  assert(M==N);
648  return sy.conformingSubtract(*this);
649  }
650 
651  // m= this * m
652  template <int N2, class E2, int CS2, int RS2>
653  typename Result<Mat<N,N2,E2,CS2,RS2> >::Mul
655  typename Result<Mat<N,N2,E2,CS2,RS2> >::Mul result;
656  for (int j=0;j<N2;++j)
657  for (int i=0;i<M;++i)
658  result(i,j) = (*this)[i].conformingMultiply(m(j)); // i.e., dot()
659  return result;
660  }
661  // m= m * this
662  template <int M2, class E2, int CS2, int RS2>
665  return m.conformingMultiply(*this);
666  }
667 
668  // m= this / m = this * m.invert()
669  template <int M2, class E2, int CS2, int RS2>
670  typename Result<Mat<M2,N,E2,CS2,RS2> >::Dvd
672  return conformingMultiply(m.invert());
673  }
674  // m= m / this = m * this.invert()
675  template <int M2, class E2, int CS2, int RS2>
678  return m.conformingMultiply((*this).invert());
679  }
680 
681  const TRow& operator[](int i) const { return row(i); }
682  TRow& operator[](int i) { return row(i); }
683  const TCol& operator()(int j) const { return col(j); }
684  TCol& operator()(int j) { return col(j); }
685 
686  const E& operator()(int i,int j) const { return elt(i,j); }
687  E& operator()(int i,int j) { return elt(i,j); }
688 
689  // This is the scalar Frobenius norm.
690  ScalarNormSq normSqr() const { return scalarNormSqr(); }
691  typename CNT<ScalarNormSq>::TSqrt
693 
694  // There is no conventional meaning for normalize() applied to a matrix. We
695  // choose to define it as follows:
696  // If the elements of this Mat are scalars, the result is what you get by
697  // dividing each element by the Frobenius norm() calculated above. If the elements are
698  // *not* scalars, then the elements are *separately* normalized. That means
699  // you will get a different answer from Mat<2,2,Mat33>::normalize() than you
700  // would from a Mat<6,6>::normalize() containing the same scalars.
701  //
702  // Normalize returns a matrix of the same dimension but in new, packed storage
703  // and with a return type that does not include negator<> even if the original
704  // Mat<> does, because we can eliminate the negation here almost for free.
705  // But we can't standardize (change conjugate to complex) for free, so we'll retain
706  // conjugates if there are any.
707  TNormalize normalize() const {
708  if (CNT<E>::IsScalar) {
710  } else {
711  TNormalize elementwiseNormalized;
712  // punt to the column Vec to deal with the elements
713  for (int j=0; j<N; ++j)
714  elementwiseNormalized(j) = (*this)(j).normalize();
715  return elementwiseNormalized;
716  }
717  }
718 
719  // Default inversion. Assume full rank if square, otherwise return
720  // pseudoinverse. (Mostly TODO)
721  TInvert invert() const;
722 
723  const Mat& operator+() const { return *this; }
724  const TNeg& operator-() const { return negate(); }
725  TNeg& operator-() { return updNegate(); }
726  const THerm& operator~() const { return transpose(); }
727  THerm& operator~() { return updTranspose(); }
728 
729  const TNeg& negate() const { return *reinterpret_cast<const TNeg*>(this); }
730  TNeg& updNegate() { return *reinterpret_cast<TNeg*>(this); }
731 
732  const THerm& transpose() const { return *reinterpret_cast<const THerm*>(this); }
733  THerm& updTranspose() { return *reinterpret_cast<THerm*>(this); }
734 
735  const TPosTrans& positionalTranspose() const
736  { return *reinterpret_cast<const TPosTrans*>(this); }
738  { return *reinterpret_cast<TPosTrans*>(this); }
739 
740  // If the underlying scalars are complex or conjugate, we can return a
741  // reference to the real part just by recasting the data to a matrix of
742  // the same dimensions but containing real elements, with the scalar
743  // spacing doubled.
744  const TReal& real() const { return *reinterpret_cast<const TReal*>(this); }
745  TReal& real() { return *reinterpret_cast< TReal*>(this); }
746 
747  // Getting the imaginary part is almost the same as real, but we have
748  // to shift the starting address of the returned object by 1 real-size
749  // ("precision") scalar so that the first element is the imaginary part
750  // of the original first element.
751  // TODO: should blow up or return a reference to a zero matrix if called
752  // on a real object.
753  // Had to contort these routines to get them through VC++ 7.net
754  const TImag& imag() const {
755  const int offs = ImagOffset;
756  const Precision* p = reinterpret_cast<const Precision*>(this);
757  return *reinterpret_cast<const TImag*>(p+offs);
758  }
759  TImag& imag() {
760  const int offs = ImagOffset;
761  Precision* p = reinterpret_cast<Precision*>(this);
762  return *reinterpret_cast<TImag*>(p+offs);
763  }
764 
765  const TWithoutNegator& castAwayNegatorIfAny() const {return *reinterpret_cast<const TWithoutNegator*>(this);}
766  TWithoutNegator& updCastAwayNegatorIfAny() {return *reinterpret_cast<TWithoutNegator*>(this);}
767 
768  const TRow& row(int i) const {
769  SimTK_INDEXCHECK(i,M, "Mat::row[i]");
770  return *reinterpret_cast<const TRow*>(&d[i*RS]);
771  }
772  TRow& row(int i) {
773  SimTK_INDEXCHECK(i,M, "Mat::row[i]");
774  return *reinterpret_cast<TRow*>(&d[i*RS]);
775  }
776 
777  const TCol& col(int j) const {
778  SimTK_INDEXCHECK(j,N, "Mat::col(j)");
779  return *reinterpret_cast<const TCol*>(&d[j*CS]);
780  }
781  TCol& col(int j) {
782  SimTK_INDEXCHECK(j,N, "Mat::col(j)");
783  return *reinterpret_cast<TCol*>(&d[j*CS]);
784  }
785 
786  const E& elt(int i, int j) const {
787  SimTK_INDEXCHECK(i,M, "Mat::elt(i,j)");
788  SimTK_INDEXCHECK(j,N, "Mat::elt(i,j)");
789  return d[i*RS+j*CS];
790  }
791  E& elt(int i, int j) {
792  SimTK_INDEXCHECK(i,M, "Mat::elt(i,j)");
793  SimTK_INDEXCHECK(j,N, "Mat::elt(i,j)");
794  return d[i*RS+j*CS];
795  }
796 
800  const TDiag& diag() const { return *reinterpret_cast<const TDiag*>(d); }
804  TDiag& updDiag() { return *reinterpret_cast<TDiag*>(d); }
807  TDiag& diag() { return *reinterpret_cast<TDiag*>(d); }
808 
809  EStandard trace() const {return diag().sum();}
810 
811  // These are elementwise binary operators, (this op ee) by default but (ee op this) if
812  // 'FromLeft' appears in the name. The result is a packed Mat<M,N> but the element type
813  // may change. These are mostly used to implement global operators.
814  // We call these "scalar" operators but actually the "scalar" can be a composite type.
815 
816  //TODO: consider converting 'e' to Standard Numbers as precalculation and changing
817  // return type appropriately.
819  scalarMultiply(const EE& e) const {
821  for (int j=0; j<N; ++j) result(j) = (*this)(j).scalarMultiply(e);
822  return result;
823  }
825  scalarMultiplyFromLeft(const EE& e) const {
827  for (int j=0; j<N; ++j) result(j) = (*this)(j).scalarMultiplyFromLeft(e);
828  return result;
829  }
830 
831  // TODO: should precalculate and store 1/e, while converting to Standard Numbers. Note
832  // that return type should change appropriately.
834  scalarDivide(const EE& e) const {
836  for (int j=0; j<N; ++j) result(j) = (*this)(j).scalarDivide(e);
837  return result;
838  }
840  scalarDivideFromLeft(const EE& e) const {
842  for (int j=0; j<N; ++j) result(j) = (*this)(j).scalarDivideFromLeft(e);
843  return result;
844  }
845 
846  // Additive operators for scalars operate only on the diagonal.
848  scalarAdd(const EE& e) const {
850  result.diag() += e;
851  return result;
852  }
853  // Add is commutative, so no 'FromLeft'.
854 
856  scalarSubtract(const EE& e) const {
858  result.diag() -= e;
859  return result;
860  }
861  // Should probably do something clever with negation here (s - m)
863  scalarSubtractFromLeft(const EE& e) const {
865  result.diag() += e; // yes, add
866  return result;
867  }
868 
869  // Generic assignments for any element type not listed explicitly, including scalars.
870  // These are done repeatedly for each element and only work if the operation can
871  // be performed leaving the original element type.
872  template <class EE> Mat& operator =(const EE& e) {return scalarEq(e);}
873  template <class EE> Mat& operator+=(const EE& e) {return scalarPlusEq(e);}
874  template <class EE> Mat& operator-=(const EE& e) {return scalarMinusEq(e);}
875  template <class EE> Mat& operator*=(const EE& e) {return scalarTimesEq(e);}
876  template <class EE> Mat& operator/=(const EE& e) {return scalarDivideEq(e);}
877 
878  // Generalized scalar assignment & computed assignment methods. These will work
879  // for any assignment-compatible element, not just scalars.
880  template <class EE> Mat& scalarEq(const EE& ee)
881  { for(int j=0; j<N; ++j) (*this)(j).scalarEq(EE(0));
882  diag().scalarEq(ee);
883  return *this; }
884 
885  template <class EE> Mat& scalarPlusEq(const EE& ee)
886  { diag().scalarPlusEq(ee); return *this; }
887 
888  template <class EE> Mat& scalarMinusEq(const EE& ee)
889  { diag().scalarMinusEq(ee); return *this; }
890  // m = s - m; negate m, then add s
891  template <class EE> Mat& scalarMinusEqFromLeft(const EE& ee)
892  { scalarTimesEq(E(-1)); diag().scalarAdd(ee); return *this; }
893 
894  template <class EE> Mat& scalarTimesEq(const EE& ee)
895  { for(int j=0; j<N; ++j) (*this)(j).scalarTimesEq(ee); return *this; }
896  template <class EE> Mat& scalarTimesEqFromLeft(const EE& ee)
897  { for(int j=0; j<N; ++j) (*this)(j).scalarTimesEqFromLeft(ee); return *this; }
898 
899  template <class EE> Mat& scalarDivideEq(const EE& ee)
900  { for(int j=0; j<N; ++j) (*this)(j).scalarDivideEq(ee); return *this; }
901  template <class EE> Mat& scalarDivideEqFromLeft(const EE& ee)
902  { for(int j=0; j<N; ++j) (*this)(j).scalarDivideEqFromLeft(ee); return *this; }
903 
904  void setToNaN() {
905  for (int j=0; j<N; ++j)
906  (*this)(j).setToNaN();
907  }
908 
909  void setToZero() {
910  for (int j=0; j<N; ++j)
911  (*this)(j).setToZero();
912  }
913 
914  // Extract a sub-Mat with size known at compile time. These have to be
915  // called with explicit template arguments, e.g. getSubMat<3,4>(i,j).
916 
917  template <int MM, int NN> struct SubMat {
919  };
920 
921  template <int MM, int NN>
922  const typename SubMat<MM,NN>::Type& getSubMat(int i, int j) const {
923  assert(0 <= i && i + MM <= M);
924  assert(0 <= j && j + NN <= N);
925  return SubMat<MM,NN>::Type::getAs(&(*this)(i,j));
926  }
927  template <int MM, int NN>
928  typename SubMat<MM,NN>::Type& updSubMat(int i, int j) {
929  assert(0 <= i && i + MM <= M);
930  assert(0 <= j && j + NN <= N);
931  return SubMat<MM,NN>::Type::updAs(&(*this)(i,j));
932  }
933  template <int MM, int NN>
934  void setSubMat(int i, int j, const typename SubMat<MM,NN>::Type& value) {
935  assert(0 <= i && i + MM <= M);
936  assert(0 <= j && j + NN <= N);
937  SubMat<MM,NN>::Type::updAs(&(*this)(i,j)) = value;
938  }
939 
942  TDropRow dropRow(int i) const {
943  assert(0 <= i && i < M);
944  TDropRow out;
945  for (int r=0, nxt=0; r<M-1; ++r, ++nxt) {
946  if (nxt==i) ++nxt; // skip the loser
947  out[r] = (*this)[nxt];
948  }
949  return out;
950  }
951 
954  TDropCol dropCol(int j) const {
955  assert(0 <= j && j < N);
956  TDropCol out;
957  for (int c=0, nxt=0; c<N-1; ++c, ++nxt) {
958  if (nxt==j) ++nxt; // skip the loser
959  out(c) = (*this)(nxt);
960  }
961  return out;
962  }
963 
967  TDropRowCol dropRowCol(int i, int j) const {
968  assert(0 <= i && i < M);
969  assert(0 <= j && j < N);
970  TDropRowCol out;
971  for (int c=0, nxtc=0; c<N-1; ++c, ++nxtc) {
972  if (nxtc==j) ++nxtc;
973  for (int r=0, nxtr=0; r<M-1; ++r, ++nxtr) {
974  if (nxtr==i) ++nxtr;
975  out(r,c) = (*this)(nxtr,nxtc);
976  }
977  }
978  return out;
979  }
980 
984  template <class EE, int SS>
985  TAppendRow appendRow(const Row<N,EE,SS>& row) const {
986  TAppendRow out;
987  out.template updSubMat<M,N>(0,0) = (*this);
988  out[M] = row;
989  return out;
990  }
991 
995  template <class EE, int SS>
996  TAppendCol appendCol(const Vec<M,EE,SS>& col) const {
997  TAppendCol out;
998  out.template updSubMat<M,N>(0,0) = (*this);
999  out(N) = col;
1000  return out;
1001  }
1002 
1009  template <class ER, int SR, class EC, int SC>
1010  TAppendRowCol appendRowCol(const Row<N+1,ER,SR>& row,
1011  const Vec<M+1,EC,SC>& col) const
1012  {
1013  TAppendRowCol out;
1014  out.template updSubMat<M,N>(0,0) = (*this);
1015  out[M].template updSubRow<N>(0) =
1016  row.template getSubRow<N>(0); // ignore last element
1017  out(N) = col;
1018  return out;
1019  }
1020 
1026  template <class EE, int SS>
1027  TAppendRow insertRow(int i, const Row<N,EE,SS>& row) const {
1028  assert(0 <= i && i <= M);
1029  if (i==M) return appendRow(row);
1030  TAppendRow out;
1031  for (int r=0, nxt=0; r<M; ++r, ++nxt) {
1032  if (nxt==i) out[nxt++] = row;
1033  out[nxt] = (*this)[r];
1034  }
1035  return out;
1036  }
1037 
1043  template <class EE, int SS>
1044  TAppendCol insertCol(int j, const Vec<M,EE,SS>& col) const {
1045  assert(0 <= j && j <= N);
1046  if (j==N) return appendCol(col);
1047  TAppendCol out;
1048  for (int c=0, nxt=0; c<N; ++c, ++nxt) {
1049  if (nxt==j) out(nxt++) = col;
1050  out(nxt) = (*this)(c);
1051  }
1052  return out;
1053  }
1054 
1062  template <class ER, int SR, class EC, int SC>
1063  TAppendRowCol insertRowCol(int i, int j, const Row<N+1,ER,SR>& row,
1064  const Vec<M+1,EC,SC>& col) const {
1065  assert(0 <= i && i <= M);
1066  assert(0 <= j && j <= N);
1067  TAppendRowCol out;
1068  for (int c=0, nxtc=0; c<N; ++c, ++nxtc) {
1069  if (nxtc==j) ++nxtc; // leave room
1070  for (int r=0, nxtr=0; r<M; ++r, ++nxtr) {
1071  if (nxtr==i) ++nxtr;
1072  out(nxtr,nxtc) = (*this)(r,c);
1073  }
1074  }
1075  out[i] = row;
1076  out(j) = col; // overwrites row's j'th element
1077  return out;
1078  }
1079 
1080  // These assume we are given a pointer to d[0] of a Mat<M,N,E,CS,RS> like this one.
1081  static const Mat& getAs(const ELT* p) {return *reinterpret_cast<const Mat*>(p);}
1082  static Mat& updAs(ELT* p) {return *reinterpret_cast<Mat*>(p);}
1083 
1084  // Note packed spacing
1086  Mat<M,N,ELT,M,1> m;
1087  m.setToNaN();
1088  return m;
1089  }
1090 
1092  bool isNaN() const {
1093  for (int j=0; j<N; ++j)
1094  if (this->col(j).isNaN())
1095  return true;
1096  return false;
1097  }
1098 
1101  bool isInf() const {
1102  bool seenInf = false;
1103  for (int j=0; j<N; ++j) {
1104  if (!this->col(j).isFinite()) {
1105  if (!this->col(j).isInf())
1106  return false; // something bad was found
1107  seenInf = true;
1108  }
1109  }
1110  return seenInf;
1111  }
1112 
1114  bool isFinite() const {
1115  for (int j=0; j<N; ++j)
1116  if (!this->col(j).isFinite())
1117  return false;
1118  return true;
1119  }
1120 
1124 
1127  template <class E2, int CS2, int RS2>
1128  bool isNumericallyEqual(const Mat<M,N,E2,CS2,RS2>& m, double tol) const {
1129  for (int j=0; j < N; ++j)
1130  if (!(*this)(j).isNumericallyEqual(m(j), tol))
1131  return false;
1132  return true;
1133  }
1134 
1138  template <class E2, int CS2, int RS2>
1140  const double tol = std::max(getDefaultTolerance(),m.getDefaultTolerance());
1141  return isNumericallyEqual(m, tol);
1142  }
1143 
1149  bool isNumericallyEqual
1150  (const ELT& e,
1151  double tol = getDefaultTolerance()) const
1152  {
1153  for (int i=0; i<M; ++i)
1154  for (int j=0; j<N; ++j) {
1155  if (i==j) {
1156  if (!CNT<ELT>::isNumericallyEqual((*this)(i,i), e, tol))
1157  return false;
1158  } else {
1159  // off-diagonals must be zero
1160  if (!CNT<ELT>::isNumericallyEqual((*this)(i,j), ELT(0), tol))
1161  return false;
1162  }
1163  }
1164  return true;
1165  }
1166 
1174  bool isNumericallySymmetric(double tol = getDefaultTolerance()) const {
1175  if (M != N) return false; // handled at compile time
1176  for (int j=0; j<M; ++j)
1177  for (int i=j; i<M; ++i)
1179  return false;
1180  return true;
1181  }
1182 
1189  bool isExactlySymmetric() const {
1190  if (M != N) return false; // handled at compile time
1191  for (int j=0; j<M; ++j)
1192  for (int i=j; i<M; ++i)
1193  if (elt(j,i) != CNT<ELT>::transpose(elt(i,j)))
1194  return false;
1195  return true;
1196  }
1197 
1199  TRow colSum() const {
1200  TRow temp;
1201  for (int j = 0; j < N; ++j)
1202  temp[j] = col(j).sum();
1203  return temp;
1204  }
1207  TRow sum() const {return colSum();}
1208 
1210  TCol rowSum() const {
1211  TCol temp;
1212  for (int i = 0; i < M; ++i)
1213  temp[i] = row(i).sum();
1214  return temp;
1215  }
1216 
1217  // Functions to be used for Scripting in MATLAB and languages that do not support operator overloading
1219  std::string toString() const {
1220  std::stringstream stream;
1221  stream << (*this) ;
1222  return stream.str();
1223  }
1225  const ELT& get(int i,int j) const { return elt(i,j); }
1227  void set(int i,int j, const ELT& value) { elt(i,j)=value; }
1228 
1229 private:
1230  E d[NActualElements];
1231 
1232  // This permits running through d as though it were stored
1233  // in row order with packed elements. Pass in the k'th value
1234  // of the row ordering and get back the index into d where
1235  // that element is stored.
1236  int rIx(int k) const {
1237  const int row = k / N;
1238  const int col = k % N; // that's modulus, not cross product!
1239  return row*RS + col*CS;
1240  }
1241 };
1242 
1244 // Global operators involving two matrices. //
1245 // m+m, m-m, m*m, m==m, m!=m //
1247 
1248 template <int M, int N, class EL, int CSL, int RSL, class ER, int CSR, int RSR> inline
1249 typename Mat<M,N,EL,CSL,RSL>::template Result<Mat<M,N,ER,CSR,RSR> >::Add
1251  return Mat<M,N,EL,CSL,RSL>::template Result<Mat<M,N,ER,CSR,RSR> >
1252  ::AddOp::perform(l,r);
1253 }
1254 
1255 template <int M, int N, class EL, int CSL, int RSL, class ER, int CSR, int RSR> inline
1256 typename Mat<M,N,EL,CSL,RSL>::template Result<Mat<M,N,ER,CSR,RSR> >::Sub
1258  return Mat<M,N,EL,CSL,RSL>::template Result<Mat<M,N,ER,CSR,RSR> >
1259  ::SubOp::perform(l,r);
1260 }
1261 
1262 // Matrix multiply of an MxN by NxP to produce a packed MxP.
1263 template <int M, int N, class EL, int CSL, int RSL, int P, class ER, int CSR, int RSR> inline
1264 typename Mat<M,N,EL,CSL,RSL>::template Result<Mat<N,P,ER,CSR,RSR> >::Mul
1266  return Mat<M,N,EL,CSL,RSL>::template Result<Mat<N,P,ER,CSR,RSR> >
1267  ::MulOp::perform(l,r);
1268 }
1269 
1270 // Non-conforming matrix multiply of an MxN by MMxNN; will be a scalar multiply if one
1271 // has scalar elements and the other has composite elements.
1272 template <int M, int N, class EL, int CSL, int RSL, int MM, int NN, class ER, int CSR, int RSR> inline
1273 typename Mat<M,N,EL,CSL,RSL>::template Result<Mat<MM,NN,ER,CSR,RSR> >::MulNon
1275  return Mat<M,N,EL,CSL,RSL>::template Result<Mat<MM,NN,ER,CSR,RSR> >
1276  ::MulOpNonConforming::perform(l,r);
1277 }
1278 
1279 template <int M, int N, class EL, int CSL, int RSL, class ER, int CSR, int RSR> inline
1281  for (int j=0; j<N; ++j)
1282  if (l(j) != r(j)) return false;
1283  return true;
1284 }
1285 template <int M, int N, class EL, int CSL, int RSL, class ER, int CSR, int RSR> inline
1287  return !(l==r);
1288 }
1289 
1290 
1292 // Global operators involving a matrix and a scalar. //
1294 
1295 // SCALAR MULTIPLY
1296 
1297 // m = m*real, real*m
1298 template <int M, int N, class E, int CS, int RS> inline
1300 operator*(const Mat<M,N,E,CS,RS>& l, const float& r)
1301  { return Mat<M,N,E,CS,RS>::template Result<float>::MulOp::perform(l,r); }
1302 template <int M, int N, class E, int CS, int RS> inline
1304 operator*(const float& l, const Mat<M,N,E,CS,RS>& r) {return r*l;}
1305 
1306 template <int M, int N, class E, int CS, int RS> inline
1308 operator*(const Mat<M,N,E,CS,RS>& l, const double& r)
1309  { return Mat<M,N,E,CS,RS>::template Result<double>::MulOp::perform(l,r); }
1310 template <int M, int N, class E, int CS, int RS> inline
1312 operator*(const double& l, const Mat<M,N,E,CS,RS>& r) {return r*l;}
1313 
1314 // m = m*int, int*m -- just convert int to m's precision float
1315 template <int M, int N, class E, int CS, int RS> inline
1316 typename Mat<M,N,E,CS,RS>::template Result<typename CNT<E>::Precision>::Mul
1317 operator*(const Mat<M,N,E,CS,RS>& l, int r) {return l * (typename CNT<E>::Precision)r;}
1318 template <int M, int N, class E, int CS, int RS> inline
1319 typename Mat<M,N,E,CS,RS>::template Result<typename CNT<E>::Precision>::Mul
1320 operator*(int l, const Mat<M,N,E,CS,RS>& r) {return r * (typename CNT<E>::Precision)l;}
1321 
1322 // Complex, conjugate, and negator are all easy to templatize.
1323 
1324 // m = m*complex, complex*m
1325 template <int M, int N, class E, int CS, int RS, class R> inline
1326 typename Mat<M,N,E,CS,RS>::template Result<std::complex<R> >::Mul
1327 operator*(const Mat<M,N,E,CS,RS>& l, const std::complex<R>& r)
1328  { return Mat<M,N,E,CS,RS>::template Result<std::complex<R> >::MulOp::perform(l,r); }
1329 template <int M, int N, class E, int CS, int RS, class R> inline
1330 typename Mat<M,N,E,CS,RS>::template Result<std::complex<R> >::Mul
1331 operator*(const std::complex<R>& l, const Mat<M,N,E,CS,RS>& r) {return r*l;}
1332 
1333 // m = m*conjugate, conjugate*m (convert conjugate->complex)
1334 template <int M, int N, class E, int CS, int RS, class R> inline
1335 typename Mat<M,N,E,CS,RS>::template Result<std::complex<R> >::Mul
1336 operator*(const Mat<M,N,E,CS,RS>& l, const conjugate<R>& r) {return l*(std::complex<R>)r;}
1337 template <int M, int N, class E, int CS, int RS, class R> inline
1338 typename Mat<M,N,E,CS,RS>::template Result<std::complex<R> >::Mul
1339 operator*(const conjugate<R>& l, const Mat<M,N,E,CS,RS>& r) {return r*(std::complex<R>)l;}
1340 
1341 // m = m*negator, negator*m: convert negator to standard number
1342 template <int M, int N, class E, int CS, int RS, class R> inline
1343 typename Mat<M,N,E,CS,RS>::template Result<typename negator<R>::StdNumber>::Mul
1344 operator*(const Mat<M,N,E,CS,RS>& l, const negator<R>& r) {return l * (typename negator<R>::StdNumber)(R)r;}
1345 template <int M, int N, class E, int CS, int RS, class R> inline
1346 typename Mat<M,N,E,CS,RS>::template Result<typename negator<R>::StdNumber>::Mul
1347 operator*(const negator<R>& l, const Mat<M,N,E,CS,RS>& r) {return r * (typename negator<R>::StdNumber)(R)l;}
1348 
1349 
1350 // SCALAR DIVIDE. This is a scalar operation when the scalar is on the right,
1351 // but when it is on the left it means scalar * pseudoInverse(mat),
1352 // which is a matrix whose type is like the matrix's Hermitian transpose.
1353 // TODO: for now it is just going to call mat.invert() which will fail on
1354 // singular matrices.
1355 
1356 // m = m/real, real/m
1357 template <int M, int N, class E, int CS, int RS> inline
1359 operator/(const Mat<M,N,E,CS,RS>& l, const float& r)
1360 { return Mat<M,N,E,CS,RS>::template Result<float>::DvdOp::perform(l,r); }
1361 
1362 template <int M, int N, class E, int CS, int RS> inline
1363 typename CNT<float>::template Result<Mat<M,N,E,CS,RS> >::Dvd
1364 operator/(const float& l, const Mat<M,N,E,CS,RS>& r)
1365 { return l * r.invert(); }
1366 
1367 template <int M, int N, class E, int CS, int RS> inline
1369 operator/(const Mat<M,N,E,CS,RS>& l, const double& r)
1370 { return Mat<M,N,E,CS,RS>::template Result<double>::DvdOp::perform(l,r); }
1371 
1372 template <int M, int N, class E, int CS, int RS> inline
1373 typename CNT<double>::template Result<Mat<M,N,E,CS,RS> >::Dvd
1374 operator/(const double& l, const Mat<M,N,E,CS,RS>& r)
1375 { return l * r.invert(); }
1376 
1377 // m = m/int, int/m -- just convert int to m's precision float
1378 template <int M, int N, class E, int CS, int RS> inline
1379 typename Mat<M,N,E,CS,RS>::template Result<typename CNT<E>::Precision>::Dvd
1380 operator/(const Mat<M,N,E,CS,RS>& l, int r)
1381 { return l / (typename CNT<E>::Precision)r; }
1382 
1383 template <int M, int N, class E, int CS, int RS> inline
1384 typename CNT<typename CNT<E>::Precision>::template Result<Mat<M,N,E,CS,RS> >::Dvd
1385 operator/(int l, const Mat<M,N,E,CS,RS>& r)
1386 { return (typename CNT<E>::Precision)l / r; }
1387 
1388 
1389 // Complex, conjugate, and negator are all easy to templatize.
1390 
1391 // m = m/complex, complex/m
1392 template <int M, int N, class E, int CS, int RS, class R> inline
1393 typename Mat<M,N,E,CS,RS>::template Result<std::complex<R> >::Dvd
1394 operator/(const Mat<M,N,E,CS,RS>& l, const std::complex<R>& r)
1395  { return Mat<M,N,E,CS,RS>::template Result<std::complex<R> >::DvdOp::perform(l,r); }
1396 template <int M, int N, class E, int CS, int RS, class R> inline
1397 typename CNT<std::complex<R> >::template Result<Mat<M,N,E,CS,RS> >::Dvd
1398 operator/(const std::complex<R>& l, const Mat<M,N,E,CS,RS>& r)
1399  { return CNT<std::complex<R> >::template Result<Mat<M,N,E,CS,RS> >::DvdOp::perform(l,r); }
1400 
1401 // m = m/conjugate, conjugate/m (convert conjugate->complex)
1402 template <int M, int N, class E, int CS, int RS, class R> inline
1403 typename Mat<M,N,E,CS,RS>::template Result<std::complex<R> >::Dvd
1404 operator/(const Mat<M,N,E,CS,RS>& l, const conjugate<R>& r) {return l/(std::complex<R>)r;}
1405 template <int M, int N, class E, int CS, int RS, class R> inline
1406 typename CNT<std::complex<R> >::template Result<Mat<M,N,E,CS,RS> >::Dvd
1407 operator/(const conjugate<R>& l, const Mat<M,N,E,CS,RS>& r) {return (std::complex<R>)l/r;}
1408 
1409 // m = m/negator, negator/m: convert negator to a standard number
1410 template <int M, int N, class E, int CS, int RS, class R> inline
1411 typename Mat<M,N,E,CS,RS>::template Result<typename negator<R>::StdNumber>::Dvd
1412 operator/(const Mat<M,N,E,CS,RS>& l, const negator<R>& r) {return l/(typename negator<R>::StdNumber)(R)r;}
1413 template <int M, int N, class E, int CS, int RS, class R> inline
1414 typename CNT<R>::template Result<Mat<M,N,E,CS,RS> >::Dvd
1415 operator/(const negator<R>& l, const Mat<M,N,E,CS,RS>& r) {return (typename negator<R>::StdNumber)(R)l/r;}
1416 
1417 
1418 // Add and subtract are odd as scalar ops. They behave as though the
1419 // scalar stands for a conforming matrix whose diagonal elements are that,
1420 // scalar and then a normal matrix add or subtract is done.
1421 
1422 // SCALAR ADD
1423 
1424 // m = m+real, real+m
1425 template <int M, int N, class E, int CS, int RS> inline
1427 operator+(const Mat<M,N,E,CS,RS>& l, const float& r)
1428  { return Mat<M,N,E,CS,RS>::template Result<float>::AddOp::perform(l,r); }
1429 template <int M, int N, class E, int CS, int RS> inline
1431 operator+(const float& l, const Mat<M,N,E,CS,RS>& r) {return r+l;}
1432 
1433 template <int M, int N, class E, int CS, int RS> inline
1435 operator+(const Mat<M,N,E,CS,RS>& l, const double& r)
1436  { return Mat<M,N,E,CS,RS>::template Result<double>::AddOp::perform(l,r); }
1437 template <int M, int N, class E, int CS, int RS> inline
1439 operator+(const double& l, const Mat<M,N,E,CS,RS>& r) {return r+l;}
1440 
1441 // m = m+int, int+m -- just convert int to m's precision float
1442 template <int M, int N, class E, int CS, int RS> inline
1443 typename Mat<M,N,E,CS,RS>::template Result<typename CNT<E>::Precision>::Add
1444 operator+(const Mat<M,N,E,CS,RS>& l, int r) {return l + (typename CNT<E>::Precision)r;}
1445 template <int M, int N, class E, int CS, int RS> inline
1446 typename Mat<M,N,E,CS,RS>::template Result<typename CNT<E>::Precision>::Add
1447 operator+(int l, const Mat<M,N,E,CS,RS>& r) {return r + (typename CNT<E>::Precision)l;}
1448 
1449 // Complex, conjugate, and negator are all easy to templatize.
1450 
1451 // m = m+complex, complex+m
1452 template <int M, int N, class E, int CS, int RS, class R> inline
1453 typename Mat<M,N,E,CS,RS>::template Result<std::complex<R> >::Add
1454 operator+(const Mat<M,N,E,CS,RS>& l, const std::complex<R>& r)
1455  { return Mat<M,N,E,CS,RS>::template Result<std::complex<R> >::AddOp::perform(l,r); }
1456 template <int M, int N, class E, int CS, int RS, class R> inline
1457 typename Mat<M,N,E,CS,RS>::template Result<std::complex<R> >::Add
1458 operator+(const std::complex<R>& l, const Mat<M,N,E,CS,RS>& r) {return r+l;}
1459 
1460 // m = m+conjugate, conjugate+m (convert conjugate->complex)
1461 template <int M, int N, class E, int CS, int RS, class R> inline
1462 typename Mat<M,N,E,CS,RS>::template Result<std::complex<R> >::Add
1463 operator+(const Mat<M,N,E,CS,RS>& l, const conjugate<R>& r) {return l+(std::complex<R>)r;}
1464 template <int M, int N, class E, int CS, int RS, class R> inline
1465 typename Mat<M,N,E,CS,RS>::template Result<std::complex<R> >::Add
1466 operator+(const conjugate<R>& l, const Mat<M,N,E,CS,RS>& r) {return r+(std::complex<R>)l;}
1467 
1468 // m = m+negator, negator+m: convert negator to standard number
1469 template <int M, int N, class E, int CS, int RS, class R> inline
1470 typename Mat<M,N,E,CS,RS>::template Result<typename negator<R>::StdNumber>::Add
1471 operator+(const Mat<M,N,E,CS,RS>& l, const negator<R>& r) {return l + (typename negator<R>::StdNumber)(R)r;}
1472 template <int M, int N, class E, int CS, int RS, class R> inline
1473 typename Mat<M,N,E,CS,RS>::template Result<typename negator<R>::StdNumber>::Add
1474 operator+(const negator<R>& l, const Mat<M,N,E,CS,RS>& r) {return r + (typename negator<R>::StdNumber)(R)l;}
1475 
1476 // SCALAR SUBTRACT -- careful, not commutative.
1477 
1478 // m = m-real, real-m
1479 template <int M, int N, class E, int CS, int RS> inline
1481 operator-(const Mat<M,N,E,CS,RS>& l, const float& r)
1482  { return Mat<M,N,E,CS,RS>::template Result<float>::SubOp::perform(l,r); }
1483 template <int M, int N, class E, int CS, int RS> inline
1484 typename CNT<float>::template Result<Mat<M,N,E,CS,RS> >::Sub
1485 operator-(const float& l, const Mat<M,N,E,CS,RS>& r)
1486  { return CNT<float>::template Result<Mat<M,N,E,CS,RS> >::SubOp::perform(l,r); }
1487 
1488 template <int M, int N, class E, int CS, int RS> inline
1490 operator-(const Mat<M,N,E,CS,RS>& l, const double& r)
1491  { return Mat<M,N,E,CS,RS>::template Result<double>::SubOp::perform(l,r); }
1492 template <int M, int N, class E, int CS, int RS> inline
1493 typename CNT<double>::template Result<Mat<M,N,E,CS,RS> >::Sub
1494 operator-(const double& l, const Mat<M,N,E,CS,RS>& r)
1495  { return CNT<double>::template Result<Mat<M,N,E,CS,RS> >::SubOp::perform(l,r); }
1496 
1497 // m = m-int, int-m // just convert int to m's precision float
1498 template <int M, int N, class E, int CS, int RS> inline
1499 typename Mat<M,N,E,CS,RS>::template Result<typename CNT<E>::Precision>::Sub
1500 operator-(const Mat<M,N,E,CS,RS>& l, int r) {return l - (typename CNT<E>::Precision)r;}
1501 template <int M, int N, class E, int CS, int RS> inline
1502 typename CNT<typename CNT<E>::Precision>::template Result<Mat<M,N,E,CS,RS> >::Sub
1503 operator-(int l, const Mat<M,N,E,CS,RS>& r) {return (typename CNT<E>::Precision)l - r;}
1504 
1505 
1506 // Complex, conjugate, and negator are all easy to templatize.
1507 
1508 // m = m-complex, complex-m
1509 template <int M, int N, class E, int CS, int RS, class R> inline
1510 typename Mat<M,N,E,CS,RS>::template Result<std::complex<R> >::Sub
1511 operator-(const Mat<M,N,E,CS,RS>& l, const std::complex<R>& r)
1512  { return Mat<M,N,E,CS,RS>::template Result<std::complex<R> >::SubOp::perform(l,r); }
1513 template <int M, int N, class E, int CS, int RS, class R> inline
1514 typename CNT<std::complex<R> >::template Result<Mat<M,N,E,CS,RS> >::Sub
1515 operator-(const std::complex<R>& l, const Mat<M,N,E,CS,RS>& r)
1516  { return CNT<std::complex<R> >::template Result<Mat<M,N,E,CS,RS> >::SubOp::perform(l,r); }
1517 
1518 // m = m-conjugate, conjugate-m (convert conjugate->complex)
1519 template <int M, int N, class E, int CS, int RS, class R> inline
1520 typename Mat<M,N,E,CS,RS>::template Result<std::complex<R> >::Sub
1521 operator-(const Mat<M,N,E,CS,RS>& l, const conjugate<R>& r) {return l-(std::complex<R>)r;}
1522 template <int M, int N, class E, int CS, int RS, class R> inline
1523 typename CNT<std::complex<R> >::template Result<Mat<M,N,E,CS,RS> >::Sub
1524 operator-(const conjugate<R>& l, const Mat<M,N,E,CS,RS>& r) {return (std::complex<R>)l-r;}
1525 
1526 // m = m-negator, negator-m: convert negator to standard number
1527 template <int M, int N, class E, int CS, int RS, class R> inline
1528 typename Mat<M,N,E,CS,RS>::template Result<typename negator<R>::StdNumber>::Sub
1529 operator-(const Mat<M,N,E,CS,RS>& l, const negator<R>& r) {return l-(typename negator<R>::StdNumber)(R)r;}
1530 template <int M, int N, class E, int CS, int RS, class R> inline
1531 typename CNT<R>::template Result<Mat<M,N,E,CS,RS> >::Sub
1532 operator-(const negator<R>& l, const Mat<M,N,E,CS,RS>& r) {return (typename negator<R>::StdNumber)(R)l-r;}
1533 
1534 
1535 // Mat I/O
1536 template <int M, int N, class E, int CS, int RS, class CHAR, class TRAITS> inline
1537 std::basic_ostream<CHAR,TRAITS>&
1538 operator<<(std::basic_ostream<CHAR,TRAITS>& o, const Mat<M,N,E,CS,RS>& m) {
1539  for (int i=0;i<M;++i) {
1540  o << std::endl << "[";
1541  for (int j=0;j<N;++j)
1542  o << (j>0?",":"") << m(i,j);
1543  o << "]";
1544  }
1545  if (M) o << std::endl;
1546  return o;
1547 }
1548 
1549 template <int M, int N, class E, int CS, int RS, class CHAR, class TRAITS> inline
1550 std::basic_istream<CHAR,TRAITS>&
1551 operator>>(std::basic_istream<CHAR,TRAITS>& is, Mat<M,N,E,CS,RS>& m) {
1552  // TODO: not sure how to do Vec input yet
1553  assert(false);
1554  return is;
1555 }
1556 
1557 } //namespace SimTK
1558 
1559 
1560 #endif //SimTK_SIMMATRIX_SMALLMATRIX_MAT_H_
EltResult< E2 >::Dvd elementwiseDivide(const Mat< M, N, E2, CS2, RS2 > &r) const
Definition: Mat.h:621
CNT< E >::TNormalize ENormalize
Definition: Mat.h:114
Mat(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, const E &e9, const E &e10)
Definition: Mat.h:408
Matrix_< E > operator/(const MatrixBase< E > &l, const typename CNT< E >::StdNumber &r)
Definition: BigMatrix.h:613
TAbs abs() const
Elementwise absolute value; that is, the return value has the same dimensions as this Mat but with ea...
Definition: Mat.h:228
Definition: Mat.h:126
bool isFinite() const
Return true if no element contains an Infinity or a NaN.
Definition: Mat.h:1114
Mat & operator=(const EE *p)
Definition: Mat.h:543
Mat & scalarDivideEq(const EE &ee)
Definition: Mat.h:899
Mat(const TRow &r0)
Definition: Mat.h:446
SymMat< M, E2, RS2 >::template Result< Mat >::Sub conformingSubtractFromLeft(const SymMat< M, E2, RS2 > &sy) const
Definition: Mat.h:646
CNT< E >::TStandard EStandard
Definition: Mat.h:112
K::ScalarNormSq ScalarNormSq
Definition: CompositeNumericalTypes.h:166
Mat(const TCol &r0, const TCol &r1, const TCol &r2, const TCol &r3)
Definition: Mat.h:495
Mat & scalarTimesEq(const EE &ee)
Definition: Mat.h:894
K::ULessScalar ULessScalar
Definition: CompositeNumericalTypes.h:161
Mat(const Row< N, EE, SS > &r0, const Row< N, EE, SS > &r1, const Row< N, EE, SS > &r2, const Row< N, EE, SS > &r3, const Row< N, EE, SS > &r4)
Definition: Mat.h:477
Mat< M, N, typename CNT< E >::template Result< EE >::Add > scalarAdd(const EE &e) const
Definition: Mat.h:848
Mat< M+1, N+1, E, M+1, 1 > TAppendRowCol
Definition: Mat.h:185
Mat< M, N, typename CNT< EE >::template Result< E >::Dvd > scalarDivideFromLeft(const EE &e) const
Definition: Mat.h:840
Mat & operator/=(const EE &e)
Definition: Mat.h:876
K::TReal TReal
Definition: CompositeNumericalTypes.h:141
EScalar Scalar
Definition: Mat.h:187
TStandard standardize() const
Definition: Mat.h:234
TDiag & diag()
This non-const version of diag() is an alternate name for updDiag() available for historical reasons...
Definition: Mat.h:807
Mat(const TRow &r0, const TRow &r1, const TRow &r2, const TRow &r3, const TRow &r4)
Definition: Mat.h:455
MulCNTs< M, N, ArgDepth, Mat, ColSpacing, RowSpacing, CNT< P >::NRows, CNT< P >::NCols, CNT< P >::ArgDepth, P, CNT< P >::ColSpacing, CNT< P >::RowSpacing > MulOp
Definition: Mat.h:256
This is a small, fixed-size symmetric or Hermitian matrix designed for no-overhead inline computation...
Definition: SimTKcommon/include/SimTKcommon/internal/common.h:621
Mat(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, const E &e9, const E &e10, const E &e11, const E &e12, const E &e13)
Definition: Mat.h:425
bool isInf() const
Return true if any element of this Mat contains a +Inf or -Inf somewhere but no element contains a Na...
Definition: Mat.h:1101
void setToZero()
Definition: Mat.h:909
CNT< E >::StdNumber EStdNumber
Definition: Mat.h:119
Vec & scalarEq(const EE &ee)
Definition: Vec.h:783
Mat(int i)
Explicit construction from an int value means we convert the int into an object of this Mat&#39;s element...
Definition: Mat.h:376
Mat< M, N, ENeg, CS, RS > TNeg
Definition: Mat.h:151
TRow sum() const
This is an alternate name for colSum(); behaves like the Matlab function of the same name...
Definition: Mat.h:1207
CNT< E >::ULessScalar EULessScalar
Definition: Mat.h:117
Mat< M, N, ENormalize, M, 1 > TNormalize
Definition: Mat.h:172
Mat< M, N+1, E, M, 1 > TAppendCol
Definition: Mat.h:184
Mat(const Row< N, EE, SS > &r0, const Row< N, EE, SS > &r1, const Row< N, EE, SS > &r2, const Row< N, EE, SS > &r3)
Definition: Mat.h:473
Mat(const EE *p)
Definition: Mat.h:532
Vec< MinDim, E, RS+CS > TDiag
Definition: Mat.h:164
This is the top-level SimTK namespace into which all SimTK names are placed to avoid collision with o...
Definition: Assembler.h:37
Mat(const Vec< M, EE, SS > &r0, const Vec< M, EE, SS > &r1, const Vec< M, EE, SS > &r2, const Vec< M, EE, SS > &r3, const Vec< M, EE, SS > &r4, const Vec< M, EE, SS > &r5)
Definition: Mat.h:525
NTraits< N >::StdNumber StdNumber
Definition: negator.h:107
SimTK::conjugate<R> should be instantiated only for float, double.
Definition: String.h:45
K::TSqrt TSqrt
Definition: CompositeNumericalTypes.h:154
TReal & real()
Definition: Mat.h:745
static TSqrt sqrt(const K &t)
Definition: CompositeNumericalTypes.h:239
Mat< M+1, N, E, M+1, 1 > TAppendRow
Definition: Mat.h:183
static int nrow()
Return the number of rows in this Mat, echoing the value supplied for the template parameter M...
Definition: Mat.h:200
CNT< E >::TNeg ENeg
Definition: Mat.h:100
Mat & operator*=(const EE &e)
Definition: Mat.h:875
Mat< M, N, typename CNT< E >::template Result< EE >::Sub > scalarSubtract(const EE &e) const
Definition: Mat.h:856
Mat & operator=(const Mat< M, N, EE, CSS, RSS > &mm)
Definition: Mat.h:538
Definition: Mat.h:131
K::Scalar Scalar
Definition: CompositeNumericalTypes.h:160
EStdNumber StdNumber
Definition: Mat.h:190
static const Mat & getAs(const ELT *p)
Definition: Mat.h:1081
Definition: Mat.h:127
K::TNormalize TNormalize
Definition: CompositeNumericalTypes.h:158
Mat(const Row< N, EE, SS > &r0)
Definition: Mat.h:465
CNT< E >::TSqHermT ESqHermT
Definition: Mat.h:107
const THerm & transpose() const
Definition: Mat.h:732
DvdOp::Type Dvd
Definition: Mat.h:267
MulOp::Type Mul
Definition: Mat.h:257
Mat(const E &e0, const E &e1)
Definition: Mat.h:380
CNT< E >::TWithoutNegator EWithoutNegator
Definition: Mat.h:101
const TCol & col(int j) const
Definition: Mat.h:777
TRow colSum() const
Returns a row vector (Row) containing the column sums of this matrix.
Definition: Mat.h:1199
Result< Mat< M, N, E2, CS2, RS2 > >::Sub conformingSubtract(const Mat< M, N, E2, CS2, RS2 > &r) const
Definition: Mat.h:596
TRow & operator[](int i)
Definition: Mat.h:682
Definition: Mat.h:130
K::TImag TImag
Definition: CompositeNumericalTypes.h:142
Definition: Mat.h:244
Mat(const SymMat< M, ELT > &src)
Explicit construction of a Mat from a SymMat (symmetric/Hermitian matrix).
Definition: Mat.h:316
Mat(const Mat< M, N, ENeg, CSS, RSS > &src)
This provides an implicit conversion from a Mat of the same dimensions and negated element type...
Definition: Mat.h:341
const TCol & operator()(int j) const
Definition: Mat.h:683
Mat()
Default construction initializes to NaN when debugging but is left uninitialized otherwise to ensure ...
Definition: Mat.h:287
E & operator()(int i, int j)
Definition: Mat.h:687
EStandard sum() const
Sum just adds up all the elements into a single return element that is the same type as this Vec&#39;s el...
Definition: Vec.h:366
AddCNTs< M, N, ArgDepth, Mat, ColSpacing, RowSpacing, CNT< P >::NRows, CNT< P >::NCols, CNT< P >::ArgDepth, P, CNT< P >::ColSpacing, CNT< P >::RowSpacing > AddOp
Definition: Mat.h:271
Mat< M, N, E, M, 1 > TPacked
Definition: Mat.h:179
std::basic_istream< CHAR, TRAITS > & operator>>(std::basic_istream< CHAR, TRAITS > &is, conjugate< R > &c)
Definition: conjugate.h:505
Definition: Mat.h:135
CNT< E >::ScalarNormSq EScalarNormSq
Definition: Mat.h:121
Mat(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, const E &e9, const E &e10, const E &e11, const E &e12)
Definition: Mat.h:419
Definition: Mat.h:281
Result< Mat< M, N, E2, CS2, RS2 > >::Add conformingAdd(const Mat< M, N, E2, CS2, RS2 > &r) const
Definition: Mat.h:588
negator<N>, where N is a number type (real, complex, conjugate), is represented in memory identically...
Definition: String.h:44
Mat(const E &e)
Explicit construction from a single element e of this Mat&#39;s element type E sets all the main diagonal...
Definition: Mat.h:360
Definition: CompositeNumericalTypes.h:120
Mat< M, N, ESqrt, M, 1 > TSqrt
Definition: Mat.h:168
TWithoutNegator & updCastAwayNegatorIfAny()
Definition: Mat.h:766
const EHerm & getEltUpper(int i, int j) const
Definition: SymMat.h:842
static double getDefaultTolerance()
Definition: CompositeNumericalTypes.h:269
Mat(const E &e0, const E &e1, const E &e2, const E &e3, const E &e4, const E &e5)
Definition: Mat.h:388
Mat< M, N, EStandard, M, 1 > TStandard
Definition: Mat.h:170
Definition: Mat.h:143
bool operator==(const PhiMatrix &p1, const PhiMatrix &p2)
Definition: SpatialAlgebra.h:791
THerm TransposeType
Definition: Mat.h:194
Mat< M, N, EImag, CS *CNT< E >::RealStrideFactor, RS *CNT< E >::RealStrideFactor > TImag
Definition: Mat.h:157
TDiag & updDiag()
Select main diagonal (of largest leading square if rectangular) and return it as a writable view (as ...
Definition: Mat.h:804
Definition: Mat.h:134
Mat< M, N, typename CNT< EE >::template Result< E >::Mul > scalarMultiplyFromLeft(const EE &e) const
Definition: Mat.h:825
Mat< M, N, P > Type
Definition: Mat.h:282
Mat< M2, N, E2, CS2, RS2 >::template Result< Mat >::Dvd conformingDivideFromLeft(const Mat< M2, N, E2, CS2, RS2 > &m) const
Definition: Mat.h:677
Mat(const Mat< M, N, EE, CSS, RSS > &mm)
Explicit construction of a Mat from a source Mat of the same dimensions and an assignment-compatible ...
Definition: Mat.h:354
const TNeg & negate() const
Definition: Mat.h:729
Mat< M, N, EReal, CS *CNT< E >::RealStrideFactor, RS *CNT< E >::RealStrideFactor > TReal
Definition: Mat.h:155
Definition: Mat.h:133
Mat< M, N, typename CNT< EE >::template Result< E >::Sub > scalarSubtractFromLeft(const EE &e) const
Definition: Mat.h:863
static int size()
Return the total number of elements M*N contained in this Mat.
Definition: Mat.h:197
const E & getEltLower(int i, int j) const
Definition: SymMat.h:838
TAppendRow insertRow(int i, const Row< N, EE, SS > &row) const
Return a matrix one row larger than this one by inserting a row before row i.
Definition: Mat.h:1027
Vec & scalarPlusEq(const EE &ee)
Definition: Vec.h:785
Mat & operator+=(const EE &e)
Definition: Mat.h:873
Mat(const ENeg &e)
Explicit construction from a single element e whose type is negator<E> (abbreviated ENeg here) where ...
Definition: Mat.h:367
Mat(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, const E &e9, const E &e10, const E &e11, const E &e12, const E &e13, const E &e14, const E &e15)
Definition: Mat.h:437
TDropCol dropCol(int j) const
Return a matrix one column smaller than this one by dropping column j.
Definition: Mat.h:954
TCol & operator()(int j)
Definition: Mat.h:684
K::TSqTHerm TSqTHerm
Definition: CompositeNumericalTypes.h:147
const TRow & operator[](int i) const
Definition: Mat.h:681
ENumber Number
Definition: Mat.h:189
Mat< M-1, N-1, E, M-1, 1 > TDropRowCol
Definition: Mat.h:182
TNormalize normalize() const
Definition: Mat.h:707
Mat(const TRow &r0, const TRow &r1, const TRow &r2, const TRow &r3, const TRow &r4, const TRow &r5)
Definition: Mat.h:459
bool isNumericallySymmetric(double tol=getDefaultTolerance()) const
A Matrix is symmetric (actually Hermitian) if it is square and each element (i,j) is the Hermitian tr...
Definition: Mat.h:1174
bool isNumericallyEqual(const Mat< M, N, E2, CS2, RS2 > &m, double tol) const
Test whether this matrix is numerically equal to some other matrix with the same shape, using a specified tolerance.
Definition: Mat.h:1128
static int ncol()
Return the number of columns in this Mat, echoing the value supplied for the template parameter N...
Definition: Mat.h:203
Definition: Mat.h:138
This is a fixed-length column vector designed for no-overhead inline computation. ...
Definition: SimTKcommon/include/SimTKcommon/internal/common.h:618
TInvert invert() const
Definition: SmallMatrixMixed.h:1004
TImag & imag()
Definition: Mat.h:759
TAppendCol appendCol(const Vec< M, EE, SS > &col) const
Return a matrix one column larger than this one by adding a column to the end.
Definition: Mat.h:996
Mat(const Vec< M, EE, SS > &r0, const Vec< M, EE, SS > &r1, const Vec< M, EE, SS > &r2, const Vec< M, EE, SS > &r3, const Vec< M, EE, SS > &r4)
Definition: Mat.h:520
TCol rowSum() const
Returns a column vector (Vec) containing the row sums of this matrix.
Definition: Mat.h:1210
Mat(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, const E &e9)
Definition: Mat.h:404
TAppendRowCol appendRowCol(const Row< N+1, ER, SR > &row, const Vec< M+1, EC, SC > &col) const
Return a matrix one row and one column larger than this one by adding a row to the bottom and a colum...
Definition: Mat.h:1010
const TRow & row(int i) const
Definition: Mat.h:768
K::Precision Precision
Definition: CompositeNumericalTypes.h:164
Definition: Mat.h:917
EltResult< E2 >::Mul elementwiseMultiply(const Mat< M, N, E2, CS2, RS2 > &r) const
Definition: Mat.h:611
CNT< E >::TReal EReal
Definition: Mat.h:102
Mat & scalarEq(const EE &ee)
Definition: Mat.h:880
Matrix_< E > operator*(const MatrixBase< E > &l, const typename CNT< E >::StdNumber &r)
Definition: BigMatrix.h:605
TDropRowCol dropRowCol(int i, int j) const
Return a matrix one row and one column smaller than this one by dropping row i and column j...
Definition: Mat.h:967
ELT E
Definition: Mat.h:99
Definition: Mat.h:136
Mat(const TCol &r0)
Definition: Mat.h:489
K::TInvert TInvert
Definition: CompositeNumericalTypes.h:157
Definition: Mat.h:132
Mat< MM, NN, ELT, CS, RS > Type
Definition: Mat.h:918
CNT< E >::TSqTHerm ESqTHerm
Definition: Mat.h:108
Mat< M, N, typename CNT< E >::template Result< P >::Sub, M, 1 > Sub
Definition: Mat.h:248
Mat(const Vec< M, EE, SS > &r0, const Vec< M, EE, SS > &r1, const Vec< M, EE, SS > &r2, const Vec< M, EE, SS > &r3)
Definition: Mat.h:516
Mat & operator*=(const Mat< N, N, EE, CSS, RSS > &mm)
Definition: Mat.h:574
THerm & operator~()
Definition: Mat.h:727
ELEM max(const VectorBase< ELEM > &v)
Definition: VectorMath.h:251
EScalarNormSq ScalarNormSq
Definition: Mat.h:192
Mat & operator-=(const Mat< M, N, negator< EE >, CSS, RSS > &mm)
Definition: Mat.h:566
const TNeg & operator-() const
Definition: Mat.h:724
const SubMat< MM, NN >::Type & getSubMat(int i, int j) const
Definition: Mat.h:922
Mat(const TCol &r0, const TCol &r1, const TCol &r2)
Definition: Mat.h:493
Result< SymMat< M, E2, RS2 > >::Sub conformingSubtract(const SymMat< M, E2, RS2 > &sy) const
Definition: Mat.h:639
SubOp::Type Sub
Definition: Mat.h:277
Definition: Mat.h:142
TAppendCol insertCol(int j, const Vec< M, EE, SS > &col) const
Return a matrix one column larger than this one by inserting a column before column j...
Definition: Mat.h:1044
Mat & scalarDivideEqFromLeft(const EE &ee)
Definition: Mat.h:901
Mat(const E &e0, const E &e1, const E &e2, const E &e3)
Definition: Mat.h:384
K::TPosTrans TPosTrans
Definition: CompositeNumericalTypes.h:145
Mat(const TRow &r0, const TRow &r1, const TRow &r2, const TRow &r3)
Definition: Mat.h:452
Vec< M, typename CNT< E >::template Result< EE >::Add > scalarAdd(const EE &e) const
Definition: Vec.h:752
Mat(const E &e0, const E &e1, const E &e2)
Definition: Mat.h:382
Definition: Mat.h:253
TAppendRow appendRow(const Row< N, EE, SS > &row) const
Return a matrix one row larger than this one by adding a row to the end.
Definition: Mat.h:985
Mat< M, N, E2, CS2, RS2 >::template Result< Mat >::Sub conformingSubtractFromLeft(const Mat< M, N, E2, CS2, RS2 > &l) const
Definition: Mat.h:604
bool isExactlySymmetric() const
A Matrix is symmetric (actually Hermitian) if it is square and each element (i,j) is the Hermitian (c...
Definition: Mat.h:1189
void setToNaN()
Definition: Mat.h:904
Mat< M, N, typename CNT< E >::template Result< EE >::Mul > scalarMultiply(const EE &e) const
Definition: Mat.h:819
Result< SymMat< M, E2, RS2 > >::Sub conformingSubtract(const SymMat< M, E2, RS2 > &r) const
Definition: SymMat.h:519
Mat(const TRow &r0, const TRow &r1)
Definition: Mat.h:448
static Mat< M, N, ELT, M, 1 > getNaN()
Definition: Mat.h:1085
Mat & operator+=(const Mat< M, N, EE, CSS, RSS > &mm)
Definition: Mat.h:550
ScalarNormSq normSqr() const
Definition: Mat.h:690
const E & elt(int i, int j) const
Definition: Mat.h:786
TSqrt sqrt() const
Elementwise square root; that is, the return value has the same dimensions as this Mat but with each ...
Definition: Mat.h:219
MulCNTsNonConforming< M, N, ArgDepth, Mat, ColSpacing, RowSpacing, CNT< P >::NRows, CNT< P >::NCols, CNT< P >::ArgDepth, P, CNT< P >::ColSpacing, CNT< P >::RowSpacing > MulOpNonConforming
Definition: Mat.h:261
bool isNumericallyEqual(const Mat< M, N, E2, CS2, RS2 > &m) const
Test whether this matrix is numerically equal to some other matrix with the same shape, using a default tolerance which is the looser of the default tolerances of the two objects being compared.
Definition: Mat.h:1139
K::StdNumber StdNumber
Definition: CompositeNumericalTypes.h:163
CNT< E >::THerm EHerm
Definition: Mat.h:105
bool isNaN() const
Return true if any element of this Mat contains a NaN anywhere.
Definition: Mat.h:1092
void setSubMat(int i, int j, const typename SubMat< MM, NN >::Type &value)
Definition: Mat.h:934
DvdCNTs< M, N, ArgDepth, Mat, ColSpacing, RowSpacing, CNT< P >::NRows, CNT< P >::NCols, CNT< P >::ArgDepth, P, CNT< P >::ColSpacing, CNT< P >::RowSpacing > DvdOp
Definition: Mat.h:266
SymMat< N, ESqHermT > TSqHermT
Definition: Mat.h:174
const TReal & real() const
Definition: Mat.h:744
Mat(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)
Definition: Mat.h:400
Specialized information about Composite Numerical Types which allows us to define appropriate templat...
Definition: CompositeNumericalTypes.h:136
CNT< E >::TInvert EInvert
Definition: Mat.h:113
Mat(const Vec< M, EE, SS > &r0)
Definition: Mat.h:508
EStandard sum() const
Definition: Row.h:254
Mat(const Vec< M, EE, SS > &r0, const Vec< M, EE, SS > &r1, const Vec< M, EE, SS > &r2)
Definition: Mat.h:513
This is a fixed-length row vector designed for no-overhead inline computation.
Definition: SimTKcommon/include/SimTKcommon/internal/common.h:619
Mat & operator=(const Mat &src)
Copy assignment copies only the elements that are present and does not touch any unused memory space ...
Definition: Mat.h:306
THerm & updTranspose()
Definition: Mat.h:733
CNT< E >::TAbs EAbs
Definition: Mat.h:111
SymMat< M, ESqTHerm > TSqTHerm
Definition: Mat.h:175
E & elt(int i, int j)
Definition: Mat.h:791
Definition: Mat.h:128
CNT< E >::TSqrt ESqrt
Definition: Mat.h:110
Definition: Mat.h:144
Mat(const Mat &src)
Copy constructor copies only the elements that are present and does not touch any unused memory space...
Definition: Mat.h:299
SubMat< MM, NN >::Type & updSubMat(int i, int j)
Definition: Mat.h:928
Mandatory first inclusion for any Simbody source or header file.
std::string toString() const
toString() returns a string representation of the Mat.
Definition: Mat.h:1219
TPosTrans & updPositionalTranspose()
Definition: Mat.h:737
Mat & operator+=(const Mat< M, N, negator< EE >, CSS, RSS > &mm)
Definition: Mat.h:555
Mat(const Row< N, EE, SS > &r0, const Row< N, EE, SS > &r1, const Row< N, EE, SS > &r2)
Definition: Mat.h:470
Mat< M2, M, E2, CS2, RS2 >::template Result< Mat >::Mul conformingMultiplyFromLeft(const Mat< M2, M, E2, CS2, RS2 > &m) const
Definition: Mat.h:664
Definition: Mat.h:146
CNT< E >::Scalar EScalar
Definition: Mat.h:116
Mat(const TCol &r0, const TCol &r1)
Definition: Mat.h:491
Mat< M, N, typename CNT< E >::template Result< EE >::Dvd > scalarDivide(const EE &e) const
Definition: Mat.h:834
TAppendRowCol insertRowCol(int i, int j, const Row< N+1, ER, SR > &row, const Vec< M+1, EC, SC > &col) const
Return a matrix one row and one column larger than this one by inserting a row before row i and a col...
Definition: Mat.h:1063
K::TNeg TNeg
Definition: CompositeNumericalTypes.h:139
Mat< M, N, E, CS, RS > T
Definition: Mat.h:150
static double getDefaultTolerance()
For approximate comparisons, the default tolerance to use for a matrix is its shortest dimension time...
Definition: Mat.h:1123
Result< SymMat< M, E2, RS2 > >::Add conformingAdd(const SymMat< M, E2, RS2 > &sy) const
Definition: Mat.h:632
K::TStandard TStandard
Definition: CompositeNumericalTypes.h:156
CNT< E >::TComplex EComplex
Definition: Mat.h:104
Mat< M, N, typename CNT< E >::template Result< P >::Mul, M, 1 > Mul
Definition: Mat.h:245
#define SimTK_INDEXCHECK(ix, ub, where)
Definition: ExceptionMacros.h:145
Mat< N, M, EInvert, N, 1 > TInvert
Definition: Mat.h:171
Mat(const E &e0, const E &e1, const E &e2, const E &e3, const E &e4, const E &e5, const E &e6)
Definition: Mat.h:392
EPrecision Precision
Definition: Mat.h:191
MulOpNonConforming::Type MulNon
Definition: Mat.h:262
K::TWithoutNegator TWithoutNegator
Definition: CompositeNumericalTypes.h:140
Mat(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, const E &e9, const E &e10, const E &e11, const E &e12, const E &e13, const E &e14)
Definition: Mat.h:431
const TImag & imag() const
Definition: Mat.h:754
Result< SymMat< M, E2, RS2 > >::Add conformingAdd(const SymMat< M, E2, RS2 > &r) const
Definition: SymMat.h:512
Vec< M, E, RS > TCol
Definition: Mat.h:163
EULessScalar ULessScalar
Definition: Mat.h:188
const TDiag & diag() const
Select main diagonal (of largest leading square if rectangular) and return it as a read-only view (as...
Definition: Mat.h:800
CNT< E >::TPosTrans EPosTrans
Definition: Mat.h:106
TDropRow dropRow(int i) const
Return a matrix one row smaller than this one by dropping row i.
Definition: Mat.h:942
Mat(const TRow &r0, const TRow &r1, const TRow &r2)
Definition: Mat.h:450
CNT< E >::TImag EImag
Definition: Mat.h:103
Mat< M, N, typename CNT< E >::template Result< P >::Dvd, M, 1 > Dvd
Definition: Mat.h:246
Row< N, E, CS > TRow
Definition: Mat.h:162
const THerm & operator~() const
Definition: Mat.h:726
const TWithoutNegator & castAwayNegatorIfAny() const
Definition: Mat.h:765
Mat< M, N, EComplex, CS, RS > TComplex
Definition: Mat.h:158
Mat & scalarMinusEqFromLeft(const EE &ee)
Definition: Mat.h:891
Mat(const Mat< M, N, E, CSS, RSS > &src)
This provides an implicit conversion from a Mat of the same dimensions and element type but with diff...
Definition: Mat.h:330
bool operator!=(const L &left, const R &right)
Definition: SimTKcommon/include/SimTKcommon/internal/common.h:641
Mat< M-1, N, E, M-1, 1 > TDropRow
Definition: Mat.h:180
This class represents a small matrix whose size is known at compile time, containing elements of any ...
Definition: SimTKcommon/include/SimTKcommon/internal/common.h:620
K::TComplex TComplex
Definition: CompositeNumericalTypes.h:143
TNeg & operator-()
Definition: Mat.h:725
Mat< M, N, typename CNT< E >::template Result< P >::Add, M, 1 > Add
Definition: Mat.h:247
Mat(const Row< N, EE, SS > &r0, const Row< N, EE, SS > &r1, const Row< N, EE, SS > &r2, const Row< N, EE, SS > &r3, const Row< N, EE, SS > &r4, const Row< N, EE, SS > &r5)
Definition: Mat.h:482
Mat & operator-=(const Mat< M, N, EE, CSS, RSS > &mm)
Definition: Mat.h:561
Mat(const Vec< M, EE, SS > &r0, const Vec< M, EE, SS > &r1)
Definition: Mat.h:510
K::Number Number
Definition: CompositeNumericalTypes.h:162
Mat< N, M, E, RS, CS > TPosTrans
Definition: Mat.h:160
Definition: Mat.h:141
AddOp::Type Add
Definition: Mat.h:272
Mat(const TCol &r0, const TCol &r1, const TCol &r2, const TCol &r3, const TCol &r4)
Definition: Mat.h:498
Mat & operator-=(const EE &e)
Definition: Mat.h:874
Mat(const E &e0, const E &e1, const E &e2, const E &e3, const E &e4)
Definition: Mat.h:386
CNT< E >::Precision EPrecision
Definition: Mat.h:120
Mat(const E &e0, const E &e1, const E &e2, const E &e3, const E &e4, const E &e5, const E &e6, const E &e7)
Definition: Mat.h:396
ScalarNormSq scalarNormSqr() const
Scalar norm square is the sum of squares of all the scalars that comprise the value of this Mat...
Definition: Mat.h:210
Mat(const TCol &r0, const TCol &r1, const TCol &r2, const TCol &r3, const TCol &r4, const TCol &r5)
Definition: Mat.h:502
Mat & scalarMinusEq(const EE &ee)
Definition: Mat.h:888
const E & operator()(int i, int j) const
Definition: Mat.h:686
TRow & row(int i)
Definition: Mat.h:772
const Mat & operator+() const
Definition: Mat.h:723
Result< Mat< M2, N, E2, CS2, RS2 > >::Dvd conformingDivide(const Mat< M2, N, E2, CS2, RS2 > &m) const
Definition: Mat.h:671
SubCNTs< M, N, ArgDepth, Mat, ColSpacing, RowSpacing, CNT< P >::NRows, CNT< P >::NCols, CNT< P >::ArgDepth, P, CNT< P >::ColSpacing, CNT< P >::RowSpacing > SubOp
Definition: Mat.h:276
Mat(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, const E &e9, const E &e10, const E &e11)
Definition: Mat.h:413
K::TSqHermT TSqHermT
Definition: CompositeNumericalTypes.h:146
TNeg & updNegate()
Definition: Mat.h:730
const TPosTrans & positionalTranspose() const
Definition: Mat.h:735
TCol & col(int j)
Definition: Mat.h:781
K::THerm THerm
Definition: CompositeNumericalTypes.h:144
EStandard trace() const
Definition: Mat.h:809
Vec & scalarMinusEq(const EE &ee)
Definition: Vec.h:787
Result< Mat< N, N2, E2, CS2, RS2 > >::Mul conformingMultiply(const Mat< N, N2, E2, CS2, RS2 > &m) const
Definition: Mat.h:654
Mat & scalarPlusEq(const EE &ee)
Definition: Mat.h:885
Definition: negator.h:64
Definition: Mat.h:145
CNT< ScalarNormSq >::TSqrt norm() const
Definition: Mat.h:692
const TDiag & diag() const
Definition: SymMat.h:822
CNT< E >::Number ENumber
Definition: Mat.h:118
Mat & scalarTimesEqFromLeft(const EE &ee)
Definition: Mat.h:896
Mat(const Row< N, EE, SS > &r0, const Row< N, EE, SS > &r1)
Definition: Mat.h:467
E TElement
Definition: Mat.h:161
Mat< M, N, EAbs, M, 1 > TAbs
Definition: Mat.h:169
K::TAbs TAbs
Definition: CompositeNumericalTypes.h:155
Mat< N, M, EHerm, RS, CS > THerm
Definition: Mat.h:159
static Mat & updAs(ELT *p)
Definition: Mat.h:1082
Mat< M, N, EWithoutNegator, CS, RS > TWithoutNegator
Definition: Mat.h:152
Mat< M, N-1, E, M, 1 > TDropCol
Definition: Mat.h:181