Posts Tagged valarray
C++ Valarray Matrix Madness
Numeric STL container valarray began life partially as an attempt to make C++ appealing to the supercomputing community. At the time the big thing in those big machines was, the ironically named, vector processing. However the valarray fell by the wayside as the people driving its development left the STL development group. Perhaps they realised that it didn’t really fit 100% with the STL, or maybe they just got sidetracked; who knows. But it is still useful, and here are a few reasons why:
- Can be used to write faster code for numeric spaces than possible with other STL types like the ubiquitous vector template.
- Offers the coder the possibility of staying within the comfort zone of the STL and familiar object oriented concepts with a small speed sacrifice over hand carved C.
- Allows library developers a way of writing optimized libraries for different environments so that coder can concentrate on the development at hand and not loose track in the complexity of the target environment.
So I’ve been playing around with valarray. It seemed that the best example to play with is that old classic the matrix. So here it is. Yes I know that there are some particularly hairy things wrong with it, but its not meant as a copy and paste solution. Its here as the results of a learning exercise and an example of what’s possible with valarray. There are one or two places which are not implemented or have not been tested, but you should be able to complete or test these by just looking at the rest of the examples. It should compile and run as is.
Have a look and have fun.
#include <fstream>
#include <iostream>
#include <iterator>
#include <string>
#include <valarray>
#include <vector>
/*------------------------------------------------------------------------------
matrix2d interface
This is the template for the matrix class. Its not that fancy and there could
be a lot more added to it but that was not the point of the exercise. It is
also limited by the number of types supported by the underlying valarray data
container. From an OO point of view its really not that good either given
that a lot of the manipulation or generator methods really don't need to have
access directly to the data part of the construct. However on a practical
basis including those methods in this class provide small benefits in speed
and allow things to be a little more obvious. There are also problems in the
way that the generators return new matrix2d objects. But I've no intention to
fix them as this was just a learning exercise. I've tossed the implementation
of the methods into seperate compilation objects to make the interface cleaner
for inspection purposes. In keeping with the valarray perspective there is no
bounds checking anywhere - you have been warned.
----------------------------------------------------------------------------*/
template<ypename T>
class matrix2d {
public:
// creates based on the rows and data size
matrix2d(std::size_t rows, std::valarray<T> data);
// creates an empty rows x size matrix
matrix2d(std::size_t rows, std::size_t cols);
// direct initialisation - beware that rows x cols must equal data.size()
matrix2d(std::size_t rows, std::size_t cols, std::valarray<T> data);
// get the number of rows in the matrix2d
std::size_t rows() const;
// get the number of columns in the matrix2d
std::size_t cols() const;
// get a copy of the data in the matrix2d
std::valarray<T> array() const;
// retrieve the data from row r of the matrix
std::valarray<T> row(std::size_t r) const;
// retrieve the data from col c of the matrix
std::valarray<T> col(std::size_t c) const;
// retrieve refernce to the data from row r of the matrix
std::slice_array<T> row(std::size_t r);
// retrieve refernce to the data from col c of the matrix
std::slice_array<T> col(std::size_t c);
// basic item reference
T & operator()(std::size_t r, std::size_t c);
// basic item retrieval
T operator()(std::size_t r, std::size_t c) const;
// generate a matrix sort each row then sort the rows - UNIMPLEMENTED
matrix2d<T> sort();
// genetate a new matrix that is the transposition of this one
matrix2d<T> transpose();
// generate a new matrix with this matrix's data and sort each row
matrix2d<T> rowItmSort();
// generate a new matrix with this matrix's data and sort each row in reverse
matrix2d<T> rowItmSortRev();
// generate a new matrix with this matrix's data and sort each col
matrix2d<T> colItmSort();
// generate a new matrix with this matrix's data and sort each col in reverse
matrix2d<T> colItmSortRev();
// generate a new matrix of this one with m appended below
matrix2d<T> appendRows(matrix2d<T> &m);
// generate a new matrix of this one with m appended to the right
matrix2d<T> appendCols(matrix2d<T> &m);
// generate a matrix of this one, upper left corner at row t col l - UNTESTED
matrix2d<T> extractMatrix2d(size_t t, size_t l, size_t w, size_t h);
protected:
std::size_t rows_;
std::size_t cols_;
std::valarray<T> data_;
};
/*------------------------------------------------------------------------------
matrix2d implementation
----------------------------------------------------------------------------*/
/*------------------------------------------------------------------------------
matrix2d constructors
----------------------------------------------------------------------------*/
template<class T>
matrix2d<T>::matrix2d(std::size_t rows, std::valarray<T> data) :
rows_(rows), cols_(data.size() / rows), data_(data) {}
template<class T>
matrix2d<T>::matrix2d(std::size_t rows, std::size_t columns) :
rows_(rows), cols_(columns), data_(rows * columns) {}
template<class T>
matrix2d<T>::matrix2d(std::size_t rows, std::size_t columns,
std::valarray<T> data) :
rows_(rows), cols_(columns), data_(data) {}
/*------------------------------------------------------------------------------
matrix2d operations
----------------------------------------------------------------------------*/
template<class T>
std::size_t matrix2d<T>::rows() const {
return rows_;
}
template<class T>
std::size_t matrix2d<T>::cols() const {
return cols_;
}
template<class T>
std::valarray<T> matrix2d<T>::array() const {
return data_;
}
template<class T>
std::valarray<T> matrix2d<T>::row(std::size_t r) const {
return data_[std::slice(r * cols(), cols(), 1)];
}
template<class T>
std::valarray<T> matrix2d<T>::col(std::size_t c) const {
return data_[std::slice(c, rows(), cols())];
}
template<class T>
std::slice_array<T> matrix2d<T>::row(std::size_t r) {
return data_[std::slice(r * cols(), cols(), 1)];
}
template<class T>
std::slice_array<T> matrix2d<T>::col(std::size_t c) {
return data_[std::slice(c, rows(), cols())];
}
template<class T>
T& matrix2d<T>::operator()(std::size_t r, std::size_t c) {
return data_[r * cols() + c];
}
template<class T>
T matrix2d<T>::operator()(std::size_t r, std::size_t c) const {
return row(r)[c];
}
/*------------------------------------------------------------------------------
matrix2d generators
----------------------------------------------------------------------------*/
template<class T>
matrix2d<T> matrix2d<T>::sort() {
matrix2d<T> result(rows_, cols_);
/* TO DO TO DO TO DO TO DO TO DO TO DO TO DO TO DO TO DO TO DO TO DO TO DO */
return result;
}
template<class T>
matrix2d<T> matrix2d<T>::transpose() {
matrix2d<T> result(cols_, rows_);
for (std::size_t i = 0; i < rows_; ++i)
result.col(i) = static_cast<std::valarray<T> > (row(i));
return result;
}
template<class T>
matrix2d<T> matrix2d<T>::rowItmSort() {
matrix2d<T> result(rows_, cols_);
for (std::size_t i = 0; i < rows_; ++i) {
std::valarray<T> x = static_cast<std::valarray<T> > (row(i));
std::sort(&x[0], &x[cols_]);
result.row(i) = x;
}
return result;
}
template<class T> bool rev (const T & a, const T & b) { return a > b; }
template<class T>
matrix2d<T> matrix2d<T>::rowItmSortRev() {
matrix2d<T> result(rows_, cols_);
for (std::size_t i = 0; i < rows_; ++i) {
std::valarray<T> x = static_cast<std::valarray<T> > (row(i));
std::sort(&x[0], &x[cols_], rev<T>);
result.row(i) = x;
}
return result;
}
template<class T>
matrix2d<T> matrix2d<T>::colItmSort() {
matrix2d<T> result(rows_, cols_);
for (std::size_t i = 0; i < cols_; ++i) {
std::valarray<T> x = static_cast<std::valarray<T> > (col(i));
std::sort(&x[0], &x[rows_]);
result.col(i) = x;
}
return result;
}
template<class T>
matrix2d<T> matrix2d<T>::colItmSortRev() {
matrix2d<T> result(rows_, cols_);
for (std::size_t i = 0; i < cols_; ++i) {
std::valarray<T> x = static_cast<std::valarray<T> > (col(i));
std::sort(&x[0], &x[rows_], rev<T>);
result.col(i) = x;
}
return result;
}
template<class T>
matrix2d<T> matrix2d<T>::appendRows(matrix2d<T> &m) {
matrix2d<T> result(rows_ + m.rows_, cols_);
result.data_[std::slice(0, rows_ * cols_, 1)] = data_;
result.data_[std::slice(rows_ * cols_, m.rows_ * m.cols_, 1)] = m.data_;
return result;
}
template<class T>
matrix2d<T> matrix2d<T>::appendCols(matrix2d<T> &m) {
matrix2d<T> result(rows_, cols_ + m.cols_);
std::size_t s1[] = {rows_,cols_}; // shape of left matrix
std::size_t p1[] = {result.cols_,1}; // position of left matrix in result
std::size_t s2[] = {m.rows_,m.cols_}; // shape of right matrix
std::size_t p2[] = {result.cols_,1}; // position or right matrix in result
std::valarray<std::size_t> sv1(s1, 2);
std::valarray<std::size_t> pv1(p1, 2);
std::valarray<std::size_t> sv2(s2, 2);
std::valarray<std::size_t> pv2(p2, 2);
result.data_[std::gslice(0, sv1, pv1)] = data_; // copy left matrix into place
result.data_[std::gslice(cols_, sv2, pv2)] = m.data_; // repeat for m
return result;
}
template<class T>
matrix2d<T> matrix2d<T>::extractMatrix2d(size_t x, size_t y, size_t w,
size_t h) {
/* TEST ME TEST ME TEST ME TEST ME TEST ME TEST ME TEST ME TEST ME TEST ME */
matrix2d<T> result(h, w);
size_t x2[] = {h, w}, s[] = {rows_, 1};
std::valarray<size_t> xa(x2, 2), sa(s, 2);
result.data_ = data_[(const std::gslice)std::gslice(y * rows_ + x, xa, sa)];
return result;
}
/*------------------------------------------------------------------------------
matrix2d general input output
This is a simple class to help collect methods of serialising the matrix2d
data forms. Really they should be done another way but once again I don't
really care as they are throw away code for the purpoese of demonstration
only.
----------------------------------------------------------------------------*/
template<typename T>
class matrix2dio {
public:
matrix2d<T> textToM2d(std::istream & is, size_t w);
matrix2d<T> fileToM2d(std::string file, size_t w);
void m2dToText(std::ostream & os, const matrix2d<T> & m);
void printValarray(std::ostream & os, const std::valarray<int> & va);
};
/*------------------------------------------------------------------------------
matrix2dio operations
----------------------------------------------------------------------------*/
typedef std::istream_iterator<int> int_istrm_iter;
template<class T>
matrix2d<T> matrix2dio<T>::textToM2d(std::istream & is, size_t w) {
std::vector<T> v;
std::valarray<T> a;
if (!is.good() || is.eof())
return;
copy(int_istrm_iter(is), int_istrm_iter(), back_inserter(v));
a.resize(v.size(), sizeof ( int));
copy(v.begin(), v.end(), &a[0]);
return new matrix2d<T > (w, a);
}
template<class T>
matrix2d<T> matrix2dio<T>::fileToM2d(std::string file, size_t w) {
std::filebuf fb;
std::istream is(&fb);
fb.open(file.c_str(), std::ios::in);
matrix2d<T> m = textToM2d(is, w);
fb.close();
return m;
}
template<class T>
void matrix2dio<T>::m2dToText(std::ostream & os, const matrix2d<T> & m) {
size_t i = 0, j = 0;
for (i = 0; i < m.rows(); i++) {
std::valarray<T> r = m.row(i);
os << r[0];
for (j = 1; j < m.cols(); j++)
std::cout << ' ' << r[j];
os << '\n';
}
os << std::flush;
}
template<class T>
void matrix2dio<T>::printValarray(std::ostream & os,
const std::valarray<int> & va) {
copy(&va[0], &va[va.size() - 1], std::ostream_iterator<T > (os, " "));
os << va[va.size() - 1] << std::endl;
}
/*------------------------------------------------------------------------------
matrix2d tests
----------------------------------------------------------------------------*/
void testConstructors() {
std::cout << "\n\n\nRunning Constructor Tests\n\n";
int a[] = { 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12 };
std::valarray<int> v( a, 12 );
std::size_t x = 3, y = 4;
matrix2dio<int> i;
std::cout <<
"\nTesting: matrix2d(std::size_t rows, std::size_t columns, "
"std::valarray<T> data);\n";
matrix2d<int> m1( x, y, v );
std::cout << "number of rows: " << m1.rows() << '\n'
<< "number of cols: " << m1.cols() << '\n'
<< "matrix content:\n";
i.m2dToText( std::cout, m1 );
std::cout << "\nTesting: matrix2d(std::size_t rows, std::size_t columns);\n";
matrix2d<int> m2( x, y );
std::cout << "number of rows: " << m2.rows() << '\n'
<< "number of cols: " << m2.cols() << '\n'
<< "matrix content:\n";
i.m2dToText( std::cout, m2 );
std::cout <<
"\nTesting: matrix2d(std::size_t rows, std::valarray<T> data);\n";
matrix2d<int> m3( x, v );
std::cout << "number of rows: " << m3.rows() << '\n'
<< "number of cols: " << m3.cols() << '\n'
<< "matrix content:\n";
i.m2dToText( std::cout, m3 );
}
void testRowsAndCols() {
std::cout << "\n\n\nRunning Row/Col Accessor Tests\n\n";
int a[] = { 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12 };
std::valarray<int> v( a, 12 );
std::size_t x = 3, y = 4;
matrix2dio<int> i;
matrix2d<int> m1( x, y, v );
std::cout << "\nTesting: std::valarray<T> row(std::size_t r) const;\n";
std::valarray<int> r1 = m1.row(0);
std::valarray<int> r2 = m1.row(1);
std::valarray<int> r3 = m1.row(2);
std::cout << "row 1:\n";
i.printValarray( std::cout, r1 );
std::cout << "row 2:\n";
i.printValarray( std::cout, r2 );
std::cout << "row 3:\n";
i.printValarray( std::cout, r3 );
std::cout << "\nTesting: std::valarray<T> col(std::size_t r) const;\n";
std::valarray<int> c1 = m1.col(0);
std::valarray<int> c2 = m1.col(1);
std::valarray<int> c3 = m1.col(2);
std::valarray<int> c4 = m1.col(3);
std::cout << "col 1:\n";
i.printValarray( std::cout, c1 );
std::cout << "col 2:\n";
i.printValarray( std::cout, c2 );
std::cout << "col 3:\n";
i.printValarray( std::cout, c3 );
std::cout << "col 4:\n";
i.printValarray( std::cout, c4 );
std::cout << "\nTesting: std::slice_array<T> row(std::size_t r);\n";
std::slice_array<int> rs1 = m1.row(0);
std::slice_array<int> rs2 = m1.row(1);
std::slice_array<int> rs3 = m1.row(2);
std::cout << "row 1:\n";
i.printValarray( std::cout, rs1 );
std::cout << "row 2:\n";
i.printValarray( std::cout, rs2 );
std::cout << "row 3:\n";
i.printValarray( std::cout, rs3 );
std::cout << "\nTesting: std::slice_array<T> col(std::size_t r);\n";
std::slice_array<int> cs1 = m1.col(0);
std::slice_array<int> cs2 = m1.col(1);
std::slice_array<int> cs3 = m1.col(2);
std::slice_array<int> cs4 = m1.col(3);
std::cout << "col 1:\n";
i.printValarray( std::cout, cs1 );
std::cout << "col 2:\n";
i.printValarray( std::cout, cs2 );
std::cout << "col 3:\n";
i.printValarray( std::cout, cs3 );
std::cout << "col 4:\n";
i.printValarray( std::cout, cs4 );
}
void testGenerators() {
std::cout << "\n\n\nRunning Generator Tests\n\n";
int a[] = { 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12 };
int b[] = { 10, 20, 30, 40, 50, 60, 70, 80, 90, 100, 110, 120 };
std::valarray<int> v1( a, 12 );
std::valarray<int> v2( b, 12 );
std::size_t x = 3, y = 4;
matrix2dio<int> i;
matrix2d<int> m1( x, y, v1 );
matrix2d<int> m2( x, y, v2 );
std::cout << "\nTesting: matrix2d<T> transpose();\noriginal:\n";
matrix2d<int> m3 = m1.transpose();
i.m2dToText( std::cout, m1 );
std::cout << "transposed:\n";
i.m2dToText( std::cout, m3 );
std::cout << "\nTesting: matrix2d<T> appendRows(matrix2d<T> &v);\n";
matrix2d<int> m4 = m2.appendRows(m1);
i.m2dToText( std::cout, m4 );
std::cout << "\nTesting: matrix2d<T> appendCols(matrix2d<T> &m);\n";
matrix2d<int> m5 = m2.appendCols(m1);
i.m2dToText( std::cout, m5 );
std::cout << "\nTesting: matrix2d<T> rowItmSort();\n";
matrix2d<int> m6 = m5.rowItmSort();
i.m2dToText( std::cout, m6 );
std::cout << "\nTesting: matrix2d<T> rowItmSortRev();\n";
matrix2d<int> m7 = m5.rowItmSortRev();
i.m2dToText( std::cout, m7 );
std::cout << "\nTesting: matrix2d<T> colItmSort();\n";
matrix2d<int> m8 = m4.colItmSort();
i.m2dToText( std::cout, m8 );
std::cout << "\nTesting: matrix2d<T> colItmSortRev();\n";
matrix2d<int> m9 = m4.colItmSortRev();
i.m2dToText( std::cout, m9 );
}
void testMatrix2d() {
testConstructors();
testRowsAndCols();
testGenerators();
}
/*------------------------------------------------------------------------------
matrix2d test entry point
----------------------------------------------------------------------------*/
int main( int argc, char** argv ) {
testMatrix2d();
return (EXIT_SUCCESS);
}
Fussing with slice_array
Some times we just make things far to complicated for ourselves: note to self (and whoever cares),do not put compiler hints into unproven code. Every now and then you overstep your confidence and do something that you look back at and think, “Now that’s just stupid!”. I’ve just spent the last several days wondering why the following does not compile.
void foo( const std::valarray<int> & va ) {
std::slice_array<int> sa = va[std::slice( 1, 1, 1 )];
}
After all the line in the middle is exactly the same as sooo much example code out there it isn’t funny. However the compiler throws an awful wobbly and kicks back an error like this.
example.cc: In function ‘void foo(const std::valarray<int>&)’:
example.cc:20: error:
conversion from
‘std::_Expr<std::_SClos<std::_ValArray, int>, int>’
to non-scalar type
‘std::slice_array<int>’
requested
Which is not quite what I expected. I’ve laid it out a bit different to the compiler for readability. Essentially the error means that the type of the expression on the right is not equivalent to the type of the variable on the left. However you can do this.
void foo( const std::valarray<int> & va ) {
std::valarray<int> sa = va[std::slice( 1, 1, 1 )];
}
Which of course will compile without issue. OK so the operator() is overloaded, what’s its problem? The solution is blindingly obvious, but for those with a sense of anticipation I’ll drag it out just a tiny wee bit more. The mystery, and partly where I went wrong, was deepened by the valarray header.
/** * @brief Return an array subset. * * Returns a new valarray containing the elements of the array * indicated by the slice argument. The new valarray has the same size * as the input slice. @see slice. * * @param s The source slice. * @return New valarray containing elements in @a s. */ _Expr<_SClos<_ValArray, _Tp>, _Tp> operator[](slice) const;/**
* @brief Return a reference to an array subset.
*
* Returns a new valarray containing the elements of the array
* indicated by the slice argument. The new valarray has the same size
* as the input slice. @see slice.
*
* @param s The source slice.
* @return New valarray containing elements in @a s.
*/
slice_array<_Tp> operator[](slice);
Now that's probably give it away if you didn't already know. But once again I wasn't looking properly. There were two real problems and none of them to do with my code: firstly I was programming by Google, secondly I was trying to be smart. The first error was as stupid as a cut and paste error. I'd looked up similar issues in Google and cut someone else's understanding of the problem out of a page that I had read and stapled it into my brain. Seems someone else had the same problem but not the nouse to get over it and was blaming the spec. Very silly, I should have known better. The second thing was that I'd tried to hint to the compiler before understanding its implications. Never hint to the compiler before you have the thing working, I'd have never wasted my time with it, although programming by debugger has its own worms.
Turns out this is the code that I wanted.
void foo( std::valarray<int> & va ) {
std::slice_array<int> sa = va[std::slice( 1, 1, 1 )];
}
See it? If you still can't see it compare it to the earlier version and have a look at the prototypes in copied from the header above, which should explain almost everything. All that time. Silly, silly silly ... proceed with self flagellation of 50 lashings with a wet noodle! I have 15 years experience programming with C++ and 20 with C (which can generate similar errors), very embarrassing.
Playing with Valarray
In the past few weeks I’ve been looking at the C++ STL numeric container valarray. Its kind of like that forgotten ternary operator : ? in C (yes I know lots of people use it but most just use it for Shiek appeal) and there are quite a few ancient posts out there that have faded into the mists of time. Prime conjecture says its because valarray, like the ternary operator, is really a niche item seemingly randomly placed in a generalised library. So lets take a look at it.
My wanderings started with a database where I was trying to do numeric analysis. Not normally a bad thing when you have a very simple task, sum this row, add that column to that thingy over there. But even a simple task like a pivot table starts taking on a whole complex life of its own; and you can just forget brute force function sequence processing using sliding windows (generating a sequence from functions over subsequences of another). Believe it or not I got this working using MySQL but it was pure nightmare territory. Hence I started looking at coding it up by hand.
First thought was to use a library that someone else had written. But it soon became obvious that most of the BLAS libraries were written by boffins dedicated to their own arcane mathematical religious systems and not really understandable to someone with only a second year University course of linear algebra under his belt. Others were too simple and lacked the power of expression that was needed. What was really painfully obvious was that levels of abstraction were, in most of the examples out there quite an abstraction from the main goals of anyone using them, speed and comprehension. There’s a rule in security that, in my paraphrase, says ‘complexity is where the issues hide’. KISS
So on to the C and C++ libraries the search went at last there was something to use, valarray. Now valarray is NOT the fastest thing in code. I read lots of complaints on the Interweb about newbies who claim to make faster code with their hand coded C. They’d be right. Its blindingly simple to write faster code. However there are a couple of things to consider before throwing in the object-oriented towel and carving your own C. Firstly, and not very helpfully, valarray is an incomplete additive to a general purpose library. Its inventors were really looking to take advantage of vector processors, not something you see much of today – or so you think. So its intentions may not fit your needs. Next the valarray is slower than your carefully carved C but it is already there for you and it works faster than a vector; the C++ programmers back stop. Lastly coding to a standard helps because there are so many people who have been there before with it its probably, little though it might be, a lot better thought out than what you may come up with in the time that you have. Still its horses for courses so don’t be afraid to ditch it if its not right for you.
Here’s a little code to get stuck into:
#include <algorithm>
#include <iostream>
#include <valarray>
#include <vector>
using namespace std;
/*
* import a valarray from a whitespace seperated stream of numbers
*/
typedef istream_iterator<int> int_istrm_iter;
void importTextArray( istream & is, valarray<int> & va ) {
vector<int> v;
if (!is.good() || is.eof())
return;
copy( int_istrm_iter( is ), int_istrm_iter(), back_inserter( v ) );
va.resize( v.size(), sizeof( int ) );
copy( v.begin(), v.end(), &va[0] );
}
Hopefully you should be able to work out that it imports a file full of white space separated numbers into a vector and then sizes a valarray up for the data and copies it in. Naturally the first parameter is a file or stringbuffer or whatever that contains the data. Once you’ve done this you’ve got something to play with. So let’s play.
Of course a single dimension like valarray isn’t really that much use. What we want to do is be able to do mad things with matricies of n dimensions and repeatedly alter them till we’re happy. A little bit of what-if analysis if you please. There are a number of ways to do this and they involve the slice and gslice classes. You can work that out for yourself, its in most references on the topic. So lets do something a little more fun, lets stick two matricies together side by side.
/* * creates a new matrix by appending one matrix to the side of the other, note * that the three matrcies have the same row count r * * for example: * * matrix a matrix b matrix c * 01 02 03 04 01 02 03 01 02 03 04 01 02 03 * 05 06 07 08 + 04 05 06 = 05 06 07 08 04 05 06 * 09 10 11 12 07 08 09 09 10 11 12 07 08 09 */
void appendMatRows( const int height, const valarray<int> &first,
const valarray<int> &second, valarray<int> &result ) {
size_t firstwidth = first.size() / height;
size_t secondwidth = second.size() / height;
size_t resultwidth = firstwidth + secondwidth;
// make the result matrix the right size
result.resize( first.size() + second.size() );
// use a gslice to get a valarray for inserting the first matrix
size_t x1[] = { height, firstwidth }; // shape of extracted array
size_t s1[] = { resultwidth, 1 }; // parent row length and row item distance
valarray<size_t> xa1( x1, 2 );
valarray<size_t> sa1( s1, 2 );
result[( const gslice )gslice( 0, xa1, sa1 )] = first;
// then repeat with the second matrix
size_t x2[] = { height, secondwidth };
size_t s2[] = { resultwidth, 1 };
valarray<size_t> xa2( x2, 2 );
valarray<size_t> sa2( s2, 2 );
result[( const gslice )gslice( firstwidth, xa2, sa2 )] = second;
}
So that’s about it for now. I guess the thing about valarray for me is that there is a good base to work from which, at a little cost, saves me some time and pain in hand coding my own. It still needs a lot of basic work to do some of the things that are needed to play with but its pretty good for what my requirements are. Incidentally MySQL goes away for a while when performing these sorts of operations; whereas the equivalent coded version is uicker to develop, clearer to debug and over a thousand row matrix seems to process with out even bothering the CPU. That’s a great place to start.