Simbody
3.5
|
00001 #ifndef SimTK_SIMMATRIX_MATRIXBASE_H_ 00002 #define SimTK_SIMMATRIX_MATRIXBASE_H_ 00003 00004 /* -------------------------------------------------------------------------- * 00005 * Simbody(tm): SimTKcommon * 00006 * -------------------------------------------------------------------------- * 00007 * This is part of the SimTK biosimulation toolkit originating from * 00008 * Simbios, the NIH National Center for Physics-Based Simulation of * 00009 * Biological Structures at Stanford, funded under the NIH Roadmap for * 00010 * Medical Research, grant U54 GM072970. See https://simtk.org/home/simbody. * 00011 * * 00012 * Portions copyright (c) 2005-13 Stanford University and the Authors. * 00013 * Authors: Michael Sherman * 00014 * Contributors: * 00015 * * 00016 * Licensed under the Apache License, Version 2.0 (the "License"); you may * 00017 * not use this file except in compliance with the License. You may obtain a * 00018 * copy of the License at http://www.apache.org/licenses/LICENSE-2.0. * 00019 * * 00020 * Unless required by applicable law or agreed to in writing, software * 00021 * distributed under the License is distributed on an "AS IS" BASIS, * 00022 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. * 00023 * See the License for the specific language governing permissions and * 00024 * limitations under the License. * 00025 * -------------------------------------------------------------------------- */ 00026 00031 namespace SimTK { 00032 00033 //============================================================================== 00034 // MATRIX BASE 00035 //============================================================================== 00067 // ---------------------------------------------------------------------------- 00068 template <class ELT> class MatrixBase { 00069 public: 00070 // These typedefs are handy, but despite appearances you cannot 00071 // treat a MatrixBase as a composite numerical type. That is, 00072 // CNT<MatrixBase> will not compile, or if it does it won't be 00073 // meaningful. 00074 00075 typedef ELT E; 00076 typedef typename CNT<E>::TNeg ENeg; 00077 typedef typename CNT<E>::TWithoutNegator EWithoutNegator; 00078 typedef typename CNT<E>::TReal EReal; 00079 typedef typename CNT<E>::TImag EImag; 00080 typedef typename CNT<E>::TComplex EComplex; 00081 typedef typename CNT<E>::THerm EHerm; 00082 typedef typename CNT<E>::TPosTrans EPosTrans; 00083 00084 typedef typename CNT<E>::TAbs EAbs; 00085 typedef typename CNT<E>::TStandard EStandard; 00086 typedef typename CNT<E>::TInvert EInvert; 00087 typedef typename CNT<E>::TNormalize ENormalize; 00088 typedef typename CNT<E>::TSqHermT ESqHermT; 00089 typedef typename CNT<E>::TSqTHerm ESqTHerm; 00090 00091 typedef typename CNT<E>::Scalar EScalar; 00092 typedef typename CNT<E>::Number ENumber; 00093 typedef typename CNT<E>::StdNumber EStdNumber; 00094 typedef typename CNT<E>::Precision EPrecision; 00095 typedef typename CNT<E>::ScalarNormSq EScalarNormSq; 00096 00097 typedef EScalar Scalar; // the underlying Scalar type 00098 typedef ENumber Number; // negator removed from Scalar 00099 typedef EStdNumber StdNumber; // conjugate goes to complex 00100 typedef EPrecision Precision; // complex removed from StdNumber 00101 typedef EScalarNormSq ScalarNormSq; // type of scalar^2 00102 00103 typedef MatrixBase<E> T; 00104 typedef MatrixBase<ENeg> TNeg; 00105 typedef MatrixBase<EWithoutNegator> TWithoutNegator; 00106 typedef MatrixBase<EReal> TReal; 00107 typedef MatrixBase<EImag> TImag; 00108 typedef MatrixBase<EComplex> TComplex; 00109 typedef MatrixBase<EHerm> THerm; 00110 typedef MatrixBase<E> TPosTrans; 00111 00112 typedef MatrixBase<EAbs> TAbs; 00113 typedef MatrixBase<EStandard> TStandard; 00114 typedef MatrixBase<EInvert> TInvert; 00115 typedef MatrixBase<ENormalize> TNormalize; 00116 typedef MatrixBase<ESqHermT> TSqHermT; // ~Mat*Mat 00117 typedef MatrixBase<ESqTHerm> TSqTHerm; // Mat*~Mat 00118 00119 const MatrixCommitment& getCharacterCommitment() const {return helper.getCharacterCommitment();} 00120 const MatrixCharacter& getMatrixCharacter() const {return helper.getMatrixCharacter();} 00121 00124 void commitTo(const MatrixCommitment& mc) 00125 { helper.commitTo(mc); } 00126 00127 // This gives the resulting matrix type when (m(i,j) op P) is applied to each element. 00128 // It will have element types which are the regular composite result of E op P. 00129 template <class P> struct EltResult { 00130 typedef MatrixBase<typename CNT<E>::template Result<P>::Mul> Mul; 00131 typedef MatrixBase<typename CNT<E>::template Result<P>::Dvd> Dvd; 00132 typedef MatrixBase<typename CNT<E>::template Result<P>::Add> Add; 00133 typedef MatrixBase<typename CNT<E>::template Result<P>::Sub> Sub; 00134 }; 00135 00137 int nrow() const {return helper.nrow();} 00139 int ncol() const {return helper.ncol();} 00140 00148 ptrdiff_t nelt() const {return helper.nelt();} 00149 00151 bool isResizeable() const {return getCharacterCommitment().isResizeable();} 00152 00153 enum { 00154 NScalarsPerElement = CNT<E>::NActualScalars, 00155 CppNScalarsPerElement = sizeof(E) / sizeof(Scalar) 00156 }; 00157 00161 MatrixBase() : helper(NScalarsPerElement,CppNScalarsPerElement) {} 00162 00165 MatrixBase(int m, int n) 00166 : helper(NScalarsPerElement,CppNScalarsPerElement,MatrixCommitment(),m,n) {} 00167 00173 explicit MatrixBase(const MatrixCommitment& commitment) 00174 : helper(NScalarsPerElement,CppNScalarsPerElement,commitment) {} 00175 00176 00181 MatrixBase(const MatrixCommitment& commitment, int m, int n) 00182 : helper(NScalarsPerElement,CppNScalarsPerElement,commitment,m,n) {} 00183 00185 MatrixBase(const MatrixBase& b) 00186 : helper(b.helper.getCharacterCommitment(), 00187 b.helper, typename MatrixHelper<Scalar>::DeepCopy()) { } 00188 00191 MatrixBase(const TNeg& b) 00192 : helper(b.helper.getCharacterCommitment(), 00193 b.helper, typename MatrixHelper<Scalar>::DeepCopy()) { } 00194 00197 MatrixBase& copyAssign(const MatrixBase& b) { 00198 helper.copyAssign(b.helper); 00199 return *this; 00200 } 00201 MatrixBase& operator=(const MatrixBase& b) { return copyAssign(b); } 00202 00203 00212 MatrixBase& viewAssign(const MatrixBase& src) { 00213 helper.writableViewAssign(const_cast<MatrixHelper<Scalar>&>(src.helper)); 00214 return *this; 00215 } 00216 00217 // default destructor 00218 00225 MatrixBase(const MatrixCommitment& commitment, int m, int n, const ELT& initialValue) 00226 : helper(NScalarsPerElement, CppNScalarsPerElement, commitment, m, n) 00227 { helper.fillWith(reinterpret_cast<const Scalar*>(&initialValue)); } 00228 00239 MatrixBase(const MatrixCommitment& commitment, int m, int n, 00240 const ELT* cppInitialValuesByRow) 00241 : helper(NScalarsPerElement, CppNScalarsPerElement, commitment, m, n) 00242 { helper.copyInByRowsFromCpp(reinterpret_cast<const Scalar*>(cppInitialValuesByRow)); } 00243 00255 00257 MatrixBase(const MatrixCommitment& commitment, 00258 const MatrixCharacter& character, 00259 int spacing, const Scalar* data) // read only data 00260 : helper(NScalarsPerElement, CppNScalarsPerElement, 00261 commitment, character, spacing, data) {} 00262 00264 MatrixBase(const MatrixCommitment& commitment, 00265 const MatrixCharacter& character, 00266 int spacing, Scalar* data) // writable data 00267 : helper(NScalarsPerElement, CppNScalarsPerElement, 00268 commitment, character, spacing, data) {} 00270 00271 // Create a new MatrixBase from an existing helper. Both shallow and deep copies are possible. 00272 MatrixBase(const MatrixCommitment& commitment, 00273 MatrixHelper<Scalar>& source, 00274 const typename MatrixHelper<Scalar>::ShallowCopy& shallow) 00275 : helper(commitment, source, shallow) {} 00276 MatrixBase(const MatrixCommitment& commitment, 00277 const MatrixHelper<Scalar>& source, 00278 const typename MatrixHelper<Scalar>::ShallowCopy& shallow) 00279 : helper(commitment, source, shallow) {} 00280 MatrixBase(const MatrixCommitment& commitment, 00281 const MatrixHelper<Scalar>& source, 00282 const typename MatrixHelper<Scalar>::DeepCopy& deep) 00283 : helper(commitment, source, deep) {} 00284 00288 void clear() {helper.clear();} 00289 00290 MatrixBase& operator*=(const StdNumber& t) { helper.scaleBy(t); return *this; } 00291 MatrixBase& operator/=(const StdNumber& t) { helper.scaleBy(StdNumber(1)/t); return *this; } 00292 MatrixBase& operator+=(const MatrixBase& r) { helper.addIn(r.helper); return *this; } 00293 MatrixBase& operator-=(const MatrixBase& r) { helper.subIn(r.helper); return *this; } 00294 00295 template <class EE> MatrixBase(const MatrixBase<EE>& b) 00296 : helper(MatrixCommitment(),b.helper, typename MatrixHelper<Scalar>::DeepCopy()) { } 00297 00298 template <class EE> MatrixBase& operator=(const MatrixBase<EE>& b) 00299 { helper = b.helper; return *this; } 00300 template <class EE> MatrixBase& operator+=(const MatrixBase<EE>& b) 00301 { helper.addIn(b.helper); return *this; } 00302 template <class EE> MatrixBase& operator-=(const MatrixBase<EE>& b) 00303 { helper.subIn(b.helper); return *this; } 00304 00315 MatrixBase& operator=(const ELT& t) { 00316 setToZero(); updDiag().setTo(t); 00317 return *this; 00318 } 00319 00325 template <class S> inline MatrixBase& 00326 scalarAssign(const S& s) { 00327 setToZero(); updDiag().setTo(s); 00328 return *this; 00329 } 00330 00334 template <class S> inline MatrixBase& 00335 scalarAddInPlace(const S& s) { 00336 updDiag().elementwiseAddScalarInPlace(s); 00337 } 00338 00339 00343 template <class S> inline MatrixBase& 00344 scalarSubtractInPlace(const S& s) { 00345 updDiag().elementwiseSubtractScalarInPlace(s); 00346 } 00347 00350 template <class S> inline MatrixBase& 00351 scalarSubtractFromLeftInPlace(const S& s) { 00352 negateInPlace(); 00353 updDiag().elementwiseAddScalarInPlace(s); // yes, add 00354 } 00355 00362 template <class S> inline MatrixBase& 00363 scalarMultiplyInPlace(const S&); 00364 00368 template <class S> inline MatrixBase& 00369 scalarMultiplyFromLeftInPlace(const S&); 00370 00377 template <class S> inline MatrixBase& 00378 scalarDivideInPlace(const S&); 00379 00383 template <class S> inline MatrixBase& 00384 scalarDivideFromLeftInPlace(const S&); 00385 00386 00389 template <class EE> inline MatrixBase& 00390 rowScaleInPlace(const VectorBase<EE>&); 00391 00394 template <class EE> inline void 00395 rowScale(const VectorBase<EE>& r, typename EltResult<EE>::Mul& out) const; 00396 00397 template <class EE> inline typename EltResult<EE>::Mul 00398 rowScale(const VectorBase<EE>& r) const { 00399 typename EltResult<EE>::Mul out(nrow(), ncol()); rowScale(r,out); return out; 00400 } 00401 00404 template <class EE> inline MatrixBase& 00405 colScaleInPlace(const VectorBase<EE>&); 00406 00407 template <class EE> inline void 00408 colScale(const VectorBase<EE>& c, typename EltResult<EE>::Mul& out) const; 00409 00410 template <class EE> inline typename EltResult<EE>::Mul 00411 colScale(const VectorBase<EE>& c) const { 00412 typename EltResult<EE>::Mul out(nrow(), ncol()); colScale(c,out); return out; 00413 } 00414 00415 00420 template <class ER, class EC> inline MatrixBase& 00421 rowAndColScaleInPlace(const VectorBase<ER>& r, const VectorBase<EC>& c); 00422 00423 template <class ER, class EC> inline void 00424 rowAndColScale(const VectorBase<ER>& r, const VectorBase<EC>& c, 00425 typename EltResult<typename VectorBase<ER>::template EltResult<EC>::Mul>::Mul& out) const; 00426 00427 template <class ER, class EC> inline typename EltResult<typename VectorBase<ER>::template EltResult<EC>::Mul>::Mul 00428 rowAndColScale(const VectorBase<ER>& r, const VectorBase<EC>& c) const { 00429 typename EltResult<typename VectorBase<ER>::template EltResult<EC>::Mul>::Mul 00430 out(nrow(), ncol()); 00431 rowAndColScale(r,c,out); return out; 00432 } 00433 00441 template <class S> inline MatrixBase& 00442 elementwiseAssign(const S& s); 00443 00445 MatrixBase& elementwiseAssign(int s) 00446 { return elementwiseAssign<Real>(Real(s)); } 00447 00449 MatrixBase& elementwiseInvertInPlace(); 00450 00451 void elementwiseInvert(MatrixBase<typename CNT<E>::TInvert>& out) const; 00452 00453 MatrixBase<typename CNT<E>::TInvert> elementwiseInvert() const { 00454 MatrixBase<typename CNT<E>::TInvert> out(nrow(), ncol()); 00455 elementwiseInvert(out); 00456 return out; 00457 } 00458 00466 template <class S> inline MatrixBase& 00467 elementwiseAddScalarInPlace(const S& s); 00468 00469 template <class S> inline void 00470 elementwiseAddScalar(const S& s, typename EltResult<S>::Add&) const; 00471 00472 template <class S> inline typename EltResult<S>::Add 00473 elementwiseAddScalar(const S& s) const { 00474 typename EltResult<S>::Add out(nrow(), ncol()); 00475 elementwiseAddScalar(s,out); 00476 return out; 00477 } 00478 00486 template <class S> inline MatrixBase& 00487 elementwiseSubtractScalarInPlace(const S& s); 00488 00489 template <class S> inline void 00490 elementwiseSubtractScalar(const S& s, typename EltResult<S>::Sub&) const; 00491 00492 template <class S> inline typename EltResult<S>::Sub 00493 elementwiseSubtractScalar(const S& s) const { 00494 typename EltResult<S>::Sub out(nrow(), ncol()); 00495 elementwiseSubtractScalar(s,out); 00496 return out; 00497 } 00498 00507 template <class S> inline MatrixBase& 00508 elementwiseSubtractFromScalarInPlace(const S& s); 00509 00510 template <class S> inline void 00511 elementwiseSubtractFromScalar( 00512 const S&, 00513 typename MatrixBase<S>::template EltResult<E>::Sub&) const; 00514 00515 template <class S> inline typename MatrixBase<S>::template EltResult<E>::Sub 00516 elementwiseSubtractFromScalar(const S& s) const { 00517 typename MatrixBase<S>::template EltResult<E>::Sub out(nrow(), ncol()); 00518 elementwiseSubtractFromScalar<S>(s,out); 00519 return out; 00520 } 00521 00523 template <class EE> inline MatrixBase& 00524 elementwiseMultiplyInPlace(const MatrixBase<EE>&); 00525 00526 template <class EE> inline void 00527 elementwiseMultiply(const MatrixBase<EE>&, typename EltResult<EE>::Mul&) const; 00528 00529 template <class EE> inline typename EltResult<EE>::Mul 00530 elementwiseMultiply(const MatrixBase<EE>& m) const { 00531 typename EltResult<EE>::Mul out(nrow(), ncol()); 00532 elementwiseMultiply<EE>(m,out); 00533 return out; 00534 } 00535 00537 template <class EE> inline MatrixBase& 00538 elementwiseMultiplyFromLeftInPlace(const MatrixBase<EE>&); 00539 00540 template <class EE> inline void 00541 elementwiseMultiplyFromLeft( 00542 const MatrixBase<EE>&, 00543 typename MatrixBase<EE>::template EltResult<E>::Mul&) const; 00544 00545 template <class EE> inline typename MatrixBase<EE>::template EltResult<E>::Mul 00546 elementwiseMultiplyFromLeft(const MatrixBase<EE>& m) const { 00547 typename EltResult<EE>::Mul out(nrow(), ncol()); 00548 elementwiseMultiplyFromLeft<EE>(m,out); 00549 return out; 00550 } 00551 00553 template <class EE> inline MatrixBase& 00554 elementwiseDivideInPlace(const MatrixBase<EE>&); 00555 00556 template <class EE> inline void 00557 elementwiseDivide(const MatrixBase<EE>&, typename EltResult<EE>::Dvd&) const; 00558 00559 template <class EE> inline typename EltResult<EE>::Dvd 00560 elementwiseDivide(const MatrixBase<EE>& m) const { 00561 typename EltResult<EE>::Dvd out(nrow(), ncol()); 00562 elementwiseDivide<EE>(m,out); 00563 return out; 00564 } 00565 00567 template <class EE> inline MatrixBase& 00568 elementwiseDivideFromLeftInPlace(const MatrixBase<EE>&); 00569 00570 template <class EE> inline void 00571 elementwiseDivideFromLeft( 00572 const MatrixBase<EE>&, 00573 typename MatrixBase<EE>::template EltResult<E>::Dvd&) const; 00574 00575 template <class EE> inline typename MatrixBase<EE>::template EltResult<EE>::Dvd 00576 elementwiseDivideFromLeft(const MatrixBase<EE>& m) const { 00577 typename MatrixBase<EE>::template EltResult<E>::Dvd out(nrow(), ncol()); 00578 elementwiseDivideFromLeft<EE>(m,out); 00579 return out; 00580 } 00581 00583 MatrixBase& setTo(const ELT& t) {helper.fillWith(reinterpret_cast<const Scalar*>(&t)); return *this;} 00584 MatrixBase& setToNaN() {helper.fillWithScalar(CNT<StdNumber>::getNaN()); return *this;} 00585 MatrixBase& setToZero() {helper.fillWithScalar(StdNumber(0)); return *this;} 00586 00587 // View creating operators. 00588 inline RowVectorView_<ELT> row(int i) const; // select a row 00589 inline RowVectorView_<ELT> updRow(int i); 00590 inline VectorView_<ELT> col(int j) const; // select a column 00591 inline VectorView_<ELT> updCol(int j); 00592 00593 RowVectorView_<ELT> operator[](int i) const {return row(i);} 00594 RowVectorView_<ELT> operator[](int i) {return updRow(i);} 00595 VectorView_<ELT> operator()(int j) const {return col(j);} 00596 VectorView_<ELT> operator()(int j) {return updCol(j);} 00597 00598 // Select a block. 00599 inline MatrixView_<ELT> block(int i, int j, int m, int n) const; 00600 inline MatrixView_<ELT> updBlock(int i, int j, int m, int n); 00601 00602 MatrixView_<ELT> operator()(int i, int j, int m, int n) const 00603 { return block(i,j,m,n); } 00604 MatrixView_<ELT> operator()(int i, int j, int m, int n) 00605 { return updBlock(i,j,m,n); } 00606 00607 // Hermitian transpose. 00608 inline MatrixView_<EHerm> transpose() const; 00609 inline MatrixView_<EHerm> updTranspose(); 00610 00611 MatrixView_<EHerm> operator~() const {return transpose();} 00612 MatrixView_<EHerm> operator~() {return updTranspose();} 00613 00616 inline VectorView_<ELT> diag() const; 00619 inline VectorView_<ELT> updDiag(); 00622 VectorView_<ELT> diag() {return updDiag();} 00623 00624 // Create a view of the real or imaginary elements. TODO 00625 //inline MatrixView_<EReal> real() const; 00626 //inline MatrixView_<EReal> updReal(); 00627 //inline MatrixView_<EImag> imag() const; 00628 //inline MatrixView_<EImag> updImag(); 00629 00630 // Overload "real" and "imag" for both read and write as a nod to convention. TODO 00631 //MatrixView_<EReal> real() {return updReal();} 00632 //MatrixView_<EReal> imag() {return updImag();} 00633 00634 // TODO: this routine seems ill-advised but I need it for the IVM port at the moment 00635 TInvert invert() const { // return a newly-allocated inverse; dump negator 00636 TInvert m(*this); 00637 m.helper.invertInPlace(); 00638 return m; // TODO - bad: makes an extra copy 00639 } 00640 00641 void invertInPlace() {helper.invertInPlace();} 00642 00644 void dump(const char* msg=0) const { 00645 helper.dump(msg); 00646 } 00647 00656 const ELT& getElt(int i, int j) const { return *reinterpret_cast<const ELT*>(helper.getElt(i,j)); } 00657 ELT& updElt(int i, int j) { return *reinterpret_cast< ELT*>(helper.updElt(i,j)); } 00658 00659 const ELT& operator()(int i, int j) const {return getElt(i,j);} 00660 ELT& operator()(int i, int j) {return updElt(i,j);} 00661 00666 void getAnyElt(int i, int j, ELT& value) const 00667 { helper.getAnyElt(i,j,reinterpret_cast<Scalar*>(&value)); } 00668 ELT getAnyElt(int i, int j) const {ELT e; getAnyElt(i,j,e); return e;} 00669 00672 // TODO: very slow! Should be optimized at least for the case 00673 // where ELT is a Scalar. 00674 ScalarNormSq scalarNormSqr() const { 00675 const int nr=nrow(), nc=ncol(); 00676 ScalarNormSq sum(0); 00677 for(int j=0;j<nc;++j) 00678 for (int i=0; i<nr; ++i) 00679 sum += CNT<E>::scalarNormSqr((*this)(i,j)); 00680 return sum; 00681 } 00682 00686 // TODO: very slow! Should be optimized at least for the case 00687 // where ELT is a Scalar. 00688 void abs(TAbs& mabs) const { 00689 const int nr=nrow(), nc=ncol(); 00690 mabs.resize(nr,nc); 00691 for(int j=0;j<nc;++j) 00692 for (int i=0; i<nr; ++i) 00693 mabs(i,j) = CNT<E>::abs((*this)(i,j)); 00694 } 00695 00699 TAbs abs() const { TAbs mabs; abs(mabs); return mabs; } 00700 00711 TStandard standardize() const { 00712 const int nr=nrow(), nc=ncol(); 00713 TStandard mstd(nr, nc); 00714 for(int j=0;j<nc;++j) 00715 for (int i=0; i<nr; ++i) 00716 mstd(i,j) = CNT<E>::standardize((*this)(i,j)); 00717 return mstd; 00718 } 00719 00723 ScalarNormSq normSqr() const { return scalarNormSqr(); } 00724 // TODO -- not good; unnecessary overflow 00725 typename CNT<ScalarNormSq>::TSqrt 00726 norm() const { return CNT<ScalarNormSq>::sqrt(scalarNormSqr()); } 00727 00731 typename CNT<ScalarNormSq>::TSqrt 00732 normRMS() const { 00733 if (!CNT<ELT>::IsScalar) 00734 SimTK_THROW1(Exception::Cant, "normRMS() only defined for scalar elements"); 00735 if (nelt() == 0) 00736 return typename CNT<ScalarNormSq>::TSqrt(0); 00737 return CNT<ScalarNormSq>::sqrt(scalarNormSqr()/nelt()); 00738 } 00739 00741 RowVector_<ELT> colSum() const { 00742 const int cols = ncol(); 00743 RowVector_<ELT> row(cols); 00744 for (int j = 0; j < cols; ++j) 00745 helper.colSum(j, reinterpret_cast<Scalar*>(&row[j])); 00746 return row; 00747 } 00749 RowVector_<ELT> sum() const {return colSum();} 00750 00752 Vector_<ELT> rowSum() const { 00753 const int rows = nrow(); 00754 Vector_<ELT> col(rows); 00755 for (int i = 0; i < rows; ++i) 00756 helper.rowSum(i, reinterpret_cast<Scalar*>(&col[i])); 00757 return col; 00758 } 00759 00760 //TODO: make unary +/- return a self-view so they won't reallocate? 00761 const MatrixBase& operator+() const {return *this; } 00762 const TNeg& negate() const {return *reinterpret_cast<const TNeg*>(this); } 00763 TNeg& updNegate() {return *reinterpret_cast<TNeg*>(this); } 00764 00765 const TNeg& operator-() const {return negate();} 00766 TNeg& operator-() {return updNegate();} 00767 00768 MatrixBase& negateInPlace() {(*this) *= EPrecision(-1); return *this;} 00769 00774 MatrixBase& resize(int m, int n) { helper.resize(m,n); return *this; } 00780 MatrixBase& resizeKeep(int m, int n) { helper.resizeKeep(m,n); return *this; } 00781 00782 // This prevents shape changes in a Matrix that would otherwise allow it. No harm if is 00783 // are called on a Matrix that is locked already; it always succeeds. 00784 void lockShape() {helper.lockShape();} 00785 00786 // This allows shape changes again for a Matrix which was constructed to allow them 00787 // but had them locked with the above routine. No harm if this is called on a Matrix 00788 // that is already unlocked, but it is not allowed to call this on a Matrix which 00789 // *never* allowed resizing. An exception will be thrown in that case. 00790 void unlockShape() {helper.unlockShape();} 00791 00792 // An assortment of handy conversions 00793 const MatrixView_<ELT>& getAsMatrixView() const { return *reinterpret_cast<const MatrixView_<ELT>*>(this); } 00794 MatrixView_<ELT>& updAsMatrixView() { return *reinterpret_cast< MatrixView_<ELT>*>(this); } 00795 const Matrix_<ELT>& getAsMatrix() const { return *reinterpret_cast<const Matrix_<ELT>*>(this); } 00796 Matrix_<ELT>& updAsMatrix() { return *reinterpret_cast< Matrix_<ELT>*>(this); } 00797 00798 const VectorView_<ELT>& getAsVectorView() const 00799 { assert(ncol()==1); return *reinterpret_cast<const VectorView_<ELT>*>(this); } 00800 VectorView_<ELT>& updAsVectorView() 00801 { assert(ncol()==1); return *reinterpret_cast< VectorView_<ELT>*>(this); } 00802 const Vector_<ELT>& getAsVector() const 00803 { assert(ncol()==1); return *reinterpret_cast<const Vector_<ELT>*>(this); } 00804 Vector_<ELT>& updAsVector() 00805 { assert(ncol()==1); return *reinterpret_cast< Vector_<ELT>*>(this); } 00806 const VectorBase<ELT>& getAsVectorBase() const 00807 { assert(ncol()==1); return *reinterpret_cast<const VectorBase<ELT>*>(this); } 00808 VectorBase<ELT>& updAsVectorBase() 00809 { assert(ncol()==1); return *reinterpret_cast< VectorBase<ELT>*>(this); } 00810 00811 const RowVectorView_<ELT>& getAsRowVectorView() const 00812 { assert(nrow()==1); return *reinterpret_cast<const RowVectorView_<ELT>*>(this); } 00813 RowVectorView_<ELT>& updAsRowVectorView() 00814 { assert(nrow()==1); return *reinterpret_cast< RowVectorView_<ELT>*>(this); } 00815 const RowVector_<ELT>& getAsRowVector() const 00816 { assert(nrow()==1); return *reinterpret_cast<const RowVector_<ELT>*>(this); } 00817 RowVector_<ELT>& updAsRowVector() 00818 { assert(nrow()==1); return *reinterpret_cast< RowVector_<ELT>*>(this); } 00819 const RowVectorBase<ELT>& getAsRowVectorBase() const 00820 { assert(nrow()==1); return *reinterpret_cast<const RowVectorBase<ELT>*>(this); } 00821 RowVectorBase<ELT>& updAsRowVectorBase() 00822 { assert(nrow()==1); return *reinterpret_cast< RowVectorBase<ELT>*>(this); } 00823 00824 // Access to raw data. We have to return the raw data 00825 // pointer as pointer-to-scalar because we may pack the elements tighter 00826 // than a C++ array would. 00827 00831 int getNScalarsPerElement() const {return NScalarsPerElement;} 00832 00836 int getPackedSizeofElement() const {return NScalarsPerElement*sizeof(Scalar);} 00837 00838 bool hasContiguousData() const {return helper.hasContiguousData();} 00839 ptrdiff_t getContiguousScalarDataLength() const { 00840 return helper.getContiguousDataLength(); 00841 } 00842 const Scalar* getContiguousScalarData() const { 00843 return helper.getContiguousData(); 00844 } 00845 Scalar* updContiguousScalarData() { 00846 return helper.updContiguousData(); 00847 } 00848 void replaceContiguousScalarData(Scalar* newData, ptrdiff_t length, bool takeOwnership) { 00849 helper.replaceContiguousData(newData,length,takeOwnership); 00850 } 00851 void replaceContiguousScalarData(const Scalar* newData, ptrdiff_t length) { 00852 helper.replaceContiguousData(newData,length); 00853 } 00854 void swapOwnedContiguousScalarData(Scalar* newData, ptrdiff_t length, Scalar*& oldData) { 00855 helper.swapOwnedContiguousData(newData,length,oldData); 00856 } 00857 00862 explicit MatrixBase(MatrixHelperRep<Scalar>* hrep) : helper(hrep) {} 00863 00864 protected: 00865 const MatrixHelper<Scalar>& getHelper() const {return helper;} 00866 MatrixHelper<Scalar>& updHelper() {return helper;} 00867 00868 private: 00869 MatrixHelper<Scalar> helper; // this is just one pointer 00870 00871 template <class EE> friend class MatrixBase; 00872 00873 // ============================= Unimplemented ============================= 00874 // This routine is useful for implementing friendlier Matrix expressions and operators. 00875 // It maps closely to the Level-3 BLAS family of pxxmm() routines like sgemm(). The 00876 // operation performed assumes that "this" is the result, and that "this" has 00877 // already been sized correctly to receive the result. We'll compute 00878 // this = beta*this + alpha*A*B 00879 // If beta is 0 then "this" can be uninitialized. If alpha is 0 we promise not 00880 // to look at A or B. The routine can work efficiently even if A and/or B are transposed 00881 // by their views, so an expression like 00882 // C += s * ~A * ~B 00883 // can be performed with the single equivalent call 00884 // C.matmul(1., s, Tr(A), Tr(B)) 00885 // where Tr(A) indicates a transposed view of the original A. 00886 // The ultimate efficiency of this call depends on the data layout and views which 00887 // are used for the three matrices. 00888 // NOTE: neither A nor B can be the same matrix as 'this', nor views of the same data 00889 // which would expose elements of 'this' that will be modified by this operation. 00890 template <class ELT_A, class ELT_B> 00891 MatrixBase& matmul(const StdNumber& beta, // applied to 'this' 00892 const StdNumber& alpha, const MatrixBase<ELT_A>& A, const MatrixBase<ELT_B>& B) 00893 { 00894 helper.matmul(beta,alpha,A.helper,B.helper); 00895 return *this; 00896 } 00897 00898 }; 00899 00900 } //namespace SimTK 00901 00902 #endif // SimTK_SIMMATRIX_MATRIXBASE_H_