SymMat.h

Go to the documentation of this file.
00001 #ifndef SimTK_SIMMATRIX_SMALLMATRIX_SYMMAT_H_
00002 #define SimTK_SIMMATRIX_SMALLMATRIX_SYMMAT_H_
00003 
00004 /* -------------------------------------------------------------------------- *
00005  *                      SimTK Core: SimTK Simmatrix(tm)                       *
00006  * -------------------------------------------------------------------------- *
00007  * This is part of the SimTK Core 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.               *
00011  *                                                                            *
00012  * Portions copyright (c) 2005-7 Stanford University and the Authors.         *
00013  * Authors: Michael Sherman                                                   *
00014  * Contributors:                                                              *
00015  *                                                                            *
00016  * Permission is hereby granted, free of charge, to any person obtaining a    *
00017  * copy of this software and associated documentation files (the "Software"), *
00018  * to deal in the Software without restriction, including without limitation  *
00019  * the rights to use, copy, modify, merge, publish, distribute, sublicense,   *
00020  * and/or sell copies of the Software, and to permit persons to whom the      *
00021  * Software is furnished to do so, subject to the following conditions:       *
00022  *                                                                            *
00023  * The above copyright notice and this permission notice shall be included in *
00024  * all copies or substantial portions of the Software.                        *
00025  *                                                                            *
00026  * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR *
00027  * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,   *
00028  * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL    *
00029  * THE AUTHORS, CONTRIBUTORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM,    *
00030  * DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR      *
00031  * OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE  *
00032  * USE OR OTHER DEALINGS IN THE SOFTWARE.                                     *
00033  * -------------------------------------------------------------------------- */
00034 
00074 #include "SimTKcommon/internal/common.h"
00075 
00076 namespace SimTK {
00077 
00079 template <int M, class ELT, int RS> class SymMat {
00080     typedef ELT                                 E;
00081     typedef typename CNT<E>::TNeg               ENeg;
00082     typedef typename CNT<E>::TWithoutNegator    EWithoutNegator;
00083     typedef typename CNT<E>::TReal              EReal;
00084     typedef typename CNT<E>::TImag              EImag;
00085     typedef typename CNT<E>::TComplex           EComplex;
00086     typedef typename CNT<E>::THerm              EHerm;
00087     typedef typename CNT<E>::TPosTrans          EPosTrans;
00088     typedef typename CNT<E>::TSqHermT           ESqHermT;
00089     typedef typename CNT<E>::TSqTHerm           ESqTHerm;
00090 
00091     typedef typename CNT<E>::TSqrt              ESqrt;
00092     typedef typename CNT<E>::TAbs               EAbs;
00093     typedef typename CNT<E>::TStandard          EStandard;
00094     typedef typename CNT<E>::TInvert            EInvert;
00095     typedef typename CNT<E>::TNormalize         ENormalize;
00096 
00097     typedef typename CNT<E>::Scalar             EScalar;
00098     typedef typename CNT<E>::ULessScalar        EULessScalar;
00099     typedef typename CNT<E>::Number             ENumber;
00100     typedef typename CNT<E>::StdNumber          EStdNumber;
00101     typedef typename CNT<E>::Precision          EPrecision;
00102     typedef typename CNT<E>::ScalarNormSq       EScalarNormSq;
00103 
00104 public:
00105 
00106     enum {
00107         NRows               = M,
00108         NCols               = M,
00109         NPackedElements     = (M*(M+1))/2,
00110         NActualElements     = RS * NPackedElements,
00111         NActualScalars      = CNT<E>::NActualScalars * NActualElements,
00112         RowSpacing          = RS,
00113         ColSpacing          = NActualElements,
00114         ImagOffset          = NTraits<ENumber>::ImagOffset,
00115         RealStrideFactor    = 1, // composite types don't change size when
00116                                  // cast from complex to real or imaginary
00117         ArgDepth            = ((int)CNT<E>::ArgDepth < (int)MAX_RESOLVED_DEPTH 
00118                                 ? CNT<E>::ArgDepth + 1 
00119                                 : MAX_RESOLVED_DEPTH),
00120         IsScalar            = 0,
00121         IsULessScalar       = 0,
00122         IsNumber            = 0,
00123         IsStdNumber         = 0,
00124         IsPrecision         = 0,
00125         SignInterpretation  = CNT<E>::SignInterpretation
00126     };
00127 
00128     typedef SymMat<M,E,RS>                  T;
00129     typedef SymMat<M,ENeg,RS>               TNeg;
00130     typedef SymMat<M,EWithoutNegator,RS>    TWithoutNegator;
00131 
00132     typedef SymMat<M,EReal,RS*CNT<E>::RealStrideFactor>          
00133                                             TReal;
00134     typedef SymMat<M,EImag,RS*CNT<E>::RealStrideFactor>          
00135                                             TImag;
00136     typedef SymMat<M,EComplex,RS>           TComplex;
00137     typedef T                               THerm;   // These two are opposite of what you might expect
00138     typedef SymMat<M,EHerm,RS>              TPosTrans;
00139     typedef E                               TElement;
00140     typedef Vec<M,E,RS>                     TDiag;   // storage type for the diagonal elements
00141     typedef Vec<(M*(M-1))/2,E,RS>           TLower;  // storage type for the below-diag elements
00142     typedef Vec<(M*(M-1))/2,EHerm,RS>       TUpper;  // cast TLower to this for upper elements
00143     typedef Vec<(M*(M+1))/2,E,RS>           TAsVec;  // the whole SymMat as a single Vec
00144 
00145     // These are the results of calculations, so are returned in new, packed
00146     // memory. Be sure to refer to element types here which are also packed.
00147     typedef Row<M,E,1>                  TRow; // packed since we have to copy
00148     typedef Vec<M,E,1>                  TCol;
00149     typedef SymMat<M,ESqrt,1>           TSqrt;
00150     typedef SymMat<M,EAbs,1>            TAbs;
00151     typedef SymMat<M,EStandard,1>       TStandard;
00152     typedef SymMat<M,EInvert,1>         TInvert;
00153     typedef SymMat<M,ENormalize,1>      TNormalize;
00154     typedef SymMat<M,ESqHermT,1>        TSqHermT; // ~Mat*Mat
00155     typedef SymMat<M,ESqTHerm,1>        TSqTHerm; // Mat*~Mat
00156     typedef SymMat<M,E,1>               TPacked;  // no extra row spacing for new data
00157     
00158     typedef EScalar                     Scalar;
00159     typedef EULessScalar                ULessScalar;
00160     typedef ENumber                     Number;
00161     typedef EStdNumber                  StdNumber;
00162     typedef EPrecision                  Precision;
00163     typedef EScalarNormSq               ScalarNormSq;
00164 
00165     int size() const { return (M*(M+1))/2; }
00166     int nrow() const { return M; }
00167     int ncol() const { return M; }
00168 
00169     // Scalar norm square is sum( squares of all scalars ). The off-diagonals
00170     // come up twice.
00171     ScalarNormSq scalarNormSqr() const { 
00172         return getDiag().scalarNormSqr() + 2.*getLower().scalarNormSqr();
00173     }
00174 
00175     // sqrt() is elementwise square root; that is, the return value has the same
00176     // dimension as this SymMat but with each element replaced by whatever it thinks
00177     // its square root is.
00178     TSqrt sqrt() const { 
00179         return TSqrt(getAsVec().sqrt());
00180     }
00181 
00182     // abs() is elementwise absolute value; that is, the return value has the same
00183     // dimension as this SymMat but with each element replaced by whatever it thinks
00184     // its absolute value is.
00185     TAbs abs() const { 
00186         return TAbs(getAsVec().abs());
00187     }
00188 
00189     TStandard standardize() const {
00190         return TStandard(getAsVec().standardize());
00191     }
00192 
00193     EStandard trace() const {return getDiag().sum();}
00194 
00195     // This gives the resulting SymMat type when (m[i] op P) is applied to each element of m.
00196     // It is a SymMat of dimension M, spacing 1, and element types which are the regular
00197     // composite result of E op P. Typically P is a scalar type but it doesn't have to be.
00198     template <class P> struct EltResult { 
00199         typedef SymMat<M, typename CNT<E>::template Result<P>::Mul, 1> Mul;
00200         typedef SymMat<M, typename CNT<E>::template Result<P>::Dvd, 1> Dvd;
00201         typedef SymMat<M, typename CNT<E>::template Result<P>::Add, 1> Add;
00202         typedef SymMat<M, typename CNT<E>::template Result<P>::Sub, 1> Sub;
00203     };
00204 
00205     // This is the composite result for m op P where P is some kind of appropriately shaped
00206     // non-scalar type.
00207     template <class P> struct Result { 
00208         typedef MulCNTs<M,M,ArgDepth,SymMat,ColSpacing,RowSpacing,
00209             CNT<P>::NRows, CNT<P>::NCols, CNT<P>::ArgDepth,
00210             P, CNT<P>::ColSpacing, CNT<P>::RowSpacing> MulOp;
00211         typedef typename MulOp::Type Mul;
00212 
00213         typedef MulCNTsNonConforming<M,M,ArgDepth,SymMat,ColSpacing,RowSpacing,
00214             CNT<P>::NRows, CNT<P>::NCols, CNT<P>::ArgDepth,
00215             P, CNT<P>::ColSpacing, CNT<P>::RowSpacing> MulOpNonConforming;
00216         typedef typename MulOpNonConforming::Type MulNon;
00217 
00218 
00219         typedef DvdCNTs<M,M,ArgDepth,SymMat,ColSpacing,RowSpacing,
00220             CNT<P>::NRows, CNT<P>::NCols, CNT<P>::ArgDepth,
00221             P, CNT<P>::ColSpacing, CNT<P>::RowSpacing> DvdOp;
00222         typedef typename DvdOp::Type Dvd;
00223 
00224         typedef AddCNTs<M,M,ArgDepth,SymMat,ColSpacing,RowSpacing,
00225             CNT<P>::NRows, CNT<P>::NCols, CNT<P>::ArgDepth,
00226             P, CNT<P>::ColSpacing, CNT<P>::RowSpacing> AddOp;
00227         typedef typename AddOp::Type Add;
00228 
00229         typedef SubCNTs<M,M,ArgDepth,SymMat,ColSpacing,RowSpacing,
00230             CNT<P>::NRows, CNT<P>::NCols, CNT<P>::ArgDepth,
00231             P, CNT<P>::ColSpacing, CNT<P>::RowSpacing> SubOp;
00232         typedef typename SubOp::Type Sub;
00233     };
00234 
00235     // Shape-preserving element substitution (always packed)
00236     template <class P> struct Substitute {
00237         typedef SymMat<M,P> Type;
00238     };
00239 
00240     // Default construction initializes to NaN when debugging but
00241     // is left uninitialized otherwise.
00242     SymMat(){ 
00243     #ifndef NDEBUG
00244         setToNaN();
00245     #endif
00246     }
00247 
00248     SymMat(const SymMat& src) {
00249         updAsVec() = src.getAsVec();
00250     }
00251 
00252     SymMat& operator=(const SymMat& src) {    // no harm if src and 'this' are the same
00253         updAsVec() = src.getAsVec();
00254         return *this;
00255     }
00256 
00257     // Allow an explicit conversion from square Mat of right size, looking only at lower
00258     // elements and real part of diagonal elements.
00259     template <class EE, int CSS, int RSS>
00260     explicit SymMat(const Mat<M,M,EE,CSS,RSS>& m) {
00261         updDiag() = m.diag().real();
00262         for (int j=0; j<M; ++j)
00263             for (int i=j+1; i<M; ++i)
00264                 updEltLower(i,j) = m(i,j);
00265     }
00266 
00267     // We want an implicit conversion from a SymMat of the same length
00268     // and element type but with different spacings.
00269     template <int RSS> SymMat(const SymMat<M,E,RSS>& src) 
00270       { updAsVec() = src.getAsVec(); }
00271 
00272     // We want an implicit conversion from a SymMat of the same length
00273     // and *negated* element type, possibly with different spacings.
00274     template <int RSS> SymMat(const SymMat<M,ENeg,RSS>& src)
00275       { updAsVec() = src.getAsVec(); }
00276 
00277     // Construct a SymMat from a SymMat of the same dimensions, with any
00278     // spacings. Works as long as the element types are assignment compatible.
00279     template <class EE, int RSS> explicit SymMat(const SymMat<M,EE,RSS>& src)
00280       { updAsVec() = src.getAsVec(); }
00281 
00282     // Construction using an element repeats that element on the diagonal
00283     // but sets the rest of the matrix to zero.
00284     // TODO: diag should just use real part
00285     explicit SymMat(const E& e)
00286       { updDiag() = e; updLower() = E(0); }
00287 
00301 
00302     SymMat(const E& e0,
00303            const E& e1,const E& e2)
00304     {   assert(M==2); TDiag& d=updDiag(); TLower& l=updLower();
00305         d[0]=CNT<E>::real(e0); d[1]=CNT<E>::real(e2); 
00306         l[0]=e1; }
00307 
00308     SymMat(const E& e0,
00309            const E& e1,const E& e2,
00310            const E& e3,const E& e4, const E& e5)
00311     {   assert(M==3); TDiag& d=updDiag(); TLower& l=updLower();
00312         d[0]=CNT<E>::real(e0);d[1]=CNT<E>::real(e2);d[2]=CNT<E>::real(e5); 
00313         l[0]=e1;l[1]=e3;
00314         l[2]=e4; }
00315 
00316     SymMat(const E& e0,
00317            const E& e1,const E& e2,
00318            const E& e3,const E& e4, const E& e5,
00319            const E& e6,const E& e7, const E& e8, const E& e9)
00320     {   assert(M==4); TDiag& d=updDiag(); TLower& l=updLower();
00321         d[0]=CNT<E>::real(e0);d[1]=CNT<E>::real(e2);d[2]=CNT<E>::real(e5);d[3]=CNT<E>::real(e9); 
00322         l[0]=e1;l[1]=e3;l[2]=e6;
00323         l[3]=e4;l[4]=e7;
00324         l[5]=e8; }
00325 
00326     SymMat(const E& e0,
00327            const E& e1, const E& e2,
00328            const E& e3, const E& e4,  const E& e5,
00329            const E& e6, const E& e7,  const E& e8,  const E& e9,
00330            const E& e10,const E& e11, const E& e12, const E& e13, const E& e14)
00331     {   assert(M==5); TDiag& d=updDiag(); TLower& l=updLower();
00332         d[0]=CNT<E>::real(e0);d[1]=CNT<E>::real(e2);d[2]=CNT<E>::real(e5);d[3]=CNT<E>::real(e9);d[4]=CNT<E>::real(e14); 
00333         l[0]=e1;l[1]=e3;l[2]=e6;l[3]=e10;
00334         l[4]=e4;l[5]=e7;l[6]=e11;
00335         l[7]=e8;l[8]=e12;
00336         l[9]=e13; }
00337 
00338     SymMat(const E& e0,
00339            const E& e1, const E& e2,
00340            const E& e3, const E& e4,  const E& e5,
00341            const E& e6, const E& e7,  const E& e8,  const E& e9,
00342            const E& e10,const E& e11, const E& e12, const E& e13, const E& e14,
00343            const E& e15,const E& e16, const E& e17, const E& e18, const E& e19, const E& e20)
00344     {   assert(M==6); TDiag& d=updDiag(); TLower& l=updLower();
00345         d[0]=CNT<E>::real(e0);d[1]=CNT<E>::real(e2);d[2]=CNT<E>::real(e5);
00346         d[3]=CNT<E>::real(e9);d[4]=CNT<E>::real(e14);d[5]=CNT<E>::real(e20); 
00347         l[0] =e1; l[1] =e3; l[2] =e6;  l[3]=e10; l[4]=e15;
00348         l[5] =e4; l[6] =e7; l[7] =e11; l[8]=e16;
00349         l[9] =e8; l[10]=e12;l[11]=e17;
00350         l[12]=e13;l[13]=e18;
00351         l[14]=e19; }
00352 
00353     // Construction from a pointer to anything assumes we're pointing
00354     // at a packed element list of the right length, providing the
00355     // lower triangle in row order, so a b c d e f means
00356     //      a
00357     //      b c
00358     //      d e f
00359     //      g h i j
00360     // This has to be mapped to our diagonal/lower layout, which in
00361     // the above example will be:
00362     //      [a c f j][b d g e h i]
00363     //
00364     // In the input layout, the i'th row begins at element i(i+1)/2,
00365     // so diagonals are at i(i+1)/2 + i, while lower
00366     // elements (i,j; i>j) are at i(i+1)/2 + j.
00367     template <class EE> explicit SymMat(const EE* p) {
00368         assert(p);
00369         for (int i=0; i<M; ++i) {
00370             const int rowStart = (i*(i+1))/2;
00371             updDiag()[i] = CNT<EE>::real(p[rowStart + i]);
00372             for (int j=0; j<i; ++j)
00373                 updEltLower(i,j) = p[rowStart + j];
00374         }
00375     }
00376 
00377     // This is the same thing except as an assignment to pointer rather
00378     // than a constructor from pointer.
00379     template <class EE> SymMat& operator=(const EE* p) {
00380         assert(p);
00381         for (int i=0; i<M; ++i) {
00382             const int rowStart = (i*(i+1))/2;
00383             updDiag()[i] = CNT<EE>::real(p[rowStart + i]);
00384             for (int j=0; j<i; ++j)
00385                 updEltLower(i,j) = p[rowStart + j];
00386         }
00387         return *this;
00388     }
00389 
00390     // Assignment works similarly to copy -- if the lengths match,
00391     // go element-by-element. Otherwise, zero and then copy to each 
00392     // diagonal element.
00393     template <class EE, int RSS> SymMat& operator=(const SymMat<M,EE,RSS>& mm) {
00394         updAsVec() = mm.getAsVec();
00395         return *this;
00396     }
00397 
00398 
00399     // Conforming assignment ops
00400     template <class EE, int RSS> SymMat& 
00401     operator+=(const SymMat<M,EE,RSS>& mm) {
00402         updAsVec() += mm.getAsVec();
00403         return *this;
00404     }
00405     template <class EE, int RSS> SymMat&
00406     operator+=(const SymMat<M,negator<EE>,RSS>& mm) {
00407         updAsVec() -= -mm.getAsVec();   // negation of negator is free
00408         return *this;
00409     }
00410 
00411     template <class EE, int RSS> SymMat&
00412     operator-=(const SymMat<M,EE,RSS>& mm) {
00413         updAsVec() -= mm.getAsVec();
00414         return *this;
00415     }
00416     template <class EE, int RSS> SymMat&
00417     operator-=(const SymMat<M,negator<EE>,RSS>& mm) {
00418         updAsVec() += -mm.getAsVec();   // negation of negator is free
00419         return *this;
00420     }
00421 
00422     // In place matrix multiply can only be done when the RHS matrix is the same
00423     // size and the elements are also *= compatible.
00424     template <class EE, int RSS> SymMat&
00425     operator*=(const SymMat<M,EE,RSS>& mm) {
00426         assert(false); // TODO
00427         return *this;
00428     }
00429 
00430     // Conforming binary ops with 'this' on left, producing new packed result.
00431     // Cases: sy=sy+-sy, m=sy+-m, m=sy*sy, m=sy*m, v=sy*v
00432 
00433     // sy= this + sy
00434     template <class E2, int RS2> 
00435     typename Result<SymMat<M,E2,RS2> >::Add
00436     conformingAdd(const SymMat<M,E2,RS2>& r) const {
00437         return typename Result<SymMat<M,E2,RS2> >::Add
00438             (getAsVec().conformingAdd(r.getAsVec()));
00439     }
00440     // m= this - m
00441     template <class E2, int RS2> 
00442     typename Result<SymMat<M,E2,RS2> >::Sub
00443     conformingSubtract(const SymMat<M,E2,RS2>& r) const {
00444         return typename Result<SymMat<M,E2,RS2> >::Sub
00445             (getAsVec().conformingSubtract(r.getAsVec()));
00446     }
00447 
00448     // TODO: need the rest of the SymMat operators
00449     
00450     // Must be i >= j.
00451     const E& operator()(int i,int j) const 
00452       { return i==j ? getDiag()[i] : getEltLower(i,j); }
00453     E& operator()(int i,int j)
00454       { return i==j ? updDiag()[i] : updEltLower(i,j); }
00455 
00456     // This is the scalar Frobenius norm.
00457     ScalarNormSq normSqr() const { return scalarNormSqr(); }
00458     typename CNT<ScalarNormSq>::TSqrt 
00459         norm() const { return CNT<ScalarNormSq>::sqrt(scalarNormSqr()); }
00460 
00461     // There is no conventional meaning for normalize() applied to a matrix. We
00462     // choose to define it as follows:
00463     // If the elements of this SymMat are scalars, the result is what you get by
00464     // dividing each element by the Frobenius norm() calculated above. If the elements are
00465     // *not* scalars, then the elements are *separately* normalized.
00466     //
00467     // Normalize returns a matrix of the same dimension but in new, packed storage
00468     // and with a return type that does not include negator<> even if the original
00469     // SymMat<> does, because we can eliminate the negation here almost for free.
00470     // But we can't standardize (change conjugate to complex) for free, so we'll retain
00471     // conjugates if there are any.
00472     TNormalize normalize() const {
00473         if (CNT<E>::IsScalar) {
00474             return castAwayNegatorIfAny() / (SignInterpretation*norm());
00475         } else {
00476             TNormalize elementwiseNormalized;
00477             // punt to the equivalent Vec to get the elements normalized
00478             elementwiseNormalized.updAsVec() = getAsVec().normalize();
00479             return elementwiseNormalized;
00480         }
00481     }
00482 
00483     TInvert invert() const {assert(false); return TInvert();} // TODO default inversion
00484 
00485     const SymMat& operator+() const { return *this; }
00486     const TNeg&   operator-() const { return negate(); }
00487     TNeg&         operator-()       { return updNegate(); }
00488     const THerm&  operator~() const { return transpose(); }
00489     THerm&        operator~()       { return updTranspose(); }
00490 
00491     const TNeg&  negate() const { return *reinterpret_cast<const TNeg*>(this); }
00492     TNeg&        updNegate()    { return *reinterpret_cast<TNeg*>(this); }
00493 
00494     const THerm& transpose()    const { return *reinterpret_cast<const THerm*>(this); }
00495     THerm&       updTranspose()       { return *reinterpret_cast<THerm*>(this); }
00496 
00497     const TPosTrans& positionalTranspose() const
00498         { return *reinterpret_cast<const TPosTrans*>(this); }
00499     TPosTrans&       updPositionalTranspose()
00500         { return *reinterpret_cast<TPosTrans*>(this); }
00501 
00502     const TReal& real() const { return *reinterpret_cast<const TReal*>(this); }
00503     TReal&       real()       { return *reinterpret_cast<      TReal*>(this); }
00504 
00505     // Had to contort these routines to get them through VC++ 7.net
00506     const TImag& imag()    const { 
00507         const int offs = ImagOffset;
00508         const EImag* p = reinterpret_cast<const EImag*>(this);
00509         return *reinterpret_cast<const TImag*>(p+offs);
00510     }
00511     TImag& imag() { 
00512         const int offs = ImagOffset;
00513         EImag* p = reinterpret_cast<EImag*>(this);
00514         return *reinterpret_cast<TImag*>(p+offs);
00515     }
00516 
00517     const TWithoutNegator& castAwayNegatorIfAny() const {return *reinterpret_cast<const TWithoutNegator*>(this);}
00518     TWithoutNegator&       updCastAwayNegatorIfAny()    {return *reinterpret_cast<TWithoutNegator*>(this);}
00519 
00520     // These are elementwise binary operators, (this op ee) by default but (ee op this) if
00521     // 'FromLeft' appears in the name. The result is a packed SymMat<M> but the element type
00522     // may change. These are mostly used to implement global operators.
00523     // We call these "scalar" operators but actually the "scalar" can be a composite type.
00524 
00525     //TODO: consider converting 'e' to Standard Numbers as precalculation and changing
00526     // return type appropriately.
00527     template <class EE> SymMat<M, typename CNT<E>::template Result<EE>::Mul>
00528     scalarMultiply(const EE& e) const {
00529         SymMat<M, typename CNT<E>::template Result<EE>::Mul> result;
00530         result.updAsVec() = getAsVec().scalarMultiply(e);
00531         return result;
00532     }
00533     template <class EE> SymMat<M, typename CNT<EE>::template Result<E>::Mul>
00534     scalarMultiplyFromLeft(const EE& e) const {
00535         SymMat<M, typename CNT<EE>::template Result<E>::Mul> result;
00536         result.updAsVec() = getAsVec().scalarMultiplyFromLeft(e);
00537         return result;
00538     }
00539 
00540     // TODO: should precalculate and store 1/e, while converting to Standard Numbers. Note
00541     // that return type should change appropriately.
00542     template <class EE> SymMat<M, typename CNT<E>::template Result<EE>::Dvd>
00543     scalarDivide(const EE& e) const {
00544         SymMat<M, typename CNT<E>::template Result<EE>::Dvd> result;
00545         result.updAsVec() = getAsVec().scalarDivide(e);
00546         return result;
00547     }
00548     template <class EE> SymMat<M, typename CNT<EE>::template Result<E>::Dvd>
00549     scalarDivideFromLeft(const EE& e) const {
00550         SymMat<M, typename CNT<EE>::template Result<E>::Dvd> result;
00551         result.updAsVec() = getAsVec().scalarDivideFromLeft(e);
00552         return result;
00553     }
00554 
00555     // Additive ops involving a scalar update only the diagonal
00556     template <class EE> SymMat<M, typename CNT<E>::template Result<EE>::Add>
00557     scalarAdd(const EE& e) const {
00558         SymMat<M, typename CNT<E>::template Result<EE>::Add> result(*this);
00559         result.updDiag() += e;
00560         return result;
00561     }
00562     // Add is commutative, so no 'FromLeft'.
00563 
00564     template <class EE> SymMat<M, typename CNT<E>::template Result<EE>::Sub>
00565     scalarSubtract(const EE& e) const {
00566         SymMat<M, typename CNT<E>::template Result<EE>::Sub> result(*this);
00567         result.diag() -= e;
00568         return result;
00569     }
00570     // This is s-m; negate m and add s to diagonal
00571     // TODO: Should do something clever with negation here.
00572     template <class EE> SymMat<M, typename CNT<EE>::template Result<E>::Sub>
00573     scalarSubtractFromLeft(const EE& e) const {
00574         SymMat<M, typename CNT<EE>::template Result<E>::Sub> result(-(*this));
00575         result.diag() += e;
00576         return result;
00577     }
00578 
00579     // Generic assignments for any element type not listed explicitly, including scalars.
00580     // These are done repeatedly for each element and only work if the operation can
00581     // be performed leaving the original element type.
00582     template <class EE> SymMat& operator =(const EE& e) {return scalarEq(e);}
00583     template <class EE> SymMat& operator+=(const EE& e) {return scalarPlusEq(e);}
00584     template <class EE> SymMat& operator-=(const EE& e) {return scalarMinusEq(e);}
00585     template <class EE> SymMat& operator*=(const EE& e) {return scalarTimesEq(e);}
00586     template <class EE> SymMat& operator/=(const EE& e) {return scalarDivideEq(e);}
00587 
00588     // Generalized scalar assignment & computed assignment methods. These will work
00589     // for any assignment-compatible element, not just scalars.
00590     template <class EE> SymMat& scalarEq(const EE& ee)
00591       { updDiag() = ee; updLower() = E(0); return *this; }
00592     template <class EE> SymMat& scalarPlusEq(const EE& ee)
00593       { updDiag() += ee; return *this; }
00594     template <class EE> SymMat& scalarMinusEq(const EE& ee)
00595       { updDiag() -= ee; return *this; }
00596 
00597     // this is m = s-m; negate off diagonal, do d=s-d for each diagonal element d
00598     template <class EE> SymMat& scalarMinusEqFromLeft(const EE& ee)
00599       { updLower() *= E(-1); updDiag().scalarMinusEqFromLeft(ee); return *this; }
00600 
00601     template <class EE> SymMat& scalarTimesEq(const EE& ee)
00602       { updAsVec() *= ee; return *this; }
00603     template <class EE> SymMat& scalarTimesEqFromLeft(const EE& ee)
00604       { updAsVec().scalarTimesEqFromLeft(ee); return *this; }
00605     template <class EE> SymMat& scalarDivideEq(const EE& ee)
00606       { updAsVec() /= ee; return *this; }
00607     template <class EE> SymMat& scalarDivideEqFromLeft(const EE& ee)
00608       { updAsVec().scalarDivideEqFromLeft(ee); return *this; } 
00609 
00610     void setToNaN() { updAsVec().setToNaN(); }
00611 
00612     // These assume we are given a pointer to d[0] of a SymMat<M,E,RS> like this one.
00613     static const SymMat& getAs(const ELT* p)  {return *reinterpret_cast<const SymMat*>(p);}
00614     static SymMat&       updAs(ELT* p)        {return *reinterpret_cast<SymMat*>(p);}
00615 
00616     // Note packed spacing
00617     static TPacked getNaN() {
00618         return TPacked(CNT<typename TPacked::TDiag>::getNaN(),
00619                        CNT<typename TPacked::TLower>::getNaN());
00620     }
00621 
00622     const TDiag&  getDiag()  const {return TDiag::getAs(d);}
00623     TDiag&        updDiag()        {return TDiag::updAs(d);}
00624 
00625     // Conventional synonym
00626     const TDiag& diag() const {return getDiag();}
00627     TDiag&       diag()       {return updDiag();}
00628 
00629     const TLower& getLower() const {return TLower::getAs(&d[M*RS]);}
00630     TLower&       updLower()       {return TLower::updAs(&d[M*RS]);}
00631 
00632     const TUpper& getUpper() const {return reinterpret_cast<const TUpper&>(getLower());}
00633     TUpper&       updUpper()       {return reinterpret_cast<      TUpper&>(updLower());}
00634 
00635     const TAsVec& getAsVec() const {return TAsVec::getAs(d);}
00636     TAsVec&       updAsVec()       {return TAsVec::updAs(d);}
00637 
00638     // must be i > j
00639     const E& getEltLower(int i, int j) const {return getLower()[lowerIx(i,j)];}
00640     E&       updEltLower(int i, int j)       {return updLower()[lowerIx(i,j)];}
00641 
00642     // must be i < j
00643     const EHerm& getEltUpper(int i, int j) const {return getUpper()[lowerIx(j,i)];}
00644     EHerm&       updEltUpper(int i, int j)       {return updUpper()[lowerIx(j,i)];}
00645     
00646     TRow sum() const {
00647         TRow temp(~getDiag());
00648         for (int i = 1; i < M; ++i)
00649             for (int j = 0; j < i; ++j) {
00650                 E value = getEltLower(i, j);;
00651                 temp[i] += value;
00652                 temp[j] += value;
00653             }
00654         return temp;
00655     }
00656 
00657 private:
00658     E d[NActualElements];
00659 
00660     SymMat(const TDiag& di, const TLower& low) {
00661         updDiag() = di; updLower() = low;
00662     }
00663 
00664     explicit SymMat(const TAsVec& v) {
00665         TAsVec::updAs(d) = v;
00666     }
00667     
00668     // Convert i,j with i>j to an index in "lower". This also gives the
00669     // right index for "upper" with i and j reversed.
00670     static int lowerIx(int i, int j) {
00671         assert(0 <= j && j < i && i < M);
00672         return (i-j-1) + j*(M-1) - (j*(j-1))/2;
00673     }
00674 
00675     template <int MM, class EE, int RSS> friend class SymMat;
00676 };
00677 
00679 // Global operators involving two symmetric matrices. //
00680 //   sy=sy+sy, sy=sy-sy, m=sy*sy, sy==sy, sy!=sy      //
00682 
00683 // m3 = m1 + m2 where all m's have the same dimension M. 
00684 template <int M, class E1, int S1, class E2, int S2> inline
00685 typename SymMat<M,E1,S1>::template Result< SymMat<M,E2,S2> >::Add
00686 operator+(const SymMat<M,E1,S1>& l, const SymMat<M,E2,S2>& r) {
00687     return SymMat<M,E1,S1>::template Result< SymMat<M,E2,S2> >
00688         ::AddOp::perform(l,r);
00689 }
00690 
00691 // m3 = m1 - m2 where all m's have the same dimension M. 
00692 template <int M, class E1, int S1, class E2, int S2> inline
00693 typename SymMat<M,E1,S1>::template Result< SymMat<M,E2,S2> >::Sub
00694 operator-(const SymMat<M,E1,S1>& l, const SymMat<M,E2,S2>& r) {
00695     return SymMat<M,E1,S1>::template Result< SymMat<M,E2,S2> >
00696         ::SubOp::perform(l,r);
00697 }
00698 
00699 // m3 = m1 * m2 where all m's have the same dimension M. 
00700 // The result will not be symmetric.
00701 template <int M, class E1, int S1, class E2, int S2> inline
00702 typename SymMat<M,E1,S1>::template Result< SymMat<M,E2,S2> >::Mul
00703 operator-(const SymMat<M,E1,S1>& l, const SymMat<M,E2,S2>& r) {
00704     return SymMat<M,E1,S1>::template Result< SymMat<M,E2,S2> >
00705         ::MulOp::perform(l,r);
00706 }
00707 
00708 // bool = m1 == m2, m1 and m2 have the same dimension M
00709 template <int M, class E1, int S1, class E2, int S2> inline bool
00710 operator==(const SymMat<M,E1,S1>& l, const SymMat<M,E2,S2>& r) {
00711     return l.getAsVec() == r.getAsVec();
00712 }
00713 
00714 // bool = m1 == m2, m1 and m2 have the same dimension M
00715 template <int M, class E1, int S1, class E2, int S2> inline bool
00716 operator!=(const SymMat<M,E1,S1>& l, const SymMat<M,E2,S2>& r) {return !(l==r);} 
00717 
00719 // Global operators involving a SymMat and a scalar. //
00721 
00722 // SCALAR MULTIPLY
00723 
00724 // m = m*real, real*m 
00725 template <int M, class E, int S> inline
00726 typename SymMat<M,E,S>::template Result<float>::Mul
00727 operator*(const SymMat<M,E,S>& l, const float& r)
00728   { return SymMat<M,E,S>::template Result<float>::MulOp::perform(l,r); }
00729 template <int M, class E, int S> inline
00730 typename SymMat<M,E,S>::template Result<float>::Mul
00731 operator*(const float& l, const SymMat<M,E,S>& r) {return r*l;}
00732 
00733 template <int M, class E, int S> inline
00734 typename SymMat<M,E,S>::template Result<double>::Mul
00735 operator*(const SymMat<M,E,S>& l, const double& r)
00736   { return SymMat<M,E,S>::template Result<double>::MulOp::perform(l,r); }
00737 template <int M, class E, int S> inline
00738 typename SymMat<M,E,S>::template Result<double>::Mul
00739 operator*(const double& l, const SymMat<M,E,S>& r) {return r*l;}
00740 
00741 template <int M, class E, int S> inline
00742 typename SymMat<M,E,S>::template Result<long double>::Mul
00743 operator*(const SymMat<M,E,S>& l, const long double& r)
00744   { return SymMat<M,E,S>::template Result<long double>::MulOp::perform(l,r); }
00745 template <int M, class E, int S> inline
00746 typename SymMat<M,E,S>::template Result<long double>::Mul
00747 operator*(const long double& l, const SymMat<M,E,S>& r) {return r*l;}
00748 
00749 // m = m*int, int*m -- just convert int to m's precision float
00750 template <int M, class E, int S> inline
00751 typename SymMat<M,E,S>::template Result<typename CNT<E>::Precision>::Mul
00752 operator*(const SymMat<M,E,S>& l, int r) {return l * (typename CNT<E>::Precision)r;}
00753 template <int M, class E, int S> inline
00754 typename SymMat<M,E,S>::template Result<typename CNT<E>::Precision>::Mul
00755 operator*(int l, const SymMat<M,E,S>& r) {return r * (typename CNT<E>::Precision)l;}
00756 
00757 // Complex, conjugate, and negator are all easy to templatize.
00758 
00759 // m = m*complex, complex*m
00760 template <int M, class E, int S, class R> inline
00761 typename SymMat<M,E,S>::template Result<std::complex<R> >::Mul
00762 operator*(const SymMat<M,E,S>& l, const std::complex<R>& r)
00763   { return SymMat<M,E,S>::template Result<std::complex<R> >::MulOp::perform(l,r); }
00764 template <int M, class E, int S, class R> inline
00765 typename SymMat<M,E,S>::template Result<std::complex<R> >::Mul
00766 operator*(const std::complex<R>& l, const SymMat<M,E,S>& r) {return r*l;}
00767 
00768 // m = m*conjugate, conjugate*m (convert conjugate->complex)
00769 template <int M, class E, int S, class R> inline
00770 typename SymMat<M,E,S>::template Result<std::complex<R> >::Mul
00771 operator*(const SymMat<M,E,S>& l, const conjugate<R>& r) {return l*(std::complex<R>)r;}
00772 template <int M, class E, int S, class R> inline
00773 typename SymMat<M,E,S>::template Result<std::complex<R> >::Mul
00774 operator*(const conjugate<R>& l, const SymMat<M,E,S>& r) {return r*(std::complex<R>)l;}
00775 
00776 // m = m*negator, negator*m: convert negator to standard number
00777 template <int M, class E, int S, class R> inline
00778 typename SymMat<M,E,S>::template Result<typename negator<R>::StdNumber>::Mul
00779 operator*(const SymMat<M,E,S>& l, const negator<R>& r) {return l * (typename negator<R>::StdNumber)(R)r;}
00780 template <int M, class E, int S, class R> inline
00781 typename SymMat<M,E,S>::template Result<typename negator<R>::StdNumber>::Mul
00782 operator*(const negator<R>& l, const SymMat<M,E,S>& r) {return r * (typename negator<R>::StdNumber)(R)l;}
00783 
00784 
00785 // SCALAR DIVIDE. This is a scalar operation when the scalar is on the right,
00786 // but when it is on the left it means scalar * pseudoInverse(mat), which is 
00787 // a similar symmetric matrix.
00788 
00789 // m = m/real, real/m 
00790 template <int M, class E, int S> inline
00791 typename SymMat<M,E,S>::template Result<float>::Dvd
00792 operator/(const SymMat<M,E,S>& l, const float& r)
00793   { return SymMat<M,E,S>::template Result<float>::DvdOp::perform(l,r); }
00794 template <int M, class E, int S> inline
00795 typename CNT<float>::template Result<SymMat<M,E,S> >::Dvd
00796 operator/(const float& l, const SymMat<M,E,S>& r)
00797   { return CNT<float>::template Result<SymMat<M,E,S> >::DvdOp::perform(l,r); }
00798 
00799 template <int M, class E, int S> inline
00800 typename SymMat<M,E,S>::template Result<double>::Dvd
00801 operator/(const SymMat<M,E,S>& l, const double& r)
00802   { return SymMat<M,E,S>::template Result<double>::DvdOp::perform(l,r); }
00803 template <int M, class E, int S> inline
00804 typename CNT<double>::template Result<SymMat<M,E,S> >::Dvd
00805 operator/(const double& l, const SymMat<M,E,S>& r)
00806   { return CNT<double>::template Result<SymMat<M,E,S> >::DvdOp::perform(l,r); }
00807 
00808 template <int M, class E, int S> inline
00809 typename SymMat<M,E,S>::template Result<long double>::Dvd
00810 operator/(const SymMat<M,E,S>& l, const long double& r)
00811   { return SymMat<M,E,S>::template Result<long double>::DvdOp::perform(l,r); }
00812 template <int M, class E, int S> inline
00813 typename CNT<long double>::template Result<SymMat<M,E,S> >::Dvd
00814 operator/(const long double& l, const SymMat<M,E,S>& r)
00815   { return CNT<long double>::template Result<SymMat<M,E,S> >::DvdOp::perform(l,r); }
00816 
00817 // m = m/int, int/m -- just convert int to m's precision float
00818 template <int M, class E, int S> inline
00819 typename SymMat<M,E,S>::template Result<typename CNT<E>::Precision>::Dvd
00820 operator/(const SymMat<M,E,S>& l, int r) {return l / (typename CNT<E>::Precision)r;}
00821 template <int M, class E, int S> inline
00822 typename CNT<typename CNT<E>::Precision>::template Result<SymMat<M,E,S> >::Dvd
00823 operator/(int l, const SymMat<M,E,S>& r) {return (typename CNT<E>::Precision)l / r;}
00824 
00825 
00826 // Complex, conjugate, and negator are all easy to templatize.
00827 
00828 // m = m/complex, complex/m
00829 template <int M, class E, int S, class R> inline
00830 typename SymMat<M,E,S>::template Result<std::complex<R> >::Dvd
00831 operator/(const SymMat<M,E,S>& l, const std::complex<R>& r)
00832   { return SymMat<M,E,S>::template Result<std::complex<R> >::DvdOp::perform(l,r); }
00833 template <int M, class E, int S, class R> inline
00834 typename CNT<std::complex<R> >::template Result<SymMat<M,E,S> >::Dvd
00835 operator/(const std::complex<R>& l, const SymMat<M,E,S>& r)
00836   { return CNT<std::complex<R> >::template Result<SymMat<M,E,S> >::DvdOp::perform(l,r); }
00837 
00838 // m = m/conjugate, conjugate/m (convert conjugate->complex)
00839 template <int M, class E, int S, class R> inline
00840 typename SymMat<M,E,S>::template Result<std::complex<R> >::Dvd
00841 operator/(const SymMat<M,E,S>& l, const conjugate<R>& r) {return l/(std::complex<R>)r;}
00842 template <int M, class E, int S, class R> inline
00843 typename CNT<std::complex<R> >::template Result<SymMat<M,E,S> >::Dvd
00844 operator/(const conjugate<R>& l, const SymMat<M,E,S>& r) {return (std::complex<R>)l/r;}
00845 
00846 // m = m/negator, negator/m: convert negator to number
00847 template <int M, class E, int S, class R> inline
00848 typename SymMat<M,E,S>::template Result<typename negator<R>::StdNumber>::Dvd
00849 operator/(const SymMat<M,E,S>& l, const negator<R>& r) {return l/(typename negator<R>::StdNumber)(R)r;}
00850 template <int M, class E, int S, class R> inline
00851 typename CNT<R>::template Result<SymMat<M,E,S> >::Dvd
00852 operator/(const negator<R>& l, const SymMat<M,E,S>& r) {return (typename negator<R>::StdNumber)(R)l/r;}
00853 
00854 
00855 // Add and subtract are odd as scalar ops. They behave as though the
00856 // scalar stands for a conforming matrix whose diagonal elements are that,
00857 // scalar and then a normal matrix add or subtract is done.
00858 
00859 // SCALAR ADD
00860 
00861 // m = m+real, real+m 
00862 template <int M, class E, int S> inline
00863 typename SymMat<M,E,S>::template Result<float>::Add
00864 operator+(const SymMat<M,E,S>& l, const float& r)
00865   { return SymMat<M,E,S>::template Result<float>::AddOp::perform(l,r); }
00866 template <int M, class E, int S> inline
00867 typename SymMat<M,E,S>::template Result<float>::Add
00868 operator+(const float& l, const SymMat<M,E,S>& r) {return r+l;}
00869 
00870 template <int M, class E, int S> inline
00871 typename SymMat<M,E,S>::template Result<double>::Add
00872 operator+(const SymMat<M,E,S>& l, const double& r)
00873   { return SymMat<M,E,S>::template Result<double>::AddOp::perform(l,r); }
00874 template <int M, class E, int S> inline
00875 typename SymMat<M,E,S>::template Result<double>::Add
00876 operator+(const double& l, const SymMat<M,E,S>& r) {return r+l;}
00877 
00878 template <int M, class E, int S> inline
00879 typename SymMat<M,E,S>::template Result<long double>::Add
00880 operator+(const SymMat<M,E,S>& l, const long double& r)
00881   { return SymMat<M,E,S>::template Result<long double>::AddOp::perform(l,r); }
00882 template <int M, class E, int S> inline
00883 typename SymMat<M,E,S>::template Result<long double>::Add
00884 operator+(const long double& l, const SymMat<M,E,S>& r) {return r+l;}
00885 
00886 // m = m+int, int+m -- just convert int to m's precision float
00887 template <int M, class E, int S> inline
00888 typename SymMat<M,E,S>::template Result<typename CNT<E>::Precision>::Add
00889 operator+(const SymMat<M,E,S>& l, int r) {return l + (typename CNT<E>::Precision)r;}
00890 template <int M, class E, int S> inline
00891 typename SymMat<M,E,S>::template Result<typename CNT<E>::Precision>::Add
00892 operator+(int l, const SymMat<M,E,S>& r) {return r + (typename CNT<E>::Precision)l;}
00893 
00894 // Complex, conjugate, and negator are all easy to templatize.
00895 
00896 // m = m+complex, complex+m
00897 template <int M, class E, int S, class R> inline
00898 typename SymMat<M,E,S>::template Result<std::complex<R> >::Add
00899 operator+(const SymMat<M,E,S>& l, const std::complex<R>& r)
00900   { return SymMat<M,E,S>::template Result<std::complex<R> >::AddOp::perform(l,r); }
00901 template <int M, class E, int S, class R> inline
00902 typename SymMat<M,E,S>::template Result<std::complex<R> >::Add
00903 operator+(const std::complex<R>& l, const SymMat<M,E,S>& r) {return r+l;}
00904 
00905 // m = m+conjugate, conjugate+m (convert conjugate->complex)
00906 template <int M, class E, int S, class R> inline
00907 typename SymMat<M,E,S>::template Result<std::complex<R> >::Add
00908 operator+(const SymMat<M,E,S>& l, const conjugate<R>& r) {return l+(std::complex<R>)r;}
00909 template <int M, class E, int S, class R> inline
00910 typename SymMat<M,E,S>::template Result<std::complex<R> >::Add
00911 operator+(const conjugate<R>& l, const SymMat<M,E,S>& r) {return r+(std::complex<R>)l;}
00912 
00913 // m = m+negator, negator+m: convert negator to standard number
00914 template <int M, class E, int S, class R> inline
00915 typename SymMat<M,E,S>::template Result<typename negator<R>::StdNumber>::Add
00916 operator+(const SymMat<M,E,S>& l, const negator<R>& r) {return l + (typename negator<R>::StdNumber)(R)r;}
00917 template <int M, class E, int S, class R> inline
00918 typename SymMat<M,E,S>::template Result<typename negator<R>::StdNumber>::Add
00919 operator+(const negator<R>& l, const SymMat<M,E,S>& r) {return r + (typename negator<R>::StdNumber)(R)l;}
00920 
00921 // SCALAR SUBTRACT -- careful, not commutative.
00922 
00923 // m = m-real, real-m 
00924 template <int M, class E, int S> inline
00925 typename SymMat<M,E,S>::template Result<float>::Sub
00926 operator-(const SymMat<M,E,S>& l, const float& r)
00927   { return SymMat<M,E,S>::template Result<float>::SubOp::perform(l,r); }
00928 template <int M, class E, int S> inline
00929 typename CNT<float>::template Result<SymMat<M,E,S> >::Sub
00930 operator-(const float& l, const SymMat<M,E,S>& r)
00931   { return CNT<float>::template Result<SymMat<M,E,S> >::SubOp::perform(l,r); }
00932 
00933 template <int M, class E, int S> inline
00934 typename SymMat<M,E,S>::template Result<double>::Sub
00935 operator-(const SymMat<M,E,S>& l, const double& r)
00936   { return SymMat<M,E,S>::template Result<double>::SubOp::perform(l,r); }
00937 template <int M, class E, int S> inline
00938 typename CNT<double>::template Result<SymMat<M,E,S> >::Sub
00939 operator-(const double& l, const SymMat<M,E,S>& r)
00940   { return CNT<double>::template Result<SymMat<M,E,S> >::SubOp::perform(l,r); }
00941 
00942 template <int M, class E, int S> inline
00943 typename SymMat<M,E,S>::template Result<long double>::Sub
00944 operator-(const SymMat<M,E,S>& l, const long double& r)
00945   { return SymMat<M,E,S>::template Result<long double>::SubOp::perform(l,r); }
00946 template <int M, class E, int S> inline
00947 typename CNT<long double>::template Result<SymMat<M,E,S> >::Sub
00948 operator-(const long double& l, const SymMat<M,E,S>& r)
00949   { return CNT<long double>::template Result<SymMat<M,E,S> >::SubOp::perform(l,r); }
00950 
00951 // m = m-int, int-m // just convert int to m's precision float
00952 template <int M, class E, int S> inline
00953 typename SymMat<M,E,S>::template Result<typename CNT<E>::Precision>::Sub
00954 operator-(const SymMat<M,E,S>& l, int r) {return l - (typename CNT<E>::Precision)r;}
00955 template <int M, class E, int S> inline
00956 typename CNT<typename CNT<E>::Precision>::template Result<SymMat<M,E,S> >::Sub
00957 operator-(int l, const SymMat<M,E,S>& r) {return (typename CNT<E>::Precision)l - r;}
00958 
00959 
00960 // Complex, conjugate, and negator are all easy to templatize.
00961 
00962 // m = m-complex, complex-m
00963 template <int M, class E, int S, class R> inline
00964 typename SymMat<M,E,S>::template Result<std::complex<R> >::Sub
00965 operator-(const SymMat<M,E,S>& l, const std::complex<R>& r)
00966   { return SymMat<M,E,S>::template Result<std::complex<R> >::SubOp::perform(l,r); }
00967 template <int M, class E, int S, class R> inline
00968 typename CNT<std::complex<R> >::template Result<SymMat<M,E,S> >::Sub
00969 operator-(const std::complex<R>& l, const SymMat<M,E,S>& r)
00970   { return CNT<std::complex<R> >::template Result<SymMat<M,E,S> >::SubOp::perform(l,r); }
00971 
00972 // m = m-conjugate, conjugate-m (convert conjugate->complex)
00973 template <int M, class E, int S, class R> inline
00974 typename SymMat<M,E,S>::template Result<std::complex<R> >::Sub
00975 operator-(const SymMat<M,E,S>& l, const conjugate<R>& r) {return l-(std::complex<R>)r;}
00976 template <int M, class E, int S, class R> inline
00977 typename CNT<std::complex<R> >::template Result<SymMat<M,E,S> >::Sub
00978 operator-(const conjugate<R>& l, const SymMat<M,E,S>& r) {return (std::complex<R>)l-r;}
00979 
00980 // m = m-negator, negator-m: convert negator to standard number
00981 template <int M, class E, int S, class R> inline
00982 typename SymMat<M,E,S>::template Result<typename negator<R>::StdNumber>::Sub
00983 operator-(const SymMat<M,E,S>& l, const negator<R>& r) {return l-(typename negator<R>::StdNumber)(R)r;}
00984 template <int M, class E, int S, class R> inline
00985 typename CNT<R>::template Result<SymMat<M,E,S> >::Sub
00986 operator-(const negator<R>& l, const SymMat<M,E,S>& r) {return (typename negator<R>::StdNumber)(R)l-r;}
00987 
00988 
00989 // SymMat I/O
00990 template <int M, class E, int RS, class CHAR, class TRAITS> inline
00991 std::basic_ostream<CHAR,TRAITS>&
00992 operator<<(std::basic_ostream<CHAR,TRAITS>& o, const SymMat<M,E,RS>& m) {
00993     for (int i=0;i<M;++i) {
00994         o << std::endl << "[";
00995         for (int j=0; j<=i; ++j)         
00996             o << (j>0?" ":"") << m(i,j);
00997         for (int j=i+1; j<M; ++j)
00998             o << " *";
00999         o << "]";
01000     }
01001     if (M) o << std::endl;
01002     return o; 
01003 }
01004 
01005 template <int M, class E, int RS, class CHAR, class TRAITS> inline
01006 std::basic_istream<CHAR,TRAITS>&
01007 operator>>(std::basic_istream<CHAR,TRAITS>& is, SymMat<M,E,RS>& m) {
01008     // TODO: not sure how to do input yet
01009     assert(false);
01010     return is;
01011 }
01012 
01013 } //namespace SimTK
01014 
01015 
01016 #endif //SimTK_SIMMATRIX_SMALLMATRIX_SYMMAT_H_

Generated by  doxygen 1.6.2