ACM Computing Seminar C++ Guide

Pointers are variables that hold the memory address for a variable of a specific type. Pointers are declared by specifying the variable type, followed by the * symbol, followed by the name of the pointer variable, e.g., double * x defines a "pointer to double" variable. The variable, x, therefore, does not hold the value of a double type, but rather, the memory address for a variable of type, double. The memory address for a variable can be obtained by the & operator.

1: double * a;
2: double b = 7;
3: 
4: // This obtains the memory address of `b`.
5: a = &b;
6: 
7: // Prints some memory address (starts with 0x)
8: std::cout << a << std::endl;

0x7ffe0d98f7b8

Similar to obtaining the memory address from a regular variable, using the & operator, you can use the * symbol before a pointer to access the variable value held at the memory location of the pointer. In this context, the * symbol is called the dereference operator. This is probably better understood with a short example:

 1: double * a;
 2: double b = 7.3;
 3: double c;
 4: 
 5: // Now `a` holds the memory address of `b`.
 6: a = &b;
 7: 
 8: // `*a` obtains the value of the variable
 9: // at the memory address held by `a`.
10: // So, `c` is 7.3.
11: c = *a;
12: 
13: std::cout << c << "\n";

7.3

A reference is a sort of like a pointer, but not quite. There are differences. A good analogy, which you can find in the previous link, is that a reference is similar to a symbolic link, or "shortcut" if you're on Windows. You can treat it more-or-less like the original variable, but it's not the original.

 1: double a = 1.1;
 2: // `b` is a reference to `a`.
 3: double & b = a;
 4: 
 5: std::cout << "a: " << a << ", b: " << b << "\n";
 6: 
 7: a = 2.1;
 8: 
 9: std::cout << "a: " << a << ", b: " << b << "\n";
10: 
11: b = 3.1;
12: 
13: std::cout << "a: " << a << ", b: " << b << "\n";
14: 
15: std::cout << "\n\n";
16: std::cout << "&a: " << &a << "\n" << "&b: " << &b << "\n";

a: 1.1, b: 1.1
a: 2.1, b: 2.1
a: 3.1, b: 3.1


&a: 0x7ffcfbe7b1e8
&b: 0x7ffcfbe7b1e8

References are useful for passing around large objects, so that the object doesn't need to be copied. References are also useful as a return type for functions (to be discussed later) because it allows to assign to assign a value to a function, which is useful if the function, for example, returns a reference to the element of an array.

Dynamic length arrays are made possible through pointers:

 1: // This allocates memory for 5 double types.
 2: double * a = new double[5];
 3: 
 4: // Afterwards, you can treat `a` like a normal array.
 5: a[0] = 1.0;
 6: // etc...
 7: 
 8: // Whenever you use the `new` keyword, you must
 9: // delete the memory allocated when you're done by hand.
10: delete [] a;
11: 
12: // We can change the size of `a`.
13: a = new double [10];
14: 
15: a[0] = 2.0;
16: // etc...
17: 
18: delete [] a;

Note that omitting the first delete statement will cause no error. However, the memory allocated by the first new statement will not be freed, and thus inaccessible. This is bad because the memory cannot be allocated to other resources. You should generally try to avoid manually memory management when possible, but a good tool for debugging memory problems is called valgrind.

Often a code block should only be executed if some condition is true. Below, we generate a random number between 0 and 1; print the number; and, print whether or not the number was greater than 0.5.

 1: #include <iostream>
 2: #include <stdlib.h>
 3: #include <time.h>
 4: 
 5: int main() {
 6:   // Seed a random number generator.
 7:   srand(123);
 8: 
 9:   // rand() produces a random integer between 0 and RAND_MAX.
10:   double num = rand() / ((double) RAND_MAX);
11: 
12:   std::cout << "num: " << num << "\n";
13: 
14:   if (num < 0.5) {
15:     std::cout << "num was less than 0.5.\n";
16:   }
17:   else {
18:     std::cout << "num was greater than 0.5.\n";
19:   }
20: 
21:   // Do it again.
22:   num = rand() / ((double) RAND_MAX);
23: 
24:   std::cout << "num: " << num << "\n";
25: 
26:   if (num < 0.5) {
27:     std::cout << "num was less than 0.5.\n";
28:   }
29:   else {
30:     std::cout << "num was greater than 0.5.\n";
31:   }
32: 
33:   return 0;
34: }

num: 0.0600514
num was less than 0.5.
num: 0.788318
num was greater than 0.5.

You can follow else immediate by another if to have multiple mutually- exclusive blocks:

 1: #include <iostream>
 2: #include <stdlib.h>
 3: #include <time.h>
 4: 
 5: int main() {
 6:   // Seed the random number generator based on the current time.
 7:   srand(time(NULL));
 8: 
 9:   // rand() produces a random integer between 0 and RAND_MAX.
10:   double num = rand() / ((double) RAND_MAX);
11: 
12:   std::cout << "num: " << num << "\n";
13: 
14:   if (num >= 0.75) {
15:     std::cout << "num was between 0.75 and 1.\n";
16:   }
17:   else if (num >= 0.5) {
18:     std::cout << "num was between 0.5 and 0.75.";
19:   }
20:   else if (num >= 0.25) {
21:     std::cout << "num was between 0.25 and 0.5.";
22:   }
23:   else {
24:     std::cout << "num was between 0 and 0.25";
25:   }
26: 
27:   return 0;
28: }

num: 0.0456405
num was between 0 and 0.25

The conditions are checked in the order that they're written. So, for example, in the second condition, we don't need to specify num >= 0.5 && num < 0.75 because we know that this condition will only be checked if the previous was false.

The for loop requires three specifications – the iteration variable initialization, the termination condition, and the update rule. The body of the loop follows these three specifications. Shown below, we declare an array; assign to its components; and, print the current component to the screen.

 1: int length = 11;
 2: double x[length];
 3: 
 4: for(int i=0; i < length; i++) {
 5:   // Assign to each array component.
 6:   x[i] = (double) i / (length - 1);
 7: 
 8:   // Print the current component.
 9:   std::cout << "x[" << i << "] = " << x[i] << std::endl;
10: }

x[0] = 0
x[1] = 0.1
x[2] = 0.2
x[3] = 0.3
x[4] = 0.4
x[5] = 0.5
x[6] = 0.6
x[7] = 0.7
x[8] = 0.8
x[9] = 0.9
x[10] = 1

A while loop iterates while a condition is true. Essentially, it is a for loop without an update variable.

The break keyword provides a mechanism for exiting the direct parent loop for which the break statement is placed. For example:

1: for(int i=0; i < 3; i++) {
2:   while(true) {
3:     std::cout << "Entering infinite loop number " << (i+1) << "\n";
4:     break;
5:   }
6:   std::cout << "We escaped the infinite loop!\n";
7: }

Entering infinite loop number 1
We escaped the infinite loop!
Entering infinite loop number 2
We escaped the infinite loop!
Entering infinite loop number 3
We escaped the infinite loop!

The previous example is contrived, but there are situations, where you might find the break statement within an infinite loop useful. Of course, you should avoid this sort of thing if there is a more straight-forward approach.

In general, how data is read in depends heavily on how the data is stored. Nevertheless, we will give an example of reading in a vector stored in a particular fashion. Suppose a text file exists in the directory, ./data/vector.txt, containing

1 2 3.14 4 5 6.28

 1: #include <iostream>
 2: #include <fstream>
 3: 
 4: int main() {
 5:     std::fstream fin("./data/vector.txt", std::ios_base::in);
 6:     double vector[6];
 7:     int i = 0;
 8:     while(fin >> vector[i]) {
 9:       std::cout << vector[i] << " ";
10:       i++;
11:     }
12:     return 0;
13: }

g++ main.cpp && ./a.out

1 2 3.14 4 5 6.28

This simply prints the data in the file back out to the console. Note, however, that the data is read into an array of type double, so it can be processed numerically thereafter.

In this example dealt with simply stored data, and it was assumed that the number of data entries was known beforehand. Parsing data can become quite complicated depending on how it is stored, and depending on the intended format of the data.

Writing to a file is similar, using the <fstream> library.

 1: #include <fstream>
 2: #include <cmath>
 3: 
 4: int main() {
 5:   std::fstream fout("./data/new_shiny_data.txt", std::ios_base::out);
 6:   double x;
 7: 
 8:   fout << "x\tsin(x)\n";
 9: 
10:   for(int i=0; i < 11; i++) {
11:     x = i / 10.0;
12:     fout << x << "\t" << sin(x) << "\n";
13:   }
14: 
15:   fout.close();
16: 
17:   return 0;
18: }

This produces a file called new_shiny_data.txt in the folder, data, containing:

x   sin(x)
0   0
0.1 0.0998334
0.2 0.198669
0.3 0.29552
0.4 0.389418
0.5 0.479426
0.6 0.564642
0.7 0.644218
0.8 0.717356
0.9 0.783327
1   0.841471

 1: #include <iostream>
 2: 
 3: // This is the function declaration.
 4: // You should describe the functions arguments
 5: // and what is returned by the function in comments
 6: // near the declaration.
 7: //
 8: // `linspace` returns an array of doubles containing
 9: // `n_points` entries which are equally-spaced, starting
10: // at `start` and ending at `stop`.
11: double * linspace(double start, double stop, int n_points);
12: 
13: // `void` is a function with no return type.
14: // `print_array` takes an array and prints it to std out.
15: void print_array(double * arr, int arr_len);
16: 
17: int main() {
18:   double * xs = linspace(-1, 1, 5);
19:   print_array(xs, 5);
20:   delete [] xs;
21: 
22:   return 0;
23: }
24: 
25: // Implementation of `linspace`.
26: double * linspace(double start, double stop, int n_points) {
27:   double * arr = new double [n_points];
28:   double dx = (stop-start) / (n_points-1.0);
29: 
30:   for(int i=0; i < n_points; i++) {
31:     arr[i] = start + i*dx;
32:   }
33: 
34:   return arr;
35: }
36: 
37: // Implementation of `print_array`.
38: void print_array(double * arr, int arr_len) {
39:   for(int i=0; i < arr_len; i++) {
40:     std::cout << arr[i] << "\n";
41:   }
42: }

-1
-0.5
0
0.5
1

Put the declarations from the into a header file, called my_library.h:

 1: #ifndef MY_LIBRARY_H
 2: #define MY_LIBRARY_H
 3: 
 4: #include <iostream>
 5: 
 6: namespace my_namespace {
 7:     // `linspace` returns an array of doubles containing
 8:     // `n_points` entries which are equally-spaced, starting
 9:     // at `start` and ending at `stop`.
10:     double * linspace(double start, double stop, int n_points);
11: 
12:     // `void` is a function with no return type.
13:     // `print_array` takes an array and prints it to std out.
14:     void print_array(double * arr, int arr_len);
15: }
16: 
17: #endif

Note the the function declarations are wrapped in conditional "macro" statments, #ifndef, #define, and #endif. You can think of this as protecting your library from being imported twice.

We have also introduced the notion of a namespace above. Namespaces help to prevent naming clashes between separate libraries. When calling a function from a particular namespace, you must write the namespace followed by :: and then the function name. This is why many standard library functions like <iostream> begin with std::.

Create a file called my_library.cpp containing the implementations as follows:

 1: #include "my_library.h"
 2: 
 3: // Implementation of `linspace`.
 4: double * my_namespace::linspace(double start, double stop, int n_points) {
 5:   double * arr = new double [n_points];
 6:   double dx = (stop-start) / (n_points-1.0);
 7: 
 8:   for(int i=0; i < n_points; i++) {
 9:     arr[i] = start + i*dx;
10:   }
11: 
12:   return arr;
13: }
14: 
15: // Implementation of `print_array`.
16: void my_namespace::print_array(double * arr, int arr_len) {
17:   for(int i=0; i < arr_len; i++) {
18:     std::cout << arr[i] << "\n";
19:   }
20: }

Note that we have to include the header file in quotations at the beginning, and the names of the functions must be prepended by the namespace that we've given in the header file.

Create a file with the main function, say main.cpp:

 1: #include <iostream>
 2: #include "my_library.h"
 3: 
 4: int main() {
 5:     double * xs = my_namespace::linspace(-1,1,5);
 6:     my_namespace::print_array(xs, 5);
 7:     delete [] xs;
 8: 
 9:     return 0;
10: }

Now the main function is very nice and clean, but now we 3 separate files we must compile into one executable. This is done as follows:

# Convert the library into an object file.
g++ -c my_library.cpp
# Compile the main to an executable.
g++ my_library.o main.cpp
# Run it.
./a.out

-1
-0.5
0
0.5
1

If successful, you will see the same output as previously.

Suppose \(f: \mathbb{R} \to \mathbb{R}\), and we'd like to find a root of \(f\). Newton's method is an iterative method for finding roots, which, starting from some initial guess, \(x_0\), iterates:

\[ x_{n+1} \leftarrow x_n - \frac{f(x_n)}{f'(x_n)} \]

For simplicity, we'll dump everything into the file containing main, but you could imagine a libary with many methods for finding roots, which would contain Newton's method.

Let's consider \(f(x) = x^2 - 2\).

 1: #include <cmath>
 2: #include <iostream>
 3: 
 4: // The function to find the root of.
 5: double func(double x);
 6: // Its derivative.
 7: double dfunc(double x);
 8: 
 9: // Find the root of `f` using Newton's method,
10: // starting from `x0` until |f(x)| < `tol` or `max_iters`
11: // is reached.
12: //
13: // Note the first and second arguments are function pointers.
14: double newton_root(double (*f)(double), double (*df)(double), double x0,
15:                    double tol, int max_iters, bool print_iters); 
16: 
17: int main() {
18:   double x = newton_root(&func, &dfunc, 1.0, 1e-6, 1000, true);
19: 
20:   return 0;
21: }
22: 
23: double func( double x) { return x*x - 2; }
24: double dfunc(double x) { return 2*x; }
25: 
26: double newton_root(double (*f)(double), double (*df)(double), double x0,
27:                    double tol, int max_iters, bool print_iters) {
28:   double x  = x0;
29:   int iter  = 0;
30: 
31:   while (std::abs(f(x)) > tol && iter < max_iters) {
32:     if (print_iters) { 
33:       std::cout << "f(" << x << ") = " << f(x) << "\n";
34:     }
35: 
36:     // Newton's method update.
37:     x -= f(x) / df(x);
38:     iter++;
39:   }
40: 
41:   // One last print if necessary.
42:   if (print_iters) { 
43:     std::cout << "f(" << x << ") = " << f(x) << "\n";
44:   }
45: 
46:   return x;
47: }

f(1) = -1
f(1.5) = 0.25
f(1.41667) = 0.00694444
f(1.41422) = 6.0073e-06
f(1.41421) = 4.51061e-12

The midpoint rule is a numerical integration method which approximates the definite integral of a specified function over a specified interval using a specified number of subintervals where on each subinterval, the area under the curve is approximated by a rectangle whose width is the width of the subinterval and whose height is the height of the function at the midpoint between the points defining the end points of the subinterval.

Specifically, if \(n\) equally-sized subintervals are used on \([a,b]\), then the midpoint rule approximation, \(M_n\), to the definite integral of \(f(x)\) on \([a,b]\) is:

\[ \int_a^b f(x) \; dx \approx \sum_{i=1}^n f\left( \frac{x_{i-1}+x_i}{2} \right) \Delta x =: M_n \]

where \(\Delta x = \frac{b-a}{n}\), and \(x_i = a + i \cdot \Delta x, \;\; i=0, 1, \ldots, n\).

Let's consider \(f(x) = \frac{1}{x}\) on \([1, e]\).

 1: #include <iostream>
 2: #include <cmath>
 3: 
 4: const double E = std::exp(1.0);
 5: 
 6: // The function to be integrated.
 7: double func(double x);
 8: 
 9: // Compute the midpoint rule approximation to
10: // the definite integral of `f` from `a` to `b`
11: // using `n` subintervals.
12: double midpoint_rule(double (*f)(double), double a, double b, int n);
13: 
14: 
15: int main() {
16:   for(int n=2; n <= 20; n += 2) {
17:     std::cout << "n = " << n << ", "
18:               << "M_n = " << midpoint_rule(&func, 1, E, n) << "\n";
19:   }
20: 
21:   return 0;
22: }
23: 
24: double func(double x) { return 1.0 / x; }
25: 
26: double midpoint_rule(double (*f)(double), double a, double b, int n) {
27:   double xi;
28:   double xi_prev = a;
29:   double dx = (b-a) / n;
30:   double sum;
31: 
32:   for(int i=1; i <= n; i++) {
33:     xi = a + i*dx;
34:     sum += f(0.5*(xi_prev + xi));
35:     xi_prev = xi;
36:   }
37: 
38:   return sum*dx;
39: }

n = 2, M_n = 0.97636
n = 4, M_n = 0.993575
n = 6, M_n = 0.997091
n = 8, M_n = 0.998353
n = 10, M_n = 0.998942
n = 12, M_n = 0.999264
n = 14, M_n = 0.999459
n = 16, M_n = 0.999585
n = 18, M_n = 0.999672
n = 20, M_n = 0.999734

Create the header file, vector.h:

 1: #ifndef VECTOR_H
 2: #define VECTOR_H
 3: 
 4: namespace vec {
 5:   class vector {
 6:   public:
 7:     // Constructor. This function is called when the object is created.
 8:     vector(unsigned len);
 9: 
10:     // Destructor. This function is called when the object is destroyed.
11:     ~vector();
12: 
13:     // length accessor.
14:     unsigned len() const;
15: 
16:     // data accessor.
17:     double & element(unsigned i) const;
18: 
19:     // Simple print function.
20:     void print() const;
21: 
22:   private:
23:     unsigned length;
24:     double * data;
25:     void check_index(unsigned i) const;
26:   };
27: }
28: #endif

First note the macro guards, #ifndef, #define, and #endif, as well as the namespace, vec, wrapping the class declaration. Macro guards and namespaces were previously discussed when we initially introduced how to write header and implementation files.

Now, within the namespace, we've declared a class, vector, which contains public and private variables and function declarations. Private functions and variables may only be accessed through the public methods. This means if you created an instance of the class, vector, you would not be able to access the private variable directly. You could only call the public member-functions, which, in turn, may manipulate the private member-variables, or call the private member-functions. Splitting variables and functions into public and private helps to ensure that other libraries and programs use your class as intended.

Thus far, this class has 5 public member-functions, 2 private member-variables, and 1 private member-function. The first two member functions are special, the constructor and destructor, respectively. The constructor is called explicitly when you declare a new instance of this class, while the destructor is usually called implicitly when the object is deleted or when it goes out of scope.

Notice that the method for accessing elements of vector is called element and its return type is a reference (discussed previously). This allows us to use this function on both the left side of assignment operators, i.e., to assign values to vector components, and on the right side of assignments, i.e., to treat it as the value of the component.

Finally, notice that some member function declarations end with the keyword, const. Functions with such a signature are not allowed to modify member variables, and they are also not allowed to call other non const member functions.

Create the implementation file, vector.cpp:

 1: #include <iostream>
 2: #include <cstdlib>
 3: #include "vector.h"
 4: 
 5: namespace vec {
 6:   vector::vector(unsigned len) {
 7:     this->length = len;
 8:     this->data = new double[len];
 9:     // Initialize data to zeros.
10:     for(int i=0; i < this->len(); i++) { this->data[i] = 0.0; }
11:   }
12: 
13:   vector::~vector() {
14:     delete [] this->data;
15:   }
16: 
17:   unsigned vector::len() const {
18:     return this->length;
19:   }
20: 
21:   double & vector::element(unsigned i) const {
22:     #ifndef NDEBUG
23:     check_index(i);
24:     #endif
25:     return this->data[i];
26:   }
27: 
28:   void vector::print() const {
29:     for(int i=0; i < this->len(); i++) {
30:       std::cout << this->data[i] << '\n';
31:     }
32:   }
33: 
34:   void vector::check_index(unsigned i) const {
35:     if (i < 0 || i >= this->length) {
36:       std::cerr << "ERROR: index, " << i << ", is out-of-bounds.\n"
37:                 << "(valid indices are 0-" << (this->length-1) << ")\n";
38:       exit(1);
39:     }
40:   }
41: }

Note that we again wrap the implementations in the same namespace as wrapped by the class declaration. Also observe how each member-function is prepended by vector::.

The keyword, this, is a pointer to the calling object. Writing, this->, is equivalent to (*this)., and in fact, can be used for any pointer. Thus, this->length is equivalent to (*this).length.

Observe how the private member function, check_index, is used in the public element accessor function. If this library is compiled with the flag, -DNEDUBG, then the check function will not be called. You could read this flag as "define no debug". Thus, when this flag is present, the debug function check_index is called whenever the element accessor is called. The check_index function simply checks if the provided index is out-of-bounds for the vector. If it is, an informative message is printed, and the program terminates prematurely by calling exit(1). Such assertions with informative messages are a good practice, and will likely save you lots of headaches in the future.

Ok. Let's see some example usage, by creating a main.cpp, containing:

 1: #include <iostream>
 2: #include "vector.h"
 3: 
 4: int main() {
 5:     vec::vector v(5);
 6: 
 7:     std::cout << "`v` has length = " << v.len() << "\n";
 8: 
 9:     v.element(0) = -1.27;
10:     v.element(3) = 3.1;
11: 
12:     v.print();
13: 
14:     v.element(5) = 1234.0;
15: 
16:     return 0;
17: }

Let's first compile with our check_index debugger function in place:

g++ -c vector.cpp
g++ vector.o main.cpp
./a.out

If successful, you should see:

`v` has length = 5
-1.27
0
0
3.1
0
ERROR: index, 5, is out-of-bounds.
(valid indices are 0-4)

Now let's run without check_index:

g++ -DNDEBUG -c vector.cpp
g++ vector.o main.cpp
./a.out

Upon running, you will likely see some extensive list of errors when the element beyond the array's length is attempted to be accessed. Again, by liberally sprinkling these sorts of assertions through your code, you will (sometimes) find debugging much easier. After you're fairly certain that your code is working, you can simply compile with -DNEDUBG.

The v.element(i) accessor is a bit clunky. We can replace this with the more natural, v[i], by overloading the [] operator. Indeed, we can overload many of the normal C++ operators, e.g. +, -, =, etc. Some of these operators require more careful consideration when implementing class that utilizes dynamic memory allocation, such as our vector class.

It may be tempting at this point to attempt to initialize w from v directly:

1: vec::vector w = v;

If you attempt this currently, you will see all sorts of errors. This is because this type of intialization does not call the assignment operator. It calls the copy constructor. The assignment operator is only called when the object has already been intialized. Writing the previous line of code is essentially equivalent to

1: vec::vector w(v);

In other words, the constructor is called with the existing vector, v, as the argument, but we have not written a constructor yet with such a call signature.

The constructor can be overloaded, i.e., we can write multiple versions of the constructor function, and the one that matches the correct call signature will be used. This function overloading behavior actually applies to all functions in C++.

Let's add the copy constructor declaration to the header file:

1: // Copy constructor.
2: vector(const vector & src);

and let's add its implementation:

1: vector::vector(const vector & src) {
2:   this->length = src.len();
3:   this->data = new double[this->len()];
4: 
5:   // Copy over the data.
6:   for(int i=0; i < this->len(); i++) {
7:     this->data[i] = src[i];
8:   }
9: }

Now we compile and run something like:

1: vec::vector w = v;
2: w.print();

we will see:

-1.27
0
0
3.1
0

Non-member functions may access private member variables and private member functions of by labeling them as friend. This is useful in situations where it is not clear that function should be "called on" an object, i.e. object.method(params). Friend functions should be declared in the class declaration, but their implementation is not prepended with class::, which is necessary for member functions.

Let us overload operator* to implement matrix multiplication. It would also make sense to overload operator* with other call signatures to implement scalar multiplication or vector broadcasting, but for now we will just stick with matrix multiplication.

Add the following two non-friend, non-member functions. They should be defined just after the class definition, but they should remain in the body of namespace, mtx.

 1: template<class T>
 2: matrix<T> operator*(const matrix<T> & L, const matrix<T> & R) {
 3:     #ifndef NDEBUG
 4:     check_sizes_for_matmul(L, R);
 5:     #endif
 6: 
 7:     matrix<T> P(L.n_rows(), R.n_cols());
 8:     for(int i=0; i < L.n_rows(); i++) {
 9:         for(int j=0; j < R.n_cols(); j++) {
10:             // Inner product between row `i` of matrix, `L`
11:             // and column `j` of matrix, `R`.
12:             for(int k=0; k < L.n_cols(); k++) {
13:                 P(i,j) += L(i,k)*R(k,j);
14:             }
15:         }
16:     }
17: 
18:     return P;
19: }
20: 
21: template<class T>
22: void check_sizes_for_matmul(const matrix<T> & L, const matrix<T> & R) {
23:     if (L.n_cols() != R.n_rows()) {
24:         std::cerr << "Size mismatch for matrix multiplication.\n"
25:                   << "Left matrix has " << L.n_cols() << " cols, but\n"
26:                   << "right matrix has " << R.n_rows() << " rows.\n";
27:         exit(1);
28:     }
29: }

Note that because we've chosen to implement these functions as non-member and non-friend, we must prefix the function definitions with template<class T>. This isn't necessary when the functions were declared inside of the class body (i.e., member or friend functions) since the class declaration already begins with template<class T>.

We've added a simple matrix multiplcation function above and a debugger functions to ensure that if two matrices are multiplied, that it is well- defined to do so. Otherwise, we print an error message and abort, as with the element accessor function. We've also wrapped the debugger function with compiler macros, so that it can be left out if we choose to do so.

Here is some example usage:

 1: #include <iostream>
 2: #include "matrix.h"
 3: 
 4: int main() {
 5:     mtx::matrix<double> A(2,2);
 6:     mtx::matrix<double> B(2,3);
 7: 
 8:     A(0,0) = A(0,1) = A(1,0) = 1; A(1,1) = -1;
 9: 
10:     B(0,0) = B(1,0) = 1;
11:     B(0,1) = B(1,1) = -1;
12:     B(0,2) = B(1,2) = 3;
13: 
14: 
15:     mtx::matrix<double> C = A*B;
16: 
17:     for(int i=0; i < C.n_rows(); i++) {
18:         for(int j=0; j < C.n_cols(); j++) {
19:             std::cout << C(i,j);
20:             if (j < C.n_cols()-1)
21:                 std::cout << ", ";
22:         }
23:         std::cout << "\n";
24:     }
25: 
26:     B*A;
27: 
28:     return 0;
29: }

Output:

2, -2, 6
0, 0, 0
Size mismatch for matrix multiplication.
Left matrix has 3 cols, but
right matrix has 2 rows.