Alrighty. So here we are. The format I was thinking of going with is: Reading through: Howard Anton's Elementary Linear Algebra 7th Edition and extending my cpp library as I learn new things.
1.1 System of Equation and Matricies
In this chapter most of the writing seems to be focused on explaining how systems of equations can be represented as a Matrix. So, here I am going to focus on the Matrix side of things and start our first class. The Matrix class.
The Matrix class ( Outline )
I want to plan out a bit how this Matrix class will look. In cpp we have the power of non-type template parameters. This means we can pass values (like constants!) into templates rather than types.
This seems to be the classic convention (see Eigen) and allows for compile time checks. These compile time checks are actually pretty cool in my opinion. We know (from reading!) that to multiply matricies A & B together, Matrix A must have the exact same number of rows as Matrix B has columns. Since we can pass values (the number of rows and columns) of our matricies in at compile time, we can enforce these constraints then (at compile time) rather than at runtime!
I also really like how Eigen provides the option to initialize matricies as follows:
Matrix3x3 mat;
mat << 1, 2, 3,
4, 5, 6,
7, 8, 9;
So let's try to throw that into our class as well. Well call that "comma initialization" (sounds cool).
So know how we want to initialize our class, we know how we want to define its size and value type. Time to choose how we want to store the data inside it.
I had given some thought to storing the values in a single array. Something like:
[1, 2, 3, 4, 5, 6, 7, 8, 9]
Rather than the classical 2d:
[
[1, 2, 3],
[4, 5, 6],
[7, 8, 9]
]
I think given a row, and column (r, c) getting a value is basically as easy in both cases.
Keeping with the classic (2d array) seems like the best approach and most straight forward. In swift there is a good reason for going the single array route (Wanting to use ublas/cblas in Accelerate) but the purpose of this class isn't for cpp as much as allowing us to implement some linear algebra so let's go with with the 2d approach.
I'll use std::vector as my array and create a type-alias so I don't have to continously type the arduous std::vector<std::vector....> each and every time.
Digression: Halfway through this post and I can tell this blog will be teaching me how to spell better.
Anyways, our type-alias will look like this:
template<typename K>
using Dataframe = std::vector<std::vector<K> >;
I'm calling this 2d vector a "Dataframe". I feel like a Dataframe is a more relaxed Matrix. It can take any type of value and it doesn't have any special methods you can call on it. I guess I think of it as the data structure behind a Matrix (??? does that make sense ???).
We use
usinghere rather thantypedefbecausetypedefdoes not work withtemplateparameters.
That is all that we need for section 1.1. Let's begin
The Matrix class
Start
First off is defining the class. We want the rows, columns, and type being stored to be passed in at compile time via non-type template arguments so lets start there.
Our template will look something like:
template<typename Type, uint32_t NumRows = 0, uint32_t NumCols = 0>
class Matrix {};
Here we've defined our template to take 3 parameters. The first of the three is the type which the matrix will hold, the second is the number of rows our matrix should have, and the last is the number of columns.
At the moment (and this is unlikely the change) our matrix will be filled stricly with numeric values. We shouldn't allow the user to pass in std::string or some other types.
Additionally, we want this check to happen at compile time again. To enable this we will use static_assert to run a check at compile time on the passed in type (Type), verifying that it in indeed a numeric value.
More specifically, we want the type to be either an integral value (short, int, long, long long..etc) or a floating-point value (float, double...).
Fortunately for us, cpp has a built in check: std:: is_arithmetic, which is exactly what we are looking for.
We end up with something like this,
template<typename Type, uint32_t NumRows = 0, uint32_t NumCols = 0>
class Matrix {
static_assert(std::is_arithmetic<Type>::value,
"Matrix value must have an integral or floating point type");
}; // Matrix
This is looking good!
Let's add in our Dataframe type-alias, and create our instance variable:
template<typename Type, uint32_t NumRows = 0, uint32_t NumCols = 0>
class Matrix {
// Make sure that the templated `Type` parameter is valid
static_assert(std::is_arithmetic<Type>::value,
"Matrix value must have an integral or floating point type");
public:
template<typename K>
using Dataframe = std::vector<std::vector<K> >;
private:
Dataframe<Type> matrix_;
}; // Matrix
Since we already know (at compile time) the number of rows and columns, we might as well initialize our matrix_ instance variable now.
// Matrix Class
template<typename Type, uint32_t NumRows = 0, uint32_t NumCols = 0>
class Matrix {
// Make sure that the templated `Type` parameter is valid
static_assert(std::is_arithmetic<Type>::value,
"Matrix value must have an integral or floating point type");
public:
template<typename K>
using Dataframe = std::vector<std::vector<K> >;
private:
Dataframe<Type> matrix_ = Dataframe<Type>(NumRows,
std::vector<Type>(NumCols));
}; // Matrix
Awesome. As narcissistic as I am, it's always good to test to make sure the code you wrote runs. Let's add in some #inlcude and build a main.cpp.
// Main.cpp
#include "linalg/matrix.h" // Import out lib
int main() {
mathlib::linalg::Matrix<float, 1, 1> matrix;
// Type Rows Cols
}
Works!
As hoped, building with:
// Main.cpp
#include "linalg/matrix.h" // Import out lib
#include <string>
int main() {
mathlib::linalg::Matrix<std::string, 1, 1> matrix;
// Type Rows Cols
}
fails with,
error: static_assert failed due to requirement 'std::is_arithmetic<std::__1::basic_string<char, std::__1::char_traits<char>, std::__1::allocator<char> > >::value' "Matrix value must have an integral or floating point type" static_assert(std::is_arithmetic<Type>::value,
Comma Init
Time to initialize. Building the comma init. is 100% unnecessary, but fun. So...
The idea is as follows:
- We want to overload the
<<operator so that when its called on our Matrix class, it returns a new class called CommaInit. - The CommaInit class has an overloaded
,operator which:- Adds the next value to our Matrix
- Returns itself (CommaInit) to pick up the next
,(if one exists))
This may seem compilicated, so let's dive into the code side of things and walk through it.
Overload <<
Overloading functions is common in many languages. In cpp we can overload the << by adding the following into our class:
CommaInit operator<<(Type v) {}
We have yet to write our CommaInit class but we already know we want our << overload to return that class type. Additionally we want to kick off this initialization when the type on the right side of the << matches the type we are storing in our Matrix.
Recall our goal:
Matrix3x3 mat;
mat << 1, 2, 3,
4, 5, 6,
7, 8, 9;
I will put into parenthesis where our operator is taking place:
(mat << 1), 2, 3,
4, 5, 6,
7, 8, 9;
After this operator occurs we will have:
commaInit, 2, 3,
4, 5, 6,
7, 8, 9;
which is exactly what we want. Now the commaInit object is next to a comma which will kick off that part of the sequence.
Before getting ahead of ourselves let's finish implementing the overload.
Notice how the first value (the 1) is gone after the overloaded function finishes. This means we need to place that value into our matrix before the overloaded function ends.
CommaInit operator<<(Type v) {
matrix[0][0] = v;
}
Now we must create the CommaInit object and return it. One thing we know is that at this point the CommaInit class is about to take over and we still somehow need a way of altering the matrix_. To achieve this we will pass the address og our matrix_ into the CommaInit constructor. This will allow CommaInit to mutate our matrix_.
CommaInit operator<<(Type v) {
matrix_[0][0] = v;
return CommaInit(&matrix_, 1);
// The 1 here denotes that we are upto the 1th index of our matrix.
// We already ingested the 0th index on the line before.
}
All that is left is to create the CommaInit class. I wont walk through this part exactly but the class is small and hopefully understandable:
struct CommaInit {
Dataframe<Type> *df;
uint32_t index;
CommaInit(Dataframe<Type> *df, uint32_t index)
: df(df), index(index) {};
// Overload the , operator
CommaInit operator,(Type x) {
(*df)[index / NumCols][index % NumCols] = x;
return CommaInit(df, index + 1); // Recursion-like
}
};
Let's add this to our Matrix class and try updating our main.cpp
// Matrix Class
template<typename Type, uint32_t NumRows = 0, uint32_t NumCols = 0>
class Matrix {
// Make sure that the templated `Type` parameter is valid
static_assert(std::is_arithmetic<Type>::value,
"Matrix value must have an integral or floating point type");
public:
template<typename K>
using Dataframe = std::vector<std::vector<K> >;
struct CommaInit {
Dataframe<Type> *df;
uint32_t index;
CommaInit(Dataframe<Type> *df, uint32_t index)
: df(df), index(index) {};
// Overload the , operator
CommaInit operator,(Type x) {
(*df)[index / NumCols][index % NumCols] = x;
return CommaInit(df, index + 1); // Recursion-like
}
};
// Everything here will be public
CommaInit operator<<(Type v) {
matrix_[0][0] = v;
return CommaInit(&matrix_, 1);
// The 1 here denotes that we are upto the 1th index of our matrix.
// We already ingested the 0th index on the line before.
}
private:
Dataframe<Type> matrix_ = Dataframe<Type>(NumRows,
std::vector<Type>(NumCols));
}; // Matrix
// main.cpp
#include "linalg/matrix.h"
int main() {
mathlib::linalg::Matrix<float, 1, 1> matrix
matrix << 1;
}
Seems to work. Let's confirm by adding a Print function which prints our matrix. Again, I'll just show the function. Nothing special here.
void Print() {
for(size_t i = 0; i < NumRows * NumCols; i++) {
std::cout << matrix_[i / NumCols][i % NumCols];
if ((i + 1) % NumCols == 0) {
std::cout << std::endl;
} else {
std::cout << " ";
}
}
}
We get,
// Matrix Class
template<typename Type, uint32_t NumRows = 0, uint32_t NumCols = 0>
class Matrix {
// Make sure that the templated `Type` parameter is valid
static_assert(std::is_arithmetic<Type>::value,
"Matrix value must have an integral or floating point type");
public:
template<typename K>
using Dataframe = std::vector<std::vector<K> >;
... See above for other code
void Print() {
for(size_t i = 0; i < NumRows * NumCols; i++) {
std::cout << matrix_[i / NumCols][i % NumCols];
if ((i + 1) % NumCols == 0) {
std::cout << std::endl;
} else {
std::cout << " ";
}
}
}
private:
Dataframe<Type> matrix_ = Dataframe<Type>(NumRows,
std::vector<Type>(NumCols));
}; // Matrix
// main.cpp
#include "linalg/matrix.h"
int main() {
mathlib::linalg::Matrix<float, 3, 3> matrix;
matrix << 1, 2, 3,
4, 5, 6,
7, 8, 9;
matrix.Print();
}
Running this give:
1 2 3
4 5 6
7 8 9
Exactly what we hoped for.
It's getting late and I feel this was a good start at our matrix class. Section 1.2 deals with Gaussian Elimination so looks like that is our next step!
-- Ehud