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