Hints to Exercises: C++
Problem Set 1

Computational Physics with Numerical Recipes
James P. Sethna
Physics 480/680, Spring 2008
http://www.physics.cornell.edu/~sethna/teaching/ComputationalPhysics/

1.1 Preliminaries

I recommend setting up Makefiles and using xmgrace for plots. 

1.1 You'll need to #include <fstream> and <cmath>, and probably <iostream>
for debugging (writing things to cout). You probably will want to use the 
standard namespace (unless you are a purist). So, the beginning of your 
C++ program should look like

	#include <iostream>
	#include <fstream>
	#include <cmath>

	using namespace std;

	int main ()
	{
	    ifstream infile;
	    ofstream outfile;
	...
	}

Your Makefile should look something like

	CPPFLAGS= -g
	# CPPFLAGS= -O3

	Preliminaries: Preliminaries.o
        	$(CXX) $(CPPFLAGS) -o Preliminaries Preliminaries.o

	Preliminaries.o: Preliminaries.cpp
        	$(CXX) $(CPPFLAGS) -c Preliminaries.cpp

Notice that tabs are not the same as spaces in Makefiles. The two lines
starting $(CXX) should have a tab indent; the other lines should not.
CPPFLAGS= -g is for debugging; CPPFLAGS = -O3 is a sometimes risky level of
optimization -- if your answers change, try -O2. Your programming session
should look like

	[Edit Preliminaries.cpp]
	make Preliminaries
	[Debug for compiler errors, repeat]
	Preliminaries
	[Debug for run-time errors, repeat]
	xmgrace OutputFileName.dat
	[Debug, repeat, print or save figures]

(a) You'll need to use outfile.open("OutputFileName.dat") and outfile.close().
M_PI is the name for pi. You can write x and y out to outfile with
	outfile << x << " " << y << endl;
unless you aren't using namespace std, in which case you need to say std::endl.

(b) For reading in NoisyYOfT.dat (notice the four capital letters), you may find
that a while loop checking for infile.eof() works well.



1.2 Condition Number

I looked up the matrix manipulation packages at Netlib, and the Lapack++ 
routines I thought I wanted suggested using TNT (which my students have been
using for a while now). You can get the source files from 
http://math.nist.gov/tnt/; get the beta version 3.0.8 (which supports matrix
and vector linear algebra). You'll want to include tnt_matrix.h and 
tnt_linalg.h, and probably want to use the namespaces TNT and Linear_Algebra.
Add a variable TNTINCLUDES to your Makefile pointing to the TNT directory, and 
add $(TNTINCLUDES) to your compilation line:
	TNTINCLUDES= -I tnt_3_0_8/tnt
	ConditionNumber: ConditionNumber.o
		$(CXX) $(CPPFLAGS) -o ConditionNumber ConditionNumber.o $(TNTINCLUDES)
	ConditionNumber.o: ConditionNumber.cpp
		$(CXX) $(CPPFLAGS) -c ConditionNumber.cpp $(TNTINCLUDES)

Setting up your vectors and matrices is slightly awkward. It would be nice to
be able to just type them in:
	Vector<double> b(2) = {0.217, 0.254};	 XXX Doesn't work
Yong Chen and I came up with a slightly more awkward method; the 
vector constructor allows us to send in an array for initialization:
	double bStart[] = {0.217, 0.254};
	Vector<double> b(2, bStart);
The array for A must be one-dimensional, in C-order:
	double AStart[] = {0.780, 0.563, 0.913, 0.659};
	Matrix<double> A(2,2, AStart);

(a) TNT defines matrix multiplication and dot products with *, so if
x is also a vector then A*x-b will be the residual. Output with
	cout << A << endl;
will print out the shape of the matrix or vector along with the values. 

(b) You'll want to construct an SVD for A and then use getSingularValues(S)
to fill a vector S with the singular values. 

(c) You'll want to construct an LU decomposition for A, and then use solve(b)
to return a vector with the solution. Don't forget to subtract (1.,-1.) to 
see what the error is, and divide by the machine precision to see how
serious it is. Notice that TNT is still unfinished: to dividing a vector V by
a scalar eps use (1.0/eps)*V. (Division by a scalar is not yet defined.)

Python: 

To type in a scipy array, you need to explicitly convert it from a list, e.g.
	b = scipy.array([0.217, 0.254])
You'll want to import scipy.linalg as well as scipy.
The functions scipy.dot, scipy.linalg.svd, and scipy.linalg.solve will
be useful. (Use help(scipy.linalg.svd), etc.)

1.3 Sherman Morrison formula

First, there is a bug in the current beta release of TNT (3.0.8), in
tnt_linalg.h. Download and replace the TNT version with the one on the
course Web site (under the link to homework exercises at the bottom of
the page). 

You'll need to write a short routine to generate the outer product of
two matrices. It should take as arguments two Vector<double>, and return
a Matrix<double>. To be perfectly general, you can get the length of the
vectors using u.size(). You can find the inverse of a matrix T by 
constructing the identity matrix I and a LU decomposition of T, and then
calling solve(I). It appears that TNT does not yet support multiplication of
a matrix by a scalar; use matrix operations to form the numerator matrix
for Delta and the denominator scalar, and then output Delta with
a double loop over the indices.


1.4 Timing Sine.

We believe there is a way to measure time in C++ to high resolution; we'd
love to hear about it if you figure it out. Here is a way that just barely
can measure the time we need.
You'll want to include <ctime> (or perhaps <time.h>, build an array of x's,
measure the start time 
	clock_t t1 = clock();
do a loop over the million sines, and then measure the stop time.
You may need to sum the sines and print the result to fool the compiler;
sometimes it will remove commands that are never used.
To convert t1 and t2 to to process time in seconds, you use 
	(double(t2)-double(t1))/CLOCKS_PER_SEC
(for Linux: CLK_TCK for Windows). You'll probably want to allocate your 
x-array as a 
	static double x[1000000];
especially if you want to try 10^7 or 10^8 spins. 

"Static" allocates the space before the program is run (on the 'heap',
whose size is known in advance); without "static" it stores x on it's
scratch space (called the 'stack'), whose size is pretty small and
depends on the machine. The stack is where your main-routine variables
are pushed when you call a subroutine. Static variables won't work if
you're doing a recursive function; in that case you need to figure out
the compiler options to increase the stack size.