Problem Set #2¶

More on Visual Studio Code¶

If you are interested in making better use of Visual Studio Code, Microsoft has a number of tutorials on-line about how to use some of its various features:

• Visual Studio Code Documentation Home

• Using Clang in Visual Studio Code

• Integrate with External Tools via Tasks

I particularly encourage you to explore using Intellisense and how to compile, find errors, and debug from within VS Code.

Extract the homework files¶

We will be using files as data input and output mechanisms for the computational capabilities we will be developing. As is customary with large software libraries, we will factor the software into multiple files. Regardless of the size of a software project, we will want to start automating the compilation process using make. (As the course progresses we will also be introducing more advanced features of make.)

Create a new subdirectory named “ps2” to hold your files for this assignment. The files for this assignment are in a “tar ball” that can be downloaded with ps2.tar.gz.

To extract the files for this assignment, download the tar ball and copy it (or move it) to your ps2 subdirectory. From a shell program, change your working directory to be the ps2 subdirectory and execute the command

$ tar zxvf ps2.tar.gz

This will extract all of the files that are contained in ps2.tar.gz. (Tar files are very much like zip files or other files built by programs that bundle multiple files together into one. The name “tar” is originally derived from “tape archive”.)

The files that should be extracted from the tape archive file are

• Makefile

• Vector.hpp – contains the declaration of the Vector class as developed in lecture 4

• Matrix.hpp – contains the declaration of the Matrix class as developed in lecture 5

• amath583.hpp – contains declarations for the functions defined in amath583.cpp

• amath583.cpp – contains assorted operations on Vectors and Matrices

• catch.hpp – header file for the “Catch2” testing framework

• CImg.h - header file for CImg image library, used by julia.cpp

• julia.cpp – main file for julia.exe

• Questions.rst - text file for short-answer questions

• test_*.cpp – files to test the programs you will be writing

• verify.bash - bash script to build and run all test programs

Make¶

The Makefile provided to you has rules for building the various executables in this assignment. As mentioned in lecture, Makefiles contain these rules, but make also supports a simple macro language. We can use this language to achieve similar goals as with any programming language – to automate repetitive manual tasks, making them repeatable and consistent.

The ability to define and use macros is one of the most basic, but most powerful, aspects of make. Consider a rule such as one we presented in lecture:

amath583.o  : amath583.cpp amath583.hpp Vector.hpp
          c++ -c amath583.cpp -o amath583.o

This rule will invoke the compilation step with the c++ compiler to build amath583.o. But note that this rule is extremely specific. The compiler must be named c++. If there are other flags that are needed, such as an include path or library, we would have to add those:

amath583.o  : amath583.cpp amath583.hpp Vector.hpp
          c++ -I/usr/include/python3.8 -c amath583.cpp -o amath583.o -lpython3.8

Again, the paths would have to be exactly the ones specified. That would be fine for building just one program in one environment. But suppose, say, you wanted to be able to build multiple executables using the same Makefile. You would need to copy the same paths and compiler name for every executable. More problematically, suppose you want to make these same programs in different environments, where the compiler might have a different name, or the paths might have different names. In that case, you would have to manually edit each line to put in the right compiler name, right paths, etc — and you would have to do this every time you changed environments. And every time you do something manually, there is the opportunity to introduce errors. Suddenly, the automation and ease promised by make doesn’t seem like such a big win.

This is where abstractions comes to the rescue (as it always does) – in this case, in the form of macros. Rather than specifying what compiler we actually want to use, or what compiler flags we pass to it, we can write our rules using placeholders:

amath583.o  : amath583.cpp amath583.hpp Vector.hpp
        $(CXX) $(CXXFLAGS) -c amath583.cpp -o amath583.o

Make will substitute whatever we define CXX to be where it sees $(CXX) and substitute whatever we define CXXFLAGS to be where it sees $(CXXFLAGS). Now when we write rules using $(CXX) we are writing generic rules that can be made specific by how we set CXX et al.

Try the following experiment in your ps2 directory. First, execute

$ make clean

This is a rule defined at the bottom of the Makefile and deletes all of the generated files in your project (.o and .exe files) as well as other detritus that might have appeared during the compilation and editing process.

Now, issue the command

$ make -n julia.exe

This will process the rules for building julia.exe – but since we are passing the -n option to make – important – make will print out what would be done, without actually doing those things.

You should see something like

c++ -std=c++11 -g -O0 -Wall -c julia.cpp -o julia.o
c++ -std=c++11 -g -O0 -Wall -c amath583.cpp -o amath583.o
c++ -std=c++11 -g -O0 -Wall julia.o -o julia.exe

If you look at the Makefile you will see thes lines

CXX             := c++

OPTS            := -g -O0
LANG            := -std=c++11
PICKY           := -Wall

CXXFLAGS        := $(LANG) $(OPTS) $(PICKY)

These define the macros CXX, OPTS, etc. Note that CXXFLAGS is defined in terms of three other macros. The convention we will be using in this course is to separately define different classes of flags that we might pass to the c++ compiler and then bundle them all together into the CXXFLAGS macro. We use OPTS for setting optimization level (and debugging level), LANG to set the C++ language level, and PICKY to set how picky we want the compiler to be when processing our code. The options here indicate to use no optimization, c++11 as the language level, and to use maximum pickiness when compiling.

Do you recognize the values of each of these macros in the text that was printed when you invoked “make -n julia.exe”?

To see how the macros work, make a temporary change to the Makefile by changing LANG to be

LANG     := -std=c++17

Invoke make -n julia.exe again. Is there a difference between what was printed the first time and what was printed now?

Be aware that the macro processing is simple text processing, it has no actual meaning relative to what flags the compiler actually supports or what the compilation process actually is. We choose CXX and LANG and so forth as aids for our own comprehension of what we are trying to accomplish. For example, you could put an invalid flag into a macro and make would still invoke the command – and you would see the resulting error from the compiler when it is invoked.

Edit your Makefile again, this time adding a nonsense string to LANG e.g.,

LANG     := -std=c++17 asdfasdfasdf

What happens when you invoke

• make -n julia.exe

• make julia.exe

• make -k julia.exe

Lastly. If you look at the lines in the Makefile, there is still alot of repetition. Not repetition in the exact text, but repetition in the following sense. For the rules that build .o files, there is a rule for the .o, the corresponding .cpp dependency and a compilation step to turn that .cpp into a .o.

Every time you see repetition you should be thinking about how to create rules and process so that the repetitions can be made by the computer rather than you. In a future assignment we will look at how these kinds of patterns can also be handled by make.

Defensive Programming and Assertions¶

Maurice Wilkes was one of the founders of modern computing and, in some sense, of debugging. One of his most poignant quotes is:

It was on one of my journeys between the EDSAC room and the punching equipment that “hesitating at the angles of stairs” the realization came over me with full force that a good part of the remainder of my life was going to be spent in finding errors in my own programs.

Over the years, defensive programming has evolved as a set of techniques that can be used to help you find your own errors. One fundamental part of the process of defensive programming is to develop your program so that it supports two different modes of compilation and execution: debug mode and release mode. In debug mode, a program is first of all, compiled to enable use of a symbolic debugger (using the -g flag to the compiler). In addition, the program itself includes self-checks, or assertions, (inserted by the programmer) to insure that necessary program invariants hold.

#include <iostream>
#include <cassert>
using namespace std;

int main() {
  assert(1 == 1 && "This is true");       // will not be triggered
  assert(1 == 0 && "This is never true"); // will be triggered
  cout << "Hello World" << endl;
  return 0;
}

#define NDEBUG

#include <iostream>
#include <cassert>
using namespace std;

int main() {
  assert(1 == 0 && "This is never true"); // will be triggered
  cout << "Hello World" << endl;
  return 0;
}

#include <iostream>
#include <cassert>
using namespace std;

#define NDEBUG


int main() {
  assert(1 == 0 && "This is never true"); // will be triggered
  cout << "Hello World" << endl;
  return 0;
}

double sqrt583(const double& y) {
  assert(y >= 0);

  double x = 0.0, dx;
  do {
    dx = - (x*x-y) / (2.0*x);
    x += dx;
  } while (abs(dx) > 1.e-9);

  return x;
}

Assertions¶

In your subdirectory there is also a file assert3.cpp that has the same assertions as assert1.cpp and assert2.cpp – but without NDEBUG being defined in the program. Modify the rule for compiling assert3.exe – do not modify assert3.cpp – so that the assertion in assert3.exe will not be triggered.

Hello World!¶

Write a new hello world program – put it in a file called “hello_world.cpp”. Rather than saying hello to the world, I want you to simply print out a number indicating how far away you are from Lewis Hall when you did this assignment (in kilometers). You don’t need to compute anything (or even be exactly correct, you won’t be graded on what the number is) – you can just write the number directly in the file and print it (either as a number or as a string). We will collect the output from these programs to test the “leaderboard” feature of our autograder.

Zeroize¶

In the file amath583.cpp are five functions, three of which you need to complete and two of which are already written as examples. You will also need to add some functions from scratch. Every function that is defined in amath583.cpp also needs to be declared in amath583.hpp.

Complete the definition and declaration of zeroize() in amath583.hpp and amath583.cpp. This function should take a Vector as an argument and fill each entry of the Vector with a zero.

This function will be tested by linking amath583.o against test_zeroize.o to create test_zeroize.exe.

Vector Norm¶

The file amath583.cpp contains an implementation of one_norm as follows

double one_norm(const Vector& x) {
  double sum = 0.0;
  for (size_t i = 0; i < x.num_rows(); ++i) {
    sum += std::abs(x(i));
  }
  return sum;
}

The corresponding declaration in amath583.hpp is:

double one_norm(const Vector& x);

Complete the bodies of two_norm() and inf_norm() which respectively compute and return the two norm and infinity norm of the Vector that is passed in as a parameter.

These functions will be tested by linking them with test_norms.o which contains a series of unit tests that check different expected behavior of the norm computations.

Inner Product¶

Add a function called dot to amath583.cpp and amath583.hpp (definitions and declaration). The function should take two vectors as arguments and return their inner product (the sum of the pairwise products of each corresponding entry in the vectors).

The tests for this function are in test_dot.cpp, which should be compiled as with the other test programs into test_dot.exe.

Question.

1. How would you use dot to implement two_norm?

Matrix Norm¶

If we consider a typical vector space of N-tuples, e.g., \(R^N\), the set of all linear operators mapping \(R^N\) to \(R^M\) (say) is the set of all \(M \times N\) matrices. But, if we think about the kinds of things we can do with a matrix we can add them, scale them, etc – all the things we can do with objects of a vector space. And, in fact, \(M\times N\) matrices over a field, form a vector space. If the field is the real numbers, we write \(R^{M\times N}\).

Just as with other kinds of vectors, it is useful to define a norm over matrices. One simple option would be to just “unravel” the matrix into an \(M \times N\) vector and use one of our vector norms. That would be a norm for the matrix in some sense, but in linear algebra we are often more interested in how a matrix behaves as an operator and we would like for matrix norm to capture something about the operator-ness of the matrix. This kind of norm is called an induced norm and is defined as follows.

If \(\| \cdot \|\) is a norm on any n-dimensional vector space, e.g., \(R^n\), then the induced matrix norm of a matrix \(A\in R^{M\times N}\) is given by

For the one norm we then have

And for the infinity norm we have

The vector norms we looked at earlier are specific cases of \(p\)-norms, defined as

While it is straightforward to compute these for vectors using addition, subtraction, multiplication, etc., the only directly computable (meaning computable in a fixed number of steps) induced \(p\)-norms for a matrix are the 1-norm and the \(\infty\)-norm.

We can implement the 1-norm for the Matrix class we defined in lecture like this:

double one_norm(const Matrix& x) {
  double sum = 0.0;
  for (size_t j = 0; j < x.num_cols(); ++j) {  // for each column
   double tmp = 0.0;
   for (size_t i = 0; i < x.num_rows(); ++i) { // compute column sum
     tmp += std::abs(x(i, j));
   }
   sum = std::max(sum, tmp);
  }
  return sum;
}

You can apply the same formatting for all of your projects (or all files within a project) by saving the settings in a .clang-format file, either in the project directory or in your home directory. The settings I use are in this file, which you are welcome to use (and modify).

Remember that in C++, we can have functions with the same name as long as the types of the arguments differ (overloading), so we can have a one_norm for Vector and a one_norm for Matrix.

One seemingly obvious extension from the vector two norm to matrices is the following

However, this is not the induced matrix two-norm, rather it is called the “Frobenius” norm. Rather, the matrix two-norm is the largest singular value of the matrix \(A\). While this is computable (and we will be computing the singular value decomposition later in the course), it is not computable in a specified finite number of steps (though the iterations involved do converge quite rapidly). No other matrix \(p\)-norms are computable at all (they require exponential work as a function of size, see Matrix p-Norms Are NP-Hard to Approximate If $pneq1,2,infty$).

(583 ONLY) Implement the \(\infty\)-norm for a Matrix. You will need to add the definition to amath583.cpp and the declaration to amath583.hpp.

The tests for this function are in test_matrix_norm.cpp, which should be compiled as with the other test programs into test_matrix_norm.exe.

Tensors¶

To keep up with the times, we would like to add a three-dimensional Tensor class to AMATH 483/583. It is fairly straightforward to extend the type of interface that we have been using for Vector and Matrix:

class Tensor {
public:
  Tensor(size_t L, size_t M, size_t N) : num_slabs_(L), num_rows_(M), num_cols_(N),
                                         storage_(num_slabs_ * num_rows_ * num_cols_) {}

  double& operator()(size_t i, size_t j, size_t k) {
    return storage_[ /* something */ ];
  }
  const double& operator()(size_t i, size_t j, size_t k) const {
    return storage_[ /* something */ ];
  }

  size_t num_slabs() const { return num_slabs_; }
  size_t num_rows()  const { return num_rows_; }
  size_t num_cols()  const { return num_cols_; }

private:
  size_t              num_slabs_, num_rows_, num_cols_;
  std::vector<double> storage_;
};

Questions:

2. What is the formula that we would use for indexing if we wanted “slab-major” ordering? Hint: What is the size of a slab? It might be helpful to draw a picture of the data in 3D and see what the size is of a slab.

3. (Extra credit.) I claim that there are six distinct indexing formulas. Show that this claim is correct or that it is incorrect.

Julia Sets¶

Julia sets are the name given to sets that exhibit complex or “chaotic” or “fractal” behavior. Interestingly, Newton’s method (described in week 1) applied to problems with complex values can exhibit chaotic behavior even for very simple functions.

In fact, computing the cube-root of a complex number using Newton’s method (or, more generally, finding the complex root of a polynomial) can generate quite interesting sets. There are different ways of visualizing the complex behavior of Newton’s method in this case. One way is to simply count the number of iterations required for the algorithm to converge, and apply the method across a range of the complex plane.

The file julia.cpp is a program for doing a simple visualization in this way.

The C++ standard library includes a complex (template) class, with overloads for all of the standard arithmetic operations – it can be used just like a built-in numeric type.

The program should be fairly easy to understand. It iterates through the pixels in a 1024 by 1024 array (the screen) and converts each pixel location to a complex number. The x coordinate is the real part of the complex number and the y coordinate is the imaginary part. It then calls Newton’s method to solve \(f(x) = 0\), using the number we pass as the initial value. The program colors the pixel based on the number of iterations required for convergence. To make some interesting colors rather than just gray scale, we scale up the red, green, and blue channels by different factors.

The area of the complex plane that is covered by the algorithm is controlled by two parameters in the program: stretch and origin:

The default values of stretch and origin are 2.5 and (0, 0), respectively.

The part of the program that is not written are the bodies of f() and fp() (the function and its derivative). In class we constructed these for computing the square root.

For this exercise you just need to fill in those two functions for computing a root of the polynomial

The function f() should compute \(f(x)\) as above, while fp() should compute its derivative.

The test for this program will compile and run it, verifying that the output image is produced.

Again, although we are using complex numbers, which have a real and imaginary parts, we can express the complex arithmetic with the complex numbers – we don’t have to work directly with the real and imaginary parts.

If your functions are correctly defined, when the program runs it will create an image file called “julia.bmp” that looks like the following when visualized.

I encourage you to experiment with some different polynomials (or even trigonometric functions) and see what kinds of images you can generate. You can also experiment with different values of stretch and origin to view different regions of the complex plane, as well as different scalings of the color planes.

If you have some particularly interesting results, please share them and I can create a small gallery on the course web site.