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_