Simbody 3.7
Loading...
Searching...
No Matches
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
45namespace SimTK {
46
47
48template <class S> class MatrixHelperRep;
49class MatrixCharacter;
51
52// ------------------------------ MatrixHelper --------------------------------
53
76
77// ----------------------------------------------------------------------------
78template <class S>
80 typedef MatrixHelper<S> This;
83public:
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).
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.)
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.
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.
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.
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.
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.
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();
317
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;
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
341
342 // Rep-stealing constructor. We're taking ownership of the supplied rep.
343 // Internal use only!
345
346private:
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.
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
360friend class MatrixHelper<typename CNT<S>::TNeg>;
361friend class MatrixHelper<typename CNT<S>::THerm>;
362};
363
364} //namespace SimTK
365
366#endif // SimTK_SIMMATRIX_MATRIX_HELPER_H_
This is a user-includable header which includes everything needed to make use of SimMatrix Scalar cod...
#define SimTK_SimTKCOMMON_EXPORT
Definition SimTKcommon/include/SimTKcommon/internal/common.h:224
The Array_<T> container class is a plug-compatible replacement for the C++ standard template library ...
Definition Array.h:1520
size_type size() const
Return the current number of elements stored in this array.
Definition Array.h:2075
const T * cbegin() const
Return a const pointer to the first element of this array if any, otherwise cend(),...
Definition Array.h:2212
Specialized information about Composite Numerical Types which allows us to define appropriate templat...
Definition CompositeNumericalTypes.h:136
K::StdNumber StdNumber
Definition CompositeNumericalTypes.h:163
K::THerm THerm
Definition CompositeNumericalTypes.h:144
K::TNeg TNeg
Definition CompositeNumericalTypes.h:139
K::Precision Precision
Definition CompositeNumericalTypes.h:164
K::Number Number
Definition CompositeNumericalTypes.h:162
A MatrixCharacter is a set containing a value for each of the matrix characteristics except element t...
Definition MatrixCharacteristics.h:597
A MatrixCommitment provides a set of acceptable matrix characteristics.
Definition MatrixCharacteristics.h:832
Definition MatrixHelper.h:48
Definition MatrixHelper.h:96
Definition MatrixHelper.h:98
Definition MatrixHelper.h:95
Definition MatrixHelper.h:97
Here we define class MatrixHelper, the scalar-type templatized helper class for the more general,...
Definition MatrixHelper.h:79
MatrixHelper(const MatrixCommitment &, MatrixHelper &, int i, int j, int nrow, int ncol)
MatrixHelper(const MatrixCommitment &mc, MatrixHelper &h, const Array_< int > &indices)
Definition MatrixHelper.h:183
void rowSum(int i, S *eltp) const
void dump(const char *msg=0) const
MatrixHelper & negatedCopyAssign(const MatrixHelper< typename CNT< S >::TNeg > &)
const MatrixCommitment & getCharacterCommitment() const
const S * getElt(int i, int j) const
CNT< S >::Number Number
Definition MatrixHelper.h:84
CNT< S >::StdNumber StdNumber
Definition MatrixHelper.h:85
MatrixHelper & writableViewAssign(MatrixHelper &source)
const MatrixHelperRep< S > & getRep() const
Definition MatrixHelper.h:333
MatrixHelper(const MatrixCommitment &, const MatrixHelper &, int i, int j, int nrow, int ncol)
const S * getContiguousData() const
MatrixHelper(int esz, int cppEsz, const MatrixCommitment &)
void getAnyElt(int i, S *value) const
const MatrixCharacter & getMatrixCharacter() const
bool hasContiguousData() const
MatrixHelper(const MatrixCommitment &, const MatrixHelper< typename CNT< S >::THerm > &, const TransposeView &)
void fillWith(const S *eltp)
void addIn(const MatrixHelper< typename CNT< S >::TNeg > &)
void resizeKeep(int m, int n)
void replaceRep(MatrixHelperRep< S > *)
void subIn(const MatrixHelper &)
MatrixHelper(const MatrixCommitment &, const MatrixHelper &, int n, const int *indices)
bool isClear() const
S * updElt(int i, int j)
MatrixHelper(const MatrixCommitment &, MatrixHelper &, const DiagonalView &)
void setRep(MatrixHelperRep< S > *hrep)
Definition MatrixHelper.h:335
void sum(S *eltp) const
MatrixHelper(const MatrixCommitment &, const MatrixHelper &source, const ShallowCopy &)
MatrixHelper(const MatrixCommitment &, const MatrixHelper< typename CNT< S >::TNeg > &source, const DeepCopy &)
MatrixHelper(int esz, int cppEsz, const MatrixCommitment &, int m, int n)
void resize(int m, int n)
MatrixHelper(const MatrixCommitment &, MatrixHelper &source, const ShallowCopy &)
MatrixHelper(MatrixHelperRep< S > *)
void commitTo(const MatrixCommitment &)
CNT< S >::Precision Precision
Definition MatrixHelper.h:86
MatrixHelper & copyAssign(const MatrixHelper &source)
MatrixHelper(const MatrixCommitment &mc, const MatrixHelper &h, const Array_< int > &indices)
Definition MatrixHelper.h:180
MatrixHelper(const MatrixCommitment &, const MatrixHelper &, const DiagonalView &)
void swapOwnedContiguousData(S *newData, ptrdiff_t length, S *&oldData)
ptrdiff_t nelt() const
MatrixHelperRep< S > * stealRep()
Definition MatrixHelper.h:336
MatrixHelper(int esz, int cppEsz, const MatrixCommitment &, const MatrixCharacter &, int spacing, const S *data)
MatrixHelper(int esz, int cppEsz)
MatrixHelper(const MatrixCommitment &, MatrixHelper< typename CNT< S >::THerm > &, const TransposeView &)
ptrdiff_t getContiguousDataLength() const
void replaceContiguousData(const S *newData, ptrdiff_t length)
void getAnyElt(int i, int j, S *value) const
~MatrixHelper()
Definition MatrixHelper.h:92
void scaleBy(const StdNumber &)
void copyInByRowsFromCpp(const S *elts)
MatrixHelper(int esz, int cppEsz, const MatrixCommitment &, const MatrixCharacter &, int spacing, S *data)
void colSum(int j, S *eltp) const
MatrixHelper(const MatrixCommitment &, MatrixHelper &, int n, const int *indices)
void subIn(const MatrixHelper< typename CNT< S >::TNeg > &)
void addIn(const MatrixHelper &)
MatrixHelperRep< S > & updRep()
Definition MatrixHelper.h:334
MatrixHelper & readOnlyViewAssign(const MatrixHelper &source)
int length() const
const S * getElt(int i) const
void fillWithScalar(const StdNumber &)
void replaceContiguousData(S *newData, ptrdiff_t length, bool takeOwnership)
MatrixHelper(const MatrixCommitment &, const MatrixHelper &source, const DeepCopy &)
This is the top-level SimTK namespace into which all SimTK names are placed to avoid collision with o...
Definition Assembler.h:37