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-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 
00061 #include "SimTKcommon/Scalar.h"
00062 
00063 #include <iostream>
00064 #include <cassert>
00065 #include <complex>
00066 #include <cstddef>
00067 
00068 namespace SimTK {
00069 
00070 
00071 template <class S> class MatrixHelperRep;
00072 
00077 class MatrixStructures {
00078     public:
00079     enum Structure {
00080         Uncommitted        = 0,
00081         Full               = 1,  // all elements are significant
00082         Diagonal           = 2,
00083         Symmetric          = 3,  // also means Hermitian
00084         Triangular         = 4,
00085         QuasiTriangular    = 5,
00086         Hessenberg         = 6,
00087         Permutation        = 7,
00088         TriDiagonal        = 8
00089     };
00090 };
00091 
00092 class MatrixShapes {
00093     public:
00094     enum Shape {
00095         Uncommitted        = 0,
00096         General            = 1,  // 2d rectangular matrix
00097         Square             = 2,  // restriction to nrow==ncol
00098         Vector             = 3,  // 1d column vector
00099         RowVector          = 4   // 1d row vector
00100     };
00101 };
00102 
00103 class MatrixSparseFormats {
00104     public:
00105     enum Sparsity {
00106         Uncommitted        = 0,
00107         Full               = 1,
00108         Banded             = 2 // needs upper & lower bandwidth
00109     };
00110 };
00111 class MatrixStorageFormats {
00112     public:
00113     enum Storage {
00114         Uncommitted        = 0,
00115         Conventional       = 1,
00116         Packed             = 2, // triangular or symmetric
00117         HouseholderProduct = 3, // orthogonal only
00118         PivotArray         = 4  // pivot only
00119     };
00120 };
00121 class MatrixConditions {
00122     public:
00123     enum Condition {
00124         Uncommitted        = 0,
00125         Unknown            = 1,
00126         Orthogonal         = 2,
00127         PositiveDefinite   = 3,
00128         WellConditioned    = 4, // implies full rank
00129         FullRank           = 5, // but might have bad conditioning
00130         Singular           = 6  // implies possible bad conditioning
00131     };
00132 };
00133 
00134 template <class S> 
00135 class SimTK_SimTKCOMMON_EXPORT MatrixHelper {
00136     typedef typename CNT<S>::Number     Number;     // strips off negator from S
00137     typedef typename CNT<S>::StdNumber  StdNumber;  // turns conjugate to complex
00138     typedef typename CNT<S>::Precision  Precision;  // strips off complex from StdNumber
00139 public:     
00140     // no default constructor, copy constructor suppressed
00141 
00142     // Destructor eliminates MatrixHelperRep object if "this" owns it.
00143     ~MatrixHelper();
00144 
00145     // Local types for directing us to the right constructor at compile time.
00146     class ShallowCopy   { };
00147     class DeepCopy      { };
00148     class TransposeView { };
00149     class DiagonalView  { };
00150     
00151     // Matrix owner constructors
00152     
00153     // 0x0, fully resizable
00154     explicit MatrixHelper(int esz, int cppEsz);
00155 
00156     // Restore helper to its just-constructed state. We forget everything except
00157     // the element size, which can never change.
00158     void clear();
00159     
00160     // (m*esz) x n, resizable with optional restrictions 
00161     MatrixHelper(int esz, int cppEsz, int m, int n, bool lockNrow=false, bool lockNcol=false);
00162     
00163     // This is a non-resizeable owner (of the data descriptor) but shares pre-existing raw data.
00164     MatrixHelper(int esz, int cppEsz, int m, int n, int leadingDim, const S*); // read only
00165     MatrixHelper(int esz, int cppEsz, int m, int n, int leadingDim, S*);       // writable
00166    
00167     // Matrix view constructors, read only and writable. Use these for submatrices,
00168     // rows, and columns. Indices are by *element* and these may consist of multiple
00169     // scalars of type template parameter S.
00170     MatrixHelper(const MatrixHelper&, int i, int j, int nrow, int ncol);
00171     MatrixHelper(MatrixHelper&,       int i, int j, int nrow, int ncol); 
00172 
00173     MatrixHelper(const MatrixHelper<typename CNT<S>::THerm>&, 
00174                  const TransposeView&);    // a read only transposed view
00175     MatrixHelper(MatrixHelper<typename CNT<S>::THerm>&,       
00176                  const TransposeView&);    // a writable transposed view
00177 
00178     // These are the constructors for making a diagonal view of a Matrix.
00179     MatrixHelper(const MatrixHelper&, const DiagonalView&); // a read only diagonal view
00180     MatrixHelper(MatrixHelper&, const DiagonalView&);       // a writable diagonal view
00181 
00182     // Copy an existing MatrixHelper using deep or shallow copy as requested. 
00183     // For shallow copy, const form loses writability, non-const retains same
00184     // writability status as source.    
00185     MatrixHelper(const MatrixHelper&, const ShallowCopy&);
00186     MatrixHelper(MatrixHelper&,       const ShallowCopy&);
00187     MatrixHelper(const MatrixHelper&, const DeepCopy&);
00188 
00189     // Construct from negated form, performing the required multiplications.
00190     MatrixHelper(const MatrixHelper<typename CNT<S>::TNeg>&, const DeepCopy&);
00191    
00192     // Behavior of copy assignment depends on whether this is an owner or view. If
00193     // it's an owner it is resized and ends up a new, dense copy of rhs. If
00194     // it's a view, sizes must match and we copy elements from rhs to corresponding
00195     // elements of lhs.  
00196     MatrixHelper& copyAssign(const MatrixHelper&);
00197     MatrixHelper& negatedCopyAssign(const MatrixHelper<typename CNT<S>::TNeg>&);
00198 
00199     // View assign always disconnects this helper from whatever view & data
00200     // it used to have (destructing as appropriate) and then makes it a new view
00201     // of the source. Writability is lost if the source is const, otherwise 
00202     // writability is inherited from the source.
00203     MatrixHelper& readOnlyViewAssign(const MatrixHelper&);
00204     MatrixHelper& writableViewAssign(MatrixHelper&);
00205 
00206     // Using *element* indices, obtain a pointer to the beginning of a particular
00207     // element. This is always a slow operation compared to raw array access;
00208     // use sparingly.    
00209     const S* getElt(int i, int j) const;
00210     S*       updElt(int i, int j);
00211 
00212     
00213     // Add up all the *elements*, returning the answer in the element given
00214     // by pointer to its first scalar.
00215     void sum(S* eltp) const;
00216     void colSum(int j, S* eltp) const;
00217     void rowSum(int i, S* eltp) const;
00218         
00219     // addition and subtraction (+= and -=)
00220     void addIn(const MatrixHelper&);   
00221     void addIn(const MatrixHelper<typename CNT<S>::TNeg>&);   
00222     void subIn(const MatrixHelper&); 
00223     void subIn(const MatrixHelper<typename CNT<S>::TNeg>&); 
00224     
00225     void fillWith(const S* eltp);
00226     void copyInByRows(const S* elts);
00227     
00228         // Scalar operations //
00229 
00230     // Fill every element with repeated copies of a single scalar value.
00231     void fillWithScalar(const StdNumber&);
00232             
00233     // Scalar multiply (and divide). This is useful no matter what the
00234     // element structure and will produce the correct result.
00235     void scaleBy(const StdNumber&);
00236     
00237     // Sums of scalars, regardless of element structure. Not much use except
00238     // when the elements are themselves scalars.
00239     S scalarSum() const;
00240     S scalarColSum(int j) const;
00241     S scalarRowSum(int i) const; 
00242 
00243     // This is only allowed for a matrix of real or complex or neg of those,
00244     // which is square, well-conditioned, and for which we have no view,
00245     // and element size 1.
00246     void invertInPlace();
00247 
00248     void dump(const char* msg=0) const; // For debugging -- comes out in a way you can feed to Matlab
00249 
00250     // See comment in MatrixBase::matmul for an explanation.
00251     template <class SA, class SB>
00252     void matmul(const StdNumber& beta,   // applied to 'this'
00253                 const StdNumber& alpha, const MatrixHelper<SA>& A, const MatrixHelper<SB>& B);
00254     
00255     // Here we use element indices but return the first scalar we find there.
00256     // This is useless except when the elements are scalars.
00257     const S& getScalar(int i, int j) const;
00258     S&       updScalar(int i, int j);      
00259     
00260         // Bookkeeping //
00261         
00262     int  nrow() const;
00263     int  ncol() const;
00264     long size() const; // nrow*ncol  
00265 
00266     void resize(int m, int n);
00267     void resizeKeep(int m, int n);
00268 
00269     void lockNRows();
00270     void lockNCols();
00271     void lockShape();
00272     void unlockNRows();
00273     void unlockNCols();
00274     void unlockShape();
00275 
00276     void setMatrixStructure(MatrixStructures::Structure );
00277     MatrixStructures::Structure getMatrixStructure() const;
00278     void setMatrixShape(MatrixShapes::Shape );
00279     MatrixShapes::Shape getMatrixShape() const;
00280     void setMatrixSparsity(MatrixSparseFormats::Sparsity );
00281     MatrixSparseFormats::Sparsity getMatrixSparsity() const;
00282     void setMatrixStorage(MatrixStorageFormats::Storage );
00283     MatrixStorageFormats::Storage getMatrixStorage() const;
00284     void setMatrixCondition(MatrixConditions::Condition );
00285     MatrixConditions::Condition getMatrixCondition() const;
00286 
00287     // Access to raw data. For now this is only allowed if there is no view
00288     // and the raw data is contiguous.
00289     bool hasContiguousData() const;
00290     long getContiguousDataLength() const;
00291     const S* getContiguousData() const;
00292     S* updContiguousData();
00293     void replaceContiguousData(S* newData, long length, bool takeOwnership);
00294     void replaceContiguousData(const S* newData, long length);
00295     void swapOwnedContiguousData(S* newData, int length, S*& oldData);
00296 
00297     // Pointer to the private implementation object.
00298     class MatrixHelperRep<S>* rep;
00299 
00300 private:
00301     // Suppress copy constructor.
00302     MatrixHelper(const MatrixHelper&);
00303 };
00304      
00305 } //namespace SimTK
00306 
00307 #endif // SimTK_SIMMATRIX_MATRIX_HELPER_H_

Generated on Thu Feb 28 01:34:30 2008 for SimTKcommon by  doxygen 1.4.7