Simbody 3.7
Loading...
Searching...
No Matches
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
33namespace SimTK {
34
97template <int M, int N, class ELT, int CS, int RS> class Mat {
98public:
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;
115
116 typedef typename CNT<E>::Scalar EScalar;
118 typedef typename CNT<E>::Number ENumber;
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
147 };
148 #endif
149
153
161 typedef E TElement;
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
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
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
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.
250
251 // This is the composite result for m op P where P is some kind of appropriately shaped
252 // non-scalar type.
279
280 // Shape-preserving element substitution (always packed)
281 template <class P> struct Substitute {
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>
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>
471 { assert(M==3);(*this)[0]=r0;(*this)[1]=r1;(*this)[2]=r2; }
472 template <class EE, int SS>
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>
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>
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>
514 { assert(N==3);(*this)(0)=r0;(*this)(1)=r1;(*this)(2)=r2; }
515 template <class EE, int SS>
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>
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>
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&
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&
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
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(); }
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.
708 if (CNT<E>::IsScalar) {
710 } else {
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();
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
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 }
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);
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);
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);
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>
987 out.template updSubMat<M,N>(0,0) = (*this);
988 out[M] = row;
989 return out;
990 }
991
995 template <class EE, int SS>
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>
1011 const Vec<M+1,EC,SC>& col) const
1012 {
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>
1028 assert(0 <= i && i <= M);
1029 if (i==M) return appendRow(row);
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>
1045 assert(0 <= j && j <= N);
1046 if (j==N) return appendCol(col);
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>
1064 const Vec<M+1,EC,SC>& col) const {
1065 assert(0 <= i && i <= M);
1066 assert(0 <= j && j <= N);
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
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
1123 static double getDefaultTolerance() {return MinDim*CNT<ELT>::getDefaultTolerance();}
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
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
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
1229private:
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
1248template <int M, int N, class EL, int CSL, int RSL, class ER, int CSR, int RSR> inline
1249typename 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
1255template <int M, int N, class EL, int CSL, int RSL, class ER, int CSR, int RSR> inline
1256typename 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.
1263template <int M, int N, class EL, int CSL, int RSL, int P, class ER, int CSR, int RSR> inline
1264typename 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.
1272template <int M, int N, class EL, int CSL, int RSL, int MM, int NN, class ER, int CSR, int RSR> inline
1273typename 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
1279template <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}
1285template <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
1298template <int M, int N, class E, int CS, int RS> inline
1299typename Mat<M,N,E,CS,RS>::template Result<float>::Mul
1300operator*(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); }
1302template <int M, int N, class E, int CS, int RS> inline
1303typename Mat<M,N,E,CS,RS>::template Result<float>::Mul
1304operator*(const float& l, const Mat<M,N,E,CS,RS>& r) {return r*l;}
1305
1306template <int M, int N, class E, int CS, int RS> inline
1307typename Mat<M,N,E,CS,RS>::template Result<double>::Mul
1308operator*(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); }
1310template <int M, int N, class E, int CS, int RS> inline
1311typename Mat<M,N,E,CS,RS>::template Result<double>::Mul
1312operator*(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
1315template <int M, int N, class E, int CS, int RS> inline
1316typename Mat<M,N,E,CS,RS>::template Result<typename CNT<E>::Precision>::Mul
1317operator*(const Mat<M,N,E,CS,RS>& l, int r) {return l * (typename CNT<E>::Precision)r;}
1318template <int M, int N, class E, int CS, int RS> inline
1319typename Mat<M,N,E,CS,RS>::template Result<typename CNT<E>::Precision>::Mul
1320operator*(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
1325template <int M, int N, class E, int CS, int RS, class R> inline
1326typename Mat<M,N,E,CS,RS>::template Result<std::complex<R> >::Mul
1327operator*(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); }
1329template <int M, int N, class E, int CS, int RS, class R> inline
1330typename Mat<M,N,E,CS,RS>::template Result<std::complex<R> >::Mul
1331operator*(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)
1334template <int M, int N, class E, int CS, int RS, class R> inline
1335typename Mat<M,N,E,CS,RS>::template Result<std::complex<R> >::Mul
1336operator*(const Mat<M,N,E,CS,RS>& l, const conjugate<R>& r) {return l*(std::complex<R>)r;}
1337template <int M, int N, class E, int CS, int RS, class R> inline
1338typename Mat<M,N,E,CS,RS>::template Result<std::complex<R> >::Mul
1339operator*(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
1342template <int M, int N, class E, int CS, int RS, class R> inline
1343typename Mat<M,N,E,CS,RS>::template Result<typename negator<R>::StdNumber>::Mul
1344operator*(const Mat<M,N,E,CS,RS>& l, const negator<R>& r) {return l * (typename negator<R>::StdNumber)(R)r;}
1345template <int M, int N, class E, int CS, int RS, class R> inline
1346typename Mat<M,N,E,CS,RS>::template Result<typename negator<R>::StdNumber>::Mul
1347operator*(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
1357template <int M, int N, class E, int CS, int RS> inline
1358typename Mat<M,N,E,CS,RS>::template Result<float>::Dvd
1359operator/(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
1362template <int M, int N, class E, int CS, int RS> inline
1363typename CNT<float>::template Result<Mat<M,N,E,CS,RS> >::Dvd
1364operator/(const float& l, const Mat<M,N,E,CS,RS>& r)
1365{ return l * r.invert(); }
1366
1367template <int M, int N, class E, int CS, int RS> inline
1368typename Mat<M,N,E,CS,RS>::template Result<double>::Dvd
1369operator/(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
1372template <int M, int N, class E, int CS, int RS> inline
1373typename CNT<double>::template Result<Mat<M,N,E,CS,RS> >::Dvd
1374operator/(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
1378template <int M, int N, class E, int CS, int RS> inline
1379typename Mat<M,N,E,CS,RS>::template Result<typename CNT<E>::Precision>::Dvd
1381{ return l / (typename CNT<E>::Precision)r; }
1382
1383template <int M, int N, class E, int CS, int RS> inline
1384typename CNT<typename CNT<E>::Precision>::template Result<Mat<M,N,E,CS,RS> >::Dvd
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
1392template <int M, int N, class E, int CS, int RS, class R> inline
1393typename Mat<M,N,E,CS,RS>::template Result<std::complex<R> >::Dvd
1394operator/(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); }
1396template <int M, int N, class E, int CS, int RS, class R> inline
1397typename CNT<std::complex<R> >::template Result<Mat<M,N,E,CS,RS> >::Dvd
1398operator/(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)
1402template <int M, int N, class E, int CS, int RS, class R> inline
1403typename Mat<M,N,E,CS,RS>::template Result<std::complex<R> >::Dvd
1404operator/(const Mat<M,N,E,CS,RS>& l, const conjugate<R>& r) {return l/(std::complex<R>)r;}
1405template <int M, int N, class E, int CS, int RS, class R> inline
1406typename CNT<std::complex<R> >::template Result<Mat<M,N,E,CS,RS> >::Dvd
1407operator/(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
1410template <int M, int N, class E, int CS, int RS, class R> inline
1411typename Mat<M,N,E,CS,RS>::template Result<typename negator<R>::StdNumber>::Dvd
1412operator/(const Mat<M,N,E,CS,RS>& l, const negator<R>& r) {return l/(typename negator<R>::StdNumber)(R)r;}
1413template <int M, int N, class E, int CS, int RS, class R> inline
1414typename CNT<R>::template Result<Mat<M,N,E,CS,RS> >::Dvd
1415operator/(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
1425template <int M, int N, class E, int CS, int RS> inline
1426typename Mat<M,N,E,CS,RS>::template Result<float>::Add
1427operator+(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); }
1429template <int M, int N, class E, int CS, int RS> inline
1430typename Mat<M,N,E,CS,RS>::template Result<float>::Add
1431operator+(const float& l, const Mat<M,N,E,CS,RS>& r) {return r+l;}
1432
1433template <int M, int N, class E, int CS, int RS> inline
1434typename Mat<M,N,E,CS,RS>::template Result<double>::Add
1435operator+(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); }
1437template <int M, int N, class E, int CS, int RS> inline
1438typename Mat<M,N,E,CS,RS>::template Result<double>::Add
1439operator+(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
1442template <int M, int N, class E, int CS, int RS> inline
1443typename Mat<M,N,E,CS,RS>::template Result<typename CNT<E>::Precision>::Add
1444operator+(const Mat<M,N,E,CS,RS>& l, int r) {return l + (typename CNT<E>::Precision)r;}
1445template <int M, int N, class E, int CS, int RS> inline
1446typename Mat<M,N,E,CS,RS>::template Result<typename CNT<E>::Precision>::Add
1447operator+(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
1452template <int M, int N, class E, int CS, int RS, class R> inline
1453typename Mat<M,N,E,CS,RS>::template Result<std::complex<R> >::Add
1454operator+(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); }
1456template <int M, int N, class E, int CS, int RS, class R> inline
1457typename Mat<M,N,E,CS,RS>::template Result<std::complex<R> >::Add
1458operator+(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)
1461template <int M, int N, class E, int CS, int RS, class R> inline
1462typename Mat<M,N,E,CS,RS>::template Result<std::complex<R> >::Add
1463operator+(const Mat<M,N,E,CS,RS>& l, const conjugate<R>& r) {return l+(std::complex<R>)r;}
1464template <int M, int N, class E, int CS, int RS, class R> inline
1465typename Mat<M,N,E,CS,RS>::template Result<std::complex<R> >::Add
1466operator+(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
1469template <int M, int N, class E, int CS, int RS, class R> inline
1470typename Mat<M,N,E,CS,RS>::template Result<typename negator<R>::StdNumber>::Add
1471operator+(const Mat<M,N,E,CS,RS>& l, const negator<R>& r) {return l + (typename negator<R>::StdNumber)(R)r;}
1472template <int M, int N, class E, int CS, int RS, class R> inline
1473typename Mat<M,N,E,CS,RS>::template Result<typename negator<R>::StdNumber>::Add
1474operator+(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
1479template <int M, int N, class E, int CS, int RS> inline
1480typename Mat<M,N,E,CS,RS>::template Result<float>::Sub
1481operator-(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); }
1483template <int M, int N, class E, int CS, int RS> inline
1484typename CNT<float>::template Result<Mat<M,N,E,CS,RS> >::Sub
1485operator-(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
1488template <int M, int N, class E, int CS, int RS> inline
1489typename Mat<M,N,E,CS,RS>::template Result<double>::Sub
1490operator-(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); }
1492template <int M, int N, class E, int CS, int RS> inline
1493typename CNT<double>::template Result<Mat<M,N,E,CS,RS> >::Sub
1494operator-(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
1498template <int M, int N, class E, int CS, int RS> inline
1499typename Mat<M,N,E,CS,RS>::template Result<typename CNT<E>::Precision>::Sub
1500operator-(const Mat<M,N,E,CS,RS>& l, int r) {return l - (typename CNT<E>::Precision)r;}
1501template <int M, int N, class E, int CS, int RS> inline
1502typename CNT<typename CNT<E>::Precision>::template Result<Mat<M,N,E,CS,RS> >::Sub
1503operator-(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
1509template <int M, int N, class E, int CS, int RS, class R> inline
1510typename Mat<M,N,E,CS,RS>::template Result<std::complex<R> >::Sub
1511operator-(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); }
1513template <int M, int N, class E, int CS, int RS, class R> inline
1514typename CNT<std::complex<R> >::template Result<Mat<M,N,E,CS,RS> >::Sub
1515operator-(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)
1519template <int M, int N, class E, int CS, int RS, class R> inline
1520typename Mat<M,N,E,CS,RS>::template Result<std::complex<R> >::Sub
1521operator-(const Mat<M,N,E,CS,RS>& l, const conjugate<R>& r) {return l-(std::complex<R>)r;}
1522template <int M, int N, class E, int CS, int RS, class R> inline
1523typename CNT<std::complex<R> >::template Result<Mat<M,N,E,CS,RS> >::Sub
1524operator-(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
1527template <int M, int N, class E, int CS, int RS, class R> inline
1528typename Mat<M,N,E,CS,RS>::template Result<typename negator<R>::StdNumber>::Sub
1529operator-(const Mat<M,N,E,CS,RS>& l, const negator<R>& r) {return l-(typename negator<R>::StdNumber)(R)r;}
1530template <int M, int N, class E, int CS, int RS, class R> inline
1531typename CNT<R>::template Result<Mat<M,N,E,CS,RS> >::Sub
1532operator-(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
1536template <int M, int N, class E, int CS, int RS, class CHAR, class TRAITS> inline
1537std::basic_ostream<CHAR,TRAITS>&
1538operator<<(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
1549template <int M, int N, class E, int CS, int RS, class CHAR, class TRAITS> inline
1550std::basic_istream<CHAR,TRAITS>&
1551operator>>(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_
#define SimTK_INDEXCHECK(ix, ub, where)
Definition ExceptionMacros.h:145
Mandatory first inclusion for any Simbody source or header file.
Specialized information about Composite Numerical Types which allows us to define appropriate templat...
Definition CompositeNumericalTypes.h:136
K::ULessScalar ULessScalar
Definition CompositeNumericalTypes.h:161
static double getDefaultTolerance()
Definition CompositeNumericalTypes.h:269
K::ScalarNormSq ScalarNormSq
Definition CompositeNumericalTypes.h:166
K::StdNumber StdNumber
Definition CompositeNumericalTypes.h:163
static TSqrt sqrt(const K &t)
Definition CompositeNumericalTypes.h:239
K::TSqHermT TSqHermT
Definition CompositeNumericalTypes.h:146
K::TSqrt TSqrt
Definition CompositeNumericalTypes.h:154
K::TInvert TInvert
Definition CompositeNumericalTypes.h:157
K::TNormalize TNormalize
Definition CompositeNumericalTypes.h:158
K::TWithoutNegator TWithoutNegator
Definition CompositeNumericalTypes.h:140
K::TReal TReal
Definition CompositeNumericalTypes.h:141
K::THerm THerm
Definition CompositeNumericalTypes.h:144
K::TPosTrans TPosTrans
Definition CompositeNumericalTypes.h:145
K::TNeg TNeg
Definition CompositeNumericalTypes.h:139
K::TStandard TStandard
Definition CompositeNumericalTypes.h:156
K::TComplex TComplex
Definition CompositeNumericalTypes.h:143
K::TSqTHerm TSqTHerm
Definition CompositeNumericalTypes.h:147
K::TImag TImag
Definition CompositeNumericalTypes.h:142
K::Precision Precision
Definition CompositeNumericalTypes.h:164
K::Scalar Scalar
Definition CompositeNumericalTypes.h:160
K::TAbs TAbs
Definition CompositeNumericalTypes.h:155
K::Number Number
Definition CompositeNumericalTypes.h:162
This class represents a small matrix whose size is known at compile time, containing elements of any ...
Definition Mat.h:97
Mat()
Default construction initializes to NaN when debugging but is left uninitialized otherwise to ensure ...
Definition Mat.h:287
CNT< E >::TImag EImag
Definition Mat.h:103
SubMat< MM, NN >::Type & updSubMat(int i, int j)
Definition Mat.h:928
Mat< M, N, E2, CS2, RS2 >::template Result< Mat >::Sub conformingSubtractFromLeft(const Mat< M, N, E2, CS2, RS2 > &l) const
Definition Mat.h:604
Mat(const SymMat< M, ELT > &src)
Explicit construction of a Mat from a SymMat (symmetric/Hermitian matrix).
Definition Mat.h:316
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
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
const E & elt(int i, int j) const
Definition Mat.h:786
void setSubMat(int i, int j, const typename SubMat< MM, NN >::Type &value)
Definition Mat.h:934
Result< SymMat< M, E2, RS2 > >::Sub conformingSubtract(const SymMat< M, E2, RS2 > &sy) const
Definition Mat.h:639
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
Result< Mat< N, N2, E2, CS2, RS2 > >::Mul conformingMultiply(const Mat< N, N2, E2, CS2, RS2 > &m) const
Definition Mat.h:654
Mat< M, N, typename CNT< EE >::template Result< E >::Dvd > scalarDivideFromLeft(const EE &e) const
Definition Mat.h:840
EltResult< E2 >::Dvd elementwiseDivide(const Mat< M, N, E2, CS2, RS2 > &r) const
Definition Mat.h:621
Mat & scalarTimesEqFromLeft(const EE &ee)
Definition Mat.h:896
TWithoutNegator & updCastAwayNegatorIfAny()
Definition Mat.h:766
Mat< M-1, N-1, E, M-1, 1 > TDropRowCol
Definition Mat.h:182
Mat< M, N+1, E, M, 1 > TAppendCol
Definition Mat.h:184
Mat(const TCol &r0, const TCol &r1, const TCol &r2, const TCol &r3)
Definition Mat.h:495
static int nrow()
Return the number of rows in this Mat, echoing the value supplied for the template parameter M.
Definition Mat.h:200
const TReal & real() const
Definition Mat.h:744
Mat(const TCol &r0)
Definition Mat.h:489
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
const TCol & operator()(int j) const
Definition Mat.h:683
CNT< E >::TComplex EComplex
Definition Mat.h:104
CNT< E >::TStandard EStandard
Definition Mat.h:112
EStdNumber StdNumber
Definition Mat.h:190
Mat< M, N, ENeg, CS, RS > TNeg
Definition Mat.h:151
void setToNaN()
Definition Mat.h:904
E & elt(int i, int j)
Definition Mat.h:791
TInvert invert() const
Definition SmallMatrixMixed.h:1004
THerm TransposeType
Definition Mat.h:194
Mat(const TRow &r0, const TRow &r1, const TRow &r2, const TRow &r3, const TRow &r4)
Definition Mat.h:455
Mat & scalarDivideEqFromLeft(const EE &ee)
Definition Mat.h:901
static double getDefaultTolerance()
For approximate comparisons, the default tolerance to use for a matrix is its shortest dimension time...
Definition Mat.h:1123
TRow & operator[](int i)
Definition Mat.h:682
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, EAbs, M, 1 > TAbs
Definition Mat.h:169
TRow sum() const
This is an alternate name for colSum(); behaves like the Matlab function of the same name.
Definition Mat.h:1207
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,...
Definition Mat.h:1128
CNT< E >::TWithoutNegator EWithoutNegator
Definition Mat.h:101
CNT< E >::TSqHermT ESqHermT
Definition Mat.h:107
Mat< M, N, typename CNT< EE >::template Result< E >::Mul > scalarMultiplyFromLeft(const EE &e) const
Definition Mat.h:825
Mat< M, N, EStandard, M, 1 > TStandard
Definition Mat.h:170
Result< Mat< M2, N, E2, CS2, RS2 > >::Dvd conformingDivide(const Mat< M2, N, E2, CS2, RS2 > &m) const
Definition Mat.h:671
Mat(const E &e0, const E &e1, const E &e2, const E &e3)
Definition Mat.h:384
TDropRow dropRow(int i) const
Return a matrix one row smaller than this one by dropping row i.
Definition Mat.h:942
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
CNT< E >::Number ENumber
Definition Mat.h:118
Mat(const Row< N, EE, SS > &r0, const Row< N, EE, SS > &r1)
Definition Mat.h:467
const THerm & transpose() const
Definition Mat.h:732
ScalarNormSq normSqr() const
Definition Mat.h:690
SymMat< N, ESqHermT > TSqHermT
Definition Mat.h:174
Mat(const Row< N, EE, SS > &r0)
Definition Mat.h:465
Mat & scalarMinusEqFromLeft(const EE &ee)
Definition Mat.h:891
EScalar Scalar
Definition Mat.h:187
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
TStandard standardize() const
Definition Mat.h:234
EScalarNormSq ScalarNormSq
Definition Mat.h:192
TReal & real()
Definition Mat.h:745
static const Mat & getAs(const ELT *p)
Definition Mat.h:1081
EStandard trace() const
Definition Mat.h:809
CNT< E >::TAbs EAbs
Definition Mat.h:111
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
CNT< ScalarNormSq >::TSqrt norm() const
Definition Mat.h:692
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, negator< EE >, CSS, RSS > &mm)
Definition Mat.h:566
Mat(const TCol &r0, const TCol &r1, const TCol &r2, const TCol &r3, const TCol &r4, const TCol &r5)
Definition Mat.h:502
Mat< M, N, E, M, 1 > TPacked
Definition Mat.h:179
Vec< M, E, RS > TCol
Definition Mat.h:163
EltResult< E2 >::Mul elementwiseMultiply(const Mat< M, N, E2, CS2, RS2 > &r) const
Definition Mat.h:611
const TRow & operator[](int i) const
Definition Mat.h:681
Mat(const Vec< M, EE, SS > &r0)
Definition Mat.h:508
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
TPosTrans & updPositionalTranspose()
Definition Mat.h:737
const E & operator()(int i, int j) const
Definition Mat.h:686
CNT< E >::ScalarNormSq EScalarNormSq
Definition Mat.h:121
Mat< N, M, E, RS, CS > TPosTrans
Definition Mat.h:160
ELT E
Definition Mat.h:99
static int ncol()
Return the number of columns in this Mat, echoing the value supplied for the template parameter N.
Definition Mat.h:203
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
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 Row< N, EE, SS > &r0, const Row< N, EE, SS > &r1, const Row< N, EE, SS > &r2)
Definition Mat.h:470
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
CNT< E >::Precision EPrecision
Definition Mat.h:120
Mat< M, N, typename CNT< EE >::template Result< E >::Sub > scalarSubtractFromLeft(const EE &e) const
Definition Mat.h:863
CNT< E >::TPosTrans EPosTrans
Definition Mat.h:106
Mat & scalarPlusEq(const EE &ee)
Definition Mat.h:885
CNT< E >::TNormalize ENormalize
Definition Mat.h:114
const TPosTrans & positionalTranspose() const
Definition Mat.h:735
SymMat< M, ESqTHerm > TSqTHerm
Definition Mat.h:175
Mat & scalarDivideEq(const EE &ee)
Definition Mat.h:899
Row< N, E, CS > TRow
Definition Mat.h:162
Mat< M, N, E, CS, RS > T
Definition Mat.h:150
Mat(const TRow &r0, const TRow &r1, const TRow &r2)
Definition Mat.h:450
Mat< M, N, ENormalize, M, 1 > TNormalize
Definition Mat.h:172
Mat & operator+=(const Mat< M, N, EE, CSS, RSS > &mm)
Definition Mat.h:550
Mat(const TCol &r0, const TCol &r1, const TCol &r2, const TCol &r3, const TCol &r4)
Definition Mat.h:498
const TRow & row(int i) const
Definition Mat.h:768
Mat(const E &e0, const E &e1)
Definition Mat.h:380
Mat(const E &e0, const E &e1, const E &e2, const E &e3, const E &e4, const E &e5)
Definition Mat.h:388
TRow colSum() const
Returns a row vector (Row) containing the column sums of this matrix.
Definition Mat.h:1199
TCol & col(int j)
Definition Mat.h:781
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
Mat(const TRow &r0, const TRow &r1, const TRow &r2, const TRow &r3, const TRow &r4, const TRow &r5)
Definition Mat.h:459
Mat & scalarEq(const EE &ee)
Definition Mat.h:880
Result< Mat< M, N, E2, CS2, RS2 > >::Add conformingAdd(const Mat< M, N, E2, CS2, RS2 > &r) const
Definition Mat.h:588
TRow & row(int i)
Definition Mat.h:772
TNeg & updNegate()
Definition Mat.h:730
CNT< E >::TSqrt ESqrt
Definition Mat.h:110
Mat< M, N, typename CNT< E >::template Result< EE >::Mul > scalarMultiply(const EE &e) const
Definition Mat.h:819
Mat & operator*=(const Mat< N, N, EE, CSS, RSS > &mm)
Definition Mat.h:574
ENumber Number
Definition Mat.h:189
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
TNormalize normalize() const
Definition Mat.h:707
static int size()
Return the total number of elements M*N contained in this Mat.
Definition Mat.h:197
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 & operator=(const Mat< M, N, EE, CSS, RSS > &mm)
Definition Mat.h:538
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
CNT< E >::TInvert EInvert
Definition Mat.h:113
const TNeg & negate() const
Definition Mat.h:729
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:1227
E & operator()(int i, int j)
Definition Mat.h:687
Mat & operator-=(const EE &e)
Definition Mat.h:874
Mat< M+1, N, E, M+1, 1 > TAppendRow
Definition Mat.h:183
CNT< E >::Scalar EScalar
Definition Mat.h:116
bool isFinite() const
Return true if no element contains an Infinity or a NaN.
Definition Mat.h:1114
Vec< MinDim, E, RS+CS > TDiag
Definition Mat.h:164
Mat< M2, M, E2, CS2, RS2 >::template Result< Mat >::Mul conformingMultiplyFromLeft(const Mat< M2, M, E2, CS2, RS2 > &m) const
Definition Mat.h:664
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
CNT< E >::TReal EReal
Definition Mat.h:102
TDiag & diag()
This non-const version of diag() is an alternate name for updDiag() available for historical reasons.
Definition Mat.h:807
bool isNumericallyEqual(const ELT &e, double tol=getDefaultTolerance()) const
Test whether this is numerically a "scalar" matrix, meaning that it is a diagonal matrix in which eac...
Definition Mat.h:1150
@ NPackedElements
Definition Mat.h:132
@ ImagOffset
Definition Mat.h:135
@ ArgDepth
Definition Mat.h:138
@ NCols
Definition Mat.h:127
@ RowSpacing
Definition Mat.h:130
@ RealStrideFactor
Definition Mat.h:136
@ NActualElements
Definition Mat.h:133
@ IsPrecision
Definition Mat.h:145
@ IsScalar
Definition Mat.h:141
@ IsNumber
Definition Mat.h:143
@ ColSpacing
Definition Mat.h:131
@ SignInterpretation
Definition Mat.h:146
@ IsStdNumber
Definition Mat.h:144
@ NActualScalars
Definition Mat.h:134
@ NRows
Definition Mat.h:126
@ MinDim
Definition Mat.h:128
@ IsULessScalar
Definition Mat.h:142
Result< SymMat< M, E2, RS2 > >::Add conformingAdd(const SymMat< M, E2, RS2 > &sy) const
Definition Mat.h:632
const SubMat< MM, NN >::Type & getSubMat(int i, int j) const
Definition Mat.h:922
TNeg & operator-()
Definition Mat.h:725
Mat(const EE *p)
Definition Mat.h:532
static Mat & updAs(ELT *p)
Definition Mat.h:1082
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
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
static Mat< M, N, ELT, M, 1 > getNaN()
Definition Mat.h:1085
Mat< M-1, N, E, M-1, 1 > TDropRow
Definition Mat.h:180
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
THerm & updTranspose()
Definition Mat.h:733
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:360
const ELT & get(int i, int j) const
Variant of indexing operator that's scripting friendly to get entry (i, j)
Definition Mat.h:1225
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 >::TNeg ENeg
Definition Mat.h:100
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:376
Mat< N, M, EInvert, N, 1 > TInvert
Definition Mat.h:171
CNT< E >::THerm EHerm
Definition Mat.h:105
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
const TCol & col(int j) const
Definition Mat.h:777
std::string toString() const
toString() returns a string representation of the Mat.
Definition Mat.h:1219
bool isNaN() const
Return true if any element of this Mat contains a NaN anywhere.
Definition Mat.h:1092
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
TCol rowSum() const
Returns a column vector (Vec) containing the row sums of this matrix.
Definition Mat.h:1210
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 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
EULessScalar ULessScalar
Definition Mat.h:188
Mat & operator+=(const EE &e)
Definition Mat.h:873
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
Mat< M, N-1, E, M, 1 > TDropCol
Definition Mat.h:181
Mat & operator*=(const EE &e)
Definition Mat.h:875
Mat< M, N, EReal, CS *CNT< E >::RealStrideFactor, RS *CNT< E >::RealStrideFactor > TReal
Definition Mat.h:155
Mat< M, N, ESqrt, M, 1 > TSqrt
Definition Mat.h:168
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
SymMat< M, E2, RS2 >::template Result< Mat >::Sub conformingSubtractFromLeft(const SymMat< M, E2, RS2 > &sy) const
Definition Mat.h:646
const THerm & operator~() const
Definition Mat.h:726
Mat & operator/=(const EE &e)
Definition Mat.h:876
TDropCol dropCol(int j) const
Return a matrix one column smaller than this one by dropping column j.
Definition Mat.h:954
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
const TWithoutNegator & castAwayNegatorIfAny() const
Definition Mat.h:765
Mat< M+1, N+1, E, M+1, 1 > TAppendRowCol
Definition Mat.h:185
Mat< M, N, EComplex, CS, RS > TComplex
Definition Mat.h:158
CNT< E >::ULessScalar EULessScalar
Definition Mat.h:117
TDiag & updDiag()
Select main diagonal (of largest leading square if rectangular) and return it as a writable view (as ...
Definition Mat.h:804
Mat(const TCol &r0, const TCol &r1, const TCol &r2)
Definition Mat.h:493
Result< Mat< M, N, E2, CS2, RS2 > >::Sub conformingSubtract(const Mat< M, N, E2, CS2, RS2 > &r) const
Definition Mat.h:596
Mat(const E &e0, const E &e1, const E &e2, const E &e3, const E &e4)
Definition Mat.h:386
const Mat & operator+() const
Definition Mat.h:723
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
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
Mat(const TRow &r0)
Definition Mat.h:446
Mat & operator-=(const Mat< M, N, EE, CSS, RSS > &mm)
Definition Mat.h:561
Mat(const TRow &r0, const TRow &r1)
Definition Mat.h:448
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
CNT< E >::StdNumber EStdNumber
Definition Mat.h:119
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< M, N, EImag, CS *CNT< E >::RealStrideFactor, RS *CNT< E >::RealStrideFactor > TImag
Definition Mat.h:157
Mat(const TCol &r0, const TCol &r1)
Definition Mat.h:491
TImag & imag()
Definition Mat.h:759
const TImag & imag() const
Definition Mat.h:754
Mat & scalarTimesEq(const EE &ee)
Definition Mat.h:894
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
Mat & scalarMinusEq(const EE &ee)
Definition Mat.h:888
Mat< M, N, EWithoutNegator, CS, RS > TWithoutNegator
Definition Mat.h:152
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
Mat & operator+=(const Mat< M, N, negator< EE >, CSS, RSS > &mm)
Definition Mat.h:555
void setToZero()
Definition Mat.h:909
Mat(const Vec< M, EE, SS > &r0, const Vec< M, EE, SS > &r1)
Definition Mat.h:510
Mat< M, N, typename CNT< E >::template Result< EE >::Add > scalarAdd(const EE &e) const
Definition Mat.h:848
Mat< N, M, EHerm, RS, CS > THerm
Definition Mat.h:159
Mat(const Vec< M, EE, SS > &r0, const Vec< M, EE, SS > &r1, const Vec< M, EE, SS > &r2)
Definition Mat.h:513
EPrecision Precision
Definition Mat.h:191
E TElement
Definition Mat.h:161
TCol & operator()(int j)
Definition Mat.h:684
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
Mat(const E &e0, const E &e1, const E &e2)
Definition Mat.h:382
Mat & operator=(const EE *p)
Definition Mat.h:543
const TNeg & operator-() const
Definition Mat.h:724
Mat< M, N, typename CNT< E >::template Result< EE >::Dvd > scalarDivide(const EE &e) const
Definition Mat.h:834
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
Mat(const TRow &r0, const TRow &r1, const TRow &r2, const TRow &r3)
Definition Mat.h:452
THerm & operator~()
Definition Mat.h:727
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
CNT< E >::TSqTHerm ESqTHerm
Definition Mat.h:108
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,...
Definition Mat.h:1139
Definition NTraits.h:436
This is a fixed-length row vector designed for no-overhead inline computation.
Definition Row.h:132
EStandard sum() const
Definition Row.h:254
This is a small, fixed-size symmetric or Hermitian matrix designed for no-overhead inline computation...
Definition SymMat.h:87
This is a fixed-length column vector designed for no-overhead inline computation.
Definition Vec.h:184
Vec & scalarMinusEq(const EE &ee)
Definition Vec.h:787
Vec & scalarPlusEq(const EE &ee)
Definition Vec.h:785
Vec< M, typename CNT< E >::template Result< EE >::Add > scalarAdd(const EE &e) const
Definition Vec.h:752
Vec & scalarEq(const EE &ee)
Definition Vec.h:783
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:366
SimTK::conjugate<R> should be instantiated only for float, double.
Definition conjugate.h:178
negator<N>, where N is a number type (real, complex, conjugate), is represented in memory identically...
Definition negator.h:75
NTraits< N >::StdNumber StdNumber
Definition negator.h:107
This is the top-level SimTK namespace into which all SimTK names are placed to avoid collision with o...
Definition Assembler.h:37
Matrix_< E > operator*(const MatrixBase< E > &l, const typename CNT< E >::StdNumber &r)
Definition BigMatrix.h:605
Matrix_< E > operator/(const MatrixBase< E > &l, const typename CNT< E >::StdNumber &r)
Definition BigMatrix.h:613
std::basic_istream< CHAR, TRAITS > & operator>>(std::basic_istream< CHAR, TRAITS > &is, conjugate< R > &c)
Definition conjugate.h:505
@ MAX_RESOLVED_DEPTH
Definition CompositeNumericalTypes.h:120
Matrix_< typename CNT< E1 >::template Result< E2 >::Add > operator+(const MatrixBase< E1 > &l, const MatrixBase< E2 > &r)
Definition BigMatrix.h:568
std::ostream & operator<<(std::ostream &o, const ContactForce &f)
Definition CompliantContactSubsystem.h:387
Matrix_< typename CNT< E1 >::template Result< E2 >::Sub > operator-(const MatrixBase< E1 > &l, const MatrixBase< E2 > &r)
Definition BigMatrix.h:584
bool operator==(const PhiMatrix &p1, const PhiMatrix &p2)
Definition SpatialAlgebra.h:791
bool operator!=(const L &left, const R &right)
Definition SimTKcommon/include/SimTKcommon/internal/common.h:641
Definition Mat.h:244
Mat< M, N, typename CNT< E >::template Result< P >::Add, M, 1 > Add
Definition Mat.h:247
Mat< M, N, typename CNT< E >::template Result< P >::Sub, M, 1 > Sub
Definition Mat.h:248
Mat< M, N, typename CNT< E >::template Result< P >::Mul, M, 1 > Mul
Definition Mat.h:245
Mat< M, N, typename CNT< E >::template Result< P >::Dvd, M, 1 > Dvd
Definition Mat.h:246
Definition Mat.h:253
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
DvdOp::Type Dvd
Definition Mat.h:267
MulOpNonConforming::Type MulNon
Definition Mat.h:262
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
MulOp::Type Mul
Definition Mat.h:257
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
SubOp::Type Sub
Definition Mat.h:277
AddOp::Type Add
Definition Mat.h:272
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
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
Definition Mat.h:917
Mat< MM, NN, ELT, CS, RS > Type
Definition Mat.h:918
Definition Mat.h:281
Mat< M, N, P > Type
Definition Mat.h:282