Simbody  3.4
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
MatrixBase.h
Go to the documentation of this file.
1 #ifndef SimTK_SIMMATRIX_MATRIXBASE_H_
2 #define SimTK_SIMMATRIX_MATRIXBASE_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-13 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 namespace SimTK {
32 
33 //==============================================================================
34 // MATRIX BASE
35 //==============================================================================
67 // ----------------------------------------------------------------------------
68 template <class ELT> class MatrixBase {
69 public:
70  // These typedefs are handy, but despite appearances you cannot
71  // treat a MatrixBase as a composite numerical type. That is,
72  // CNT<MatrixBase> will not compile, or if it does it won't be
73  // meaningful.
74 
75  typedef ELT E;
76  typedef typename CNT<E>::TNeg ENeg;
77  typedef typename CNT<E>::TWithoutNegator EWithoutNegator;
78  typedef typename CNT<E>::TReal EReal;
79  typedef typename CNT<E>::TImag EImag;
80  typedef typename CNT<E>::TComplex EComplex;
81  typedef typename CNT<E>::THerm EHerm;
82  typedef typename CNT<E>::TPosTrans EPosTrans;
83 
84  typedef typename CNT<E>::TAbs EAbs;
85  typedef typename CNT<E>::TStandard EStandard;
86  typedef typename CNT<E>::TInvert EInvert;
87  typedef typename CNT<E>::TNormalize ENormalize;
88  typedef typename CNT<E>::TSqHermT ESqHermT;
89  typedef typename CNT<E>::TSqTHerm ESqTHerm;
90 
91  typedef typename CNT<E>::Scalar EScalar;
92  typedef typename CNT<E>::Number ENumber;
93  typedef typename CNT<E>::StdNumber EStdNumber;
94  typedef typename CNT<E>::Precision EPrecision;
95  typedef typename CNT<E>::ScalarNormSq EScalarNormSq;
96 
97  typedef EScalar Scalar; // the underlying Scalar type
98  typedef ENumber Number; // negator removed from Scalar
99  typedef EStdNumber StdNumber; // conjugate goes to complex
100  typedef EPrecision Precision; // complex removed from StdNumber
101  typedef EScalarNormSq ScalarNormSq; // type of scalar^2
102 
103  typedef MatrixBase<E> T;
104  typedef MatrixBase<ENeg> TNeg;
105  typedef MatrixBase<EWithoutNegator> TWithoutNegator;
106  typedef MatrixBase<EReal> TReal;
107  typedef MatrixBase<EImag> TImag;
108  typedef MatrixBase<EComplex> TComplex;
109  typedef MatrixBase<EHerm> THerm;
110  typedef MatrixBase<E> TPosTrans;
111 
112  typedef MatrixBase<EAbs> TAbs;
113  typedef MatrixBase<EStandard> TStandard;
114  typedef MatrixBase<EInvert> TInvert;
115  typedef MatrixBase<ENormalize> TNormalize;
116  typedef MatrixBase<ESqHermT> TSqHermT; // ~Mat*Mat
117  typedef MatrixBase<ESqTHerm> TSqTHerm; // Mat*~Mat
118 
119  const MatrixCommitment& getCharacterCommitment() const {return helper.getCharacterCommitment();}
120  const MatrixCharacter& getMatrixCharacter() const {return helper.getMatrixCharacter();}
121 
124  void commitTo(const MatrixCommitment& mc)
125  { helper.commitTo(mc); }
126 
127  // This gives the resulting matrix type when (m(i,j) op P) is applied to each element.
128  // It will have element types which are the regular composite result of E op P.
129  template <class P> struct EltResult {
130  typedef MatrixBase<typename CNT<E>::template Result<P>::Mul> Mul;
131  typedef MatrixBase<typename CNT<E>::template Result<P>::Dvd> Dvd;
132  typedef MatrixBase<typename CNT<E>::template Result<P>::Add> Add;
133  typedef MatrixBase<typename CNT<E>::template Result<P>::Sub> Sub;
134  };
135 
137  int nrow() const {return helper.nrow();}
139  int ncol() const {return helper.ncol();}
140 
148  ptrdiff_t nelt() const {return helper.nelt();}
149 
151  bool isResizeable() const {return getCharacterCommitment().isResizeable();}
152 
153  enum {
154  NScalarsPerElement = CNT<E>::NActualScalars,
155  CppNScalarsPerElement = sizeof(E) / sizeof(Scalar)
156  };
157 
161  MatrixBase() : helper(NScalarsPerElement,CppNScalarsPerElement) {}
162 
165  MatrixBase(int m, int n)
166  : helper(NScalarsPerElement,CppNScalarsPerElement,MatrixCommitment(),m,n) {}
167 
173  explicit MatrixBase(const MatrixCommitment& commitment)
174  : helper(NScalarsPerElement,CppNScalarsPerElement,commitment) {}
175 
176 
181  MatrixBase(const MatrixCommitment& commitment, int m, int n)
182  : helper(NScalarsPerElement,CppNScalarsPerElement,commitment,m,n) {}
183 
185  MatrixBase(const MatrixBase& b)
186  : helper(b.helper.getCharacterCommitment(),
187  b.helper, typename MatrixHelper<Scalar>::DeepCopy()) { }
188 
191  MatrixBase(const TNeg& b)
192  : helper(b.helper.getCharacterCommitment(),
193  b.helper, typename MatrixHelper<Scalar>::DeepCopy()) { }
194 
197  MatrixBase& copyAssign(const MatrixBase& b) {
198  helper.copyAssign(b.helper);
199  return *this;
200  }
201  MatrixBase& operator=(const MatrixBase& b) { return copyAssign(b); }
202 
203 
212  MatrixBase& viewAssign(const MatrixBase& src) {
213  helper.writableViewAssign(const_cast<MatrixHelper<Scalar>&>(src.helper));
214  return *this;
215  }
216 
217  // default destructor
218 
225  MatrixBase(const MatrixCommitment& commitment, int m, int n, const ELT& initialValue)
226  : helper(NScalarsPerElement, CppNScalarsPerElement, commitment, m, n)
227  { helper.fillWith(reinterpret_cast<const Scalar*>(&initialValue)); }
228 
239  MatrixBase(const MatrixCommitment& commitment, int m, int n,
240  const ELT* cppInitialValuesByRow)
241  : helper(NScalarsPerElement, CppNScalarsPerElement, commitment, m, n)
242  { helper.copyInByRowsFromCpp(reinterpret_cast<const Scalar*>(cppInitialValuesByRow)); }
243 
255 
257  MatrixBase(const MatrixCommitment& commitment,
258  const MatrixCharacter& character,
259  int spacing, const Scalar* data) // read only data
260  : helper(NScalarsPerElement, CppNScalarsPerElement,
261  commitment, character, spacing, data) {}
262 
264  MatrixBase(const MatrixCommitment& commitment,
265  const MatrixCharacter& character,
266  int spacing, Scalar* data) // writable data
267  : helper(NScalarsPerElement, CppNScalarsPerElement,
268  commitment, character, spacing, data) {}
270 
271  // Create a new MatrixBase from an existing helper. Both shallow and deep copies are possible.
272  MatrixBase(const MatrixCommitment& commitment,
273  MatrixHelper<Scalar>& source,
274  const typename MatrixHelper<Scalar>::ShallowCopy& shallow)
275  : helper(commitment, source, shallow) {}
276  MatrixBase(const MatrixCommitment& commitment,
277  const MatrixHelper<Scalar>& source,
278  const typename MatrixHelper<Scalar>::ShallowCopy& shallow)
279  : helper(commitment, source, shallow) {}
280  MatrixBase(const MatrixCommitment& commitment,
281  const MatrixHelper<Scalar>& source,
282  const typename MatrixHelper<Scalar>::DeepCopy& deep)
283  : helper(commitment, source, deep) {}
284 
288  void clear() {helper.clear();}
289 
290  MatrixBase& operator*=(const StdNumber& t) { helper.scaleBy(t); return *this; }
291  MatrixBase& operator/=(const StdNumber& t) { helper.scaleBy(StdNumber(1)/t); return *this; }
292  MatrixBase& operator+=(const MatrixBase& r) { helper.addIn(r.helper); return *this; }
293  MatrixBase& operator-=(const MatrixBase& r) { helper.subIn(r.helper); return *this; }
294 
295  template <class EE> MatrixBase(const MatrixBase<EE>& b)
296  : helper(MatrixCommitment(),b.helper, typename MatrixHelper<Scalar>::DeepCopy()) { }
297 
298  template <class EE> MatrixBase& operator=(const MatrixBase<EE>& b)
299  { helper = b.helper; return *this; }
300  template <class EE> MatrixBase& operator+=(const MatrixBase<EE>& b)
301  { helper.addIn(b.helper); return *this; }
302  template <class EE> MatrixBase& operator-=(const MatrixBase<EE>& b)
303  { helper.subIn(b.helper); return *this; }
304 
315  MatrixBase& operator=(const ELT& t) {
316  setToZero(); updDiag().setTo(t);
317  return *this;
318  }
319 
325  template <class S> inline MatrixBase&
326  scalarAssign(const S& s) {
327  setToZero(); updDiag().setTo(s);
328  return *this;
329  }
330 
334  template <class S> inline MatrixBase&
335  scalarAddInPlace(const S& s) {
336  updDiag().elementwiseAddScalarInPlace(s);
337  }
338 
339 
343  template <class S> inline MatrixBase&
344  scalarSubtractInPlace(const S& s) {
345  updDiag().elementwiseSubtractScalarInPlace(s);
346  }
347 
350  template <class S> inline MatrixBase&
351  scalarSubtractFromLeftInPlace(const S& s) {
352  negateInPlace();
353  updDiag().elementwiseAddScalarInPlace(s); // yes, add
354  }
355 
362  template <class S> inline MatrixBase&
363  scalarMultiplyInPlace(const S&);
364 
368  template <class S> inline MatrixBase&
369  scalarMultiplyFromLeftInPlace(const S&);
370 
377  template <class S> inline MatrixBase&
378  scalarDivideInPlace(const S&);
379 
383  template <class S> inline MatrixBase&
384  scalarDivideFromLeftInPlace(const S&);
385 
386 
389  template <class EE> inline MatrixBase&
390  rowScaleInPlace(const VectorBase<EE>&);
391 
394  template <class EE> inline void
395  rowScale(const VectorBase<EE>& r, typename EltResult<EE>::Mul& out) const;
396 
397  template <class EE> inline typename EltResult<EE>::Mul
398  rowScale(const VectorBase<EE>& r) const {
399  typename EltResult<EE>::Mul out(nrow(), ncol()); rowScale(r,out); return out;
400  }
401 
404  template <class EE> inline MatrixBase&
405  colScaleInPlace(const VectorBase<EE>&);
406 
407  template <class EE> inline void
408  colScale(const VectorBase<EE>& c, typename EltResult<EE>::Mul& out) const;
409 
410  template <class EE> inline typename EltResult<EE>::Mul
411  colScale(const VectorBase<EE>& c) const {
412  typename EltResult<EE>::Mul out(nrow(), ncol()); colScale(c,out); return out;
413  }
414 
415 
420  template <class ER, class EC> inline MatrixBase&
421  rowAndColScaleInPlace(const VectorBase<ER>& r, const VectorBase<EC>& c);
422 
423  template <class ER, class EC> inline void
424  rowAndColScale(const VectorBase<ER>& r, const VectorBase<EC>& c,
425  typename EltResult<typename VectorBase<ER>::template EltResult<EC>::Mul>::Mul& out) const;
426 
427  template <class ER, class EC> inline typename EltResult<typename VectorBase<ER>::template EltResult<EC>::Mul>::Mul
428  rowAndColScale(const VectorBase<ER>& r, const VectorBase<EC>& c) const {
429  typename EltResult<typename VectorBase<ER>::template EltResult<EC>::Mul>::Mul
430  out(nrow(), ncol());
431  rowAndColScale(r,c,out); return out;
432  }
433 
441  template <class S> inline MatrixBase&
442  elementwiseAssign(const S& s);
443 
445  MatrixBase& elementwiseAssign(int s)
446  { return elementwiseAssign<Real>(Real(s)); }
447 
449  MatrixBase& elementwiseInvertInPlace();
450 
451  void elementwiseInvert(MatrixBase<typename CNT<E>::TInvert>& out) const;
452 
453  MatrixBase<typename CNT<E>::TInvert> elementwiseInvert() const {
454  MatrixBase<typename CNT<E>::TInvert> out(nrow(), ncol());
455  elementwiseInvert(out);
456  return out;
457  }
458 
466  template <class S> inline MatrixBase&
467  elementwiseAddScalarInPlace(const S& s);
468 
469  template <class S> inline void
470  elementwiseAddScalar(const S& s, typename EltResult<S>::Add&) const;
471 
472  template <class S> inline typename EltResult<S>::Add
473  elementwiseAddScalar(const S& s) const {
474  typename EltResult<S>::Add out(nrow(), ncol());
475  elementwiseAddScalar(s,out);
476  return out;
477  }
478 
486  template <class S> inline MatrixBase&
487  elementwiseSubtractScalarInPlace(const S& s);
488 
489  template <class S> inline void
490  elementwiseSubtractScalar(const S& s, typename EltResult<S>::Sub&) const;
491 
492  template <class S> inline typename EltResult<S>::Sub
493  elementwiseSubtractScalar(const S& s) const {
494  typename EltResult<S>::Sub out(nrow(), ncol());
495  elementwiseSubtractScalar(s,out);
496  return out;
497  }
498 
507  template <class S> inline MatrixBase&
508  elementwiseSubtractFromScalarInPlace(const S& s);
509 
510  template <class S> inline void
511  elementwiseSubtractFromScalar(
512  const S&,
513  typename MatrixBase<S>::template EltResult<E>::Sub&) const;
514 
515  template <class S> inline typename MatrixBase<S>::template EltResult<E>::Sub
516  elementwiseSubtractFromScalar(const S& s) const {
517  typename MatrixBase<S>::template EltResult<E>::Sub out(nrow(), ncol());
518  elementwiseSubtractFromScalar<S>(s,out);
519  return out;
520  }
521 
523  template <class EE> inline MatrixBase&
524  elementwiseMultiplyInPlace(const MatrixBase<EE>&);
525 
526  template <class EE> inline void
527  elementwiseMultiply(const MatrixBase<EE>&, typename EltResult<EE>::Mul&) const;
528 
529  template <class EE> inline typename EltResult<EE>::Mul
530  elementwiseMultiply(const MatrixBase<EE>& m) const {
531  typename EltResult<EE>::Mul out(nrow(), ncol());
532  elementwiseMultiply<EE>(m,out);
533  return out;
534  }
535 
537  template <class EE> inline MatrixBase&
538  elementwiseMultiplyFromLeftInPlace(const MatrixBase<EE>&);
539 
540  template <class EE> inline void
541  elementwiseMultiplyFromLeft(
542  const MatrixBase<EE>&,
543  typename MatrixBase<EE>::template EltResult<E>::Mul&) const;
544 
545  template <class EE> inline typename MatrixBase<EE>::template EltResult<E>::Mul
546  elementwiseMultiplyFromLeft(const MatrixBase<EE>& m) const {
547  typename EltResult<EE>::Mul out(nrow(), ncol());
548  elementwiseMultiplyFromLeft<EE>(m,out);
549  return out;
550  }
551 
553  template <class EE> inline MatrixBase&
554  elementwiseDivideInPlace(const MatrixBase<EE>&);
555 
556  template <class EE> inline void
557  elementwiseDivide(const MatrixBase<EE>&, typename EltResult<EE>::Dvd&) const;
558 
559  template <class EE> inline typename EltResult<EE>::Dvd
560  elementwiseDivide(const MatrixBase<EE>& m) const {
561  typename EltResult<EE>::Dvd out(nrow(), ncol());
562  elementwiseDivide<EE>(m,out);
563  return out;
564  }
565 
567  template <class EE> inline MatrixBase&
568  elementwiseDivideFromLeftInPlace(const MatrixBase<EE>&);
569 
570  template <class EE> inline void
571  elementwiseDivideFromLeft(
572  const MatrixBase<EE>&,
573  typename MatrixBase<EE>::template EltResult<E>::Dvd&) const;
574 
575  template <class EE> inline typename MatrixBase<EE>::template EltResult<EE>::Dvd
576  elementwiseDivideFromLeft(const MatrixBase<EE>& m) const {
577  typename MatrixBase<EE>::template EltResult<E>::Dvd out(nrow(), ncol());
578  elementwiseDivideFromLeft<EE>(m,out);
579  return out;
580  }
581 
583  MatrixBase& setTo(const ELT& t) {helper.fillWith(reinterpret_cast<const Scalar*>(&t)); return *this;}
584  MatrixBase& setToNaN() {helper.fillWithScalar(CNT<StdNumber>::getNaN()); return *this;}
585  MatrixBase& setToZero() {helper.fillWithScalar(StdNumber(0)); return *this;}
586 
587  // View creating operators.
588  inline RowVectorView_<ELT> row(int i) const; // select a row
589  inline RowVectorView_<ELT> updRow(int i);
590  inline VectorView_<ELT> col(int j) const; // select a column
591  inline VectorView_<ELT> updCol(int j);
592 
593  RowVectorView_<ELT> operator[](int i) const {return row(i);}
594  RowVectorView_<ELT> operator[](int i) {return updRow(i);}
595  VectorView_<ELT> operator()(int j) const {return col(j);}
596  VectorView_<ELT> operator()(int j) {return updCol(j);}
597 
598  // Select a block.
599  inline MatrixView_<ELT> block(int i, int j, int m, int n) const;
600  inline MatrixView_<ELT> updBlock(int i, int j, int m, int n);
601 
602  MatrixView_<ELT> operator()(int i, int j, int m, int n) const
603  { return block(i,j,m,n); }
604  MatrixView_<ELT> operator()(int i, int j, int m, int n)
605  { return updBlock(i,j,m,n); }
606 
607  // Hermitian transpose.
608  inline MatrixView_<EHerm> transpose() const;
609  inline MatrixView_<EHerm> updTranspose();
610 
611  MatrixView_<EHerm> operator~() const {return transpose();}
612  MatrixView_<EHerm> operator~() {return updTranspose();}
613 
616  inline VectorView_<ELT> diag() const;
619  inline VectorView_<ELT> updDiag();
622  VectorView_<ELT> diag() {return updDiag();}
623 
624  // Create a view of the real or imaginary elements. TODO
625  //inline MatrixView_<EReal> real() const;
626  //inline MatrixView_<EReal> updReal();
627  //inline MatrixView_<EImag> imag() const;
628  //inline MatrixView_<EImag> updImag();
629 
630  // Overload "real" and "imag" for both read and write as a nod to convention. TODO
631  //MatrixView_<EReal> real() {return updReal();}
632  //MatrixView_<EReal> imag() {return updImag();}
633 
634  // TODO: this routine seems ill-advised but I need it for the IVM port at the moment
635  TInvert invert() const { // return a newly-allocated inverse; dump negator
636  TInvert m(*this);
637  m.helper.invertInPlace();
638  return m; // TODO - bad: makes an extra copy
639  }
640 
641  void invertInPlace() {helper.invertInPlace();}
642 
644  void dump(const char* msg=0) const {
645  helper.dump(msg);
646  }
647 
648 
649  // This routine is useful for implementing friendlier Matrix expressions and operators.
650  // It maps closely to the Level-3 BLAS family of pxxmm() routines like sgemm(). The
651  // operation performed assumes that "this" is the result, and that "this" has
652  // already been sized correctly to receive the result. We'll compute
653  // this = beta*this + alpha*A*B
654  // If beta is 0 then "this" can be uninitialized. If alpha is 0 we promise not
655  // to look at A or B. The routine can work efficiently even if A and/or B are transposed
656  // by their views, so an expression like
657  // C += s * ~A * ~B
658  // can be performed with the single equivalent call
659  // C.matmul(1., s, Tr(A), Tr(B))
660  // where Tr(A) indicates a transposed view of the original A.
661  // The ultimate efficiency of this call depends on the data layout and views which
662  // are used for the three matrices.
663  // NOTE: neither A nor B can be the same matrix as 'this', nor views of the same data
664  // which would expose elements of 'this' that will be modified by this operation.
665  template <class ELT_A, class ELT_B>
666  MatrixBase& matmul(const StdNumber& beta, // applied to 'this'
667  const StdNumber& alpha, const MatrixBase<ELT_A>& A, const MatrixBase<ELT_B>& B)
668  {
669  helper.matmul(beta,alpha,A.helper,B.helper);
670  return *this;
671  }
672 
681  const ELT& getElt(int i, int j) const { return *reinterpret_cast<const ELT*>(helper.getElt(i,j)); }
682  ELT& updElt(int i, int j) { return *reinterpret_cast< ELT*>(helper.updElt(i,j)); }
683 
684  const ELT& operator()(int i, int j) const {return getElt(i,j);}
685  ELT& operator()(int i, int j) {return updElt(i,j);}
686 
691  void getAnyElt(int i, int j, ELT& value) const
692  { helper.getAnyElt(i,j,reinterpret_cast<Scalar*>(&value)); }
693  ELT getAnyElt(int i, int j) const {ELT e; getAnyElt(i,j,e); return e;}
694 
697  // TODO: very slow! Should be optimized at least for the case
698  // where ELT is a Scalar.
699  ScalarNormSq scalarNormSqr() const {
700  const int nr=nrow(), nc=ncol();
701  ScalarNormSq sum(0);
702  for(int j=0;j<nc;++j)
703  for (int i=0; i<nr; ++i)
704  sum += CNT<E>::scalarNormSqr((*this)(i,j));
705  return sum;
706  }
707 
711  // TODO: very slow! Should be optimized at least for the case
712  // where ELT is a Scalar.
713  void abs(TAbs& mabs) const {
714  const int nr=nrow(), nc=ncol();
715  mabs.resize(nr,nc);
716  for(int j=0;j<nc;++j)
717  for (int i=0; i<nr; ++i)
718  mabs(i,j) = CNT<E>::abs((*this)(i,j));
719  }
720 
724  TAbs abs() const { TAbs mabs; abs(mabs); return mabs; }
725 
736  TStandard standardize() const {
737  const int nr=nrow(), nc=ncol();
738  TStandard mstd(nr, nc);
739  for(int j=0;j<nc;++j)
740  for (int i=0; i<nr; ++i)
741  mstd(i,j) = CNT<E>::standardize((*this)(i,j));
742  return mstd;
743  }
744 
748  ScalarNormSq normSqr() const { return scalarNormSqr(); }
749  // TODO -- not good; unnecessary overflow
750  typename CNT<ScalarNormSq>::TSqrt
751  norm() const { return CNT<ScalarNormSq>::sqrt(scalarNormSqr()); }
752 
756  typename CNT<ScalarNormSq>::TSqrt
757  normRMS() const {
758  if (!CNT<ELT>::IsScalar)
759  SimTK_THROW1(Exception::Cant, "normRMS() only defined for scalar elements");
760  if (nelt() == 0)
761  return typename CNT<ScalarNormSq>::TSqrt(0);
762  return CNT<ScalarNormSq>::sqrt(scalarNormSqr()/nelt());
763  }
764 
766  RowVector_<ELT> colSum() const {
767  const int cols = ncol();
768  RowVector_<ELT> row(cols);
769  for (int j = 0; j < cols; ++j)
770  helper.colSum(j, reinterpret_cast<Scalar*>(&row[j]));
771  return row;
772  }
774  RowVector_<ELT> sum() const {return colSum();}
775 
777  Vector_<ELT> rowSum() const {
778  const int rows = nrow();
779  Vector_<ELT> col(rows);
780  for (int i = 0; i < rows; ++i)
781  helper.rowSum(i, reinterpret_cast<Scalar*>(&col[i]));
782  return col;
783  }
784 
785  //TODO: make unary +/- return a self-view so they won't reallocate?
786  const MatrixBase& operator+() const {return *this; }
787  const TNeg& negate() const {return *reinterpret_cast<const TNeg*>(this); }
788  TNeg& updNegate() {return *reinterpret_cast<TNeg*>(this); }
789 
790  const TNeg& operator-() const {return negate();}
791  TNeg& operator-() {return updNegate();}
792 
793  MatrixBase& negateInPlace() {(*this) *= EPrecision(-1); return *this;}
794 
799  MatrixBase& resize(int m, int n) { helper.resize(m,n); return *this; }
805  MatrixBase& resizeKeep(int m, int n) { helper.resizeKeep(m,n); return *this; }
806 
807  // This prevents shape changes in a Matrix that would otherwise allow it. No harm if is
808  // are called on a Matrix that is locked already; it always succeeds.
809  void lockShape() {helper.lockShape();}
810 
811  // This allows shape changes again for a Matrix which was constructed to allow them
812  // but had them locked with the above routine. No harm if this is called on a Matrix
813  // that is already unlocked, but it is not allowed to call this on a Matrix which
814  // *never* allowed resizing. An exception will be thrown in that case.
815  void unlockShape() {helper.unlockShape();}
816 
817  // An assortment of handy conversions
818  const MatrixView_<ELT>& getAsMatrixView() const { return *reinterpret_cast<const MatrixView_<ELT>*>(this); }
819  MatrixView_<ELT>& updAsMatrixView() { return *reinterpret_cast< MatrixView_<ELT>*>(this); }
820  const Matrix_<ELT>& getAsMatrix() const { return *reinterpret_cast<const Matrix_<ELT>*>(this); }
821  Matrix_<ELT>& updAsMatrix() { return *reinterpret_cast< Matrix_<ELT>*>(this); }
822 
823  const VectorView_<ELT>& getAsVectorView() const
824  { assert(ncol()==1); return *reinterpret_cast<const VectorView_<ELT>*>(this); }
825  VectorView_<ELT>& updAsVectorView()
826  { assert(ncol()==1); return *reinterpret_cast< VectorView_<ELT>*>(this); }
827  const Vector_<ELT>& getAsVector() const
828  { assert(ncol()==1); return *reinterpret_cast<const Vector_<ELT>*>(this); }
829  Vector_<ELT>& updAsVector()
830  { assert(ncol()==1); return *reinterpret_cast< Vector_<ELT>*>(this); }
831  const VectorBase<ELT>& getAsVectorBase() const
832  { assert(ncol()==1); return *reinterpret_cast<const VectorBase<ELT>*>(this); }
833  VectorBase<ELT>& updAsVectorBase()
834  { assert(ncol()==1); return *reinterpret_cast< VectorBase<ELT>*>(this); }
835 
836  const RowVectorView_<ELT>& getAsRowVectorView() const
837  { assert(nrow()==1); return *reinterpret_cast<const RowVectorView_<ELT>*>(this); }
838  RowVectorView_<ELT>& updAsRowVectorView()
839  { assert(nrow()==1); return *reinterpret_cast< RowVectorView_<ELT>*>(this); }
840  const RowVector_<ELT>& getAsRowVector() const
841  { assert(nrow()==1); return *reinterpret_cast<const RowVector_<ELT>*>(this); }
842  RowVector_<ELT>& updAsRowVector()
843  { assert(nrow()==1); return *reinterpret_cast< RowVector_<ELT>*>(this); }
844  const RowVectorBase<ELT>& getAsRowVectorBase() const
845  { assert(nrow()==1); return *reinterpret_cast<const RowVectorBase<ELT>*>(this); }
846  RowVectorBase<ELT>& updAsRowVectorBase()
847  { assert(nrow()==1); return *reinterpret_cast< RowVectorBase<ELT>*>(this); }
848 
849  // Access to raw data. We have to return the raw data
850  // pointer as pointer-to-scalar because we may pack the elements tighter
851  // than a C++ array would.
852 
856  int getNScalarsPerElement() const {return NScalarsPerElement;}
857 
861  int getPackedSizeofElement() const {return NScalarsPerElement*sizeof(Scalar);}
862 
863  bool hasContiguousData() const {return helper.hasContiguousData();}
864  ptrdiff_t getContiguousScalarDataLength() const {
865  return helper.getContiguousDataLength();
866  }
867  const Scalar* getContiguousScalarData() const {
868  return helper.getContiguousData();
869  }
870  Scalar* updContiguousScalarData() {
871  return helper.updContiguousData();
872  }
873  void replaceContiguousScalarData(Scalar* newData, ptrdiff_t length, bool takeOwnership) {
874  helper.replaceContiguousData(newData,length,takeOwnership);
875  }
876  void replaceContiguousScalarData(const Scalar* newData, ptrdiff_t length) {
877  helper.replaceContiguousData(newData,length);
878  }
879  void swapOwnedContiguousScalarData(Scalar* newData, ptrdiff_t length, Scalar*& oldData) {
880  helper.swapOwnedContiguousData(newData,length,oldData);
881  }
882 
887  explicit MatrixBase(MatrixHelperRep<Scalar>* hrep) : helper(hrep) {}
888 
889 protected:
890  const MatrixHelper<Scalar>& getHelper() const {return helper;}
891  MatrixHelper<Scalar>& updHelper() {return helper;}
892 
893 private:
894  MatrixHelper<Scalar> helper; // this is just one pointer
895 
896  template <class EE> friend class MatrixBase;
897 };
898 
899 } //namespace SimTK
900 
901 #endif // SimTK_SIMMATRIX_MATRIXBASE_H_