Simbody  3.5
MatrixHelper.h
Go to the documentation of this file.
1 #ifndef SimTK_SIMMATRIX_MATRIX_HELPER_H_
2 #define SimTK_SIMMATRIX_MATRIX_HELPER_H_
3 
4 /* -------------------------------------------------------------------------- *
5  * Simbody(tm): SimTKcommon *
6  * -------------------------------------------------------------------------- *
7  * This is part of the SimTK biosimulation toolkit originating from *
8  * Simbios, the NIH National Center for Physics-Based Simulation of *
9  * Biological Structures at Stanford, funded under the NIH Roadmap for *
10  * Medical Research, grant U54 GM072970. See https://simtk.org/home/simbody. *
11  * *
12  * Portions copyright (c) 2005-12 Stanford University and the Authors. *
13  * Authors: Michael Sherman *
14  * Contributors: *
15  * *
16  * Licensed under the Apache License, Version 2.0 (the "License"); you may *
17  * not use this file except in compliance with the License. You may obtain a *
18  * copy of the License at http://www.apache.org/licenses/LICENSE-2.0. *
19  * *
20  * Unless required by applicable law or agreed to in writing, software *
21  * distributed under the License is distributed on an "AS IS" BASIS, *
22  * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. *
23  * See the License for the specific language governing permissions and *
24  * limitations under the License. *
25  * -------------------------------------------------------------------------- */
26 
37 #include "SimTKcommon/Scalar.h"
38 
39 #include <iostream>
40 #include <cassert>
41 #include <complex>
42 #include <cstddef>
43 #include <utility> // for std::pair
44 
45 namespace SimTK {
46 
47 
48 template <class S> class MatrixHelperRep;
49 class MatrixCharacter;
50 class MatrixCommitment;
51 
52 // ------------------------------ MatrixHelper --------------------------------
53 
76 
77 // ----------------------------------------------------------------------------
78 template <class S>
80  typedef MatrixHelper<S> This;
83 public:
84  typedef typename CNT<S>::Number Number; // strips off negator from S
85  typedef typename CNT<S>::StdNumber StdNumber; // turns conjugate to complex
86  typedef typename CNT<S>::Precision Precision; // strips off complex from StdNumber
87 
88  // no default constructor
89  // copy constructor suppressed
90 
91  // Destructor eliminates MatrixHelperRep object if "this" owns it.
92  ~MatrixHelper() {deleteRepIfOwner();}
93 
94  // Local types for directing us to the right constructor at compile time.
95  class ShallowCopy { };
96  class DeepCopy { };
97  class TransposeView { };
98  class DiagonalView { };
99 
101  // "Owner" constructors //
103 
104  // 0x0, fully resizable, fully uncommitted.
105  MatrixHelper(int esz, int cppEsz);
106 
107  // Default allocation for the given commitment.
108  MatrixHelper(int esz, int cppEsz, const MatrixCommitment&);
109 
110  // Allocate a matrix that satisfies a given commitment and has a
111  // particular initial size.
112  MatrixHelper(int esz, int cppEsz, const MatrixCommitment&, int m, int n);
113 
114  // Copy constructor that produces a new owner whose logical shape and contents are
115  // the same as the source, but with a possibly better storage layout. Data will
116  // be contiguous in the copy regardless of how spread out it was in the source.
117  // The second argument is just to disambiguate this constructor from similar ones.
118  MatrixHelper(const MatrixCommitment&, const MatrixHelper& source, const DeepCopy&);
119 
120  // This has the same semantics as the previous DeepCopy constructor, except that
121  // the source has negated elements with respect to S. The resulting logical shape
122  // and logical values are identical to the source, meaning that the negation must
123  // actually be performed here, using one flop for every meaningful source scalar.
124  MatrixHelper(const MatrixCommitment&, const MatrixHelper<typename CNT<S>::TNeg>& source, const DeepCopy&);
125 
126 
128  // External "View" constructors //
130 
131  // These constructors produce a matrix which provides a view of externally-allocated
132  // data, which is known to have the given MatrixCharacter. There is also provision
133  // for a "spacing" parameter which defines gaps in the supplied data, although
134  // the interpretation of that parameter depends on the MatrixCharacter (typically
135  // it is the leading dimension for a matrix or the stride for a vector). Note
136  // that spacing is interpreted as "number of scalars between adjacent elements"
137  // which has the usual Lapack interpretation if the elements are scalars but
138  // can be used for C++ vs. Simmatrix packing differences for composite elements.
139  // The resulting handle is *not* the owner of the data, so destruction of the
140  // handle does not delete the data.
141 
142  // Create a read-only view into existing data.
143  MatrixHelper(int esz, int cppEsz, const MatrixCommitment&,
144  const MatrixCharacter&, int spacing, const S* data);
145  // Create a writable view into existing data.
146  MatrixHelper(int esz, int cppEsz, const MatrixCommitment&,
147  const MatrixCharacter&, int spacing, S* data);
148 
149 
151  // Matrix "View" constructors //
153 
154  // Matrix view constructors, read only and writable. Use these for submatrices,
155  // rows, and columns. Indices are by *element* and these may consist of multiple
156  // scalars of type template parameter S.
157 
158  // "Block" views
159  MatrixHelper(const MatrixCommitment&, const MatrixHelper&, int i, int j, int nrow, int ncol);
160  MatrixHelper(const MatrixCommitment&, MatrixHelper&, int i, int j, int nrow, int ncol);
161 
162  // "Transpose" views (note that this is Hermitian transpose; i.e., element
163  // type is conjugated).
164  MatrixHelper(const MatrixCommitment&, const MatrixHelper<typename CNT<S>::THerm>&,
165  const TransposeView&); // a read only transposed view
167  const TransposeView&); // a writable transposed view
168 
169  // "Diagonal" views.
170  MatrixHelper(const MatrixCommitment&, const MatrixHelper&, const DiagonalView&); // a read only diagonal view
171  MatrixHelper(const MatrixCommitment&, MatrixHelper&, const DiagonalView&); // a writable diagonal view
172 
173  // "Indexed" view of a 1d matrix.
175  int n, const int* indices);
177  int n, const int* indices);
178 
179  // These invoke the previous constructors but with friendlier index source.
181  const Array_<int>& indices)
182  { new (this) MatrixHelper(mc, h, (int)indices.size(), indices.cbegin()); }
184  const Array_<int>& indices)
185  { new (this) MatrixHelper(mc, h, (int)indices.size(), indices.cbegin()); }
186 
187  // "Copy" an existing MatrixHelper by making a new view into the same data.
188  // The const form loses writability, non-const retains same writability status
189  // as source. If the source is already a view then the destination will have
190  // an identical element filter so that the logical shape and values are the
191  // same for both source and copy. (The second argument is used just to
192  // disambiguate this constructor from similar ones.)
193  MatrixHelper(const MatrixCommitment&, const MatrixHelper& source, const ShallowCopy&);
194  MatrixHelper(const MatrixCommitment&, MatrixHelper& source, const ShallowCopy&);
195 
197  // Copy assignment //
199 
200  // Behavior of copy assignment depends on whether "this" is an owner or view. If
201  // it's an owner it is resized and ends up a new, dense copy of "source" just as
202  // for the DeepCopy constructor above. If "this" is a writable view, sizes must match
203  // and we copy elements from "source" to corresponding elements of "this". If
204  // "this" is not writable then the operation will fail.
205  MatrixHelper& copyAssign(const MatrixHelper& source);
206 
207  // Same as copyAssign() except the source has element types which are logically
208  // negated from the element types of "this". Since the result must have the
209  // same value as the source, this requires that all the copied elements are
210  // actually negated at a cost of one flop per scalar.
211  MatrixHelper& negatedCopyAssign(const MatrixHelper<typename CNT<S>::TNeg>&);
212 
214  // View assignment //
216 
217  // View assign always disconnects this helper from whatever view & data
218  // it used to have (destructing as appropriate) and then makes it a new view
219  // of the source. Writability is lost if the source is const, otherwise
220  // writability is inherited from the source.
221  MatrixHelper& readOnlyViewAssign(const MatrixHelper& source);
222  MatrixHelper& writableViewAssign(MatrixHelper& source);
223 
224  // Restore helper to its just-constructed state. We forget everything except
225  // the element size and handle commitment, which can never change. The
226  // allocated helper will be the same as if we had just default-constructed
227  // using the current commitment.
228  void clear();
229 
230  // Return true if there is currently no data associated with this handle.
231  bool isClear() const;
232 
233  // Using *element* indices, obtain a pointer to the beginning of a particular
234  // element. This is always a slow operation compared to raw array access;
235  // use sparingly. These will fail if the indices are outside the stored
236  // portion of the matrix. Use getAnyElt() if you want something that always
237  // works.
238  const S* getElt(int i, int j) const;
239  S* updElt(int i, int j);
240 
241  // Faster for 1-d matrices (vectors) if you know you have one.
242  const S* getElt(int i) const;
243  S* updElt(int i);
244 
245  // This returns the correct value for any element within the logical dimensions
246  // of the matrix. In some cases it has to compute the value; in all cases
247  // it has to copy it.
248  void getAnyElt(int i, int j, S* value) const;
249 
250  // Faster for 1-d matrices (vectors) if you know you have one.
251  void getAnyElt(int i, S* value) const;
252 
253  // Add up all the *elements*, returning the answer in the element given
254  // by pointer to its first scalar.
255  void sum(S* eltp) const;
256  void colSum(int j, S* eltp) const;
257  void rowSum(int i, S* eltp) const;
258 
259  // addition and subtraction (+= and -=)
260  void addIn(const MatrixHelper&);
261  void addIn(const MatrixHelper<typename CNT<S>::TNeg>&);
262  void subIn(const MatrixHelper&);
263  void subIn(const MatrixHelper<typename CNT<S>::TNeg>&);
264 
265  // Fill all our stored data with copies of the same supplied element.
266  void fillWith(const S* eltp);
267 
268  // We're copying data from a C++ row-oriented matrix into our general
269  // Matrix. In addition to the row ordering, C++ may use different spacing
270  // for elements than Simmatrix does.
271  void copyInByRowsFromCpp(const S* elts);
272 
273  // Scalar operations //
274 
275  // Fill every element with repeated copies of a single scalar value.
276  void fillWithScalar(const StdNumber&);
277 
278  // Scalar multiply (and divide). This is useful no matter what the
279  // element structure and will produce the correct result.
280  void scaleBy(const StdNumber&);
281 
282  // This is only allowed for a matrix of real or complex or neg of those,
283  // which is square, well-conditioned, and for which we have no view,
284  // and element size 1.
285  void invertInPlace();
286 
287  void dump(const char* msg=0) const; // For debugging -- comes out in a way you can feed to Matlab
288 
289  // Bookkeeping //
290 
291  // This is the number of logical *elements* in each column of this matrix; i.e., m.
292  int nrow() const;
293  // This is the number of *elements* in each row of this matrix; i.e., n.
294  int ncol() const;
295  // This is the total number of *elements* in the logical shape of this matrix, i.e. m*n.
296  ptrdiff_t nelt() const; // nrow*ncol
297  // This is the number of elements if this is a 1d matrix (vector or rowvector). That is,
298  // it is the same as one of nrow() or ncol(); the other must be 1. It is also the
299  // same as nelt() but limited to fit in 32 bits.
300  int length() const;
301 
302  // Change the logical size of the underlying memory for this Matrix to m X n, forgetting
303  // everything that used to be there. This will fail if it would have to resize any
304  // non-resizable dimension. However, it will succeed even on non-resizable matrices and
305  // views provided the existing dimensions are already correct. If no size change is made,
306  // no action is taken and the original data is still accessible.
307  void resize(int m, int n);
308 
309  // Same as resize() except as much of the original data as will fit remains in place at
310  // the same (i,j) location it had before. This may require copying the elements after
311  // allocating new space. As for resize(), this is allowed for any Matrix whose dimensions
312  // are already right, even if that Matrix is not resizable.
313  void resizeKeep(int m, int n);
314 
315  void lockShape();
316  void unlockShape();
317 
318  const MatrixCommitment& getCharacterCommitment() const;
319  const MatrixCharacter& getMatrixCharacter() const;
320  void commitTo(const MatrixCommitment&);
321 
322  // Access to raw data. For now this is only allowed if there is no view
323  // and the raw data is contiguous.
324  bool hasContiguousData() const;
325  ptrdiff_t getContiguousDataLength() const;
326  const S* getContiguousData() const;
327  S* updContiguousData();
328 
329  void replaceContiguousData(S* newData, ptrdiff_t length, bool takeOwnership);
330  void replaceContiguousData(const S* newData, ptrdiff_t length);
331  void swapOwnedContiguousData(S* newData, ptrdiff_t length, S*& oldData);
332 
333  const MatrixHelperRep<S>& getRep() const {assert(rep); return *rep;}
334  MatrixHelperRep<S>& updRep() {assert(rep); return *rep;}
335  void setRep(MatrixHelperRep<S>* hrep) {assert(!rep); rep = hrep;}
337  { assert(rep); MatrixHelperRep<S>* stolen=rep; rep=0; return stolen; }
338 
339  void deleteRepIfOwner();
340  void replaceRep(MatrixHelperRep<S>*);
341 
342  // Rep-stealing constructor. We're taking ownership of the supplied rep.
343  // Internal use only!
344  explicit MatrixHelper(MatrixHelperRep<S>*);
345 
346 private:
347  // Pointer to the private implementation object. This is the only
348  // allowable data member in this class.
349  class MatrixHelperRep<S>* rep;
350 
351  // Suppress copy constructor.
352  MatrixHelper(const MatrixHelper&);
353 
354  // ============================= Unimplemented =============================
355  // See comment in MatrixBase::matmul for an explanation.
356  template <class SA, class SB>
357  void matmul(const StdNumber& beta, // applied to 'this'
358  const StdNumber& alpha, const MatrixHelper<SA>& A, const MatrixHelper<SB>& B);
359 
360 friend class MatrixHelper<typename CNT<S>::TNeg>;
361 friend class MatrixHelper<typename CNT<S>::THerm>;
362 };
363 
364 } //namespace SimTK
365 
366 #endif // SimTK_SIMMATRIX_MATRIX_HELPER_H_
Here we define class MatrixHelper<S>, the scalar-type templatized helper class for the more general...
Definition: MatrixHelper.h:79
MatrixHelper(const MatrixCommitment &mc, const MatrixHelper &h, const Array_< int > &indices)
Definition: MatrixHelper.h:180
#define SimTK_SimTKCOMMON_EXPORT
Definition: SimTKcommon/include/SimTKcommon/internal/common.h:218
Definition: MatrixHelper.h:48
Apache License January AND DISTRIBUTION Definitions License shall mean the terms and conditions for and distribution as defined by Sections through of this document Licensor shall mean the copyright owner or entity authorized by the copyright owner that is granting the License Legal Entity shall mean the union of the acting entity and all other entities that control are controlled by or are under common control with that entity For the purposes of this definition control direct or to cause the direction or management of such whether by contract or including but not limited to software source documentation source
Definition: LICENSE.txt:26
This is the top-level SimTK namespace into which all SimTK names are placed to avoid collision with o...
Definition: Assembler.h:37
A MatrixCharacter is a set containing a value for each of the matrix characteristics except element t...
Definition: MatrixCharacteristics.h:597
Definition: MatrixHelper.h:95
size_type size() const
Return the current number of elements stored in this array.
Definition: Array.h:2037
const T * cbegin() const
Return a const pointer to the first element of this array if any, otherwise cend(), which may be null (0) in that case but does not have to be.
Definition: Array.h:2172
m
Definition: CMakeCache.txt:469
ELEM sum(const VectorBase< ELEM > &v)
Definition: VectorMath.h:147
MatrixHelper(const MatrixCommitment &mc, MatrixHelper &h, const Array_< int > &indices)
Definition: MatrixHelper.h:183
K::Precision Precision
Definition: CompositeNumericalTypes.h:164
CNT< S >::Precision Precision
Definition: MatrixHelper.h:86
CNT< S >::StdNumber StdNumber
Definition: MatrixHelper.h:85
This is a user-includable header which includes everything needed to make use of SimMatrix Scalar cod...
MatrixHelperRep< S > * stealRep()
Definition: MatrixHelper.h:336
CNT< S >::Number Number
Definition: MatrixHelper.h:84
~MatrixHelper()
Definition: MatrixHelper.h:92
K::StdNumber StdNumber
Definition: CompositeNumericalTypes.h:163
MatrixHelperRep< S > & updRep()
Definition: MatrixHelper.h:334
gikDdMV wfaIJt A٩t1 JcA nr S q is3 ֧ VK C 9Z D q Fxn n T Y < ['jd< K JvTMH"sw>}o_o? z'z:mV$yng͖i۸J{ Ta*dE|lzbX@!^Ooi_=O}&ŲQUVWTsh!P_7DRAVfʿbOԹɫt0Y!|'x'óݥ:/ V[,}-B֞/܂;:;;Iޘ[nK4#-='Gf\lb41۩> Os7x f pZzB I g n
Definition: SimmathUserGuide.doc:2262
const MatrixHelperRep< S > & getRep() const
Definition: MatrixHelper.h:333
Specialized information about Composite Numerical Types which allows us to define appropriate templat...
Definition: CompositeNumericalTypes.h:136
Definition: MatrixHelper.h:98
K::TNeg TNeg
Definition: CompositeNumericalTypes.h:139
void setRep(MatrixHelperRep< S > *hrep)
Definition: MatrixHelper.h:335
A MatrixCommitment provides a set of acceptable matrix characteristics.
Definition: MatrixCharacteristics.h:832
Definition: MatrixHelper.h:97
K::Number Number
Definition: CompositeNumericalTypes.h:162
K::THerm THerm
Definition: CompositeNumericalTypes.h:144
Definition: MatrixHelper.h:96