Discovering Modern C++. An Intensive Course for Scientists, Engineers, and Programmers(2016)

---

Speaking of abstractions, the CG implementation in Listing A–2 does not commit to any technical detail. In no place is the function restricted to a numerical type like double. It works as well for float, GNU’s multiprecision, complex, interval arithmetic, quaternions, . . .

The matrix A can be stored in any internal format, as long as A can be multiplied with a vector. In fact, it does not even need to be a matrix but can be any linear operator. For instance, an object that performs a Fast Fourier Transformation (FFT) on a vector can be used as A when the FFT is expressed by a product of A with the vector. Similarly, the vectors do not need to be represented by finite-dimensional arrays but can be elements of any vector space that is somehow computer-representable as long as all operations in the algorithm can be performed.

We are also open to other computer architectures. If the matrix and the vectors are distributed over the nodes of a parallel supercomputer and corresponding parallel operations are available, the function runs in parallel without changing any single line. (GP)GPU acceleration can also be realized without changing the algorithm’s implementation. In general, any existing or new platform for which we can implement matrix and vector types and their corresponding operations is supported by our Generic conjugate gradient function.

As a consequence, sophisticated scientific applications of several thousand lines on top of such abstractions can be ported to new platforms without code modification.

A.2 Basics in Detail

This section accumulates the outsourced details from Chapter 1.

A.2.1 More about Qualifying Literals

Here we want to expand the examples for literal qualifications from Section 1.2.2 a bit.

Availability: The standard library provides a (template) type for complex numbers where the type for the real and imaginary parts can be parameterized by the user:

std::complex<float> z(1.3, 2.4), z2;

These complex numbers of course provide the common operations like addition and multiplication. For unknown reasons, these operations are only provided between the type itself and the underlying real type. Thus, when we write:

z2= 2 * z; // Error: no int * complex<float>
z2= 2.0 * z; // Error: no double * complex<float>

we will get an error message that the multiplication is not available. More specifically, the compiler will tell us that there is no operator*() for int and std::complex<float> respectively for double and std::complex<float>.¹ To get the simple expression above compiled, we must make sure that our value 2 has the type float:

1. Since the multiplication is implemented as a template function, the compiler does not convert the int and double to float.

z2= 2.0f * z;

Ambiguity: At some point in this book, we introduced function overloading: a function with different implementations for different argument types (§1.5.4). The compiler selects that function overload which fits best. Sometimes the best fit is not clear, for instance, if function f accepts anunsigned or a pointer and we call:

f(0);

The literal 0 is considered an int and can be implicitly converted into unsigned or any pointer type. None of the conversions is prioritized. As before, we can address the issue by a literal of the desired type:

f(0u);

Accuracy: The accuracy issue comes up when we work with long double. On the author’s computer, the format can handle at least 19 digits. Let us define 1/3 with 20 digits and print out 19 of them:

long double third= 0.3333333333333333333;
cout.precision(19);
cout "One third is " third ".\n";

The result is

One third is 0.3333333333333333148.

The program behavior is more satisfying if we append an l to the number:

long double third= 0.3333333333333333333l;

yielding the printout that we hoped for:

One third is 0.3333333333333333333.

A.2.2 static Variables

In contrast to the scoped variables from Section 1.2.4 which die at the end of the scope, static variables live till the end of the program. Thus, declaring a local variable as static has an effect only when the contained block is executed more than once: in a loop or a function. Within a function, we could implement a counter for how often the function was called:

void self_observing_function()
{
static int counter= 0; // only executed once
++counter;
cout "I was called " counter " times.\n";
...
}

To reuse the static data, the initialization is only performed once. The major motivation for static variables is to reuse helper data like lookup tables or caches in the next function call. However, if the management of the helper data reaches a certain complexity, a class-based solution (Chapter 2) probably leads to a cleaner design.

The effect of the keyword static depends on the context but the common denominator is:

• Persistence: data of static variables remains in memory for the rest of the program execution.

• File scope: static variables and functions are only visible in the compilation of the actual file and do not collide when multiple compiled programs are linked together. This is illustrated in Section 7.2.3.2.

Thus, the impact on global variables is to limit their visibility as they already live until the program ends. Conversely, the effect on local variables is the lifetime extension since their visibility is already limited. There is also a class-related meaning of static which is discussed in Chapter 2.

A.2.3 More about if

The condition of an if must be a bool expression (or something convertible to bool). Thus, the following is allowed:

int i= ...
if (i) // bad style
do_something();

This relies on the implicit conversion from int to bool. In other words, we test whether i is different from 0. It is clearer to say this instead:

if (i != 0) // better
do_something();

An if-statement can contain other if-statements:

if (weight > 100.0) {
if (weight > 200.0) {
cout "This is extremely heavy.\n";
} else {
cout "This is quite heavy.\n";
}
} else {
if (weight < 50.0) {
cout "A child can carry this.\n";
} else {
cout "I can carry this.\n";
}
}

In the above example, the parentheses could be omitted without changing the behavior but it is clearer to have them. The example is more readable when we reorganize the nesting:

if (weight < 50.0) {
cout "A child can carry this.\n";
} else if (weight <= 100.0) {
cout "I can carry this.\n";
} else if (weight <= 200.0) {
cout "This is quite heavy.\n";
} else {
cout "This is extremely heavy.\n";
}

The parentheses can be omitted here as well, and it requires less effort to figure out what is going on.

At the end of our discussion on if-then-else, we want to do something sophisticated. Let us take the then-branch of the second-last example without braces:

if (weight > 100.0)
if (weight > 200.0)
cout "This is extremely heavy.\n";
else
cout "This is quite heavy.\n";

It looks like the last line is executed when weight is between 100 and 200 assuming the first if has no else-branch. But we could also assume the second if comes without an else-branch and the last line is executed when weight is less than or equal to 100. Fortunately, the C++ standard specifies that an else-branch always belongs to the innermost possible if. So, we can count on our first interpretation. In case the else-branch belongs to the first if we need braces:

if (weight > 100.0) {
if (weight > 200.0)
cout "This is extremely heavy.\n";
} else
cout "This is not so heavy.\n";

Maybe these examples will convince you that it is more productive to set more braces and save the time guessing what the branches belong to.

A.2.4 Duff’s Device

The continued execution in switch cases allows us also to implement short loops without the termination test after each iteration. Say we have vectors with dimension ≤ 5. Then we could implement a vector addition without a loop:

assert(size(v) <= 5);
int i= 0;
switch (size(v)) {
case 5: v[i]= w[i] + x[i]; ++i; // keep going
case 4: v[i]= w[i] + x[i]; ++i; // keep going
case 3: v[i]= w[i] + x[i]; ++i; // keep going
case 2: v[i]= w[i] + x[i]; ++i; // keep going
case 1: v[i]= w[i] + x[i]; // keep going
case 0: ;
}

This technique is called Duff’s device. It is usually not used stand-alone as above (but as a cleanup of unrolled loops). Such techniques should not be used in the main development of projects but only as final tuning of performance-critical kernels.

A.2.5 More about main

Arguments containing spaces must be quoted. In the first argument, we can also see when the program is called with path information; e.g.:

../c++11/argc_argv_test first "second third "fourth

prints out:

../c++11/argc_argv_test
first
second third
fourth

Some compilers also support a vector of strings as main arguments. This is more convenient but not portable.

For calculating with the arguments, they must be converted first:

cout argv[1] " times " argv[2] " is "
stof(argv[1]) * stof(argv[2]) ".\n";

which could provide such impressive knowledge to us:

argc_argv_test 3.9 2.8
3.9 times 2.8 is 10.92.

Unfortunately, the conversion from strings to numbers does not tell us when the complete string is not convertible. As long as the string starts with a number or plus or minus, the reading is stopped if a character is found that cannot belong to a number and the sub-string read up to this point is converted to a number.

In Unix-like systems, the exit code of the last command can be accessed by $? in the shell. We can also use the exit code to execute multiple commands in a row under the condition that the preceding succeeded:

do_this && do_that && finish_it

In contrast to C and C++, the command shell interprets an exit code of 0 as true in the sense of okay. However, the handling of && is similar to that of C and C++: only when the first sub-expression is true do we need to evaluate the second. Likewise, a command is only executed when the preceding one was successful. Dually, || can be used for error handling because the command after an || is only performed when the preceding one failed.

A.2.6 Assertion or Exception?

Without going into detail, exceptions are more expensive than assertions since C++ cleans up the run-time environment when an exception is thrown. Turning off exception handling can accelerate applications perceivably. Assertions, on the other hand, immediately kill the program and do not require the cleanup. In addition, they are usually disabled in release mode anyway.

As we said before, unexpected or inconsistent values that originate from programming errors should be handled with assertions and exceptional states with exceptions. Unfortunately, this distinction is not always obvious when we encounter a problem. Consider our example of a file that cannot be opened. The reason can be either that an incorrect name was typed by a user or found in a configuration file. Then an exception would be best. The wrong file can also be a literal in the sources or the result of an incorrect string concatenation. Such program errors cannot be handled and we prefer the program to terminate with an assertion.

This dilemma is a conflict between avoiding redundancies and immediate sanity checks. At the location where the file name is entered or composed, we do not know whether it is a programming error or erroneous input. Implementing the error handling at these points could require repeating the opening test many times. This causes extra programming effort for repeated check implementation and bears the danger that the tests are not consistent with each other. Thus, it is more productive and less error-prone to test only once and not to know what caused our current trouble. In this case, we shall be prudent and throw an exception so that a cure is possible in some situations at least.

Corrupt data is usually better handled by an exception. Assume the salaries of your company are computed and the data of a newbie is not yet fully set. Raising an assertion would mean that the entire company (including you) would not be paid that month or at least until the problematic data set is fixed. If an exception is thrown during the data evaluation, the application could report the error in some way and continue for the remaining employees.

Speaking of program robustness, functions in universally used libraries should never abort. If the function is used, for instance, to implement an autopilot, we would rather turn off the autopilot than have the entire program terminated. In other words, when we do not know all application domains of a library, we cannot tell the consequences of a program abort.

Sometimes, the cause of a problem is not 100 % certain in theory but sufficiently clear in practice. The access operator of a vector or matrix should check whether the index is in a valid range. In principle, an out-of-range index could originate from user input or a configuration file, but in practically all cases it comes from a program bug. Raising an assertion is appropriate here. To comply with the robustness issue, it might be necessary to allow the user to decide between assertion and exception with compile flags (see §1.9.2.3).

A.2.7 Binary I/O

Conversions from and to strings can be quite expensive. Therefore, it is often more efficient to write data directly in their respective binary representations into files. Nonetheless, before doing so it is advisable to check with tools whether file I/O is really a significant bottleneck of the application.

When we decide to go binary, we should set the flag std::ios::binary for impeding implicit conversions like adapting line breaks to Windows, Unix, or Mac OS. It is not true that the flag makes the difference between text and binary files: binary data can be written without the flag and text with it. However, to prevent the before-mentioned surprises it is better to set the flag appropriately.

The binary output is performed with the ostream’s member function write and input with istream::read. The functions take a char address and a size as arguments. Thus, all other types have to be casted to pointers of char:

int main (int argc, char* argv[])
{
std::ofstream outfile;
with_io_exceptions(outfile);
outfile.open("fb.txt", ios::binary);

double o1= 5.2, o2= 6.2;
outfile.write(reinterpret_cast<const char *>(&o1), sizeof(o1));
outfile.write(reinterpret_cast<const char *>(&o2), sizeof(o2));
outfile.close();

std::ifstream infile;
with_io_exceptions(infile);
infile.open("fb.txt", ios::binary);

double i1, i2;
infile.read(reinterpret_cast<char *>(&i1), sizeof(i1));
infile.read(reinterpret_cast<char *>(&i2), sizeof(i2));
std::cout "i1 = " i1 ", i2 = " i2 "\n";
}

An advantage of binary I/O is that we do not need to worry about how the stream is parsed. On the other hand, non-matching types in the read and write commands result in completely unusable data. In particular, we must be utterly careful when the files are not read on the same platform as they were created: a long variable can contain 32 bits on one machine and 64 bits on another. For this purpose, library <cstdint> provides what sizes are identical on each platform. The type int32_t, for instance, is a 32-bit signed int on every platform. uint32_t is the corresponding unsigned.

The binary I/O works in the same manner for classes when they are self-contained, that is, when all data is stored within the object as opposed to external data referred to by pointers or references. Writing structures containing memory addresses—like trees or graphs—to files requires special representations since the addresses are obviously invalid in a new program run. In Section A.6.4, we will show a convenience function that allows us to write or read many objects in a single call.

A.2.8 C-Style I/O

The old-style I/O from C is also available in C++:

#include <cstdio>

int main ()
{
double x= 3.6;
printf("The square of %f is %f\n", x, x*x);
}

The command printf stands not so surprisingly for print-formatted. Its corresponding input is scanf. File I/O is realized with fprintf and fscanf.

The advantage of these functions is that the formatting is quite compact; printing the first number in 6 characters with 2 decimal places and the second number in 14 characters with 9 decimal places is expressed by the following format string:

printf("The square of %6.2f is %14.9f\n", x, x*x);

The problem with the format strings is that they are not Type-Safe. If the argument does not match its format, strange things can happen, e.g.:

int i= 7;
printf("i is %s\n", i);

The argument is an int and will be printed as a C string. A C string is passed by a pointer to its first character. Thus, the 7 is interpreted as an address and in most cases the program crashes. More recent compilers check the format string if it is provided as a literal in the printf call. But the string can be set before:

int i= 7;
char s[]= "i is %s\n";
printf(s, i);

or be the result of string operations. In this case, the compiler will not warn us.

Another disadvantage is that it cannot be extended for user types. C-style I/O can be convenient in log-based debugging, but streams are much less prone to errors and should be preferred in high-quality production software.

A.2.9 Garbarge Collection

Garbage Collection (GC) is understood as the automatic release of unused memory. Several languages (e.g., Java) from time to time discharge memory that is not referred to in the program any longer. Memory handling in C++ is designed more explicitly: the programmer controls in one way or another when memory is released. Nonetheless, there is interest among C++ programmers in GC either to make software more reliable—especially when old leaky components are contained which nobody is willing or able to fix—or when interoperability to other languages with managed memory handling is desired. An example of the latter is Managed C++ in .NET.

Given that interest, the standard (starting with C++11) defines an interface for garbage collectors. However, garbage collection is not a mandatory feature, and so far we do not know a compiler that supports it. In turn, applications relying on GC do not run with common compilers. Garbage collection should only be the last resort. Memory management should be primarily encapsulated in classes and tight to their creation and destruction, i.e., RAII (§2.4.2.1). If this is not feasible, one should consider unique_ptr (§1.8.3.1) and shared_ptr (§1.8.3.2) which automatically release memory that is not referred to. Actually, the reference counting by a shared_ptr is already a simple form of GC (although probably not everybody will agree with this definition). Only when all those techniques are non-viable due to some form of tricky cyclic dependencies, and portability is not an issue, should one resort to garbage collection.

A.2.10 Trouble with Macros

For instance, a function with the signature

double compute_something(double fnm1, double scr1*}, double scr2)

compiles on most compilers but yields a bizarre error message on Visual Studio. The reason is that scr1 is a macro defining a hexadecimal number which is substituted for our second argument name. Obviously, this is not legal C++ code any longer, but the error will contain our original source before the substitution took place. Thus, we cannot see anything suspicious. The only way to resolve such problems is to run the preprocessor and check the expanded sources, e.g.:

g++ my_code.cpp -E -o my_code.ii.cpp

This might take some effort since the expanded version contains the sources of all directly and indirectly included files. Eventually, we will find what happened with the source code line rejected by the compiler:

double compute_something(double fnm1, double 0x0490, double scr2)

We can fix this simply by changing the argument name, once we know that it is used as a macro somewhere.

Constants used in computations should not be defined by macros:

#define pi 3.1415926535897932384626433832795028841 // Do Not!!!

but as true constants:

const long double pi= 3.1415926535897932384626433832795028841L;

Otherwise, we create exactly the kind of trouble that scr1 caused in our preceding example. In C++11, we can also use constexpr to ascertain that the value is available during compilation, and C++14 also offers template constants.

Function-like macros offer yet a different spectrum of funny traps. The main problem is that macro arguments only behave like function arguments in simple use cases. For instance, when we write a macro max_square:

#define max_square(x, y) x*x >= y*y ? x*x : y*y

the implementation of the expression looks straightforward and we probably will not expect any trouble. But we will get into some when we use it with a sum or difference:

int max_result= max_square(a+b, a-b);

Then it will be evaluated as

int max_result= a+b*a+b >= a-b*a-b ? a+b*a+b : a-b*a-b;

yielding evidently wrong results. This problem can be fixed with some parentheses:

#define max_square(x, y) ((x)*(x) >= (y)*(y) ? (x)*(x) : (y)*(y))

To protect against high-priority operators, we also surrounded the entire expression with a pair of parentheses. Thus, macro expressions need parentheses around each argument and around the entire expression. Note that this is a necessary, not a sufficient, condition for correctness.

Another serious problem is the replication of arguments in expressions. If we call max_square as in

int max_result= max_square(++a, ++b);

the variables a and b are incremented four times.

Macros are a quite simple language feature, but implementing and using them is more complicated and dangerous than it seems at first glance. The hazard is that they interact with the entire program. Therefore, new software should not contain macros at all. Unfortunately, existing software already contains many of them and we have to deal with them.

Regrettably, there is no general cure for trouble with macros. Here are some tips that work in most cases:

• Avoid names of popular macros. Most prominently, assert is a macro in the standard library and giving that name to a function is asking for trouble.

• Un-define the macro with #undef.

• Include macro-heavy libraries after all others. Then the macros still pollute your application but not the other included header files.

• Impressively, some libraries offer protection against their own macros: one can define a macro that disables or renames macros’ dangerously short names.²

2. We have seen a library that defined a single underscore as a macro which created a lot of problems.

A.3 Real-World Example: Matrix Inversion

“The difference between theory and practice is larger in practice than in theory.”

—Tilmar König

To round up the basic features, we apply them to demonstrating how we can easily create new functionality. We want to give you an impression of how our ideas can evolve naturally into reliable and efficient C++ programs. Particular attention is paid to clarity and reusability. Our programs should be well-structured inside and intuitive to use from outside.

To simplify the realization of our study case, we use (the open-source part of) the author’s library Matrix Template Library 4—see http://www.mtl4.org. It already contains a lot of linear-algebra functionality that we need here.³ We hope that future C++ standards will provide such functionality out of the box. Maybe some of the book’s readers will contribute to it.

3. It actually already contains the inversion function inv we are going for, but let us pretend it does not.

As a software development approach, we will use a principle from Extreme Programming: writing tests first and implementing the functionality afterward—known as Test-Driven Development (TDD). This has two significant advantages:

• It protects us as programmers (to some extent) from Featurism—the obsession to add more and more features instead of finishing one thing after another. If we write down what we want to achieve, we work more directly toward this goal and usually accomplish it much earlier. When writing the function call, we already specify the interface of the function that we plan to implement. When we set up expected values for the tests, we say something about the semantics of our function. Thus, tests are compilable documentation. The tests might not tell everything about the functions and classes we are going to implement, but what it says does it very precisely. Documentation in text can be much more detailed and comprehensible but also much vaguer than tests (and procrastinated eternally).

• If we start writing tests after we finally finish the implementation—say, on a late Friday afternoon—we do not want to see it fail. We will write the test with nice data (whatever this means for the program in question) and minimize the risk of failure. Or we might decide to go home and swear that we will test it on Monday.

For those reasons, we will be more honest if we write our tests first. Of course, we can modify our tests later if we realize that something does not work out as we imagined or we evolved the interface design out of the experience we gained so far. Or maybe we just want to test more. It goes without saying that verifying partial implementations requires commenting out parts of our test—temporarily.

Before we start implementing our inverse function and its tests, we have to choose an algorithm. We can choose among different direct solvers: determinants of sub-matrices, block algorithms, Gauß-Jordan, or LU decomposition with or without pivoting. Let’s say we prefer LU factorization with column pivoting so that we have

LU = PA,

with a unit lower triangular matrix L, an upper triangular matrix U, and a permutation matrix P. Thus it is

A = P^–1LU

and

We use the LU factorization from MTL4, implement the inversion of the lower and upper triangular matrices, and compose it appropriately.

c++11/inverse.cpp

Now, we start with our test by defining an invertible matrix and printing it out:

int main(int argc, char* argv[])
{
const unsigned size= 3;
using Matrix= mtl::dense2D<double>; // type from MTL4
Matrix A(size, size);
A= 4, 1, 2,
1, 5, 3,
2, 6, 9;

cout "A is:\n" A;

For later abstraction, we define the type Matrix and the constant size. Using C++11, we could set up the matrix with uniform initialization:

Matrix A= {{4, 1, 2}, {1, 5, 3}, {2, 6, 9}};

However, for our implementation we get along with C++03.

The LU factorization in MTL4 is performed in place. So as not to alter our original matrix, we copy it into a new one.

Matrix LU(A);

We also define a vector for the permutation computed in the factorization:

mtl::dense_vector<unsigned> Pv(size);

These are the two arguments for the LU factorization:

lu(LU, Pv);

For our purpose, it is more convenient to represent the permutation as a matrix:

Matrix P(permutation(Pv));
cout "Permutation vector is " Pv "\nPermutation matrix is\n" P;

This allows us to express the row permutation by a matrix product:⁴

4. You might wonder why a matrix is created from the product P * A although the product is already a matrix. Well, it is not, technically. For efficiency reasons, it is an expression template (see Section 5.3) that evaluates the product in certain expressions. Output is not one of those expressions—at least not in the public version.

cout "Permuted A is \n" Matrix(P * A);

We now define an identity matrix of appropriate size and extract L and U from our in-place factorization:

Matrix I(matrix::identity(size, size)), L(I + strict_lower(LU)),
U(upper(LU));

Note that the unit diagonal of L is not stored and needs to be added for the test. It could also be treated implicitly, but we refrain from it for the sake of simplicity. We have now finished the preliminaries and come to our first test. Once we have computed the inverse of U, named UI, their product must be the identity matrix (approximately). Likewise for the inverse of L:

constexpr double eps= 0.1;

Matrix UI(inverse_upper(U));
cout "inverse(U) [permuted] is:\n" UI "UI * U
is:\n" Matrix(UI * U);
assert(one_norm(Matrix(UI * U - I)) < eps);

Testing results of non-trivial numeric calculation for equality is quite certain to fail. Therefore, we used the norm of the matrix difference as the criterion. Likewise, the inversion of L (with a different function) is tested.

Matrix LI(inverse_lower(L));
cout "inverse(L) [permuted] is:\n" LI "LI * L
is:\n" Matrix(LI * L);
assert(one_norm(Matrix(LI * L - I)) < eps);

This enables us to calculate the inverse of A itself and test its correctness:

Matrix AI(UI * LI * P);
cout "inverse(A) [UI * LI * P] is \n" AI "A * AI
is \n" Matrix(AI * A);
assert(one_norm(Matrix(AI * A - I)) < eps);

We also check our inverse function with the same criterion:

Matrix A_inverse(inverse(A));
cout "inverse(A) is \n" A_inverse "A * AI
is \n" Matrix(A_inverse * A);
assert(one_norm(Matrix(A_inverse * A - I)) < eps);

After establishing tests for all components of our calculation, we start with their implementations.

The first function that we program is the inversion of an upper triangular matrix. This function takes a dense matrix as an argument and returns a dense matrix as well:

Matrix inverse_upper(const Matrix & A) {

}

Since we do not need another copy of the input matrix, we pass it as a reference. The argument should not be changed so we can pass it as const. The constancy has several advantages:

• We improve the reliability of our program. Arguments passed as const are guaranteed not to change; if we accidentally modify them, the compiler will tell us and abort the compilation. There is a way to remove the constancy, but this should only be used as a last resort, e.g., for interfacing to obsolete libraries written by others. Everything you write yourself can be realized without eliminating the constancy of arguments.

• Compilers can better optimize when the objects are guaranteed not to be altered.

• In case of references, the function can be called with expressions. Non-const references are required to store the expression into a variable and pass the variable to the function.

Another comment: people might tell you that it is too expensive to return containers as results and it is more efficient to use references. This is true—in principle. For the moment, we accept this extra cost and pay more attention to clarity and convenience. Later in this book, we will introduce techniques for how to minimize the cost of returning containers from functions.

So much for the function signature; let us now turn our attention to the function body. The first thing we do is verify that our argument is valid. Obviously the matrix must be square:

const unsigned n= num_rows(A);
if (num_cols(A) != n)
throw "Matrix must be square";

The number of rows is needed several times in this function and is therefore stored in a variable, well, constant. Another prerequisite is that the matrix has no zero entries in the diagonal. We leave this test to the triangular solver.

Speaking of which, we can get our inverse triangular matrix with a triangular solver of a linear system, which we find in MTL4; more precisely, the k-th vector of U^–1 is the solution of

Ux = e_k

where e_k is the k-th unit vector. First we define a temporary variable for the result:

Matrix Inv(n, n);

Then we iterate over the columns of Inv:

for (unsigned k= 0; k < n; ++k) {

}

In each iteration we need the k-th unit vector:

dense_vector <double> e_k(n);
for (unsigned i= 0; i < n; ++i)
if (i == k)
e_k[i]= 1.0;
else
e_k[i]= 0.0;

The triangular solver returns a column vector. We could assign the entries of this vector directly to entries of the target matrix:

for (unsigned i= 0; i < n; ++i)
Inv[i][k]= upper_trisolve(A, e_k)[i];

This is nice and short but we would compute upper_trisolve n times! Although we said that performance is not our primary goal at this point, the increase of overall complexity from order 3 to 4 is too wasteful of resources. Many programmers make the mistake of optimizing too early, but this does not mean that we should accept (without serious reasons) an implementation with a higher order of complexity. To avoid the superfluous recomputations, we store the triangle solver result and copy the entries from there.

dense_vector<double> res_k(n);
res_k= upper_trisolve(A, e_k);

for (unsigned i= 0; i < n; ++i)
Inv[i][k]= res_k[i];

Finally, we return the temporary matrix. The function in its complete form looks as follows:

Matrix inverse_upper(Matrix const& A)
{
const unsigned n= num_rows(A);
assert(num_cols(A) == n); // Matrix must be square

Matrix Inv(n, n);

for (unsigned k= 0; k < n; ++k) {
dense_vector<double> e_k(n);
for (unsigned i= 0; i < n; ++i)
if (i == k)
e_k[i]= 1.0;
else
e_k[i]= 0.0;

dense_vector<double> res_k(n);
res_k= upper_trisolve(A, e_k);

for (unsigned i= 0; i < n; ++i)
Inv[i][k]= res_k[i];
}
return Inv;
}

Now that the function is complete, we first run our test. Evidently, we have to comment out part of the test because we have only implemented one function so far. Nonetheless, it is good to know early whether this first function behaves as expected. It does, and we could be happy with it now and turn our attention to the next task; there are still many. But we will not.

Well, at least we can be glad to have a function running correctly. Nevertheless, it is still worth spending some time to improve it. Such improvements are called Refactoring. Experience has shown that it takes much less time to refactor immediately after implementation than later when bugs are discovered or the software is ported to other platforms. Obviously, it is much easier to simplify and structure our software immediately when we still know what is going on than after some weeks/months/years.

The first thing we might dislike is that something so simple as the initialization of a unit vector takes five lines. This is rather verbose:

for (unsigned i= 0; i < n; ++i)
if (i == k)
e_k[i]= 1.0;
else
e_k[i]= 0.0;

We can write this more compactly with the conditional operator:

for (unsigned i= 0; i < n; ++i)
e_k[i]= i == k ? 1.0 : 0.0;

The conditional operator ?: usually takes some time to get used to but it results in a more concise representation.

Although we have not changed anything semantically in the program and it seems obvious that the result will still be the same, it cannot hurt to run our test again. You will see that often you are sure that your program modifications could never change the behavior—but still they do. The sooner we find an unexpected behavior, the less work it is to fix it. With the test(s) that we have already written, it only takes a few seconds and makes us feel more confident.

If we would like to be really cool, we could explore some insider knowhow. The expression i == k returns a boolean, and we know that bool can be converted implicitly into int and subsequently to double. In this conversion, false turns into 0 and true into 1 according to the standard. These are precisely the values that we want as double:

e_k[i]= static_cast<double>(i == k);

In fact, the conversion from int to double is performed implicitly and can be omitted as well:

e_k[i]= i == k;

As cute as this looks, it is some stretch to assign a logical value to a floating-point number. It is well defined by the implicit conversion chain bool → int → double, but it will confuse potential readers and you might end up explaining to them what is happening on a mailing list or adding a comment to the program. In both cases you end up writing more text for the explication than you saved in the program.

Another thought that might occur to us is that it is probably not the last time that we need a unit vector. So, why not write a function for it?

dense_vector<double> unit_vector(unsigned k, unsigned n)
{
dense_vector<double> e_k(n, 0.0);
e_k[k]= 1;
return e_k;
}

Since the function returns the unit vector, we can just take it as an argument for the triangular solver:

res_k= upper_trisolve(A, unit_vector(k, n));

For a dense matrix, MTL4 allows us to access a matrix column as a column vector (instead of a sub-matrix). Then, we can assign the result vector directly without a loop:

Inv[irange(0, n)][k]= res_k;

As a short explanation, the bracket operator is implemented in such a manner that integer indices for rows and columns return the matrix entry while ranges for rows and columns return a sub-matrix. Likewise, a range of rows and a single column gives us a column of the corresponding matrix—or part of this column. Conversely, a row vector can be extracted from a matrix with an integer as a row index and a range for the columns.

This is an interesting example of how to deal with the limitations as well as the possibilities of C++. Other languages have ranges as part of their intrinsic notation; e.g., Python has a symbol: for expressing ranges of indices. C++ does not provide such a symbol, but we can introduce a new type—like MTL4’s irange—and define the behavior of operator[] for this type. This leads to an extremely powerful mechanism!