I calls 'em as I sees 'em. — proverb often attributed to baseball umpires.
Thanks to David Serrano for review.
|
This is a container for use with the C++ programming language:
You may download the source code (~100kb) and a test program (~40kb) and modify them to suit your needs. Sections 6, 7, and 8 of this report contain links to many short but complete sample programs. |
§ 1A. This report describes a container to be used with the C++ programming language: a multi-dimensional array, its distinctive feature being that the number of dimensions need not be known at compile time. In this regard, it resembles the arrays of APL. Also, subscripts need not be zero-based or one-based, but can vary over whatever integer range the programmer specifies; this is similar to what is found in Ada.
Its namespace is called mat_gen_dim, which abbreviates matrix of generalized dimensionality; all pertinent definitions lie within this namespace (file mat_gen_dim.h). The container itself comes in two varieties, namely var_matrix (for variable) and con_matrix (for constant). In the former, components can be altered once the matrix is constructed; in the latter, they cannot. Both var_matrix and con_matrix are derived from a sub_matrix class of limited functionality.
The word matrix is used instead of array because the latter already has extensive meaning in the C++ language. The (arguably ugly) name mat_gen_dim is unlikely to clash with identifiers in any other programs.
Other endeavors in this area are few, among them Bjoern Andres' marray whose source code is readily downloadable, and Josef Weinbub et alii's approach the source code for which the present author could not find. Plentiful by contrast is software to handle the case where the dimensionality is known at compile time; for instance Bjarne Stroustrup devotes a chapter of The C++ Programming Language, fourth edition to that subject.
§ 1B. As with arrays in the C or Fortran languages, storage of a matrix is dense rather than sparse. Also in common with those languages, a matrix is rectangular in shape; here is what that means. Suppose for example there is some seven-dimensional array and all three of the following subscripts (a few numbers underlined for emphasis) are not out of range:
( +2, -6, +3, +8, -1, 0, +9)
( +2, -6, +3, +8, +7, 0, +9)
( +2, -5, +3, +8, -1, 0, +9)
then if the array is rectangular, this next subscript will also be in range:
( +2, -5, +3, +8, +7, 0, +9)
Not supported by mat_gen_dim are non-rectangular arrays, the best-known of which is triangular. Such a data structure in two dimensions might have the constraint that both subscripts are nonnegative, and the second subscript must be no greater than the first; it is capable of being efficient in both memory use and component access. Counterparts of higher dimensionality immediately follow.
The constraint of rectangularity means that mat_gen_dim cannot directly effect the "ragged-edge" array possible with nesting of the C++ Standard Library's std::vector. Still, nesting is possible with mat_gen_dim, and if a matrix has components that are themselves matrices, those inner matrices need not have the same dimensionality or subscript ranges as one another. Nonetheless, they must all have the same component type.
§ 1C. Whatever the dimensionality of a matrix, it is mapped into a one-dimensional array for internal storage; contiguity of the components maximizes the chance for cache hits. This container uses a generalization of row-major order (as in C++), as opposed to column-major (as in Fortran). Informally, "last subscript varies fastest". This influences what happens when there is a change in the matrix's number of dimensions, or in the ranges of its subscripts.
Besides the usual arithmetic, the mat_gen_dim supports the outer product, as well as generalizations of contraction and the inner product. Thus it can be used for many types of tensor calculation; note however that a matrix does not distinguish superscripts versus subscripts, nor does it limit contraction to two subscripts, so the tensor programmer will have to look after these matters — an adapter class might fit the bill.
There are no direct facilities to change the number of components in a matrix; this is because in a multi-dimensional array there is no one obvious place where old components should be deleted or new ones inserted. However, there is support to copy rectangular blocks of components from one matrix to another, and thus to build a new matrix of any size from bits and pieces of other matrices.
§ 1D. This container relies heavily on the C++ Standard Library, and employs its std::vector for many kinds of internal storage. Despite great effort, however, the author could not find a satisfactory way to give mat_gen_dim an interface similar to that typically found for the Standard Library's containers, and a much different approach ultimately emerged.
Another departure from C or C++ practice is the manner of writing subscripts. If a C-style array is defined thus:
double da[7][9][8];
then subarray expressions such as
da[3]
da[4][5]
are legitimate, as demonstrated by their incorporation in the following definitions:
double (* dap3)[7][9][8] = &(da); // a 3-dimensional subarray
double (* dap2)[9][8] = &(da[3]); // a 2-dimensional subarray
double (* dap1)[8] = &(da[4][5]); // a 1-dimensional subarray
double (* dap0) = &(da[1][6][4]); // a 0-dimensional subarray
However, subarray selection is limited to those subarrays that happen to be stored in contiguous memory. A different syntax for a C-style subarray might have been:
da[3] ↔ da[3][all][all]
da[4][5] ↔ da[4][ 5 ][all]
What C-style arrays lack is an equivalent to these:
da[all][ 2 ][all]
da[all][all][ 3 ]
da[all][ 1 ][ 4 ]
da[ 6 ][all][ 7 ]
In other words, C-style arrays require the "all" to be in trailing positions. The mat_gen_dim namespace has two lenses, var_lens_less_dim and con_lens_less_dim, that lift this constraint, bringing forth a greater symmetry among subscripts. To make this work well, a much different syntax became necessary, as will be seen later.
§ 1E. The arrays found throughout computer science might be classified four ways, according to the information available at when the compiler is run, or the interpreter is started:
| subscript ranges known | subscript ranges unknown | |
|---|---|---|
| dimensionality known | C-style arrays | Ada unconstrained arrays |
| dimensionality unknown | rare, but might appear in tensor calculations | APL, mat_gen_dim |
Compared to the arrays in the other categories, the matrix of mat_gen_dim can be expected to run slower, use somewhat more memory, and require more effort on the part of the user. However, if key specifications of the data structure are simply not available at compile time, there may be little choice. This report details the implementation more than minimally, in order to show what steps have been taken to mitigate overhead.
Plausible also is the case where a complicated program manipulates numerous arrays of many different dimensionalities and subscript ranges; and although all the details could in principle be worked out at compile time, great labor would be required in doing so. Then the programmer might find it more convenient to use mat_gen_dim and let the machine do the housekeeping.
§ 1F. There is a boolean constexpr named mat_gen_dim::check (file mat_gen_dim.h). When true, the program will verify, among many other things, that subscripts have the right number of components and that those components in range. Some exceptions will be thrown by mat_gen_dim's own facilities, and others by the Standard Library, particularly std::vector. When false, the checking is skipped and the program will run faster, but an erroneous program will lead to unpredictable behavior.
The code contains many passages such as this:
if (check)
return some_vector.at(n);
else
return some_vector[n];
Since check is a compile-time constant, the compiler will presumably optimize away the test and the unexecutable code so that there is no run-time cost of time or space associated with the test.
A programmer who wants finer control of checking can within any class, function, or other scope declare a more local constexpr named check and give it a value different from mat_gen_dim::check.
§ 2A. Every matrix contains exactly two pieces of data: a form and a body:
Each of var_matrix and con_matrix is a template class, and takes one or two template parameters:
Given a matrix, the programmer can examine its form, or indeed remove the current form and put in a different one at little cost if the counts match; neither the number of dimensions nor the subscript ranges need be the same. Also, a constant reference to the body's contents can be extracted as a std::vector<compo,alloc>.
§ 2B. For efficiency, the form (file 1E_form.h) is little more than a pointer to a subform (file 1D_subform.h) which holds the actual dimensional data. The subform is not visible to the programmer, all access being through the form. Although a form can be constant or variable, a subform is always constant. Often, several forms will point to the same subform. Every matrix has exactly one form, but several matrices might share a subform. Meanwhile, there can be free-standing forms, unaffiliated with any matrix.
An internally-used reference counting scheme means that the form's copy constructor and copy assignment operator will incur little overhead and thus may be used liberally, because the subform, which is of nontrivial size, will quite rarely be copied or assigned. Nearly all copying and assigning of forms will thus be shallow rather than deep, fast rather than slow.
§ 2C. By the same token, the body (file 1G_body.h) is little more than a pointer to a subbody (file 1F_subbody.h) which holds the actual component data; every matrix has exactly one body. Neither the body nor the subbody is visible to the programmer, all access being through the matrix. In contrast to the subform which is always constant, the subbody is variable in the case of a var_matrix, but not a con_matrix; and that is the germ of distinction of the two categories of matrix types. As with the form, the body uses reference counting, with the body's copy constructor and copy assignment operator incuring little overhead.
Because many subbodys are variable, copy-on-write semantics are employed. If, as often occurs, multiple bodys are pointing at the same subbody, and one of the bodys wants to change the value of a component in the subbody, the subbody is deep-copied (cloned). The modifying body gets the new copy and makes the change, while the other bodys continue to point to the unmodified original subbody. Cloning of the subbody is an expensive operation for a large matrix, and this scheme is designed to perform the copy only when inevitable. In the contrasting case where only one body is pointing to the subbody in question, no cloning is necessary.
§ 2D. Copy construction and copy assignment are low-overhead operations for a matrix, because they are low-overhead operations for its two elements, the form and the body. Depending on the implementation, when a program passes a matrix as a parameter (or argument) to a function, it may be more efficient to pass (or call) by copy rather than to pass by constant reference. Of course, a compiler might optimize either into something better.
Some policy on parameter passing had to be adopted for the inner workings of mat_gen_dim; among several feasible options, the following was chosen:
Extensive copying of objects of the component type is sometimes unavoidable, so mat_gen_dim assumes that the programmer has supplied a component type that is reasonably efficient to copy. For the built-in C++ types, this is not a problem; on the other hand for complicated types the programmer might benefit from using smart pointers or other indirection.
As the source code for mat_gen_dim is public, the programmer can make changes to increase efficiency in particular applications.
§ 3A1. A form behaves something like an array whose component type is blt (file 1C_blt.h), which stands for bottom-length-top. Within a form is one instance of a blt for each dimension of the matrix. A form containing zero blts describes a zero-dimensional matrix, which has exactly one component.
A blt constructor takes two integer parameters:
The most negative integer (for instance, −2,147,483,648 in a typical 32-bit implementation) is not a legal subscript. The mat_gen_dim namespace gives this number a name: null_int. Perhaps surprisingly, this number comes in handy for some special purposes to be seen later.
The length is calculated as top minus bottom. A negative length is illegal, but a zero length is permitted, and means that the entire matrix will have zero components, no matter what subscript ranges the other dimensions might have. Constructor examples (note the C++11 syntax):
blt a { 4, 9}; // the legal subscripts will be 4, 5, 6, 7, 8
blt b {-4, 3}; // the legal subscripts will be -4, -3, -2, -1, 0, +1, +2
blt c { 5, 5}; // this is a legal blt, but no legal subscripts are available
blt d { 5, 3}; // this is illegal, as the length would be −2
A non-C++ notation for blt a is 4 =< 9, for b is −4 =< 3, et cetera. The rationale for this symbol is that in x =< y, the lowest subscript is equal to x, and every subscript is less than y. This symbol =< is unrelated to the C++ operators <= and >=.
If x and y are two blts, then x is said to contain y if y's bottom is not less than x's bottom, and y's top is not greater than x's top. This is akin to a subset relationship.
§ 3A2. Subscripts are of type mat_gen_dim::sign_int, which via a using declaration turns into int_least32_t, which is an integer type guaranteed to be signed and to have at least 32 bits. For simplicity, the same type is used for all other integer purposes throughout the implementation, in particular the form's count function. A 16-bit integer would not be able to handle a matrix with more than about 32,000 elements, which could be a problem for a large but realistic data structure. The 32-bit integer supports over 2,000 million components, which is enough for almost any purpose. (The present author is reluctant to write 2 billion, as the definition of that number varies.) For utterly massive data, the programmer can substitute int_least64_t.
If the underlying storage type for a matrix is std::vector<compo,alloc>, then why was not std::vector<compo,alloc>::size_type chosen for the subscript type? First of all, mat_gen_dim subscripts can be negative, but the vector's size_type is almost certainly unsigned. Second, a more subtle problem arises — the Standard Library does not guarantee the following to be the same type:
std::vector<compo1,alloc1>::size_type std::vector<compo2,alloc2>::size_typeIn for instance a function with inputs of types var_matrix<compo1,alloc1> and var_matrix<compo2,alloc2>, which integer type should be chosen for the internals of the calculation? This ambiguity resulted in the flat choice of int_least32_t for subscripts, and ultimately everything else integer.
§ 3B1. If the information is available at compile time, a form can be constructed with an initializer list:
form e {blt{2,6},blt{3,8},blt{1,5}};
// 3 dimensions, count == 4 × 5 × 4 == 100
form f {blt{1,5}};
// 1 dimension, count == 4
form g { };
// 0 dimensions, count == 1
A non-C++ symbol for form e is ( 2 =< 6, 3 =< 8, 1 =< 5 ).
Here is the critical step for creating a matrix whose dimensionality is not known at compile time:
std::vector<blt> svblt;
// manipulate svb as desired, making it any size, loading it with any values
form h {svblt};
// obtain first() and last() from any source
form i (first(), last());
No matter how the form is obtained, it can be passed directly to the matrix constructor.
§ 3B2. The form class offers these functions to find out how many blts a form contains, and to retrieve any one of them:
sign_int dimen () const;
blt get (sign_int const n) const;
and the blt class has these inspectors:
sign_int get_bot () const;
sign_int get_len () const;
sign_int get_top () const;
Via the get_form function of the matrix classes, described later, the programmer can learn the dimensionality and subscript ranges of a matrix whose characteristics are unknown; such might be the case when a matrix is received as a parameter to a function.
If x and y are two forms, then x is said to contain y if each blt of x contains the blt in the corresponding position of y. In other words, x contains y if every subscript that is valid for y is also valid for x. An application of continence is this: if a matrix's form is x, then y may be used to specify a submatrix to be copied from it.
§ 3C. What constitutes a valid subscript? Some examples will explain. Assume this:
form i {blt{2,6},blt{-3,8},blt{0,5}};
using vec_si = std::vector<sign_int>;
Here are some candidate subscripts:§ 3D. A body behaves something like a one-dimensional array of the component type. In fact, the matrix components are stored in a std::vector<compo,alloc>. This vector, which resides in a subbody, is kept private so that the user's program cannot change its size; that action would imply changing the size of the matrix, and that action requires considerably more detail if it is to be fully specified.
§ 4A. Classes var_matrix and con_matrix are template types derived from base type sub_matrix (file 1I_matrix.h). The programmer will have little need to create a sub_matrix object directly, although a reference to a sub_matrix might be convenient as the formal parameter of a function when either a var_matrix or con_matrix could be the actual parameter; it is serendipitous that a sub_matrix will not violate any of the constancy requirements of a con_matrix. Perhaps surprisingly, there is not a class named matrix.
// EXCERPTS FROM ...
template <typename compo, typename alloc = std::allocator<compo> >
class sub_matrix : public virtual sub_matrix_base<compo,alloc> {
form f;
body<compo,alloc> mutable b;
public:
sub_matrix ();
sub_matrix (form const &, compo const c = compo {});
sub_matrix (form const &, vector<compo,alloc> const &);
template <typename input_iter>
sub_matrix (form const &, input_iter, input_iter);
sub_matrix (form const &, initializer_list<compo>)
void assign_form (form const &);
std::vector <compo,alloc> const & get_vector () const;
sign_int count () const;
form const & get_form () const override;
};
// EXCERPTS FROM ...
template <typename compo, typename alloc = std::allocator<compo> >
struct var_matrix :
public var_matrix_base<compo,alloc>,
public sub_matrix<compo,alloc>
{
var_matrix (form const &, compo const c = compo {});
var_matrix (form const &, vector<compo,alloc> const & v);
template <typename input_iter>
var_matrix (form const &, input_iter, input_iter);
var_matrix (form const &, initializer_list<compo>);
var_matrix (sub_matrix_base<compo,alloc> const &);
var_matrix (sub_matrix<compo,alloc> const &);
compo get (vec_si const &) const override;
compo get (sign_int const) const override;
compo const & con_ref (vec_si const &) const override;
compo const & con_ref (sign_int const) const override;
compo & var_ref (vec_si const &) const override;
compo & var_ref (sign_int const) const override;
void set (vec_si const &, compo const) const override;
void set (sign_int const, compo const) const override;
template <typename input_iter>
void set_all (input_iter begin, input_iter end) const;
};
// EXCERPTS FROM ...
template <typename compo, typename alloc = std::allocator<compo> >
struct con_matrix :
public con_matrix_base<compo,alloc>,
public sub_matrix <compo,alloc>
{
con_matrix (form const &, compo const c = compo {});
con_matrix (form const &, vector<compo,alloc> const &);
template <typename input_iter>
con_matrix (form const &, input_iter, input_iter);
con_matrix (form const &, initializer_list<compo>);
con_matrix (sub_matrix_base<compo,alloc> const &);
con_matrix (sub_matrix<compo,alloc> const &);
con_matrix (var_matrix<compo,alloc> const &);
compo get (vec_si const &) const override;
compo get (sign_int const) const override;
compo const & con_ref (vec_si const &) const override;
compo const & con_ref (sign_int const) const override;
};
Here are some practical ways to build either a sub_, var_, or con_matrix from scratch:
form const f {blt{-4,-1},blt{+5,+8}};
var_matrix<double> m (f, { 0.0,1.1,2.2, 3.3,4.4,5.5, 6.6,7.7,8.8 }); // LEGAL
var_matrix<double> n (f, {{0.0,1.1,2.2},{3.3,4.4,5.5},{6.6,7.7,8.8}}); // ILLEGAL
When a matrix (whether sub_, var_, or con_) is constructed without parameters, it will have zero dimensions and its sole component will be initialized to the component type's default value. The following are some valid parameterless definitions that happen to illustrate nesting of matrices:
var_matrix<double> md;
var_matrix<var_matrix<double>> mdd;
var_matrix<var_matrix<var_matrix<double>>> mddd;
Differing from some array treatments is mat_gen_dim's policy that the following are not the same data type:
§ 4B. Within sub_matrix, and by inheritance var_matrix and con_matrix:
Within all three classes:
Access to components is by way of:
The set and set_all functions will look to see if multiple matrices are sharing a subbody, and if so will clone the subbody so that this matrix has an unshared copy. The var_ref functions must do the same, because they return a reference to a non-constant component, and there is risk that it will be changed. In contrast, the get and con_ref functions need never clone because they can never change the matrix.
Why not the customary functions at and operator[] as in the Standard Library? Because in their typical implementation, they return a reference to a non-constant component when called on a non-constant container. In mat_gen_dim, however, the functions when called on a matrix would have to clone if its subbody is shared, because of the mere possibility that the component will later be changed. As a result, there would likely be unnecessary clonings, which are often expensive. This is addressed in mat_gen_dim, where a con_ref can be obtained on a non-constant matrix, eliminating the risk of cloning.
§ 4C. The mat_gen_dim notion of what "constant" means may not be intuitive, but it bears a resemblance to Standard Library iterators.
§ 4D. There is a non-member function (file 1J_matrix_manip.h) to copy a rectangular block of components from one matrix to another, without changing any forms. The source can be a sub_, var_, or con_ matrix, but the destination must be a var_matrix. Here is a schematic of its definition:
template <typename compo, typename alloc>
void copy_block (
sub_matrix<compo,alloc> const & srce_mat,
form const & srce_inner_frm,
var_matrix<compo,alloc> & dest_mat,
form const & dest_inner_frm
) {
form const & srce_outer_frm {srce_mat.get_form()};
form const & dest_outer_frm {dest_mat.get_form()};
// more code ...
}
There are six objects to consider:§ 4E. There are functions (file 1J_matrix_manip.h) to create a sub_matrix of a floating-point type conveniently loaded with numbers for study or testing. The functions' main purpose is to initialize var_matrix and con_matrix.
template <typename compo = double, typename alloc = std::allocator<compo>>
sub_matrix<compo,alloc> sub_matrix_linear (form const &, double const);
template <typename compo = double, typename alloc = std::allocator<compo>>
sub_matrix<compo,alloc> sub_matrix_random (form const &, rand_double const &);
With sub_matrix_linear, each blt in the matrix's form should have a bottom value no less than zero, and a top no greater than ten; otherwise the pattern of numbers is obscured. If check is true, this becomes a requirement. The double parameter is added as an offset to every component of the matrix; if zero, the components will lie between zero and one. When a test program has several linear matrices, this offset makes it easier to tell them apart when components are listed.
With sub_matrix_random, the rand_double (file 1A_miscellany.h) constructor takes two parameters indicating what the range of random numbers should be. When a test program has several random matrices, they can have different ranges for easier distinguishing.
§ 4F. There are two functions (file 1K_matrix_join.h) to help manage those matrices whose components are themselves matrices.
discussion, sample program, and output
§ 5A. The reader may have noticed classes sub_matrix_base, var_matrix_base, and con_matrix_base (file 1H_matrix_base.h). Their principal purpose is to establish much of the interface of sub_matrix, var_matrix, and con_matrix, and they house only functions, mostly virtual and abstract. Classes X_matrix_base have no storage for forms or bodys; that responsibility is handled by sub_matrix for the benefit of var_matrix and con_matrix.
Classes X_matrix_base greatly aid the lenses explained in section 6, as a lens merely provides the illusion of a modified matrix, and it is returned as an object of a X_matrix_base type.
For functions that require a matrix input parameter, it is generally sufficient that the parameter be of type sub_matrix_base & when the parameter is not to be modified, and var_matrix_base & when it will be changed. An X_matrix_base parameter is matched by the result of a lens (other than X_lens_less_final); an X_matrix parameter is not.
Here is the relationship of the six classes:
template <typename compo, typename alloc=std::allocator<compo>>
struct sub_matrix_base;
template <typename compo, typename alloc=std::allocator<compo>>
struct var_matrix_base :
public virtual sub_matrix_base<compo,alloc>;
template <typename compo, typename alloc=std::allocator<compo>>
struct con_matrix_base :
public virtual sub_matrix_base<compo,alloc>;
template <typename compo, typename alloc=std::allocator<compo>>
class sub_matrix :
public virtual sub_matrix_base<compo,alloc>;
template <typename compo, typename alloc=std::allocator<compo>>
class var_matrix :
public var_matrix_base<compo,alloc>,
public sub_matrix<compo,alloc>;
template <typename compo, typename alloc=std::allocator<compo>>
class con_matrix :
public con_matrix_base<compo,alloc>,
public sub_matrix<compo,alloc>;
§ 5B. There is a critical difference between function copy_block above and copy_block_base, the schematic of which is here:
template <typename compo, typename alloc>
void copy_block_base (
sub_matrix_base<compo,alloc> const & srce_mat,
form const & srce_inner_frm,
var_matrix_base<compo,alloc> & dest_mat,
form const & dest_inner_frm
) {
form const & srce_outer_frm {srce_mat.get_form()};
form const & dest_outer_frm {dest_mat.get_form()};
// more code ...
}
Although they look similar, copy_block_base cannot determine whether it is copying from a matrix into itself. If it is, the components of the matrix will have unpredictable values because some may be overwritten before they are copied from. Still, the structure of the underlying matrix and intermediate lenses will remain sound.
There is no problem if the two matrices are distinct, but share a subbody. Also, if the function is copying into a subbody that has several bodys pointing into it, cloning will be performed correctly.
§ 6A. A lens causes a matrix to appear somehow different from what it really is. This section explains fourteen lenses, namely:
| var_ | con_ | file | behaves as |
|---|---|---|---|
| var_lens_other_form | con_lens_other_form | 2A_lens_other_form.h | X_matrix_base |
| var_lens_permu_form | con_lens_permu_form | 2B_lens_permu_form.h | |
| var_lens_permu_body | con_lens_permu_body | 2C_lens_permu_body.h | |
| var_lens_caten | con_lens_caten | 2D_lens_caten.h | |
| var_lens_less_range | con_lens_less_range | 2E_lens_less_range.h | |
| var_lens_less_dimen | con_lens_less_dimen | 2F_lens_less_dimen.h | mixture of X_matrix_base and iterator |
| var_lens_less_final | con_lens_less_final | 2G_lens_less_final.h | iterator |
These lenses can be combined sequentially for a composite effect, but no lens can follow an X_lens_less_final. Too, other lenses with completely different kinds of behavior can be built from scratch.
Component access is through functions get and con_ref which can be used freely with either var_lens_X and con_lens_X; and through functions var_ref and set which can be used only with var_lens_X. This parallels the functions' usage with var_matrix and con_matrix. Constructing a lens does not modify the underlying matrix, but var_ref and set can change the underlying matrix's component values.
Lenses are fragile, meaning that if a matrix is destructed or if its form is changed, any lenses relying on it are invalidated. Unfortunately, the lenses do not receive any notice of the changes, and any attempt to use them will quietly (and perniciously) yield unpredictable results. Something very similar happens in the Standard Library when an iterator is invalidated. Also, if lens x is relying on lens y, and then lens y is destructed or significantly changed, lens x becomes invalidated even if the underlying matrix is stable. In contrast, matrices are robust in part because of their reference-counting mechanism. Recommendation: in the midst of complicated calculations, important data should not be left in its apparent state in some lens, but rather copied into a matrix. While lenses do behave deterministically, they are error-prone because they are aliases.
A word on the relationships among var_matrix, con_matrix, var_lens_X, and con_lens_X is in order. With the convention that the hash mark stands for suitable parameters, this table classifes certain sequences of two statements:
| legal | illegal | ||
|---|---|---|---|
| var_matrix< # > M1 { # }; var_lens_X< # > A1 {M1, # }; | var_matrix< # > M2 { # }; con_lens_X< # > A2 {M2, # }; | con_matrix< # > M3 { # }; con_lens_X< # > A3 {M3, # }; | con_matrix< # > M4 { # }; var_lens_X< # > A4 {M4, # }; |
| var_lens_X< # > B1 { # }; var_lens_X< # > C1 {B1, # }; | var_lens_X< # > B2 { # }; con_lens_X< # > C2 {B2, # }; | con_lens_X< # > B3 { # }; con_lens_X< # > C3 {B3, # }; | con_lens_X< # > B4 { # }; var_lens_X< # > C4 {B4, # }; |
By contrast, the following are all legal because the apparent matrix within a fragile lens is being copied into a robust matrix:
| legal | |||
|---|---|---|---|
| var_lens_X< # > D1 { # }; var_matrix< # > N1 {D1}; | var_lens_X< # > D2 { # }; con_matrix< # > N2 {D2}; | con_lens_X< # > D3 { # }; con_matrix< # > N3 {D3}; | con_lens_X< # > D4 { # }; var_matrix< # > N4 {D4}; |
Other ways to copy lens values into a matrix are the assignment operator; and when the destination is a var_matrix, the copy_block function.
Meanwhile, the following are legal because the copy-on-write mechanism will guarantee copying the instant it is needed:
| legal | |||
|---|---|---|---|
| var_matrix< # > P1 { # }; var_matrix< # > Q1 {P1}; | var_matrix< # > P2 { # }; con_matrix< # > Q2 {P2}; | con_matrix< # > P3 { # }; con_matrix< # > Q3 {P3}; | con_matrix< # > P4 { # }; var_matrix< # > Q4 {P4}; |
Finally, some lenses require a utility class named permu (file 1B_permu.h) for housing permutations of the first n nonnegative integers. To create a permu object, the programmer supplies either a std::initializer_list<sign_int> or a std::vector<sign_int>. When mat_gen_dim::check is true, the sequence of numbers will be verified as a permutation of the first n nonnegative integers.
§ 6B: X_lens_other_form. These, the simplest lenses, cause a matrix to appear to have a different form. However, the sequence of components does not appear to be changed.
If the programmer needs an identity lens, in other words a lens that makes no apparent change of any kind, the most efficient way is to use X_lens_other_form with the same form as the underlying lens or matrix.
§ 6C: X_lens_permu_form. These appear to permute the form and hence the subscripts of the underlying matrix, with an apparent change in the sequence of components.
§ 6D: X_lens_permu_body. These appear to rearrange the components of the underlying matrix according to a permutation of one subscript. The form does not appear to be changed.
§ 6E: X_lens_caten. Two matrices that have abutting subscripts appear as one. A major reason for the matrices' assign_form function is for adjusting the subscripts of the ingredient matrices to be in phase prior to catenation. A shallower subscript adjustment could also be effected by X_lens_other_form.
§ 6F: X_lens_less_range. The subscript ranges appear to be reduced, but the number of dimensions does not.
§ 6G: X_lens_less_dimen. The number of dimensions appears to be reduced, but the subscript ranges in the remaining dimensions appear unchanged. There is a state, similar to that of a Standard Library iterator, to select which components are visible.
discussion, sample program, and output
§ 6H: X_lens_less_final. At any time, exactly one component is selected. There is a state, similar to that of a Standard Library iterator, that establishes which component that is.
discussion, sample program, and output
§ 7A. Nine functions are provided to assist in calculations that access all the elements of a matrix: some modify an input matrix and return nothing; some modify an input matrix and return its new value; and some modify no input matrix but return a new matrix. In all cases the programmer supplies a function object that specifies how the matrix components are to be manipulated. A function object is used instead of a plain old function to allow the programmer to employ a "smart function" if desired, perhaps one that remembers some values from one invocation to the next.
For the functions in this category, all matrices must have the same form, because calculations are performed in parallel. In other words, if two components (whether of the same matrix or different matrices) have different subscripts, the value of one does not influence calculations involving the other. (However, this principle might be impeached if the programmer supplies a function object that has side effects.)
Many combinations of types are supported, but the default is that all matrices involved have the same component type, and any scalar is also of that type; thus the rather formidable template parameter lists will usually be of little bother.
Similar functions with greater numbers of matrix or scalar parameters can be written if needed.
§ 7B1. These three functions (file 3A_alone.h) take for input an a matrix and a function object.
§ 7B2. These three functions (file 3B_scalar.h) take for input a matrix, a scalar, and a function object.
§ 7B3. These three functions (file 3C_parallel.h) take for input two matrices and a function object.
§ 7C1. The functions described in sections 7B are very generic; section 7C responds by presenting some practical implementations. This subsection covers unary operations (file 3A_alone.h).
§ 7C2. These functions (file 3B_scalar.h) handle ordinary multiplication or division between a matrix and scalar of a numerical type.
§ 7C3. These functions (file 3C_parallel.h) handle ordinary addition or subtraction of two numerical matrices.
§ 8A. Also supported are three operations (file 3D_algebra.h) of matrix and tensor algebra: outer product, contraction, and inner product; the last two in a very broad generalization. As stated above, these functions do not distinguish the superscripts and subscripts of tensor algebra. Also, contraction is not limited to being performed on exactly two subscripts.
These functions always return a new matrix without modifying any input matrix, because input and output matrices will typically have different forms.
If contraction is performed on the outer product of two matrices, what results is precisely their inner product. However, to calculate the inner product by way of contracting an outer product can be inefficient by orders of magnitude. Thus an inner product function is provided to allow smaller and faster direct calculation when the other results are not needed. The customary matrix multiplication of linear algebra is handled with the inner product function.
§ 8B. This function produces the outer product of two matrices.
discussion, sample program, and output
§ 8C1. This function contracts a matrix according to a std::vector<bool> supplied as a parameter.
discussion, sample program, and output
§ 8C2. Repeated contractions can be used to find the sum of a matrix's components.
§ 8D. This function produces the inner product of two matrices according to a std::vector<bool> supplied as a parameter.
discussion, sample program, and output
§ 8E. This example shows the outer product, contraction, and inner product together, using mixed component types for the matrices and suitable function objects for the arithmetic.
discussion, sample program, and output
§ 8F. This example shows how to implement the most common kind of matrix multiplication.
§ 8G. This example does matrix multiplication with block matrices, illustrating nesting of the matrix data structure.
§ 9. Here is this implementation's policy on allocators.
A matrix of practical size can easily contain thousands of components, so it makes sense to give the programmer the option of specifying what allocator to use, as in some cases great efficiency can be gained thereby. Whenever an allocator is mentioned in mat_gen_dim, it is intended as the allocator for components of a matrix.
Different matrices can have different component types, so multiple allocators must be supported. This is particularly evident with those functions that take multiple matrices as inputs, because each matrix can have a different component type, and thus a different allocator. Moreover, if the function returns a matrix as a result, it might have yet another component type and another allocator. As a result, some functions have lengthy template parameter lists; in mitigation, allocator defaults are provided with any class or function that the programmer is anticipated to use directly.
In order to keep things simple, mat_gen_dim uses the Standard Library's default allocator for all other allocations. This includes:
Can a programmer install custom allocators for these situations? Certainly. An alternative approach, however, is to observe that the length of a form (and the length of a vec_si, usually) is equal to the dimensionality of a nearby matrix. If the programmer predicts that most matrices in a program will not have more than (for instance) ten dimensions, a data structure could be employed that carries enough room for ten blts and a pointer to an overflow area. This avoids dynamic memory allocation for subforms holding fewer than eleven blts, enhancing speed at the cost of space. The cut-over point (eleven in this case) would naturally be an easily-changed compile-time constant.
§ 10. Because mat_gen_dim is intended to be generic across dimensions, code to support eigenvalues or other calculations specific to two dimensions is not likely to be added. On the other hand, the author is investigating a multi-dimensional generalization of matrix inversion, which will be explained here using a notation similar to that customarily encountered in linear algebra, rather than the computer-style notation above.
Matrices (two-dimensional in this case) A and B are said to be inverses of each other if they are square, and if their product is the identity matrix:
AB = BA = I
Given A in such an equation, we may be able to find B, although for some A there is no B that will succeed. More general is the following equation, where A and C are given, and B is sought. Again, a suitable B may or may not exist:
AB = C
Function mat_gen_dim::inner_product, supported by the parameter which, is capable of performing this customary matrix multiplication. In a mathematical-style notation, this inner_product might be represented by a circled × in prefix position, and the which parameter by w:
⊗ (A, B, w) = C
The question becomes: Given any A, w, and C, is there a B to satisfy the equation? Are there multiple solutions? Necessary but not sufficient for a solution is that the variables be conformable:
The dimensionality of B will turn out to equal the number of components in w minus the dimensionality of A.
§ 11. Miscellany.
Simple program to demonstrate that the dimension of a matrix need not be determined until run time.
Table of the major classes.
Table of the major non-member functions.
Purity: Unused in this implementation are the C++ operators const_cast, dynamic_cast, static_cast and reinterpret_cast.