MatrixHelper.h

Go to the documentation of this file.
00001 #ifndef SimTK_SIMMATRIX_MATRIX_HELPER_H_
00002 #define SimTK_SIMMATRIX_MATRIX_HELPER_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-9 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 
00045 #include "SimTKcommon/Scalar.h"
00046 
00047 #include <iostream>
00048 #include <cassert>
00049 #include <complex>
00050 #include <cstddef>
00051 #include <utility> // for std::pair
00052 #include <vector>
00053 
00054 namespace SimTK {
00055 
00056 
00057 template <class S> class MatrixHelperRep;
00058 class MatrixCharacter;
00059 class MatrixCommitment;
00060 
00061 //  ------------------------------ MatrixHelper --------------------------------
00062 
00085 
00086 //  ----------------------------------------------------------------------------
00087 template <class S> 
00088 class SimTK_SimTKCOMMON_EXPORT MatrixHelper {
00089     typedef MatrixHelper<S>                         This;
00090     typedef MatrixHelper<typename CNT<S>::TNeg>     ThisNeg;
00091     typedef MatrixHelper<typename CNT<S>::THerm>    ThisHerm;
00092 public: 
00093     typedef typename CNT<S>::Number     Number;     // strips off negator from S
00094     typedef typename CNT<S>::StdNumber  StdNumber;  // turns conjugate to complex
00095     typedef typename CNT<S>::Precision  Precision;  // strips off complex from StdNumber
00096 
00097     // no default constructor
00098     // copy constructor suppressed
00099 
00100     // Destructor eliminates MatrixHelperRep object if "this" owns it.
00101     ~MatrixHelper() {deleteRepIfOwner();}
00102 
00103     // Local types for directing us to the right constructor at compile time.
00104     class ShallowCopy   { };
00105     class DeepCopy      { };
00106     class TransposeView { };
00107     class DiagonalView  { };
00108     
00110         // "Owner" constructors //
00112     
00113     // 0x0, fully resizable, fully uncommitted.
00114     MatrixHelper(int esz, int cppEsz);
00115 
00116     // Default allocation for the given commitment.
00117     MatrixHelper(int esz, int cppEsz, const MatrixCommitment&);
00118 
00119     // Allocate a matrix that satisfies a given commitment and has a
00120     // particular initial size.
00121     MatrixHelper(int esz, int cppEsz, const MatrixCommitment&, int m, int n);
00122 
00123     // Copy constructor that produces a new owner whose logical shape and contents are
00124     // the same as the source, but with a possibly better storage layout. Data will
00125     // be contiguous in the copy regardless of how spread out it was in the source.
00126     // The second argument is just to disambiguate this constructor from similar ones.
00127     MatrixHelper(const MatrixCommitment&, const MatrixHelper& source, const DeepCopy&);
00128 
00129     // This has the same semantics as the previous DeepCopy constructor, except that
00130     // the source has negated elements with respect to S. The resulting logical shape
00131     // and logical values are identical to the source, meaning that the negation must
00132     // actually be performed here, using one flop for every meaningful source scalar.
00133     MatrixHelper(const MatrixCommitment&, const MatrixHelper<typename CNT<S>::TNeg>& source, const DeepCopy&);
00134    
00135     
00137         // External "View" constructors //
00139 
00140     // These constructors produce a matrix which provides a view of externally-allocated
00141     // data, which is known to have the given MatrixCharacter. There is also provision
00142     // for a "spacing" parameter which defines gaps in the supplied data, although
00143     // the interpretation of that parameter depends on the MatrixCharacter (typically
00144     // it is the leading dimension for a matrix or the stride for a vector). Note
00145     // that spacing is interpreted as "number of scalars between adjacent elements"
00146     // which has the usual Lapack interpretation if the elements are scalars but
00147     // can be used for C++ vs. Simmatrix packing differences for composite elements.
00148     // The resulting handle is *not* the owner of the data, so destruction of the 
00149     // handle does not delete the data.
00150 
00151     // Create a read-only view into existing data.
00152     MatrixHelper(int esz, int cppEsz, const MatrixCommitment&,
00153                  const MatrixCharacter&, int spacing, const S* data);
00154     // Create a writable view into existing data.
00155     MatrixHelper(int esz, int cppEsz, const MatrixCommitment&,
00156                  const MatrixCharacter&, int spacing, S* data);
00157 
00158 
00160         // Matrix "View" constructors //
00162 
00163     // Matrix view constructors, read only and writable. Use these for submatrices,
00164     // rows, and columns. Indices are by *element* and these may consist of multiple
00165     // scalars of type template parameter S.
00166 
00167     // "Block" views
00168     MatrixHelper(const MatrixCommitment&, const MatrixHelper&, int i, int j, int nrow, int ncol);
00169     MatrixHelper(const MatrixCommitment&, MatrixHelper&,       int i, int j, int nrow, int ncol); 
00170 
00171     // "Transpose" views (note that this is Hermitian transpose; i.e., element
00172     // type is conjugated).
00173     MatrixHelper(const MatrixCommitment&, const MatrixHelper<typename CNT<S>::THerm>&, 
00174                  const TransposeView&);    // a read only transposed view
00175     MatrixHelper(const MatrixCommitment&, MatrixHelper<typename CNT<S>::THerm>&,       
00176                  const TransposeView&);    // a writable transposed view
00177 
00178     // "Diagonal" views.
00179     MatrixHelper(const MatrixCommitment&, const MatrixHelper&, const DiagonalView&); // a read only diagonal view
00180     MatrixHelper(const MatrixCommitment&, MatrixHelper&,       const DiagonalView&); // a writable diagonal view
00181 
00182     // "Indexed" view of a 1d matrix.
00183     MatrixHelper(const MatrixCommitment&, const MatrixHelper&, int n, const int* indices);
00184     MatrixHelper(const MatrixCommitment&, MatrixHelper&, int n, const int* indices);
00185 
00186     // These invoke the previous constructors but with friendlier index source.
00187     MatrixHelper(const MatrixCommitment& mc, const MatrixHelper& h, const std::vector<int>& indices)
00188     {   new (this) MatrixHelper(mc, h, (int)indices.size(), &indices[0]); }
00189     MatrixHelper(const MatrixCommitment& mc, MatrixHelper& h, const std::vector<int>& indices)
00190     {   new (this) MatrixHelper(mc, h, (int)indices.size(), &indices[0]); }
00191 
00192     // "Copy" an existing MatrixHelper by making a new view into the same data. 
00193     // The const form loses writability, non-const retains same writability status 
00194     // as source. If the source is already a view then the destination will have
00195     // an identical element filter so that the logical shape and values are the
00196     // same for both source and copy. (The second argument is used just to 
00197     // disambiguate this constructor from similar ones.)
00198     MatrixHelper(const MatrixCommitment&, const MatrixHelper& source, const ShallowCopy&);
00199     MatrixHelper(const MatrixCommitment&, MatrixHelper&       source, const ShallowCopy&);
00200 
00202         // Copy assignment //
00204 
00205     // Behavior of copy assignment depends on whether "this" is an owner or view. If
00206     // it's an owner it is resized and ends up a new, dense copy of "source" just as
00207     // for the DeepCopy constructor above. If "this" is a writable view, sizes must match 
00208     // and we copy elements from "source" to corresponding elements of "this". If
00209     // "this" is not writable then the operation will fail.
00210     MatrixHelper& copyAssign(const MatrixHelper& source);
00211 
00212     // Same as copyAssign() except the source has element types which are logically
00213     // negated from the element types of "this". Since the result must have the
00214     // same value as the source, this requires that all the copied elements are
00215     // actually negated at a cost of one flop per scalar.
00216     MatrixHelper& negatedCopyAssign(const MatrixHelper<typename CNT<S>::TNeg>&);
00217 
00219         // View assignment //
00221 
00222     // View assign always disconnects this helper from whatever view & data
00223     // it used to have (destructing as appropriate) and then makes it a new view
00224     // of the source. Writability is lost if the source is const, otherwise 
00225     // writability is inherited from the source.
00226     MatrixHelper& readOnlyViewAssign(const MatrixHelper& source);
00227     MatrixHelper& writableViewAssign(MatrixHelper&       source);
00228 
00229     // Restore helper to its just-constructed state. We forget everything except
00230     // the element size and handle commitment, which can never change. The
00231     // allocated helper will be the same as if we had just default-constructed
00232     // using the current commitment.
00233     void clear();
00234 
00235     // Return true if there is currently no data associated with this handle.
00236     bool isClear() const;
00237 
00238     // Using *element* indices, obtain a pointer to the beginning of a particular
00239     // element. This is always a slow operation compared to raw array access;
00240     // use sparingly. These will fail if the indices are outside the stored
00241     // portion of the matrix. Use getAnyElt() if you want something that always
00242     // works.
00243     const S* getElt(int i, int j) const;
00244     S*       updElt(int i, int j);
00245 
00246     // Faster for 1-d matrices (vectors) if you know you have one.
00247     const S* getElt(int i) const;
00248     S*       updElt(int i);
00249 
00250     // This returns the correct value for any element within the logical dimensions
00251     // of the matrix. In some cases it has to compute the value; in all cases
00252     // it has to copy it.
00253     void getAnyElt(int i, int j, S* value) const;
00254 
00255     // Faster for 1-d matrices (vectors) if you know you have one.
00256     void getAnyElt(int i, S* value) const;
00257 
00258     // Add up all the *elements*, returning the answer in the element given
00259     // by pointer to its first scalar.
00260     void sum(S* eltp) const;
00261     void colSum(int j, S* eltp) const;
00262     void rowSum(int i, S* eltp) const;
00263         
00264     // addition and subtraction (+= and -=)
00265     void addIn(const MatrixHelper&);   
00266     void addIn(const MatrixHelper<typename CNT<S>::TNeg>&);   
00267     void subIn(const MatrixHelper&); 
00268     void subIn(const MatrixHelper<typename CNT<S>::TNeg>&); 
00269 
00270     // Fill all our stored data with copies of the same supplied element.
00271     void fillWith(const S* eltp);
00272 
00273     // We're copying data from a C++ row-oriented matrix into our general
00274     // Matrix. In addition to the row ordering, C++ may use different spacing
00275     // for elements than Simmatrix does.
00276     void copyInByRowsFromCpp(const S* elts);
00277     
00278         // Scalar operations //
00279 
00280     // Fill every element with repeated copies of a single scalar value.
00281     void fillWithScalar(const StdNumber&);
00282             
00283     // Scalar multiply (and divide). This is useful no matter what the
00284     // element structure and will produce the correct result.
00285     void scaleBy(const StdNumber&);
00286 
00287     // This is only allowed for a matrix of real or complex or neg of those,
00288     // which is square, well-conditioned, and for which we have no view,
00289     // and element size 1.
00290     void invertInPlace();
00291 
00292     void dump(const char* msg=0) const; // For debugging -- comes out in a way you can feed to Matlab
00293 
00294     // See comment in MatrixBase::matmul for an explanation.
00295     template <class SA, class SB>
00296     void matmul(const StdNumber& beta,   // applied to 'this'
00297                 const StdNumber& alpha, const MatrixHelper<SA>& A, const MatrixHelper<SB>& B);    
00298     
00299         // Bookkeeping //
00300     
00301     // This is the number of logical *elements* in each column of this matrix; i.e., m.
00302     int    nrow() const;
00303     // This is the number of *elements* in each row of this matrix; i.e., n.
00304     int    ncol() const;
00305     // This is the total number of *elements* in the logical shape of this matrix, i.e. m*n.
00306     ptrdiff_t nelt() const; // nrow*ncol  
00307     // This is the number of elements if this is a 1d matrix (vector or rowvector). That is,
00308     // it is the same as one of nrow() or ncol(); the other must be 1. It is also the
00309     // same as nelt() but limited to fit in 32 bits.
00310     int length() const;
00311 
00312     // Change the logical size of the underlying memory for this Matrix to m X n, forgetting
00313     // everything that used to be there. This will fail if it would have to resize any
00314     // non-resizable dimension. However, it will succeed even on non-resizable matrices and
00315     // views provided the existing dimensions are already correct. If no size change is made,
00316     // no action is taken and the original data is still accessible.
00317     void resize(int m, int n);
00318 
00319     // Same as resize() except as much of the original data as will fit remains in place at
00320     // the same (i,j) location it had before. This may require copying the elements after
00321     // allocating new space. As for resize(), this is allowed for any Matrix whose dimensions
00322     // are already right, even if that Matrix is not resizable.
00323     void resizeKeep(int m, int n);
00324 
00325     void lockShape();
00326     void unlockShape();
00327 
00328     const MatrixCommitment& getCharacterCommitment() const;
00329     const MatrixCharacter&  getMatrixCharacter() const;
00330     void commitTo(const MatrixCommitment&);
00331 
00332     // Access to raw data. For now this is only allowed if there is no view
00333     // and the raw data is contiguous.
00334     bool      hasContiguousData() const;
00335     ptrdiff_t getContiguousDataLength() const;
00336     const S*  getContiguousData() const;
00337     S*        updContiguousData();
00338 
00339     void replaceContiguousData(S* newData, ptrdiff_t length, bool takeOwnership);
00340     void replaceContiguousData(const S* newData, ptrdiff_t length);
00341     void swapOwnedContiguousData(S* newData, ptrdiff_t length, S*& oldData);
00342 
00343     const MatrixHelperRep<S>& getRep() const {assert(rep); return *rep;}
00344     MatrixHelperRep<S>&       updRep()       {assert(rep); return *rep;}
00345     void setRep(MatrixHelperRep<S>* hrep)    {assert(!rep); rep = hrep;}
00346     MatrixHelperRep<S>* stealRep() 
00347     {   assert(rep); MatrixHelperRep<S>* stolen=rep; rep=0; return stolen; }
00348 
00349     void deleteRepIfOwner();
00350     void replaceRep(MatrixHelperRep<S>*);
00351 
00352     // Rep-stealing constructor. We're taking ownership of the supplied rep.
00353     // Internal use only!
00354     explicit MatrixHelper(MatrixHelperRep<S>*);
00355 
00356 private:
00357     // Pointer to the private implementation object. This is the only
00358     // allowable data member in this class.
00359     class MatrixHelperRep<S>* rep;
00360 
00361     // Suppress copy constructor.
00362     MatrixHelper(const MatrixHelper&);
00363 
00364 friend class MatrixHelper<typename CNT<S>::TNeg>;
00365 friend class MatrixHelper<typename CNT<S>::THerm>;
00366 };
00367      
00368 } //namespace SimTK
00369 
00370 #endif // SimTK_SIMMATRIX_MATRIX_HELPER_H_

Generated by  doxygen 1.6.2