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
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
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
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