thomas-solver.cpp

/* Fills solution into x. Warning: will modify c and d! */
// this is the old signature: 
//void TridiagonalSolve(const double *a, const double *b, double *c, double *d, double *x, unsigned int n){

// this is the new signature signature: 
// - we don't use arrays, instead we use the matrix and vector objects introduced in earlier pechstein exercises
// - the "unsigned int n", which denotes the size of the matrix/vectors, can be dropped because it is 
//   stored inside the matrix/vector objects and is redundant information. 
void TridiagonalSolve(const Matrix &Kh, const Vector &fh, const Vector &uh){
	int i;
	int n = fh.size(); // you could also use uh.size(), should obviously be the same! 
 
	/* Modify the coefficients. */
	// Old code: 
	// c[0] /= b[0];	/* Division by zero risk. */
	// d[0] /= b[0];	/* Division by zero would imply a singular matrix. */

	// New code: 
	// Notice that we don't have the arrays a, b, c (the diagonals of A, or "Kh", as Pechstein calls it)
	// Instead we have only one Matrix A and we know everything outside the tridiagonal 
	// structure is zero. So: 
	// - a[0] becomes Kh( 1, 0 )  (or below main diagonal)
	// - b[0] becomes Kh( 0, 0 )  (main diagonal)
	// - c[0] becomes Kh( 0, 1 )  (or above main diagonal)
	// Also note that Christoph uses parentheses "()" for matrix-entries, e.g. Kh(i, j)
	// and brackets "[]" for vector-entries, e.g. uh[i]
	// Thus two lines of code from above translate into: 
	Kh(0,1) = Kh(0,1)/Kh(0, 0);	/* Division by zero risk. */
	fh[0] = fh[0]/Kh(0,0);		/* Division by zero would imply a singular matrix. */
	
	// The rest of the code translates the same way: 
	for(i = 1; i < n; i++){
		double id = (Kh(i,i) - Kh(i,i+1) * Kh(i-1,i));  /* Division by zero risk. */
		Kh(i,i+1) = Kh(i,i+1) / id;	        /* Last value calculated is redundant. */
		fh[i] = (fh[i] - fh[i-1] * Kh(i+1,i))/id;
	}
 
	/* Now back substitute. */
	uh[n - 1] = fh[n - 1];
	for(i = n - 2; i >= 0; i--)
		uh[i] = fh[i] - A(i,i+1) * uh[i + 1];

	// done :) 
}