commit 13bf8126e1c48f251437e9bb2b1701ecd525efdb Author: Anders Damsgaard Christensen <adc@geo.au.dk> Date: Sun, 14 Oct 2012 11:07:41 +0200 First commit Diffstat:
| A | ODEs/Makefile | | | 24 | ++++++++++++++++++++++++ |
| A | ODEs/main.cpp | | | 53 | +++++++++++++++++++++++++++++++++++++++++++++++++++++ |
| A | ODEs/plot.gp | | | 9 | +++++++++ |
| A | ODEs/rkdriver.cpp | | | 46 | ++++++++++++++++++++++++++++++++++++++++++++++ |
| A | ODEs/rkstep.cpp | | | 20 | ++++++++++++++++++++ |
| A | QR-dec/Makefile | | | 58 | ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ |
| A | QR-dec/header.h | | | 17 | +++++++++++++++++ |
| A | QR-dec/main.cpp | | | 173 | +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ |
| A | QR-dec/plotall.gp | | | 9 | +++++++++ |
| A | QR-dec/qrfunc.cpp | | | 90 | +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ |
| A | QR-dec/qrfunc.h | | | 49 | +++++++++++++++++++++++++++++++++++++++++++++++++ |
| A | README.rst | | | 4 | ++++ |
| A | eigen/Makefile | | | 69 | +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ |
| A | eigen/header.h | | | 20 | ++++++++++++++++++++ |
| A | eigen/jacobi.cpp | | | 103 | +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ |
| A | eigen/jacobi.h | | | 39 | +++++++++++++++++++++++++++++++++++++++ |
| A | eigen/main.cpp | | | 88 | +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ |
| A | exam/Makefile | | | 155 | +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ |
| A | exam/README.rst | | | 89 | +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ |
| A | exam/check.cpp | | | 11 | +++++++++++ |
| A | exam/check.h | | | 7 | +++++++ |
| A | exam/functions.h | | | 25 | +++++++++++++++++++++++++ |
| A | exam/html/Makefile.html | | | 182 | +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ |
| A | exam/html/README.html | | | 406 | +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ |
| A | exam/html/check.cpp.html | | | 39 | +++++++++++++++++++++++++++++++++++++++ |
| A | exam/html/check.h.html | | | 32 | ++++++++++++++++++++++++++++++++ |
| A | exam/html/functions.h.html | | | 52 | ++++++++++++++++++++++++++++++++++++++++++++++++++++ |
| A | exam/html/mainA.cpp.html | | | 87 | +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ |
| A | exam/html/mainB.cpp.html | | | 78 | ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ |
| A | exam/html/mainC.cpp.html | | | 78 | ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ |
| A | exam/html/mainD.cpp.html | | | 96 | +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ |
| A | exam/html/ode.cpp.html | | | 191 | +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ |
| A | exam/html/ode.h.html | | | 114 | +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ |
| A | exam/html/plotA.gp.html | | | 32 | ++++++++++++++++++++++++++++++++ |
| A | exam/html/plotB.gp.html | | | 32 | ++++++++++++++++++++++++++++++++ |
| A | exam/html/plotC.gp.html | | | 33 | +++++++++++++++++++++++++++++++++ |
| A | exam/html/plotD.gp.html | | | 33 | +++++++++++++++++++++++++++++++++ |
| A | exam/html/typedefs.h.html | | | 41 | +++++++++++++++++++++++++++++++++++++++++ |
| A | exam/html/vector_arithmetic.h.html | | | 116 | +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ |
| A | exam/mainA.cpp | | | 59 | +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ |
| A | exam/mainB.cpp | | | 50 | ++++++++++++++++++++++++++++++++++++++++++++++++++ |
| A | exam/mainC.cpp | | | 50 | ++++++++++++++++++++++++++++++++++++++++++++++++++ |
| A | exam/mainD.cpp | | | 68 | ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ |
| A | exam/ode.cpp | | | 163 | +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ |
| A | exam/ode.h | | | 87 | +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ |
| A | exam/plotA.gp | | | 8 | ++++++++ |
| A | exam/plotB.gp | | | 8 | ++++++++ |
| A | exam/plotC.gp | | | 9 | +++++++++ |
| A | exam/plotD.gp | | | 9 | +++++++++ |
| A | exam/typedefs.h | | | 15 | +++++++++++++++ |
| A | exam/vector_arithmetic.h | | | 89 | +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ |
| A | hello/Makefile | | | 71 | +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ |
| A | hello/hello.B.cpp | | | 6 | ++++++ |
| A | hello/hello_user.A.cpp | | | 19 | +++++++++++++++++++ |
| A | hello/main.B.cpp | | | 14 | ++++++++++++++ |
| A | hello/main.B.h | | | 8 | ++++++++ |
| A | hello/your_username.B.cpp | | | 14 | ++++++++++++++ |
| A | integration/Makefile | | | 71 | +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ |
| A | integration/functions.h | | | 27 | +++++++++++++++++++++++++++ |
| A | integration/header.h | | | 18 | ++++++++++++++++++ |
| A | integration/integrator.cpp | | | 120 | +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ |
| A | integration/integrator.h | | | 68 | ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ |
| A | integration/main.cpp | | | 57 | +++++++++++++++++++++++++++++++++++++++++++++++++++++++++ |
| A | leastsq/Makefile | | | 67 | +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ |
| A | leastsq/data.A.txt | | | 7 | +++++++ |
| A | leastsq/data.B.txt | | | 10 | ++++++++++ |
| A | leastsq/file_o.cpp | | | 25 | +++++++++++++++++++++++++ |
| A | leastsq/header.h | | | 21 | +++++++++++++++++++++ |
| A | leastsq/lsfit.cpp | | | 61 | +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ |
| A | leastsq/lsfit.h | | | 41 | +++++++++++++++++++++++++++++++++++++++++ |
| A | leastsq/main.cpp | | | 89 | +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ |
| A | leastsq/plotall.gp | | | 17 | +++++++++++++++++ |
| A | leastsq/qrfunc.cpp | | | 90 | +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ |
| A | leastsq/qrfunc.h | | | 49 | +++++++++++++++++++++++++++++++++++++++++++++++++ |
| A | matrixmul/Makefile | | | 150 | +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ |
| A | matrixmul/c-arrofarrs.c | | | 72 | ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ |
| A | matrixmul/c-linarr.c | | | 58 | ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ |
| A | matrixmul/cpp-linvectors.cpp | | | 49 | +++++++++++++++++++++++++++++++++++++++++++++++++ |
| A | matrixmul/cpp-vectorofvectors.cpp | | | 49 | +++++++++++++++++++++++++++++++++++++++++++++++++ |
| A | matrixmul/julia.jl | | | 13 | +++++++++++++ |
| A | matrixmul/lua-arrofarrs.lua | | | 30 | ++++++++++++++++++++++++++++++ |
| A | matrixmul/lua-linarr.lua | | | 25 | +++++++++++++++++++++++++ |
| A | matrixmul/octave.m | | | 14 | ++++++++++++++ |
| A | matrixmul/plot.gp | | | 9 | +++++++++ |
| A | matrixmul/python-numpy.py | | | 16 | ++++++++++++++++ |
| A | montecarlo/Makefile | | | 89 | +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ |
| A | montecarlo/functions.h | | | 20 | ++++++++++++++++++++ |
| A | montecarlo/header.h | | | 18 | ++++++++++++++++++ |
| A | montecarlo/main.cpp | | | 62 | ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ |
| A | montecarlo/montecarlo.cpp | | | 88 | +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ |
| A | montecarlo/montecarlo.h | | | 58 | ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ |
| A | montecarlo/plot.gp | | | 17 | +++++++++++++++++ |
| A | notes.rst | | | 398 | +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ |
| A | notes/book.pdf | | | 0 | |
| A | notes/eigen.pdf | | | 0 | |
| A | notes/fft.pdf | | | 0 | |
| A | notes/integration.pdf | | | 0 | |
| A | notes/interp.pdf | | | 0 | |
| A | notes/krylov.pdf | | | 0 | |
| A | notes/leastsq.pdf | | | 0 | |
| A | notes/lineq.pdf | | | 0 | |
| A | notes/montecarlo.pdf | | | 0 | |
| A | notes/odes.pdf | | | 0 | |
| A | notes/roots.pdf | | | 0 | |
| A | notes/sfun.pdf | | | 0 | |
| A | optimization/Makefile | | | 69 | +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ |
| A | optimization/downhill_simplex.cpp | | | 173 | +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ |
| A | optimization/downhill_simplex.h | | | 61 | +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ |
| A | optimization/functions.h | | | 24 | ++++++++++++++++++++++++ |
| A | optimization/header.h | | | 12 | ++++++++++++ |
| A | optimization/main.B.cpp | | | 93 | +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ |
| A | optimization/main.cpp | | | 68 | ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ |
| A | roots/Makefile | | | 74 | ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ |
| A | roots/functions.h | | | 84 | +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ |
| A | roots/header.h | | | 20 | ++++++++++++++++++++ |
| A | roots/main.cpp | | | 79 | +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ |
| A | roots/newton.cpp | | | 89 | +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ |
| A | roots/qrfunc.cpp | | | 90 | +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ |
| A | roots/qrfunc.h | | | 49 | +++++++++++++++++++++++++++++++++++++++++++++++++ |
119 files changed, 6605 insertions(+), 0 deletions(-)
diff --git a/ODEs/Makefile b/ODEs/Makefile @@ -0,0 +1,24 @@ +#!/bin/env make + +SHELL = /bin/sh +CC = g++ +CFLAGS = -Wall -O3 +OBJ = rkdriver.o rkstep.o +BIN = ode + +all: $(BIN) ODE.output ode.png + +ODE.output: $(BIN) + ./$(BIN) > ODE.output + +$(BIN): main.o rkdriver.o rkstep.o + $(CC) $(CFLAGS) -o $(BIN) main.o $(OBJ) + +clean: + rm -f $(BIN) *.txt *~ *.output *.png *.o + +%.png: plot.gp + gnuplot < $< + +edit: + vim -p Makefile *.cpp *.gp diff --git a/ODEs/main.cpp b/ODEs/main.cpp @@ -0,0 +1,53 @@ +#include <iomanip> +#include <cstdlib> // Exit function +#include <string.h> +#include <cmath> +#include <vector> + +// File I/O +#include <iostream> +#include <fstream> + +using namespace std; + +// this function gives the sign +bool sign(double n) +{ return n >= 0; }; + +void rkdriver(void f(int, double, vector<double>*, vector<double>*),int n, vector<double>* tlist, vector<vector<double> >* ylist, double b, double h, double acc, double eps, int max ); + +// this function simulates a cannonball +void func1(int n, double x, vector<double>* y, vector<double>* dydx){ + int A = sign((*y)[2]); + int B = sign((*y)[3]); + double k = 0.02; + double g = 9.82; + (*dydx)[0]=(*y)[2]; + (*dydx)[1]=(*y)[3]; + (*dydx)[2]=sqrt((*y)[2]*(*y)[2]+(*y)[3]*(*y)[3])*-A*k; + (*dydx)[3]=-g-sqrt((*y)[2]*(*y)[2]+(*y)[3]*(*y)[3])*B*k; +}; + +int main () { + + int n=4; // Number of equations + int max=10000; // Maximum "iterations" + double a=0, b=15, h=0.01, acc=0.001, eps=0.001; + vector<double> tlist(max); // + vector<vector<double> > ylist(max,vector<double>(n)); + tlist[0]=a; ylist[0][0]=0; ylist[0][1]=0; ylist[0][2]=100; + ylist[0][3]=100; + + rkdriver(&func1,n,&tlist,&ylist,b,h,acc,eps,max); + + // Printing the output + int it = 0; + while(tlist[it]<b) { + cout << tlist[it] << " "; + for (int j = 0; j<n;j++) cout << ylist[it][j] << " "; + cout << endl; + it++; + }; + +}; + diff --git a/ODEs/plot.gp b/ODEs/plot.gp @@ -0,0 +1,9 @@ +set term png enhanced +set output "ode.png" +set size 1.0,1.0 +set xlabel 'x' +set ylabel 'y' +set title "Projectile trajectory" +plot [0:1500] [-100:500] \ + "ODE.output" u 2:3 w p title 'data' + diff --git a/ODEs/rkdriver.cpp b/ODEs/rkdriver.cpp @@ -0,0 +1,46 @@ +#include <math.h> +#include <stdlib.h> +#include <stdio.h> +#include <iostream> +#include <fstream> +#include <vector> +using namespace std; + +void rkstep(void f(int,double,vector<double>*,vector<double>*),int n, double t,vector<double>* y,double h,vector<double>* y1,vector<double>* dy); + +// Main driver function +void rkdriver(void f(int, double, vector<double>*, vector<double>*),int n, vector<double>* tlist, vector<vector<double> >* ylist, double b, double h, double acc, double eps, int max ) { + + int i=0; //iterator + double t = (*tlist)[0]; + double a = t; + vector<double> dy(n); + vector<double> y1(n); + + while((*tlist)[i]<b) { + double t=(*tlist)[i]; + vector<double> y=(*ylist)[i]; + if(t+h>b) h=b-t; + rkstep(f,n,t,&y,h,&y1,&dy); + double err=0; for(int j=0;j<n;j++)err+=dy[j]*dy[j]; + err=sqrt(err); + double normy=0; for(int j=0;j<n;j++)normy+=y1[j]*y1[j]; + normy=sqrt(normy); + double tol=(normy*eps+acc)*sqrt(h/(b-a)); + + if(tol>err){ // accept step and go on + i++; + if(i>max-1) { + cout << "Reached max step \nIncrease max step to go further \n"; + exit(1); + }; + + (*tlist)[i]=t+h; + + for(int j=0;j<n;j++) (*ylist)[i][j]=y1[j]; + }; + + if(err>0) h = h*pow(tol/err,0.25)*0.95; + else h = 2*h; + }//end while +}; diff --git a/ODEs/rkstep.cpp b/ODEs/rkstep.cpp @@ -0,0 +1,20 @@ +#include <stdlib.h> +#include <cmath> +#include <vector> + +using namespace std; + +// Stepper function +void rkstep( + void f(int, double, vector<double>*, vector<double>*), + int n, double t, vector<double>* y, double h, vector<double>* y1, vector<double>* dy) +{ + vector<double> k0(n); + vector<double> y12(n); + vector<double> k12(n); + (*f)(n,t,y,&k0); + for(int i = 0;i<n;i++) y12[i] = (*y)[i]+k0[i]*h/2; + (*f)(n,(t+h/2),&y12,&k12); + for(int i = 0;i<n;i++) (*y1)[i] = (*y)[i]+k12[i]*h; + for(int i = 0;i<n;i++) (*dy)[i] = (k0[i]-k12[i])*1*(h)/2; +}; diff --git a/QR-dec/Makefile b/QR-dec/Makefile @@ -0,0 +1,58 @@ +# Define compiler +CC=g++ + +# Define compiler flags (show all warnings) +CPPFLAGS=-Wall + +# Define linker flags +LDFLAGS= + +# Define extra libraries to be dynamically linked +LDLIBS+=-larmadillo + +# Compile optimized code +CPPFLAGS+=-O2 + +# Compile debuggable code +#CPPFLAGS+=-g + +# Compile profilable code +#CPPFLAGS+=-pg +#LDFLAGS+=-pg + +# Define linker +LD=g++ + +# Filenames of source code +SRC=$(shell ls *.cpp) + +# Filenames of object files +OBJ=$(SRC:.cpp=.o) + +# Remove file type extension for binary filename +BIN=qr + +# The default "all" depends on A and B + +plot.png: performance.dat + gnuplot plotall.gp + +%.dat: $(BIN) + ./$(BIN) + +$(BIN): $(OBJ) + @# Link object files together + $(LD) $(LDFLAGS) $(OBJ) -o $(BIN) $(LDLIBS) + @# Execute program and redirect stdout to file + @#./$(BIN) > out.txt + +clean: + @# Remove object files + rm -f $(OBJ) + @# Remove binary + rm -f $(BIN) + @# Remove datafiles and plot + rm -f *.dat *.png +edit: + vim -p Makefile *.cpp *.h *.gp + diff --git a/QR-dec/header.h b/QR-dec/header.h @@ -0,0 +1,17 @@ +// Make sure header is only included once +#ifndef header_h_ +#define header_h_ + +// Define whether the program should output values of matrices +const bool verbose = false; + +// Choose vector length variable type +typedef int Lengthtype; + +// Choose floating-point precision +typedef double Floattype; + +// Prototype functions +void check(const bool statement); + +#endif diff --git a/QR-dec/main.cpp b/QR-dec/main.cpp @@ -0,0 +1,173 @@ +#include <iostream> +#include <fstream> +#include <ctime> +#include <armadillo> +#include "header.h" +#include "qrfunc.h" + +int main(int argc, char* argv[]) +{ + // Namespace declarations + using std::cout; + + // Timer variables + double tic1, toc1, tic2, toc2; + + // Create 2D matrices from Armadillo templates + Lengthtype width; + if (argc == 1) // If no command line argument is given... + width = 4; // Matrices are 4x4. + else + width = atoi(argv[1]); // Else the specified size + + const Lengthtype height = width; + cout << "\nInitializing " << width << " by " << height << " matrices.\n"; + arma::mat A = arma::randu<arma::mat>(width, height); + + // QR decomposition is performed upon initialization of QR class + tic1 = clock(); // Start clock1 + QR qr(A); + toc1 = clock(); // Stop clock1 + + //// QR decomposition check + cout << "\n\033[1;33m--- QR decomposition check ---\033[0m\n"; + if (verbose == true) { + // Display values to stdout + qr.A.print("Original A:"); + } + + // Check QR decomposition + if (verbose == true) { + qr.Q.print("Q, after QR dec:"); + qr.Q.t().print("Q^T:"); + qr.R.print("R, after QR dec:"); + cout << '\n'; + } + + // Check that matrix is orthogonal + arma::mat QTQ = qr.Q.t()*qr.Q; + Floattype checksum = arma::sum(arma::sum(QTQ)); + if (verbose == true) { + QTQ.print("Q^T Q:"); + cout << "sum = " << checksum << '\n'; + } + cout << "Check: Q^T Q = 1: \t"; + check(checksum-(Floattype)height < 1e-12f); + + cout << "Check: QR = A by ||QR-A|| = 0: "; + checksum = sum(sum((qr.Q*qr.R)-qr.A)); + check(checksum < 1e-12f); + + + //// Solving linear equations + cout << "\n\033[1;33m--- Solving linear equations: Ax=b ---\033[0m\n"; + cout << "Solving QRx=b.\n"; + arma::vec b = arma::randu<arma::vec>(qr.size()); + + // Perform back-substitution of QR system + if (verbose == true) { + b.print("Vector b:"); + } + + arma::vec x = qr.backsub(b); + + if (verbose == true) { + x.print("Solution, x:"); + } + + cout << "Check: Ax = b by |Ax-b| = 0: "; + checksum = arma::sum(arma::sum(qr.A*x - b)); + check(checksum < 1e-12f); + + //// Calculating the determinant + cout << "\n\033[1;33m--- Determinant of A ---\033[0m\n"; + cout << "|det(A)| = " << qr.det() << '\n'; + + //// Calculating the inverse + cout << "\n\033[1;33m--- Inverse of A ---\033[0m\n"; + arma::mat Ainv = qr.inverse(); + if (verbose == true) + Ainv.print("A^(-1):"); + cout << "Check: (A^(-1))^(-1) = A: "; + QR qrinv(Ainv); + arma::mat Ainvinv = qrinv.inverse(); + bool equal = true; // Elementwise comparison of values + for (Lengthtype i=0; i<width; ++i) { + for (Lengthtype j=0; j<height; ++j) { + if (fabs(A(i,j)-Ainvinv(i,j)) > 1e12f) { + equal = false; + cout << "At (" << i << "," << j << ") = " + << "A = " << A(i,j) + << ", (A^(-1))^(-1) = " << Ainvinv(i,j) << '\n'; + } + } + } + check(equal); + + //// Use the Armadillo built-in QR decomposition + tic2 = clock(); // Start clock2 + arma::mat Q, R; + arma::qr(Q, R, A); + toc2 = clock(); // Stop clock2 + + //// Statistics + // Print resulting time of homegrown function and Armadillo function + cout << "\n\033[1;33m--- Performance comparison ---\033[0m\n"; + double t1 = (toc1 - tic1)/(CLOCKS_PER_SEC); + double t2 = (toc2 - tic2)/(CLOCKS_PER_SEC); + cout << "Homegrown implementation: \t" << t1 << " s \n" + << "Armadillo implementation: \t" << t2 << " s \n"; + + cout << "Benchmarking performance across a range of sizes...\n"; + + // Write results to file + std::ofstream outstream; + // Open outfile for write + outstream.open("performance.dat"); + double homegrown_time, armadillo_time; + + // Define sizes + Lengthtype dims[] = {32, 64, 128, 256, 512, 1024, 2048}; + Lengthtype ndims = sizeof(dims)/sizeof(int); // Number of entries in dims + + // Loop through sizes and record performance + for (Lengthtype i=0; i<ndims; ++i) { + cout << " " << dims[i] << std::flush; + // Generate input matrix + arma::mat An = arma::randu<arma::mat>(dims[i],dims[i]); + + // Homegrown implementation + tic1 = clock(); + QR qrn = QR(An); + toc1 = clock(); + + // Armadillo implementation + tic2 = clock(); + arma::mat Qn, Rn; + arma::qr(Qn, Rn, An); + toc2 = clock(); + + // Record time spent + homegrown_time = (toc1 - tic1)/(CLOCKS_PER_SEC); + armadillo_time = (toc2 - tic2)/(CLOCKS_PER_SEC); + + // Write time to file in three columns + outstream << dims[i] << '\t' << homegrown_time << '\t' << armadillo_time << '\n'; + } + + cout << '\n'; + // Close file + outstream.close(); + + // Return successfully + return 0; +} + +void check(const bool statement) +{ + using std::cout; + if (statement == true) + cout << "\t\033[0;32mPassed\033[0m\n"; + else + cout << "\t\033[1;31mFail!!\033[0m\n"; +} diff --git a/QR-dec/plotall.gp b/QR-dec/plotall.gp @@ -0,0 +1,9 @@ +set terminal png # Set output file format +set output "plot.png" # Set output filename +set key top left +set xlabel "Matrix width and height" +set ylabel "Execution time [s]" +set title "Performance comparison of QR decomposition" +set log xy +set grid +plot "performance.dat" u 1:2 title "Homegrown" w lp, "performance.dat" u 1:3 title "Armadillo" w lp diff --git a/QR-dec/qrfunc.cpp b/QR-dec/qrfunc.cpp @@ -0,0 +1,90 @@ +#include <iostream> +#include <armadillo> +#include "header.h" +#include "qrfunc.h" + +// QR decomposition constructor +QR::QR(arma::mat &A) + : n(A.n_cols), + A(A), + Q(A) +{ + // Initialize output structures + R = arma::zeros<arma::mat> (n,n); + + // Perform QR decomposition straight away + decomp(); +} + +// Class deconstructor (equivalent to compiler destructor) +QR::~QR() { }; + +// Return system size +Lengthtype QR::size() +{ + return n; +} + +// QR decomposition function of Armadillo matrix. +// Returns right matrix R, and modifies A into Q. +// Uses Gram-Schmidt orthogonalization +void QR::decomp() +{ + Floattype r, s; + Lengthtype j; + for (Lengthtype i=0; i<n; ++i) { + r = dot(Q.col(i), Q.col(i)); + R(i,i) = sqrt(r); + Q.col(i) /= sqrt(r); // Normalization + for (j=i+1; j<n; ++j) { + s = dot(Q.col(i), Q.col(j)); + Q.col(j) -= s*Q.col(i); // Orthogonalization + R(i,j) = s; + } + } +} + +// Solve the square system of linear equations with back +// substitution. T is an upper triangular matrix. +arma::vec QR::backsub(arma::vec &b) +{ + Floattype tmpsum; + arma::vec x = Q.t() * b; + for (Lengthtype i=n-1; i>=0; --i) { + tmpsum = 0.0f; + for (Lengthtype k=i+1; k<n; ++k) + tmpsum += R(i,k) * x(k); + + x(i) = 1.0f/R(i,i) * (x(i) - tmpsum); + } + return x; +} + +// Calculate the (absolute value of the) determinant of +// matrix A given the Q and R components. +// det(A) = det(Q) * det(R), |det(Q) = 1| +// => |det(A)| = |det(R)| +Floattype QR::det() +{ + Floattype determinant = 1.0f; + for (Lengthtype i=0; i<n; ++i) + determinant *= R(i,i); + return fabs(determinant); +} + +// Calculate the inverse of matrix A given the Q and R +// components. +arma::mat QR::inverse() +{ + arma::mat inv = arma::zeros<arma::mat> (n, n); + // In vector z, all elements are equal to 0, except z(i) = 1 + arma::vec z = arma::zeros<arma::mat> (n); + + for (Lengthtype i=0; i<n; ++i) { + z(i) = 1.0f; // Element i changed to 1 + inv.col(i) = backsub(z); + z(i) = 0.0f; // Element i changed back to 0 + } + + return inv; +} diff --git a/QR-dec/qrfunc.h b/QR-dec/qrfunc.h @@ -0,0 +1,49 @@ +// Make sure header is only included once +#ifndef qrfunc_h_ +#define qrfunc_h_ + +#include <armadillo> +#include "header.h" + +// QR structure +class QR { + private: + // System size + const Lengthtype n; + + public: + //// Data + + // Input data + arma::mat A; + + // QR decomposition matrices + arma::mat Q; + arma::mat R; + + //// Prototype functions + + // Constructor prototype + QR(arma::mat &A); + + // Destructor + ~QR(); + + // Return system size + Lengthtype size(); + + // QR decomposition of Armadillo matrix A, returning R + // and modified A (=Q) + void decomp(); + + // Backsubstitution of triangular system + arma::vec backsub(arma::vec &b); + + // Absolute value of the determinant of matrix R + Floattype det(); + + // Inverse of matrix A + arma::mat inverse(); +}; + +#endif diff --git a/README.rst b/README.rst @@ -0,0 +1,4 @@ +numeric +======= + +Routines written in C++ to handle various numerical routines. Code is written in 2011 for the course "Numerical methods", held by Dmitri Fedorov, at the department of Physics and Astronomy at Aarhus University, Denmark. diff --git a/eigen/Makefile b/eigen/Makefile @@ -0,0 +1,69 @@ +# Define compiler +CC=g++ + +# Define compiler flags (show all warnings) +CPPFLAGS=-Wall + +# Define linker flags +LDFLAGS= + +# Define extra libraries to be dynamically linked +LDLIBS+=-larmadillo + +# Compile optimized code +#CPPFLAGS+=-O2 + +# Compile debuggable code +#CPPFLAGS+=-g + +# Compile profilable code +#CPPFLAGS+=-pg +#LDFLAGS+=-pg + +# Define linker +LD=g++ + +# Filenames of source code +SRC=$(shell ls *.cpp) + +# Filenames of object files +OBJ=$(SRC:.cpp=.o) + +# Remove file type extension for binary filename +BIN=eigen + +# The default "all" depends on A and B + +#all: A B +A: + ./eigen 11 + +#A: plot.A.png + +#B: plot.B.png + +#plot.%.png: fit.A.dat fit.B.dat data.A.txt data.B.txt +# gnuplot plotall.gp + +#fit.A.dat: $(BIN) +# ./$(BIN) data.A.txt fit.A.dat + +#fit.B.dat: $(BIN) +# ./$(BIN) data.B.txt fit.B.dat + +$(BIN): $(OBJ) + @# Link object files together + $(LD) $(LDFLAGS) $(OBJ) -o $(BIN) $(LDLIBS) + @# Execute program and redirect stdout to file + @#./$(BIN) > out.txt + +clean: + @# Remove object files + rm -f $(OBJ) + @# Remove binary + rm -f $(BIN) + @# Remove datafiles and plot + #rm -f *.dat *.png +edit: + vim -p Makefile *.cpp *.h *.gp + diff --git a/eigen/header.h b/eigen/header.h @@ -0,0 +1,20 @@ +// Make sure header is only included once +#ifndef HEADER_H_ +#define HEADER_H_ + +// Define whether the program should output values of matrices +const bool verbose = false; +//const bool verbose = true; + +// Choose vector length variable type +typedef int Lengthtype; + +// Choose floating-point precision +//typedef float Floattype; +typedef double Floattype; +//typedef long double Floattype; + +// Prototype for checking function +void check(const bool statement); + +#endif diff --git a/eigen/jacobi.cpp b/eigen/jacobi.cpp @@ -0,0 +1,103 @@ +#include <iostream> +#include <cmath> +#include <armadillo> +#include "header.h" +#include "jacobi.h" + +// Constructor: Perform rotation straight away +Jacobi::Jacobi(const arma::Mat<Floattype> &A) + : n(A.n_rows), At(A), e(n), V(n,n) +{ + // Initialize diagonal vector to be filled with eigenvalues + e = A.diag(); + + // Set main diagonal values to 1, off-diagonal values to 0 + V.eye(); + + Lengthtype p, q, i; // Iterator vars + int changed; + Lengthtype rotations = 0; // Number of rotations performed + + // Cell value variables, used as local storage for mat/vec vals. + Floattype app, aqq, apq, phi, c, s, app1, aqq1; + Floattype aip, api, aiq, aqi, vip, viq; + + do { + changed = 0; + for (p=0; p<n; ++p) { // Row iteration + for (q=p+1; q<n; ++q) { // Cols right of diagonal + + // Initialize cell-relevant data + app = e(p); + aqq = e(q); + apq = At(p,q); + phi = 0.5f * atan2(2.0f * apq, aqq-app); + c = cos(phi); + s = sin(phi); + app1 = c*c * app - 2.0f * s * c * apq + s*s * aqq; + aqq1 = s*s * app + 2.0f * s * c * apq + c*c * aqq; + + if (app1 != app || aqq1 != aqq) { + changed = 1; + ++rotations; + e(p) = app1; + e(q) = aqq1; + At(p,q) = 0.0f; + + for (i = 0; i<p; ++i) { + aip = At(i,p); + aiq = At(i,q); + At(i,p) = c * aip - s * aiq; + At(i,q) = c * aiq + s * aip; + } + + for (i=p+1; i<q; ++i) { + api = At(p,i); + aiq = At(i,q); + At(p,i) = c * api - s * aiq; + At(i,q) = c * aiq + s * api; + } + + for (i=q+1; i<n; ++i) { + api = At(p,i); + aqi = At(q,i); + At(p,i) = c * api - s * aqi; + At(q,i) = c * aqi + s * api; + } + + for (i=0; i<n; ++i) { + vip = V(i,p); + viq = V(i,q); + V(i,p) = c * vip - s * viq; + V(i,q) = c * viq + s * vip; + } + + } // if end + } // q loop end + } // p loop end + } while (changed != 0); // do-while end + + // Make transformed matrix symmetric + At = arma::symmatu(At); + + if (verbose == true) + std::cout << "\nPerformed " << rotations << " Jacobi rotations.\n"; +} + +// Return transformed matrix +arma::Mat<Floattype> Jacobi::trans() +{ + return At; +} + +// Return matrix of eigenvectors +arma::Mat<Floattype> Jacobi::eigenvectors() +{ + return V; +} + +// Return vector of eigenvalues +arma::Col<Floattype> Jacobi::eigenvalues() +{ + return e; +} diff --git a/eigen/jacobi.h b/eigen/jacobi.h @@ -0,0 +1,39 @@ +// Make sure header is only included once per object +#ifndef JACOBI_H_ +#define JACOBI_H_ + +#include <armadillo> +#include "header.h" + +// lsfit structure +class Jacobi { + private: + const Lengthtype n; // Matrix width and height + + // Transformed matrix A: At + arma::Mat<Floattype> At; + + // Diagonal + arma::Col<Floattype> e; + + // Matrix of eigenvectors + arma::Mat<Floattype> V; + + public: + + // Constructor. Arguments: input matrix + Jacobi(const arma::Mat<Floattype> &A); + + // Destructor + //~Jacobi(); + + // Return transformed matrix + arma::Mat<Floattype> trans(); + + // Return matrix of eigenvectors + arma::Mat<Floattype> eigenvectors(); + + // Return vector of eigenvalues + arma::Col<Floattype> eigenvalues(); +}; +#endif diff --git a/eigen/main.cpp b/eigen/main.cpp @@ -0,0 +1,88 @@ +#include <iostream> +#include <armadillo> +#include "header.h" +#include "jacobi.h" + +int main(int argc, char* argv[]) +{ + // Namespace declarations + using std::cout; + + // Define matrix size + Lengthtype msize = 6; + if (argc == 2) { + if ((strcmp(argv[1], "-h") == 0) || (strcmp(argv[1], "--help") == 0)) { + cout << "Usage: " << argv[0] << " [matrix size]\n" + << "If matrix size is not specified, " + << "the matrix width and length will be " + << msize << ".\n"; + return 1; + } + msize = atoi(argv[1]); // Use specified matrix size + } + + // Calculate machine precision + Floattype eps = 1.0f; + while (1.0f + eps != 1.0f) + eps /= 2.0f; + //cout << "Machine precision of '" << typeid(eps).name() + // << "' type is: eps = " << eps << '\n'; + + Floattype checksum; + // Generate input matrix A, which is symmetric + cout << "\n\033[1;33m--- Input data check ---\033[0m\n"; + arma::Mat<Floattype> A = symmatu(arma::randu< arma::Mat<Floattype> > (msize, msize)); + checksum = arma::sum(arma::sum(A - A.t())); + cout << "Symmetry check: A = A^T: "; + check(checksum < eps); + if (verbose == true) { + A.print("Original matrix:"); + } + + // Perform Jacobi diagonalization of matrix A + Jacobi Diag = Jacobi(A); + + cout << "\n\033[1;33m--- Diagonalization check ---\033[0m\n"; + if (verbose == true) + Diag.trans().print("Transformed matrix (At):"); + cout << "Check: V V^T = 1: \t"; + checksum = arma::sum(arma::sum(Diag.eigenvectors().t() * Diag.eigenvectors())); + check(fabs(checksum - (Floattype)msize) < eps*msize*msize); + if (verbose == true) { + Diag.eigenvalues().print("Eigenvalues (e):"); + Diag.eigenvectors().print("Eigenvectors (V):"); + Diag.eigenvectors().t().print("V^T"); + (Diag.eigenvectors() * Diag.eigenvectors().t()).print("V V^T"); + (Diag.eigenvectors().t() * A * Diag.eigenvectors()).print("V^T A V"); + (Diag.eigenvectors().t() * Diag.trans() * Diag.eigenvectors()).print("V^T At V"); + } + + // Armadillo implementation + arma::Mat<Floattype> V_a (msize, msize); + arma::Col<Floattype> e_a (msize); + arma::eig_sym(e_a, V_a, A); + if (verbose == true) { + e_a.print("Armadillo eigenvalues:"); + V_a.print("Armadillo eigenvectors:"); + } + cout << "\n\033[1;33m--- Armadillo comparison ---\033[0m\n"; + checksum = arma::sum(arma::sum(V_a - Diag.eigenvectors())); + cout << "Eigenvectors identical:\t"; + check(checksum < eps*msize*msize*2); + checksum = arma::sum(arma::sum(e_a - Diag.eigenvalues())); + cout << "Eigenvalues identical: \t"; + check(checksum < eps*msize*2); + + + // Return successfully + return 0; +} + +void check(const bool statement) +{ + using std::cout; + if (statement == true) + cout << "\t\033[0;32mPassed\033[0m\n"; + else + cout << "\t\033[1;31mFail!!\033[0m\n"; +} diff --git a/exam/Makefile b/exam/Makefile @@ -0,0 +1,155 @@ +# Define compiler +#CXX=g++ + +# Define compiler flags (show all warnings) +CXXFLAGS=-Wall +#CXXFLAGS=-std=c++0x + +# Define linker flags +#LDFLAGS=-fopenmp + +# Compile optimized code +CXXFLAGS+=-O2 + +# Compile debuggable code +#CXXFLAGS+=-g + +# Compile profilable code +#CXXFLAGS+=-pg +#LDFLAGS+=-pg + +# Define linker +LD=g++ + +# All source code files +SRC=$(shell ls *.cpp) + + +# Filenames of source code +SHARED_SRC=ode.cpp check.cpp +SHARED_HEADERS=typedefs.h ode.h functions.h check.h +SRC_A=mainA.cpp $(SHARED_SRC) +HEAD_A=$(SHARED_HEADERS) +SRC_B=mainB.cpp $(SHARED_SRC) +HEAD_B=$(SHARED_HEADERS) +SRC_C=mainC.cpp $(SHARED_SRC) +HEAD_C=$(SHARED_HEADERS) +SRC_D=mainD.cpp $(SHARED_SRC) +HEAD_D=$(SHARED_HEADERS) + +# Filenames of object files +OBJ_A=$(SRC_A:.cpp=.o) +OBJ_B=$(SRC_B:.cpp=.o) +OBJ_C=$(SRC_C:.cpp=.o) +OBJ_D=$(SRC_D:.cpp=.o) + +# Remove file type extension for binary filename +BIN_A=odeA +BIN_B=odeB +BIN_C=odeC +BIN_D=odeD + +# Define editor and options for `make edit` +EDITOR=vim -p + + +# The default "all" depends on A and B +all: A B C D + +A: plotA.png + +B: plotB.png + +C: plotC.png + +D: plotD.png + +plotA.png: funcA.dat plotA.gp + # Gnuplot: plotA.png + @gnuplot plotA.gp + +plotB.png: funcB.dat plotB.gp + # Gnuplot: plotB.png + @gnuplot plotB.gp + +plotC.png: funcC.dat plotC.gp + # Gnuplot: plotC.png + @gnuplot plotC.gp + +plotD.png: funcD.dat plotD.gp + # Gnuplot: plotD.png + @gnuplot plotD.gp + +funcA.dat: $(BIN_A) + @./$(BIN_A) + +funcB.dat: $(BIN_B) + @./$(BIN_B) + +funcC.dat: $(BIN_C) + @./$(BIN_C) + +funcD.dat: $(BIN_D) + @./$(BIN_D) + +$(BIN_A): $(OBJ_A) $(HEAD_A) + @# Link object files together + $(LD) $(LDFLAGS) $(OBJ_A) -o $@ $(LDLIBS) + +$(BIN_B): $(OBJ_B) $(HEAD_B) + @# Link object files together + $(LD) $(LDFLAGS) $(OBJ_B) -o $@ $(LDLIBS) + +$(BIN_C): $(OBJ_C) $(HEAD_C) + @# Link object files together + $(LD) $(LDFLAGS) $(OBJ_C) -o $@ $(LDLIBS) + +$(BIN_D): $(OBJ_D) $(HEAD_D) + @# Link object files together + $(LD) $(LDFLAGS) $(OBJ_D) -o $@ $(LDLIBS) + +clean: cleanA cleanB cleanC cleanD + +cleanA: + @# Remove object files + rm -f $(OBJ_A) + @# Remove binaries + rm -f $(BIN_A) + @# Remove datafiles and plot + rm -f funcA.dat plotA.png + +cleanB: + @# Remove object files + rm -f $(OBJ_B) + @# Remove binaries + rm -f $(BIN_B) + @# Remove datafiles and plot + rm -f funcB.dat plotB.png + +cleanC: + @# Remove object files + rm -f $(OBJ_C) + @# Remove binaries + rm -f $(BIN_C) + @# Remove datafiles and plot + rm -f funcC.dat plotC.png + +cleanD: + @# Remove object files + rm -f $(OBJ_D) + @# Remove binaries + rm -f $(BIN_D) + # Removing datafile and plot + @rm -f funcD.dat plotD.png + +htmlfiles: html/mainA.cpp.html html/mainB.cpp.html html/mainC.cpp.html html/mainD.cpp.html html/ode.cpp.html html/check.cpp.html html/check.h.html html/functions.h.html html/ode.h.html html/typedefs.h.html html/vector_arithmetic.h.html html/plotA.gp.html html/plotB.gp.html html/plotC.gp.html html/plotD.gp.html html/Makefile.html + # Generating HTML files + rst2html2 README.rst > html/README.html + +html/%.html: % + vim $< +TOhtml +"w $@" +"qall!" + + +edit: + @$(EDITOR) Makefile README.rst *.cpp *.h *.gp + diff --git a/exam/README.rst b/exam/README.rst @@ -0,0 +1,89 @@ +============================================ +README: ODE integration with complex numbers +============================================ +Exam exercise for *Numerical Methods* by Anders D. Christensen (mail_) + +File description +---------------- +- ``Makefile``: Description for GNU Make, handles compilation and execution. +- ``README.rst`` (this file): Description of numeric implementation and usage. + Written with reStructuredText syntax. +- ``check.cpp``: Function for displaying the state of a condition to stdout. +- ``check.h``: Prototype for the check-function. +- ``functions.h``: Input functions to be evaluated. +- ``mainA.cpp``: Main source code file for part A. +- ``ode.cpp``: Constructor and functions for the ODE class, including Runge-Kutta + stepper and driver. +- ``ode.h``: Header file with the ODE class. This file must be included in all + programs that want to utilize the ODE functionality. +- ``plot.gp``: Script for plotting all graphs with Gnuplot. +- ``typedefs.h``: Header file containing definitions of two main types, + ``Inttype``, a whole-number type, and ``Floattype``, a floating point number + type. The type definitions can be changed to different lengths and precisions. + The program can be compiled for verbose output by changing the ``verbose`` + variable. +- ``vector_arithmetic.h``: Operator overloading functions for the ``std::vector`` + class. + +Problem descriptions +-------------------- +The four generated executables each demonstrate the ODE solvers functionality by +performing the following tasks. The results consist of the console output and +the corresponding plot with filename ``plot<Character>.png``. +- *A*: Construct an ODE solver that can handle functions with complex values. +Demonstrate that it solves the real component correctly, by stepping along +a path in the real range. +- *B*: Demonstrate that the ODE solver can solve the imaginary component by +stepping along a path in the imaginary range. +- *C*: Demonstrate the solution of a set of complex equations by stepping +through the complex plane. +- *D*: For an integration path in the complex plane, visualize how the +requirements of absolute- and relative precision are related to the number of +integration steps, for a given floating point precision. + +Implementation +-------------- +This exercise was written in object-oriented C++ for easy reuse. For +portability, all included classes are from the standard template library. + +The necessary ``std::vector`` arithmetic operations where overloaded to support +element-wise operations, such as vector-scalar multiplication, vector-vector +addition, etc. This approach was preferred over using ``std::valarray``, since it +is not dynamically expandable. + +When creating a new ODE object, the user specifies the end-points of the linear +range, where the specified system of ordinary differential equations with +complex values will be solved. The range end-points are complex numbers +themselves, and the user can thus specify whether the integrator steps through a +range of real values, imaginary values, or both components in the complex plane. + +The solver steps through the specified range by an adaptive step size, which is +also a complex number. The user specifies the fraction of the range to be used +as a start value for the step. The default value is 0.01. + +The ODE class contains functions for writing the ODE solution to stdout +(``ODE::print``) or to a text file (``ODE::write``). The output format is the +following; the first column is the real part of x, second column the imaginary +part. The subsequent columns do in turn consist of real- and imaginary parts of +the variables in the ODE. + +The program requires a modern C++ compiler, GNU Make and Gnuplot. It has been +tested with GCC, Clang and llvm. + +Compiliation and execution +-------------------------- +To make and execute the program, go to the root folder and type `make`. This +will compile and execute the programs for part A-D, and plot output graphs. If +desired, individual parts can be compiled and executed using `make <Character>`. + +To view the source code in a browser with vim's syntax highlighting, type `make +html`, and view the files in the `html` folder. The generation of HTML files +requires a newer vim for the source code files, and Docutils for the readme. + +All output and objects can be removed using `make clean`. + + + +.. _mail: mailto:adc@geo.au.dk + +#vim: set tw=80 diff --git a/exam/check.cpp b/exam/check.cpp @@ -0,0 +1,11 @@ +// Function used for reporting the condition of a +// statement to stdout using ANSI colors. +#include <iostream> + +void check(const bool statement) +{ + if (statement == true) + std::cout << "\t\033[0;32mPassed\033[0m\n"; + else + std::cout << "\t\033[1;31mFail!!\033[0m\n"; +} diff --git a/exam/check.h b/exam/check.h @@ -0,0 +1,7 @@ +#ifndef CHECK_H_ +#define CHECK_H_ + +// Prototype for checking function +void check(const bool statement); + +#endif diff --git a/exam/functions.h b/exam/functions.h @@ -0,0 +1,25 @@ +// Make sure file is only included once per object +#ifndef FUNCTIONS_H_ +#define FUNCTIONS_H_ + +#include <vector> +#include <complex> +#include "typedefs.h" + + +//// ODEs with real+complex parts. +//// Return the derivatives at the point x,vec(y) + +std::vector<std::complex<Floattype> > + func1(const std::complex<Floattype> z, + const std::vector<std::complex<Floattype> > &y) +{ + std::vector<std::complex<Floattype> > dydz(2); + dydz[0].real() = y[1].real(); + dydz[0].imag() = y[1].imag(); + dydz[1].real() = -y[0].real(); + dydz[1].imag() = 0.5f*y[0].imag(); + return dydz; +} + +#endif diff --git a/exam/html/Makefile.html b/exam/html/Makefile.html @@ -0,0 +1,182 @@ +<!DOCTYPE HTML PUBLIC "-//W3C//DTD HTML 4.01//EN" "http://www.w3.org/TR/html4/strict.dtd"> +<html> +<head> +<meta http-equiv="content-type" content="text/html; charset=UTF-8"> +<title>~/code/numeric/exam/Makefile.html</title> +<meta name="Generator" content="Vim/7.3"> +<meta name="plugin-version" content="vim7.3_v10"> +<meta name="syntax" content="make"> +<meta name="settings" content="number_lines,use_css,pre_wrap,expand_tabs"> +<style type="text/css"> +<!-- +pre { white-space: pre-wrap; font-family: monospace; color: #ffffff; background-color: #000000; } +body { font-family: monospace; color: #ffffff; background-color: #000000; } +.lnr { color: #ffff00; } +.Special { color: #ff40ff; } +.Constant { color: #ffff00; } +.Statement { color: #ffff00; } +.Identifier { color: #00ffff; } +.Comment { color: #00ffff; } +--> +</style> +</head> +<body> +<pre> +<span class="lnr"> 1 </span><span class="Comment"># Define compiler</span> +<span class="lnr"> 2 </span><span class="Comment">#CXX=g++</span> +<span class="lnr"> 3 </span> +<span class="lnr"> 4 </span><span class="Comment"># Define compiler flags (show all warnings)</span> +<span class="lnr"> 5 </span><span class="Identifier">CXXFLAGS</span>=-Wall +<span class="lnr"> 6 </span><span class="Comment">#CXXFLAGS=-std=c++0x</span> +<span class="lnr"> 7 </span> +<span class="lnr"> 8 </span><span class="Comment"># Define linker flags</span> +<span class="lnr"> 9 </span><span class="Comment">#LDFLAGS=-fopenmp</span> +<span class="lnr"> 10 </span> +<span class="lnr"> 11 </span><span class="Comment"># Compile optimized code</span> +<span class="lnr"> 12 </span><span class="Identifier">CXXFLAGS</span>+=-O2 +<span class="lnr"> 13 </span> +<span class="lnr"> 14 </span><span class="Comment"># Compile debuggable code</span> +<span class="lnr"> 15 </span><span class="Comment">#CXXFLAGS+=-g</span> +<span class="lnr"> 16 </span> +<span class="lnr"> 17 </span><span class="Comment"># Compile profilable code</span> +<span class="lnr"> 18 </span><span class="Comment">#CXXFLAGS+=-pg</span> +<span class="lnr"> 19 </span><span class="Comment">#LDFLAGS+=-pg</span> +<span class="lnr"> 20 </span> +<span class="lnr"> 21 </span><span class="Comment"># Define linker</span> +<span class="lnr"> 22 </span><span class="Identifier">LD</span>=g++ +<span class="lnr"> 23 </span> +<span class="lnr"> 24 </span><span class="Comment"># All source code files</span> +<span class="lnr"> 25 </span><span class="Identifier">SRC</span>=<span class="Identifier">$(</span><span class="Statement">shell</span><span class="Identifier"> ls *.cpp)</span> +<span class="lnr"> 26 </span> +<span class="lnr"> 27 </span> +<span class="lnr"> 28 </span><span class="Comment"># Filenames of source code</span> +<span class="lnr"> 29 </span><span class="Identifier">SHARED_SRC</span>=ode.cpp check.cpp +<span class="lnr"> 30 </span><span class="Identifier">SHARED_HEADERS</span>=typedefs.h ode.h functions.h check.h +<span class="lnr"> 31 </span><span class="Identifier">SRC_A</span>=mainA.cpp <span class="Identifier">$(SHARED_SRC)</span> +<span class="lnr"> 32 </span><span class="Identifier">HEAD_A</span>=<span class="Identifier">$(SHARED_HEADERS)</span> +<span class="lnr"> 33 </span><span class="Identifier">SRC_B</span>=mainB.cpp <span class="Identifier">$(SHARED_SRC)</span> +<span class="lnr"> 34 </span><span class="Identifier">HEAD_B</span>=<span class="Identifier">$(SHARED_HEADERS)</span> +<span class="lnr"> 35 </span><span class="Identifier">SRC_C</span>=mainC.cpp <span class="Identifier">$(SHARED_SRC)</span> +<span class="lnr"> 36 </span><span class="Identifier">HEAD_C</span>=<span class="Identifier">$(SHARED_HEADERS)</span> +<span class="lnr"> 37 </span><span class="Identifier">SRC_D</span>=mainD.cpp <span class="Identifier">$(SHARED_SRC)</span> +<span class="lnr"> 38 </span><span class="Identifier">HEAD_D</span>=<span class="Identifier">$(SHARED_HEADERS)</span> +<span class="lnr"> 39 </span> +<span class="lnr"> 40 </span><span class="Comment"># Filenames of object files</span> +<span class="lnr"> 41 </span><span class="Identifier">OBJ_A</span>=<span class="Identifier">$(SRC_A:.cpp=.o)</span> +<span class="lnr"> 42 </span><span class="Identifier">OBJ_B</span>=<span class="Identifier">$(SRC_B:.cpp=.o)</span> +<span class="lnr"> 43 </span><span class="Identifier">OBJ_C</span>=<span class="Identifier">$(SRC_C:.cpp=.o)</span> +<span class="lnr"> 44 </span><span class="Identifier">OBJ_D</span>=<span class="Identifier">$(SRC_D:.cpp=.o)</span> +<span class="lnr"> 45 </span> +<span class="lnr"> 46 </span><span class="Comment"># Remove file type extension for binary filename</span> +<span class="lnr"> 47 </span><span class="Identifier">BIN_A</span>=odeA +<span class="lnr"> 48 </span><span class="Identifier">BIN_B</span>=odeB +<span class="lnr"> 49 </span><span class="Identifier">BIN_C</span>=odeC +<span class="lnr"> 50 </span><span class="Identifier">BIN_D</span>=odeD +<span class="lnr"> 51 </span> +<span class="lnr"> 52 </span><span class="Comment"># Define editor and options for `make edit`</span> +<span class="lnr"> 53 </span><span class="Identifier">EDITOR</span>=vim -p +<span class="lnr"> 54 </span> +<span class="lnr"> 55 </span> +<span class="lnr"> 56 </span><span class="Comment"># The default "all" depends on A and B</span> +<span class="lnr"> 57 </span><span class="Identifier">all:</span> A B C D +<span class="lnr"> 58 </span> +<span class="lnr"> 59 </span><span class="Identifier">A:</span> plotA.png +<span class="lnr"> 60 </span> +<span class="lnr"> 61 </span><span class="Identifier">B:</span> plotB.png +<span class="lnr"> 62 </span> +<span class="lnr"> 63 </span><span class="Identifier">C:</span> plotC.png +<span class="lnr"> 64 </span> +<span class="lnr"> 65 </span><span class="Identifier">D:</span> plotD.png +<span class="lnr"> 66 </span> +<span class="lnr"> 67 </span><span class="Identifier">plotA.png:</span> funcA.dat plotA.gp +<span class="lnr"> 68 </span><span class="Constant"> </span><span class="Comment"># Gnuplot: plotA.png</span> +<span class="lnr"> 69 </span><span class="Special"> @</span><span class="Constant">gnuplot plotA.gp</span> +<span class="lnr"> 70 </span> +<span class="lnr"> 71 </span><span class="Identifier">plotB.png:</span> funcB.dat plotB.gp +<span class="lnr"> 72 </span><span class="Constant"> </span><span class="Comment"># Gnuplot: plotB.png</span> +<span class="lnr"> 73 </span><span class="Special"> @</span><span class="Constant">gnuplot plotB.gp</span> +<span class="lnr"> 74 </span> +<span class="lnr"> 75 </span><span class="Identifier">plotC.png:</span> funcC.dat plotC.gp +<span class="lnr"> 76 </span><span class="Constant"> </span><span class="Comment"># Gnuplot: plotC.png</span> +<span class="lnr"> 77 </span><span class="Special"> @</span><span class="Constant">gnuplot plotC.gp</span> +<span class="lnr"> 78 </span> +<span class="lnr"> 79 </span><span class="Identifier">plotD.png:</span> funcD.dat plotD.gp +<span class="lnr"> 80 </span><span class="Constant"> </span><span class="Comment"># Gnuplot: plotD.png</span> +<span class="lnr"> 81 </span><span class="Special"> @</span><span class="Constant">gnuplot plotD.gp</span> +<span class="lnr"> 82 </span> +<span class="lnr"> 83 </span><span class="Identifier">funcA.dat:</span> <span class="Identifier">$(BIN_A)</span> +<span class="lnr"> 84 </span><span class="Special"> @</span><span class="Constant">./</span><span class="Identifier">$(BIN_A)</span> +<span class="lnr"> 85 </span><span class="Constant"> </span> +<span class="lnr"> 86 </span><span class="Identifier">funcB.dat:</span> <span class="Identifier">$(BIN_B)</span> +<span class="lnr"> 87 </span><span class="Special"> @</span><span class="Constant">./</span><span class="Identifier">$(BIN_B)</span> +<span class="lnr"> 88 </span> +<span class="lnr"> 89 </span><span class="Identifier">funcC.dat:</span> <span class="Identifier">$(BIN_C)</span> +<span class="lnr"> 90 </span><span class="Special"> @</span><span class="Constant">./</span><span class="Identifier">$(BIN_C)</span> +<span class="lnr"> 91 </span><span class="Constant"> </span> +<span class="lnr"> 92 </span><span class="Identifier">funcD.dat:</span> <span class="Identifier">$(BIN_D)</span> +<span class="lnr"> 93 </span><span class="Special"> @</span><span class="Constant">./</span><span class="Identifier">$(BIN_D)</span> +<span class="lnr"> 94 </span><span class="Constant"> </span> +<span class="lnr"> 95 </span><span class="Identifier">$(BIN_A):</span> <span class="Identifier">$(OBJ_A)</span> <span class="Identifier">$(HEAD_A)</span> +<span class="lnr"> 96 </span><span class="Special"> @</span><span class="Comment"># Link object files together</span> +<span class="lnr"> 97 </span><span class="Constant"> </span><span class="Identifier">$(LD)</span><span class="Constant"> </span><span class="Identifier">$(LDFLAGS)</span><span class="Constant"> </span><span class="Identifier">$(OBJ_A)</span><span class="Constant"> -o </span><span class="Identifier">$@</span><span class="Constant"> </span><span class="Identifier">$(LDLIBS)</span> +<span class="lnr"> 98 </span><span class="Constant"> </span> +<span class="lnr"> 99 </span><span class="Identifier">$(BIN_B):</span> <span class="Identifier">$(OBJ_B)</span> <span class="Identifier">$(HEAD_B)</span> +<span class="lnr">100 </span><span class="Special"> @</span><span class="Comment"># Link object files together</span> +<span class="lnr">101 </span><span class="Constant"> </span><span class="Identifier">$(LD)</span><span class="Constant"> </span><span class="Identifier">$(LDFLAGS)</span><span class="Constant"> </span><span class="Identifier">$(OBJ_B)</span><span class="Constant"> -o </span><span class="Identifier">$@</span><span class="Constant"> </span><span class="Identifier">$(LDLIBS)</span> +<span class="lnr">102 </span> +<span class="lnr">103 </span><span class="Identifier">$(BIN_C):</span> <span class="Identifier">$(OBJ_C)</span> <span class="Identifier">$(HEAD_C)</span> +<span class="lnr">104 </span><span class="Special"> @</span><span class="Comment"># Link object files together</span> +<span class="lnr">105 </span><span class="Constant"> </span><span class="Identifier">$(LD)</span><span class="Constant"> </span><span class="Identifier">$(LDFLAGS)</span><span class="Constant"> </span><span class="Identifier">$(OBJ_C)</span><span class="Constant"> -o </span><span class="Identifier">$@</span><span class="Constant"> </span><span class="Identifier">$(LDLIBS)</span> +<span class="lnr">106 </span> +<span class="lnr">107 </span><span class="Identifier">$(BIN_D):</span> <span class="Identifier">$(OBJ_D)</span> <span class="Identifier">$(HEAD_D)</span> +<span class="lnr">108 </span><span class="Special"> @</span><span class="Comment"># Link object files together</span> +<span class="lnr">109 </span><span class="Constant"> </span><span class="Identifier">$(LD)</span><span class="Constant"> </span><span class="Identifier">$(LDFLAGS)</span><span class="Constant"> </span><span class="Identifier">$(OBJ_D)</span><span class="Constant"> -o </span><span class="Identifier">$@</span><span class="Constant"> </span><span class="Identifier">$(LDLIBS)</span> +<span class="lnr">110 </span> +<span class="lnr">111 </span><span class="Identifier">clean:</span> cleanA cleanB cleanC cleanD +<span class="lnr">112 </span> +<span class="lnr">113 </span><span class="Identifier">cleanA: </span> +<span class="lnr">114 </span><span class="Special"> @</span><span class="Comment"># Remove object files</span> +<span class="lnr">115 </span><span class="Constant"> rm -f </span><span class="Identifier">$(OBJ_A)</span><span class="Constant"> </span> +<span class="lnr">116 </span><span class="Special"> @</span><span class="Comment"># Remove binaries</span> +<span class="lnr">117 </span><span class="Constant"> rm -f </span><span class="Identifier">$(BIN_A)</span> +<span class="lnr">118 </span><span class="Special"> @</span><span class="Comment"># Remove datafiles and plot</span> +<span class="lnr">119 </span><span class="Constant"> rm -f funcA.dat plotA.png</span> +<span class="lnr">120 </span> +<span class="lnr">121 </span><span class="Identifier">cleanB: </span> +<span class="lnr">122 </span><span class="Special"> @</span><span class="Comment"># Remove object files</span> +<span class="lnr">123 </span><span class="Constant"> rm -f </span><span class="Identifier">$(OBJ_B)</span><span class="Constant"> </span> +<span class="lnr">124 </span><span class="Special"> @</span><span class="Comment"># Remove binaries</span> +<span class="lnr">125 </span><span class="Constant"> rm -f </span><span class="Identifier">$(BIN_B)</span> +<span class="lnr">126 </span><span class="Special"> @</span><span class="Comment"># Remove datafiles and plot</span> +<span class="lnr">127 </span><span class="Constant"> rm -f funcB.dat plotB.png</span> +<span class="lnr">128 </span> +<span class="lnr">129 </span><span class="Identifier">cleanC: </span> +<span class="lnr">130 </span><span class="Special"> @</span><span class="Comment"># Remove object files</span> +<span class="lnr">131 </span><span class="Constant"> rm -f </span><span class="Identifier">$(OBJ_C)</span><span class="Constant"> </span> +<span class="lnr">132 </span><span class="Special"> @</span><span class="Comment"># Remove binaries</span> +<span class="lnr">133 </span><span class="Constant"> rm -f </span><span class="Identifier">$(BIN_C)</span> +<span class="lnr">134 </span><span class="Special"> @</span><span class="Comment"># Remove datafiles and plot</span> +<span class="lnr">135 </span><span class="Constant"> rm -f funcC.dat plotC.png</span> +<span class="lnr">136 </span> +<span class="lnr">137 </span><span class="Identifier">cleanD: </span> +<span class="lnr">138 </span><span class="Special"> @</span><span class="Comment"># Remove object files</span> +<span class="lnr">139 </span><span class="Constant"> rm -f </span><span class="Identifier">$(OBJ_D)</span><span class="Constant"> </span> +<span class="lnr">140 </span><span class="Special"> @</span><span class="Comment"># Remove binaries</span> +<span class="lnr">141 </span><span class="Constant"> rm -f </span><span class="Identifier">$(BIN_D)</span> +<span class="lnr">142 </span><span class="Constant"> </span><span class="Comment"># Removing datafile and plot</span> +<span class="lnr">143 </span><span class="Special"> @</span><span class="Constant">rm -f funcD.dat plotD.png</span> +<span class="lnr">144 </span> +<span class="lnr">145 </span><span class="Identifier">htmlfiles:</span> html/mainA.cpp.html html/mainB.cpp.html html/mainC.cpp.html html/mainD.cpp.html html/ode.cpp.html html/check.cpp.html html/check.h.html html/functions.h.html html/ode.h.html html/typedefs.h.html html/vector_arithmetic.h.html html/plotA.gp.html html/plotB.gp.html html/plotC.gp.html html/plotD.gp.html html/Makefile.html +<span class="lnr">146 </span><span class="Constant"> </span><span class="Comment"># Generating HTML files</span> +<span class="lnr">147 </span><span class="Constant"> rst2html2 README.rst > html/README.html</span> +<span class="lnr">148 </span> +<span class="lnr">149 </span><span class="Identifier">html/%.html:</span> <span class="Identifier">%</span> +<span class="lnr">150 </span><span class="Constant"> vim </span><span class="Identifier">$<</span><span class="Constant"> +TOhtml +</span><span class="Constant">"w </span><span class="Identifier">$@</span><span class="Constant">"</span><span class="Constant"> +</span><span class="Constant">"qall!"</span> +<span class="lnr">151 </span> +<span class="lnr">152 </span> +<span class="lnr">153 </span><span class="Identifier">edit:</span> +<span class="lnr">154 </span><span class="Special"> @</span><span class="Identifier">$(EDITOR)</span><span class="Constant"> Makefile README.rst *.cpp *.h *.gp</span> +<span class="lnr">155 </span> +</pre> +</body> +</html> diff --git a/exam/html/README.html b/exam/html/README.html @@ -0,0 +1,406 @@ +<?xml version="1.0" encoding="utf-8" ?> +<!DOCTYPE html PUBLIC "-//W3C//DTD XHTML 1.0 Transitional//EN" "http://www.w3.org/TR/xhtml1/DTD/xhtml1-transitional.dtd"> +<html xmlns="http://www.w3.org/1999/xhtml" xml:lang="en" lang="en"> +<head> +<meta http-equiv="Content-Type" content="text/html; charset=utf-8" /> +<meta name="generator" content="Docutils 0.9.1: http://docutils.sourceforge.net/" /> +<title>README: ODE integration with complex numbers</title> +<style type="text/css"> + +/* +:Author: David Goodger (goodger@python.org) +:Id: $Id: html4css1.css 7434 2012-05-11 21:06:27Z milde $ +:Copyright: This stylesheet has been placed in the public domain. + +Default cascading style sheet for the HTML output of Docutils. + +See http://docutils.sf.net/docs/howto/html-stylesheets.html for how to +customize this style sheet. +*/ + +/* used to remove borders from tables and images */ +.borderless, table.borderless td, table.borderless th { + border: 0 } + +table.borderless td, table.borderless th { + /* Override padding for "table.docutils td" with "! important". + The right padding separates the table cells. */ + padding: 0 0.5em 0 0 ! important } + +.first { + /* Override more specific margin styles with "! important". */ + margin-top: 0 ! important } + +.last, .with-subtitle { + margin-bottom: 0 ! important } + +.hidden { + display: none } + +a.toc-backref { + text-decoration: none ; + color: black } + +blockquote.epigraph { + margin: 2em 5em ; } + +dl.docutils dd { + margin-bottom: 0.5em } + +object[type="image/svg+xml"], object[type="application/x-shockwave-flash"] { + overflow: hidden; +} + +/* Uncomment (and remove this text!) to get bold-faced definition list terms +dl.docutils dt { + font-weight: bold } +*/ + +div.abstract { + margin: 2em 5em } + +div.abstract p.topic-title { + font-weight: bold ; + text-align: center } + +div.admonition, div.attention, div.caution, div.danger, div.error, +div.hint, div.important, div.note, div.tip, div.warning { + margin: 2em ; + border: medium outset ; + padding: 1em } + +div.admonition p.admonition-title, div.hint p.admonition-title, +div.important p.admonition-title, div.note p.admonition-title, +div.tip p.admonition-title { + font-weight: bold ; + font-family: sans-serif } + +div.attention p.admonition-title, div.caution p.admonition-title, +div.danger p.admonition-title, div.error p.admonition-title, +div.warning p.admonition-title { + color: red ; + font-weight: bold ; + font-family: sans-serif } + +/* Uncomment (and remove this text!) to get reduced vertical space in + compound paragraphs. +div.compound .compound-first, div.compound .compound-middle { + margin-bottom: 0.5em } + +div.compound .compound-last, div.compound .compound-middle { + margin-top: 0.5em } +*/ + +div.dedication { + margin: 2em 5em ; + text-align: center ; + font-style: italic } + +div.dedication p.topic-title { + font-weight: bold ; + font-style: normal } + +div.figure { + margin-left: 2em ; + margin-right: 2em } + +div.footer, div.header { + clear: both; + font-size: smaller } + +div.line-block { + display: block ; + margin-top: 1em ; + margin-bottom: 1em } + +div.line-block div.line-block { + margin-top: 0 ; + margin-bottom: 0 ; + margin-left: 1.5em } + +div.sidebar { + margin: 0 0 0.5em 1em ; + border: medium outset ; + padding: 1em ; + background-color: #ffffee ; + width: 40% ; + float: right ; + clear: right } + +div.sidebar p.rubric { + font-family: sans-serif ; + font-size: medium } + +div.system-messages { + margin: 5em } + +div.system-messages h1 { + color: red } + +div.system-message { + border: medium outset ; + padding: 1em } + +div.system-message p.system-message-title { + color: red ; + font-weight: bold } + +div.topic { + margin: 2em } + +h1.section-subtitle, h2.section-subtitle, h3.section-subtitle, +h4.section-subtitle, h5.section-subtitle, h6.section-subtitle { + margin-top: 0.4em } + +h1.title { + text-align: center } + +h2.subtitle { + text-align: center } + +hr.docutils { + width: 75% } + +img.align-left, .figure.align-left, object.align-left { + clear: left ; + float: left ; + margin-right: 1em } + +img.align-right, .figure.align-right, object.align-right { + clear: right ; + float: right ; + margin-left: 1em } + +img.align-center, .figure.align-center, object.align-center { + display: block; + margin-left: auto; + margin-right: auto; +} + +.align-left { + text-align: left } + +.align-center { + clear: both ; + text-align: center } + +.align-right { + text-align: right } + +/* reset inner alignment in figures */ +div.align-right { + text-align: inherit } + +/* div.align-center * { */ +/* text-align: left } */ + +ol.simple, ul.simple { + margin-bottom: 1em } + +ol.arabic { + list-style: decimal } + +ol.loweralpha { + list-style: lower-alpha } + +ol.upperalpha { + list-style: upper-alpha } + +ol.lowerroman { + list-style: lower-roman } + +ol.upperroman { + list-style: upper-roman } + +p.attribution { + text-align: right ; + margin-left: 50% } + +p.caption { + font-style: italic } + +p.credits { + font-style: italic ; + font-size: smaller } + +p.label { + white-space: nowrap } + +p.rubric { + font-weight: bold ; + font-size: larger ; + color: maroon ; + text-align: center } + +p.sidebar-title { + font-family: sans-serif ; + font-weight: bold ; + font-size: larger } + +p.sidebar-subtitle { + font-family: sans-serif ; + font-weight: bold } + +p.topic-title { + font-weight: bold } + +pre.address { + margin-bottom: 0 ; + margin-top: 0 ; + font: inherit } + +pre.literal-block, pre.doctest-block, pre.math, pre.code { + margin-left: 2em ; + margin-right: 2em } + +pre.code .ln { /* line numbers */ + color: grey; +} + +.code { + background-color: #eeeeee +} + +span.classifier { + font-family: sans-serif ; + font-style: oblique } + +span.classifier-delimiter { + font-family: sans-serif ; + font-weight: bold } + +span.interpreted { + font-family: sans-serif } + +span.option { + white-space: nowrap } + +span.pre { + white-space: pre } + +span.problematic { + color: red } + +span.section-subtitle { + /* font-size relative to parent (h1..h6 element) */ + font-size: 80% } + +table.citation { + border-left: solid 1px gray; + margin-left: 1px } + +table.docinfo { + margin: 2em 4em } + +table.docutils { + margin-top: 0.5em ; + margin-bottom: 0.5em } + +table.footnote { + border-left: solid 1px black; + margin-left: 1px } + +table.docutils td, table.docutils th, +table.docinfo td, table.docinfo th { + padding-left: 0.5em ; + padding-right: 0.5em ; + vertical-align: top } + +table.docutils th.field-name, table.docinfo th.docinfo-name { + font-weight: bold ; + text-align: left ; + white-space: nowrap ; + padding-left: 0 } + +h1 tt.docutils, h2 tt.docutils, h3 tt.docutils, +h4 tt.docutils, h5 tt.docutils, h6 tt.docutils { + font-size: 100% } + +ul.auto-toc { + list-style-type: none } + +</style> +</head> +<body> +<div class="document" id="readme-ode-integration-with-complex-numbers"> +<h1 class="title">README: ODE integration with complex numbers</h1> + +<p>Exam exercise for <em>Numerical Methods</em> by Anders D. Christensen (<a class="reference external" href="mailto:adc@geo.au.dk">mail</a>)</p> +<div class="section" id="file-description"> +<h1>File description</h1> +<ul class="simple"> +<li><tt class="docutils literal">Makefile</tt>: Description for GNU Make, handles compilation and execution.</li> +<li><tt class="docutils literal">README.rst</tt> (this file): Description of numeric implementation and usage. +Written with reStructuredText syntax.</li> +<li><tt class="docutils literal">check.cpp</tt>: Function for displaying the state of a condition to stdout.</li> +<li><tt class="docutils literal">check.h</tt>: Prototype for the check-function.</li> +<li><tt class="docutils literal">functions.h</tt>: Input functions to be evaluated.</li> +<li><tt class="docutils literal">mainA.cpp</tt>: Main source code file for part A.</li> +<li><tt class="docutils literal">ode.cpp</tt>: Constructor and functions for the ODE class, including Runge-Kutta +stepper and driver.</li> +<li><tt class="docutils literal">ode.h</tt>: Header file with the ODE class. This file must be included in all +programs that want to utilize the ODE functionality.</li> +<li><tt class="docutils literal">plot.gp</tt>: Script for plotting all graphs with Gnuplot.</li> +<li><tt class="docutils literal">typedefs.h</tt>: Header file containing definitions of two main types, +<tt class="docutils literal">Inttype</tt>, a whole-number type, and <tt class="docutils literal">Floattype</tt>, a floating point number +type. The type definitions can be changed to different lengths and precisions. +The program can be compiled for verbose output by changing the <tt class="docutils literal">verbose</tt> +variable.</li> +<li><tt class="docutils literal">vector_arithmetic.h</tt>: Operator overloading functions for the <tt class="docutils literal"><span class="pre">std::vector</span></tt> +class.</li> +</ul> +</div> +<div class="section" id="problem-descriptions"> +<h1>Problem descriptions</h1> +<p>The four generated executables each demonstrate the ODE solvers functionality by +performing the following tasks. The results consist of the console output and +the corresponding plot with filename <tt class="docutils literal"><span class="pre">plot<Character>.png</span></tt>. +- <em>A</em>: Construct an ODE solver that can handle functions with complex values. +Demonstrate that it solves the real component correctly, by stepping along +a path in the real range. +- <em>B</em>: Demonstrate that the ODE solver can solve the imaginary component by +stepping along a path in the imaginary range. +- <em>C</em>: Demonstrate the solution of a set of complex equations by stepping +through the complex plane. +- <em>D</em>: For an integration path in the complex plane, visualize how the +requirements of absolute- and relative precision are related to the number of +integration steps, for a given floating point precision.</p> +</div> +<div class="section" id="implementation"> +<h1>Implementation</h1> +<p>This exercise was written in object-oriented C++ for easy reuse. For +portability, all included classes are from the standard template library.</p> +<p>The necessary <tt class="docutils literal"><span class="pre">std::vector</span></tt> arithmetic operations where overloaded to support +element-wise operations, such as vector-scalar multiplication, vector-vector +addition, etc. This approach was preferred over using <tt class="docutils literal"><span class="pre">std::valarray</span></tt>, since it +is not dynamically expandable.</p> +<p>When creating a new ODE object, the user specifies the end-points of the linear +range, where the specified system of ordinary differential equations with +complex values will be solved. The range end-points are complex numbers +themselves, and the user can thus specify whether the integrator steps through a +range of real values, imaginary values, or both components in the complex plane.</p> +<p>The solver steps through the specified range by an adaptive step size, which is +also a complex number. The user specifies the fraction of the range to be used +as a start value for the step. The default value is 0.01.</p> +<p>The ODE class contains functions for writing the ODE solution to stdout +(<tt class="docutils literal"><span class="pre">ODE::print</span></tt>) or to a text file (<tt class="docutils literal"><span class="pre">ODE::write</span></tt>). The output format is the +following; the first column is the real part of x, second column the imaginary +part. The subsequent columns do in turn consist of real- and imaginary parts of +the variables in the ODE.</p> +<p>The program requires a modern C++ compiler, GNU Make and Gnuplot. It has been +tested with GCC, Clang and llvm.</p> +</div> +<div class="section" id="compiliation-and-execution"> +<h1>Compiliation and execution</h1> +<p>To make and execute the program, go to the root folder and type <cite>make</cite>. This +will compile and execute the programs for part A-D, and plot output graphs. If +desired, individual parts can be compiled and executed using <cite>make <Character></cite>.</p> +<p>To view the source code in a browser with vim's syntax highlighting, type <cite>make +html</cite>, and view the files in the <cite>html</cite> folder. The generation of HTML files +requires a newer vim for the source code files, and Docutils for the readme.</p> +<p>All output and objects can be removed using <cite>make clean</cite>.</p> +<p>#vim: set tw=80</p> +</div> +</div> +</body> +</html> diff --git a/exam/html/check.cpp.html b/exam/html/check.cpp.html @@ -0,0 +1,39 @@ +<!DOCTYPE HTML PUBLIC "-//W3C//DTD HTML 4.01//EN" "http://www.w3.org/TR/html4/strict.dtd"> +<html> +<head> +<meta http-equiv="content-type" content="text/html; charset=UTF-8"> +<title>~/code/numeric/exam/check.cpp.html</title> +<meta name="Generator" content="Vim/7.3"> +<meta name="plugin-version" content="vim7.3_v10"> +<meta name="syntax" content="cpp"> +<meta name="settings" content="number_lines,use_css,pre_wrap,expand_tabs"> +<style type="text/css"> +<!-- +pre { white-space: pre-wrap; font-family: monospace; color: #ffffff; background-color: #000000; } +body { font-family: monospace; color: #ffffff; background-color: #000000; } +.lnr { color: #ffff00; } +.Special { color: #ff40ff; } +.Statement { color: #ffff00; } +.Type { color: #00ff00; } +.Constant { color: #ffff00; } +.PreProc { color: #ff40ff; } +.Comment { color: #00ffff; } +--> +</style> +</head> +<body> +<pre> +<span class="lnr"> 1 </span><span class="Comment">// Function used for reporting the condition of a </span> +<span class="lnr"> 2 </span><span class="Comment">// statement to stdout using ANSI colors.</span> +<span class="lnr"> 3 </span><span class="PreProc">#include </span><span class="Constant"><iostream></span> +<span class="lnr"> 4 </span> +<span class="lnr"> 5 </span><span class="Type">void</span> check(<span class="Type">const</span> <span class="Type">bool</span> statement) +<span class="lnr"> 6 </span>{ +<span class="lnr"> 7 </span> <span class="Statement">if</span> (statement == <span class="Constant">true</span>) +<span class="lnr"> 8 </span> std::cout << <span class="Constant">"</span><span class="Special">\t\033</span><span class="Constant">[0;32mPassed</span><span class="Special">\033</span><span class="Constant">[0m</span><span class="Special">\n</span><span class="Constant">"</span>; +<span class="lnr"> 9 </span> <span class="Statement">else</span> +<span class="lnr">10 </span> std::cout << <span class="Constant">"</span><span class="Special">\t\033</span><span class="Constant">[1;31mFail!!</span><span class="Special">\033</span><span class="Constant">[0m</span><span class="Special">\n</span><span class="Constant">"</span>; +<span class="lnr">11 </span>} +</pre> +</body> +</html> diff --git a/exam/html/check.h.html b/exam/html/check.h.html @@ -0,0 +1,32 @@ +<!DOCTYPE HTML PUBLIC "-//W3C//DTD HTML 4.01//EN" "http://www.w3.org/TR/html4/strict.dtd"> +<html> +<head> +<meta http-equiv="content-type" content="text/html; charset=UTF-8"> +<title>~/code/numeric/exam/check.h.html</title> +<meta name="Generator" content="Vim/7.3"> +<meta name="plugin-version" content="vim7.3_v10"> +<meta name="syntax" content="cpp"> +<meta name="settings" content="number_lines,use_css,pre_wrap,expand_tabs"> +<style type="text/css"> +<!-- +pre { white-space: pre-wrap; font-family: monospace; color: #ffffff; background-color: #000000; } +body { font-family: monospace; color: #ffffff; background-color: #000000; } +.lnr { color: #ffff00; } +.Type { color: #00ff00; } +.Comment { color: #00ffff; } +.PreProc { color: #ff40ff; } +--> +</style> +</head> +<body> +<pre> +<span class="lnr">1 </span><span class="PreProc">#ifndef CHECK_H_</span> +<span class="lnr">2 </span><span class="PreProc">#define CHECK_H_</span> +<span class="lnr">3 </span> +<span class="lnr">4 </span><span class="Comment">// Prototype for checking function</span> +<span class="lnr">5 </span><span class="Type">void</span> check(<span class="Type">const</span> <span class="Type">bool</span> statement); +<span class="lnr">6 </span> +<span class="lnr">7 </span><span class="PreProc">#endif</span> +</pre> +</body> +</html> diff --git a/exam/html/functions.h.html b/exam/html/functions.h.html @@ -0,0 +1,52 @@ +<!DOCTYPE HTML PUBLIC "-//W3C//DTD HTML 4.01//EN" "http://www.w3.org/TR/html4/strict.dtd"> +<html> +<head> +<meta http-equiv="content-type" content="text/html; charset=UTF-8"> +<title>~/code/numeric/exam/functions.h.html</title> +<meta name="Generator" content="Vim/7.3"> +<meta name="plugin-version" content="vim7.3_v10"> +<meta name="syntax" content="cpp"> +<meta name="settings" content="number_lines,use_css,pre_wrap,expand_tabs"> +<style type="text/css"> +<!-- +pre { white-space: pre-wrap; font-family: monospace; color: #ffffff; background-color: #000000; } +body { font-family: monospace; color: #ffffff; background-color: #000000; } +.lnr { color: #ffff00; } +.Statement { color: #ffff00; } +.Type { color: #00ff00; } +.Constant { color: #ffff00; } +.PreProc { color: #ff40ff; } +.Comment { color: #00ffff; } +--> +</style> +</head> +<body> +<pre> +<span class="lnr"> 1 </span><span class="Comment">// Make sure file is only included once per object</span> +<span class="lnr"> 2 </span><span class="PreProc">#ifndef FUNCTIONS_H_</span> +<span class="lnr"> 3 </span><span class="PreProc">#define FUNCTIONS_H_</span> +<span class="lnr"> 4 </span> +<span class="lnr"> 5 </span><span class="PreProc">#include </span><span class="Constant"><vector></span> +<span class="lnr"> 6 </span><span class="PreProc">#include </span><span class="Constant"><complex></span> +<span class="lnr"> 7 </span><span class="PreProc">#include </span><span class="Constant">"typedefs.h"</span> +<span class="lnr"> 8 </span> +<span class="lnr"> 9 </span> +<span class="lnr">10 </span><span class="Comment">//// ODEs with real+complex parts.</span> +<span class="lnr">11 </span><span class="Comment">//// Return the derivatives at the point x,vec(y)</span> +<span class="lnr">12 </span> +<span class="lnr">13 </span>std::vector<std::complex<Floattype> > +<span class="lnr">14 </span> func1(<span class="Type">const</span> std::complex<Floattype> z, +<span class="lnr">15 </span> <span class="Type">const</span> std::vector<std::complex<Floattype> > &y) +<span class="lnr">16 </span>{ +<span class="lnr">17 </span> std::vector<std::complex<Floattype> > dydz(<span class="Constant">2</span>); +<span class="lnr">18 </span> dydz[<span class="Constant">0</span>].real() = y[<span class="Constant">1</span>].real(); +<span class="lnr">19 </span> dydz[<span class="Constant">0</span>].imag() = y[<span class="Constant">1</span>].imag(); +<span class="lnr">20 </span> dydz[<span class="Constant">1</span>].real() = -y[<span class="Constant">0</span>].real(); +<span class="lnr">21 </span> dydz[<span class="Constant">1</span>].imag() = <span class="Constant">0.5f</span>*y[<span class="Constant">0</span>].imag(); +<span class="lnr">22 </span> <span class="Statement">return</span> dydz; +<span class="lnr">23 </span>} +<span class="lnr">24 </span> +<span class="lnr">25 </span><span class="PreProc">#endif</span> +</pre> +</body> +</html> diff --git a/exam/html/mainA.cpp.html b/exam/html/mainA.cpp.html @@ -0,0 +1,87 @@ +<!DOCTYPE HTML PUBLIC "-//W3C//DTD HTML 4.01//EN" "http://www.w3.org/TR/html4/strict.dtd"> +<html> +<head> +<meta http-equiv="content-type" content="text/html; charset=UTF-8"> +<title>~/code/numeric/exam/mainA.cpp.html</title> +<meta name="Generator" content="Vim/7.3"> +<meta name="plugin-version" content="vim7.3_v10"> +<meta name="syntax" content="cpp"> +<meta name="settings" content="number_lines,use_css,pre_wrap,expand_tabs"> +<style type="text/css"> +<!-- +pre { white-space: pre-wrap; font-family: monospace; color: #ffffff; background-color: #000000; } +body { font-family: monospace; color: #ffffff; background-color: #000000; } +.lnr { color: #ffff00; } +.Special { color: #ff40ff; } +.Statement { color: #ffff00; } +.Comment { color: #00ffff; } +.Type { color: #00ff00; } +.Constant { color: #ffff00; } +.PreProc { color: #ff40ff; } +--> +</style> +</head> +<body> +<pre> +<span class="lnr"> 1 </span><span class="PreProc">#include </span><span class="Constant"><iostream></span> +<span class="lnr"> 2 </span><span class="PreProc">#include </span><span class="Constant"><vector></span> +<span class="lnr"> 3 </span><span class="PreProc">#include </span><span class="Constant"><complex></span> +<span class="lnr"> 4 </span><span class="PreProc">#include </span><span class="Constant"><cmath></span> +<span class="lnr"> 5 </span><span class="PreProc">#include </span><span class="Constant">"typedefs.h"</span> +<span class="lnr"> 6 </span><span class="PreProc">#include </span><span class="Constant">"check.h"</span> +<span class="lnr"> 7 </span><span class="PreProc">#include </span><span class="Constant">"ode.h"</span> +<span class="lnr"> 8 </span><span class="PreProc">#include </span><span class="Constant">"functions.h"</span> +<span class="lnr"> 9 </span> +<span class="lnr">10 </span> +<span class="lnr">11 </span><span class="Type">int</span> main() +<span class="lnr">12 </span>{ +<span class="lnr">13 </span> <span class="Comment">// Namespace declarations</span> +<span class="lnr">14 </span> <span class="Statement">using</span> std::cout; +<span class="lnr">15 </span> <span class="Statement">using</span> std::vector; +<span class="lnr">16 </span> <span class="Statement">using</span> std::complex; +<span class="lnr">17 </span> +<span class="lnr">18 </span> <span class="Comment">// Calculate machine precision</span> +<span class="lnr">19 </span> Floattype eps_machine = <span class="Constant">1.0f</span>; +<span class="lnr">20 </span> <span class="Statement">while</span> (<span class="Constant">1.0f</span> + eps_machine != <span class="Constant">1.0f</span>) +<span class="lnr">21 </span> eps_machine /= <span class="Constant">2.0f</span>; +<span class="lnr">22 </span> +<span class="lnr">23 </span> <span class="Type">const</span> <span class="Type">int</span> id = <span class="Constant">20062213</span>; +<span class="lnr">24 </span> <span class="Type">const</span> <span class="Type">char</span> n = <span class="Constant">10</span>; +<span class="lnr">25 </span> cout << <span class="Constant">"</span><span class="Special">\n</span><span class="Constant">My student id is </span><span class="Special">\033</span><span class="Constant">[1;37m"</span> << id +<span class="lnr">26 </span> << <span class="Constant">"</span><span class="Special">\033</span><span class="Constant">[0m, resulting in exam exercise: </span><span class="Special">\033</span><span class="Constant">[1;31m"</span> +<span class="lnr">27 </span> << id%n << <span class="Constant">"</span><span class="Special">\033</span><span class="Constant">[0m</span><span class="Special">\n</span><span class="Constant">"</span>; +<span class="lnr">28 </span> cout << <span class="Constant">"Examination project:</span><span class="Special">\033</span><span class="Constant">[1;37m ODE integration "</span> +<span class="lnr">29 </span> << <span class="Constant">"with complex numbers</span><span class="Special">\033</span><span class="Constant">[0m</span><span class="Special">\n\n</span><span class="Constant">"</span>; +<span class="lnr">30 </span> +<span class="lnr">31 </span> cout << <span class="Constant">"</span><span class="Special">\033</span><span class="Constant">[1;33m--- Part A: Solving along a real path ---</span><span class="Special">\033</span><span class="Constant">[0m</span><span class="Special">\n</span><span class="Constant">"</span>; +<span class="lnr">32 </span> <span class="Type">complex</span><Floattype> a(<span class="Constant">0.0f</span>, <span class="Constant">0.0f</span>); <span class="Comment">// Lower limit</span> +<span class="lnr">33 </span> <span class="Type">complex</span><Floattype> b(<span class="Constant">2.0f</span>*<span class="Constant">M_PI</span>, <span class="Constant">0.0f</span>); <span class="Comment">// Upper limit</span> +<span class="lnr">34 </span> cout << <span class="Constant">"Integration path: b-a = "</span> << b-a << <span class="Special">'\n'</span>; +<span class="lnr">35 </span> Inttype n_eqs = <span class="Constant">2</span>; <span class="Comment">// Number of equations in ODE system</span> +<span class="lnr">36 </span> vector<<span class="Type">complex</span><Floattype> > y_start(n_eqs); +<span class="lnr">37 </span> <span class="Type">complex</span><Floattype> y0(<span class="Constant">0.0f</span>, <span class="Constant">0.0f</span>); +<span class="lnr">38 </span> <span class="Type">complex</span><Floattype> y1(<span class="Constant">1.0f</span>, <span class="Constant">1.0f</span>); +<span class="lnr">39 </span> y_start[<span class="Constant">0</span>] = y0; +<span class="lnr">40 </span> y_start[<span class="Constant">1</span>] = y1; +<span class="lnr">41 </span> Floattype h_start = <span class="Constant">0.01f</span>; +<span class="lnr">42 </span> ODE realode(func1, <span class="Comment">// ODE system</span> +<span class="lnr">43 </span> y_start, <span class="Comment">// Initial values</span> +<span class="lnr">44 </span> a, <span class="Comment">// Lower limit</span> +<span class="lnr">45 </span> b, <span class="Comment">// Upper limit</span> +<span class="lnr">46 </span> h_start, <span class="Comment">// Start value of step size</span> +<span class="lnr">47 </span> <span class="Constant">10000</span>, <span class="Comment">// Max. number of steps</span> +<span class="lnr">48 </span> eps_machine*<span class="Constant">1e12f</span>, <span class="Comment">// Absolute precision</span> +<span class="lnr">49 </span> eps_machine*<span class="Constant">1e12f</span>); <span class="Comment">// Relative precision</span> +<span class="lnr">50 </span> realode.write(<span class="Constant">"funcA.dat"</span>); <span class="Comment">// Write solutions to data file</span> +<span class="lnr">51 </span> +<span class="lnr">52 </span> <span class="Comment">// Report to stdout</span> +<span class="lnr">53 </span> cout << <span class="Constant">"ODE system solved in "</span> +<span class="lnr">54 </span> << realode.steps() << <span class="Constant">" steps.</span><span class="Special">\n\n</span><span class="Constant">"</span>; +<span class="lnr">55 </span> +<span class="lnr">56 </span> <span class="Comment">// Return successfully</span> +<span class="lnr">57 </span> <span class="Statement">return</span> <span class="Constant">0</span>; +<span class="lnr">58 </span>} +<span class="lnr">59 </span> +</pre> +</body> +</html> diff --git a/exam/html/mainB.cpp.html b/exam/html/mainB.cpp.html @@ -0,0 +1,78 @@ +<!DOCTYPE HTML PUBLIC "-//W3C//DTD HTML 4.01//EN" "http://www.w3.org/TR/html4/strict.dtd"> +<html> +<head> +<meta http-equiv="content-type" content="text/html; charset=UTF-8"> +<title>~/code/numeric/exam/mainB.cpp.html</title> +<meta name="Generator" content="Vim/7.3"> +<meta name="plugin-version" content="vim7.3_v10"> +<meta name="syntax" content="cpp"> +<meta name="settings" content="number_lines,use_css,pre_wrap,expand_tabs"> +<style type="text/css"> +<!-- +pre { white-space: pre-wrap; font-family: monospace; color: #ffffff; background-color: #000000; } +body { font-family: monospace; color: #ffffff; background-color: #000000; } +.lnr { color: #ffff00; } +.Special { color: #ff40ff; } +.Statement { color: #ffff00; } +.Comment { color: #00ffff; } +.Type { color: #00ff00; } +.Constant { color: #ffff00; } +.PreProc { color: #ff40ff; } +--> +</style> +</head> +<body> +<pre> +<span class="lnr"> 1 </span><span class="PreProc">#include </span><span class="Constant"><iostream></span> +<span class="lnr"> 2 </span><span class="PreProc">#include </span><span class="Constant"><vector></span> +<span class="lnr"> 3 </span><span class="PreProc">#include </span><span class="Constant"><complex></span> +<span class="lnr"> 4 </span><span class="PreProc">#include </span><span class="Constant"><cmath></span> +<span class="lnr"> 5 </span><span class="PreProc">#include </span><span class="Constant">"typedefs.h"</span> +<span class="lnr"> 6 </span><span class="PreProc">#include </span><span class="Constant">"check.h"</span> +<span class="lnr"> 7 </span><span class="PreProc">#include </span><span class="Constant">"ode.h"</span> +<span class="lnr"> 8 </span><span class="PreProc">#include </span><span class="Constant">"functions.h"</span> +<span class="lnr"> 9 </span> +<span class="lnr">10 </span><span class="Type">int</span> main() +<span class="lnr">11 </span>{ +<span class="lnr">12 </span> <span class="Comment">// Namespace declarations</span> +<span class="lnr">13 </span> <span class="Statement">using</span> std::cout; +<span class="lnr">14 </span> <span class="Statement">using</span> std::vector; +<span class="lnr">15 </span> <span class="Statement">using</span> std::complex; +<span class="lnr">16 </span> +<span class="lnr">17 </span> <span class="Comment">// Calculate machine precision</span> +<span class="lnr">18 </span> Floattype eps_machine = <span class="Constant">1.0f</span>; +<span class="lnr">19 </span> <span class="Statement">while</span> (<span class="Constant">1.0f</span> + eps_machine != <span class="Constant">1.0f</span>) +<span class="lnr">20 </span> eps_machine /= <span class="Constant">2.0f</span>; +<span class="lnr">21 </span> +<span class="lnr">22 </span> cout << <span class="Constant">"</span><span class="Special">\n\033</span><span class="Constant">[1;33m--- Part B: Solving along an imaginary path ---</span><span class="Special">\033</span><span class="Constant">[0m</span><span class="Special">\n</span><span class="Constant">"</span>; +<span class="lnr">23 </span> <span class="Type">complex</span><Floattype> a(<span class="Constant">0.0f</span>, <span class="Constant">0.0f</span>); <span class="Comment">// Lower limit</span> +<span class="lnr">24 </span> <span class="Type">complex</span><Floattype> b(<span class="Constant">0.0f</span>, <span class="Constant">2.0f</span>*<span class="Constant">M_PI</span>); <span class="Comment">// Upper limit</span> +<span class="lnr">25 </span> cout << <span class="Constant">"Integration path: b-a = "</span> << b-a << <span class="Special">'\n'</span>; +<span class="lnr">26 </span> Inttype n_eqs = <span class="Constant">2</span>; <span class="Comment">// Number of equations in ODE system</span> +<span class="lnr">27 </span> vector<<span class="Type">complex</span><Floattype> > y_start(n_eqs); +<span class="lnr">28 </span> <span class="Type">complex</span><Floattype> y0(<span class="Constant">0.0f</span>, <span class="Constant">0.0f</span>); +<span class="lnr">29 </span> <span class="Type">complex</span><Floattype> y1(<span class="Constant">1.0f</span>, <span class="Constant">1.0f</span>); +<span class="lnr">30 </span> y_start[<span class="Constant">0</span>] = y0; +<span class="lnr">31 </span> y_start[<span class="Constant">1</span>] = y1; +<span class="lnr">32 </span> Floattype h_start = <span class="Constant">0.01f</span>; +<span class="lnr">33 </span> ODE imagode(func1, <span class="Comment">// ODE system</span> +<span class="lnr">34 </span> y_start, <span class="Comment">// Initial values</span> +<span class="lnr">35 </span> a, <span class="Comment">// Lower limit</span> +<span class="lnr">36 </span> b, <span class="Comment">// Upper limit</span> +<span class="lnr">37 </span> h_start, <span class="Comment">// Start value of step size</span> +<span class="lnr">38 </span> <span class="Constant">10000</span>, <span class="Comment">// Max. number of steps</span> +<span class="lnr">39 </span> eps_machine*<span class="Constant">1e12f</span>, <span class="Comment">// Absolute precision</span> +<span class="lnr">40 </span> eps_machine*<span class="Constant">1e12f</span>); <span class="Comment">// Relative precision</span> +<span class="lnr">41 </span> imagode.write(<span class="Constant">"funcB.dat"</span>); <span class="Comment">// Write solutions to data file</span> +<span class="lnr">42 </span> +<span class="lnr">43 </span> <span class="Comment">// Report to stdout</span> +<span class="lnr">44 </span> cout << <span class="Constant">"ODE system solved in "</span> +<span class="lnr">45 </span> << imagode.steps() << <span class="Constant">" steps.</span><span class="Special">\n\n</span><span class="Constant">"</span>; +<span class="lnr">46 </span> +<span class="lnr">47 </span> <span class="Comment">// Return successfully</span> +<span class="lnr">48 </span> <span class="Statement">return</span> <span class="Constant">0</span>; +<span class="lnr">49 </span>} +<span class="lnr">50 </span> +</pre> +</body> +</html> diff --git a/exam/html/mainC.cpp.html b/exam/html/mainC.cpp.html @@ -0,0 +1,78 @@ +<!DOCTYPE HTML PUBLIC "-//W3C//DTD HTML 4.01//EN" "http://www.w3.org/TR/html4/strict.dtd"> +<html> +<head> +<meta http-equiv="content-type" content="text/html; charset=UTF-8"> +<title>~/code/numeric/exam/mainC.cpp.html</title> +<meta name="Generator" content="Vim/7.3"> +<meta name="plugin-version" content="vim7.3_v10"> +<meta name="syntax" content="cpp"> +<meta name="settings" content="number_lines,use_css,pre_wrap,expand_tabs"> +<style type="text/css"> +<!-- +pre { white-space: pre-wrap; font-family: monospace; color: #ffffff; background-color: #000000; } +body { font-family: monospace; color: #ffffff; background-color: #000000; } +.lnr { color: #ffff00; } +.Special { color: #ff40ff; } +.Statement { color: #ffff00; } +.Comment { color: #00ffff; } +.Type { color: #00ff00; } +.Constant { color: #ffff00; } +.PreProc { color: #ff40ff; } +--> +</style> +</head> +<body> +<pre> +<span class="lnr"> 1 </span><span class="PreProc">#include </span><span class="Constant"><iostream></span> +<span class="lnr"> 2 </span><span class="PreProc">#include </span><span class="Constant"><vector></span> +<span class="lnr"> 3 </span><span class="PreProc">#include </span><span class="Constant"><complex></span> +<span class="lnr"> 4 </span><span class="PreProc">#include </span><span class="Constant"><cmath></span> +<span class="lnr"> 5 </span><span class="PreProc">#include </span><span class="Constant">"typedefs.h"</span> +<span class="lnr"> 6 </span><span class="PreProc">#include </span><span class="Constant">"check.h"</span> +<span class="lnr"> 7 </span><span class="PreProc">#include </span><span class="Constant">"ode.h"</span> +<span class="lnr"> 8 </span><span class="PreProc">#include </span><span class="Constant">"functions.h"</span> +<span class="lnr"> 9 </span> +<span class="lnr">10 </span><span class="Type">int</span> main() +<span class="lnr">11 </span>{ +<span class="lnr">12 </span> <span class="Comment">// Namespace declarations</span> +<span class="lnr">13 </span> <span class="Statement">using</span> std::cout; +<span class="lnr">14 </span> <span class="Statement">using</span> std::vector; +<span class="lnr">15 </span> <span class="Statement">using</span> std::complex; +<span class="lnr">16 </span> +<span class="lnr">17 </span> <span class="Comment">// Calculate machine precision</span> +<span class="lnr">18 </span> Floattype eps_machine = <span class="Constant">1.0f</span>; +<span class="lnr">19 </span> <span class="Statement">while</span> (<span class="Constant">1.0f</span> + eps_machine != <span class="Constant">1.0f</span>) +<span class="lnr">20 </span> eps_machine /= <span class="Constant">2.0f</span>; +<span class="lnr">21 </span> +<span class="lnr">22 </span> cout << <span class="Constant">"</span><span class="Special">\n\033</span><span class="Constant">[1;33m--- Part C: Solving along a complex path ---</span><span class="Special">\033</span><span class="Constant">[0m</span><span class="Special">\n</span><span class="Constant">"</span>; +<span class="lnr">23 </span> <span class="Type">complex</span><Floattype> a(<span class="Constant">0.0f</span>, <span class="Constant">0.0f</span>); <span class="Comment">// Lower limit</span> +<span class="lnr">24 </span> <span class="Type">complex</span><Floattype> b(<span class="Constant">2.0f</span>*<span class="Constant">M_PI</span>, <span class="Constant">2.0f</span>*<span class="Constant">M_PI</span>); <span class="Comment">// Upper limit</span> +<span class="lnr">25 </span> cout << <span class="Constant">"Integration path: b-a = "</span> << b-a << <span class="Special">'\n'</span>; +<span class="lnr">26 </span> Inttype n_eqs = <span class="Constant">2</span>; <span class="Comment">// Number of equations in ODE system</span> +<span class="lnr">27 </span> vector<<span class="Type">complex</span><Floattype> > y_start(n_eqs); +<span class="lnr">28 </span> <span class="Type">complex</span><Floattype> y0(<span class="Constant">0.0f</span>, <span class="Constant">0.0f</span>); +<span class="lnr">29 </span> <span class="Type">complex</span><Floattype> y1(<span class="Constant">1.0f</span>, <span class="Constant">1.0f</span>); +<span class="lnr">30 </span> y_start[<span class="Constant">0</span>] = y0; +<span class="lnr">31 </span> y_start[<span class="Constant">1</span>] = y1; +<span class="lnr">32 </span> Floattype h_start = <span class="Constant">0.01f</span>; +<span class="lnr">33 </span> ODE ode(func1, <span class="Comment">// ODE system</span> +<span class="lnr">34 </span> y_start, <span class="Comment">// Initial values</span> +<span class="lnr">35 </span> a, <span class="Comment">// Lower limit</span> +<span class="lnr">36 </span> b, <span class="Comment">// Upper limit</span> +<span class="lnr">37 </span> h_start, <span class="Comment">// Start value of step size</span> +<span class="lnr">38 </span> <span class="Constant">10000</span>, <span class="Comment">// Max. number of steps</span> +<span class="lnr">39 </span> eps_machine*<span class="Constant">1e12f</span>, <span class="Comment">// Absolute precision</span> +<span class="lnr">40 </span> eps_machine*<span class="Constant">1e12f</span>); <span class="Comment">// Relative precision</span> +<span class="lnr">41 </span> ode.write(<span class="Constant">"funcC.dat"</span>); <span class="Comment">// Write solutions to data file</span> +<span class="lnr">42 </span> +<span class="lnr">43 </span> <span class="Comment">// Report to stdout</span> +<span class="lnr">44 </span> cout << <span class="Constant">"ODE system solved in "</span> +<span class="lnr">45 </span> << ode.steps() << <span class="Constant">" steps.</span><span class="Special">\n\n</span><span class="Constant">"</span>; +<span class="lnr">46 </span> +<span class="lnr">47 </span> <span class="Comment">// Return successfully</span> +<span class="lnr">48 </span> <span class="Statement">return</span> <span class="Constant">0</span>; +<span class="lnr">49 </span>} +<span class="lnr">50 </span> +</pre> +</body> +</html> diff --git a/exam/html/mainD.cpp.html b/exam/html/mainD.cpp.html @@ -0,0 +1,96 @@ +<!DOCTYPE HTML PUBLIC "-//W3C//DTD HTML 4.01//EN" "http://www.w3.org/TR/html4/strict.dtd"> +<html> +<head> +<meta http-equiv="content-type" content="text/html; charset=UTF-8"> +<title>~/code/numeric/exam/mainD.cpp.html</title> +<meta name="Generator" content="Vim/7.3"> +<meta name="plugin-version" content="vim7.3_v10"> +<meta name="syntax" content="cpp"> +<meta name="settings" content="number_lines,use_css,pre_wrap,expand_tabs"> +<style type="text/css"> +<!-- +pre { white-space: pre-wrap; font-family: monospace; color: #ffffff; background-color: #000000; } +body { font-family: monospace; color: #ffffff; background-color: #000000; } +.lnr { color: #ffff00; } +.Special { color: #ff40ff; } +.Statement { color: #ffff00; } +.Comment { color: #00ffff; } +.Type { color: #00ff00; } +.Constant { color: #ffff00; } +.PreProc { color: #ff40ff; } +--> +</style> +</head> +<body> +<pre> +<span class="lnr"> 1 </span><span class="PreProc">#include </span><span class="Constant"><iostream></span> +<span class="lnr"> 2 </span><span class="PreProc">#include </span><span class="Constant"><vector></span> +<span class="lnr"> 3 </span><span class="PreProc">#include </span><span class="Constant"><complex></span> +<span class="lnr"> 4 </span><span class="PreProc">#include </span><span class="Constant"><cmath></span> +<span class="lnr"> 5 </span><span class="PreProc">#include </span><span class="Constant"><fstream></span> +<span class="lnr"> 6 </span><span class="PreProc">#include </span><span class="Constant">"typedefs.h"</span> +<span class="lnr"> 7 </span><span class="PreProc">#include </span><span class="Constant">"check.h"</span> +<span class="lnr"> 8 </span><span class="PreProc">#include </span><span class="Constant">"ode.h"</span> +<span class="lnr"> 9 </span><span class="PreProc">#include </span><span class="Constant">"functions.h"</span> +<span class="lnr">10 </span> +<span class="lnr">11 </span> +<span class="lnr">12 </span><span class="Type">int</span> main() +<span class="lnr">13 </span>{ +<span class="lnr">14 </span> <span class="Comment">// Namespace declarations</span> +<span class="lnr">15 </span> <span class="Statement">using</span> std::cout; +<span class="lnr">16 </span> <span class="Statement">using</span> std::vector; +<span class="lnr">17 </span> <span class="Statement">using</span> std::complex; +<span class="lnr">18 </span> +<span class="lnr">19 </span> <span class="Comment">// Calculate machine precision</span> +<span class="lnr">20 </span> Floattype eps_machine = <span class="Constant">1.0f</span>; +<span class="lnr">21 </span> <span class="Statement">while</span> (<span class="Constant">1.0f</span> + eps_machine != <span class="Constant">1.0f</span>) +<span class="lnr">22 </span> eps_machine /= <span class="Constant">2.0f</span>; +<span class="lnr">23 </span> +<span class="lnr">24 </span> cout << <span class="Constant">"</span><span class="Special">\n\033</span><span class="Constant">[1;33m--- Part D: Precision analysis ---</span><span class="Special">\033</span><span class="Constant">[0m</span><span class="Special">\n</span><span class="Constant">"</span>; +<span class="lnr">25 </span> <span class="Type">complex</span><Floattype> a(<span class="Constant">0.0f</span>, <span class="Constant">0.0f</span>); <span class="Comment">// Lower limit</span> +<span class="lnr">26 </span> <span class="Type">complex</span><Floattype> b(<span class="Constant">2.0f</span>*<span class="Constant">M_PI</span>, <span class="Constant">2.0f</span>*<span class="Constant">M_PI</span>); <span class="Comment">// Upper limit</span> +<span class="lnr">27 </span> cout << <span class="Constant">"Integration path: b-a = "</span> << b-a << <span class="Special">'\n'</span>; +<span class="lnr">28 </span> Inttype n_eqs = <span class="Constant">2</span>; <span class="Comment">// Number of equations in ODE system</span> +<span class="lnr">29 </span> vector<<span class="Type">complex</span><Floattype> > y_start(n_eqs); +<span class="lnr">30 </span> <span class="Type">complex</span><Floattype> y0(<span class="Constant">0.0f</span>, <span class="Constant">0.0f</span>); +<span class="lnr">31 </span> <span class="Type">complex</span><Floattype> y1(<span class="Constant">1.0f</span>, <span class="Constant">1.0f</span>); +<span class="lnr">32 </span> y_start[<span class="Constant">0</span>] = y0; +<span class="lnr">33 </span> y_start[<span class="Constant">1</span>] = y1; +<span class="lnr">34 </span> Floattype h_start = <span class="Constant">0.01f</span>; +<span class="lnr">35 </span> +<span class="lnr">36 </span> vector<Floattype> precs; <span class="Comment">// Vector containing precision values</span> +<span class="lnr">37 </span> vector<Inttype> steps; <span class="Comment">// Vector containing number of steps required</span> +<span class="lnr">38 </span> +<span class="lnr">39 </span> <span class="Statement">for</span> (Floattype prec=eps_machine*<span class="Constant">10.0f</span>; prec<<span class="Constant">0.1f</span>; prec *= <span class="Constant">10.0f</span>) { +<span class="lnr">40 </span> ODE ode(func1, <span class="Comment">// ODE system</span> +<span class="lnr">41 </span> y_start, <span class="Comment">// Initial values</span> +<span class="lnr">42 </span> a, <span class="Comment">// Lower limit</span> +<span class="lnr">43 </span> b, <span class="Comment">// Upper limit</span> +<span class="lnr">44 </span> h_start, <span class="Comment">// Start value of step size</span> +<span class="lnr">45 </span> <span class="Constant">100000</span>, <span class="Comment">// Max. number of steps</span> +<span class="lnr">46 </span> prec, <span class="Comment">// Absolute precision</span> +<span class="lnr">47 </span> prec); <span class="Comment">// Relative precision</span> +<span class="lnr">48 </span> precs.push_back(prec); <span class="Comment">// Save precision</span> +<span class="lnr">49 </span> steps.push_back(ode.steps()); <span class="Comment">// Save number of steps taken</span> +<span class="lnr">50 </span> } +<span class="lnr">51 </span> +<span class="lnr">52 </span> <span class="Comment">// Save results to text file</span> +<span class="lnr">53 </span> std::ofstream ost; <span class="Comment">// Out stream object</span> +<span class="lnr">54 </span> ost.open(<span class="Constant">"funcD.dat"</span>); <span class="Comment">// Open outfile for write</span> +<span class="lnr">55 </span> <span class="Statement">if</span> (!ost) { +<span class="lnr">56 </span> std::cerr << <span class="Constant">"Error, can't open output file!</span><span class="Special">\n</span><span class="Constant">"</span>; +<span class="lnr">57 </span> <span class="Statement">return</span> <span class="Constant">1</span>; +<span class="lnr">58 </span> } +<span class="lnr">59 </span> <span class="Comment">// Write output values</span> +<span class="lnr">60 </span> <span class="Statement">for</span> (Inttype i=<span class="Constant">0</span>; i<precs.size(); ++i) +<span class="lnr">61 </span> ost << precs[i] << <span class="Special">'\t'</span> << steps[i] << <span class="Special">'\n'</span>; +<span class="lnr">62 </span> +<span class="lnr">63 </span> <span class="Comment">// Close file</span> +<span class="lnr">64 </span> ost.close(); +<span class="lnr">65 </span> +<span class="lnr">66 </span> <span class="Comment">// Return successfully</span> +<span class="lnr">67 </span> <span class="Statement">return</span> <span class="Constant">0</span>; +<span class="lnr">68 </span>} +</pre> +</body> +</html> diff --git a/exam/html/ode.cpp.html b/exam/html/ode.cpp.html @@ -0,0 +1,191 @@ +<!DOCTYPE HTML PUBLIC "-//W3C//DTD HTML 4.01//EN" "http://www.w3.org/TR/html4/strict.dtd"> +<html> +<head> +<meta http-equiv="content-type" content="text/html; charset=UTF-8"> +<title>~/code/numeric/exam/ode.cpp.html</title> +<meta name="Generator" content="Vim/7.3"> +<meta name="plugin-version" content="vim7.3_v10"> +<meta name="syntax" content="cpp"> +<meta name="settings" content="number_lines,use_css,pre_wrap,expand_tabs"> +<style type="text/css"> +<!-- +pre { white-space: pre-wrap; font-family: monospace; color: #ffffff; background-color: #000000; } +body { font-family: monospace; color: #ffffff; background-color: #000000; } +.lnr { color: #ffff00; } +.Special { color: #ff40ff; } +.Statement { color: #ffff00; } +.Type { color: #00ff00; } +.Comment { color: #00ffff; } +.Constant { color: #ffff00; } +.PreProc { color: #ff40ff; } +--> +</style> +</head> +<body> +<pre> +<span class="lnr"> 1 </span><span class="PreProc">#include </span><span class="Constant"><iostream></span> +<span class="lnr"> 2 </span><span class="PreProc">#include </span><span class="Constant"><vector></span> +<span class="lnr"> 3 </span><span class="PreProc">#include </span><span class="Constant"><cmath></span> <span class="Comment">// for sqrt and pow</span> +<span class="lnr"> 4 </span><span class="PreProc">#include </span><span class="Constant"><fstream></span> +<span class="lnr"> 5 </span><span class="PreProc">#include </span><span class="Constant">"typedefs.h"</span> +<span class="lnr"> 6 </span><span class="PreProc">#include </span><span class="Constant">"ode.h"</span> +<span class="lnr"> 7 </span><span class="PreProc">#include </span><span class="Constant">"vector_arithmetic.h"</span> +<span class="lnr"> 8 </span> +<span class="lnr"> 9 </span><span class="Comment">// Constructor</span> +<span class="lnr"> 10 </span>ODE::ODE(std::vector<std::complex<Floattype> > +<span class="lnr"> 11 </span> (*f_in)(<span class="Type">const</span> std::complex<Floattype> x, +<span class="lnr"> 12 </span> <span class="Type">const</span> std::vector<std::complex<Floattype> > &y), +<span class="lnr"> 13 </span> <span class="Type">const</span> std::vector<std::complex<Floattype> > y_start, +<span class="lnr"> 14 </span> <span class="Type">const</span> std::complex<Floattype> a_in, +<span class="lnr"> 15 </span> <span class="Type">const</span> std::complex<Floattype> b_in, +<span class="lnr"> 16 </span> <span class="Type">const</span> Floattype h_start, +<span class="lnr"> 17 </span> <span class="Type">const</span> Inttype max_steps, +<span class="lnr"> 18 </span> <span class="Type">const</span> Floattype delta_in, +<span class="lnr"> 19 </span> <span class="Type">const</span> Floattype epsilon_in, +<span class="lnr"> 20 </span> <span class="Type">const</span> Floattype power_in, +<span class="lnr"> 21 </span> <span class="Type">const</span> Floattype safety_in) +<span class="lnr"> 22 </span> : f(f_in), +<span class="lnr"> 23 </span> a(a_in), b(b_in), +<span class="lnr"> 24 </span> h((b_in-a_in)*h_start), +<span class="lnr"> 25 </span> n_max(max_steps), +<span class="lnr"> 26 </span> delta(delta_in), epsilon(epsilon_in), +<span class="lnr"> 27 </span> power(power_in), safety(safety_in) +<span class="lnr"> 28 </span>{ +<span class="lnr"> 29 </span> x_list.push_back(a); +<span class="lnr"> 30 </span> y_list.push_back(y_start); +<span class="lnr"> 31 </span> +<span class="lnr"> 32 </span> <span class="Comment">// Perform integration</span> +<span class="lnr"> 33 </span> rkdriver(); +<span class="lnr"> 34 </span>} +<span class="lnr"> 35 </span> +<span class="lnr"> 36 </span> +<span class="lnr"> 37 </span><span class="Comment">// Estimate tolerance</span> +<span class="lnr"> 38 </span>Floattype ODE::tau(<span class="Type">const</span> std::vector<std::complex<Floattype> > &y, +<span class="lnr"> 39 </span> <span class="Type">const</span> std::complex<Floattype> h_i) +<span class="lnr"> 40 </span>{ +<span class="lnr"> 41 </span> <span class="Statement">return</span> abs((epsilon * cnorm(y) + delta) * sqrt(h_i/(b-a))); +<span class="lnr"> 42 </span>} +<span class="lnr"> 43 </span> +<span class="lnr"> 44 </span><span class="Comment">// Runge-Kutta mid-point stepper</span> +<span class="lnr"> 45 </span><span class="Type">void</span> ODE::rkstep12(<span class="Type">const</span> std::complex<Floattype> x0, +<span class="lnr"> 46 </span> <span class="Type">const</span> std::vector<std::complex<Floattype> > &y0, +<span class="lnr"> 47 </span> std::vector<std::complex<Floattype> > &y1, +<span class="lnr"> 48 </span> std::vector<std::complex<Floattype> > &dy) +<span class="lnr"> 49 </span>{ +<span class="lnr"> 50 </span> <span class="Comment">// Denominator 2 used in arithmetic operations</span> +<span class="lnr"> 51 </span> <span class="Type">const</span> std::complex<Floattype> den2 (<span class="Constant">2.0f</span>,<span class="Constant">2.0f</span>); +<span class="lnr"> 52 </span> +<span class="lnr"> 53 </span> <span class="Comment">// Evaluate function at two points</span> +<span class="lnr"> 54 </span> (<span class="Type">void</span>)f(x0,y0); +<span class="lnr"> 55 </span> <span class="Type">const</span> std::vector<std::complex<Floattype> > k0 = f(x0,y0); +<span class="lnr"> 56 </span> <span class="Type">const</span> std::vector<std::complex<Floattype> > k12 = f(x0 + h/den2, y0 + k0*h/den2); +<span class="lnr"> 57 </span> +<span class="lnr"> 58 </span> <span class="Comment">// Write results to output vectors</span> +<span class="lnr"> 59 </span> y1 = y0 + k12*h; +<span class="lnr"> 60 </span> dy = (k0 - k12) * h/den2; +<span class="lnr"> 61 </span>} +<span class="lnr"> 62 </span> +<span class="lnr"> 63 </span> +<span class="lnr"> 64 </span><span class="Comment">// ODE driver with adaptive step size control</span> +<span class="lnr"> 65 </span><span class="Type">void</span> ODE::rkdriver() +<span class="lnr"> 66 </span>{ +<span class="lnr"> 67 </span> std::vector<std::complex<Floattype> > dy(n_max); +<span class="lnr"> 68 </span> std::vector<std::complex<Floattype> > y1(n_max); +<span class="lnr"> 69 </span> +<span class="lnr"> 70 </span> std::complex<Floattype> x; +<span class="lnr"> 71 </span> Floattype err, tol; +<span class="lnr"> 72 </span> std::vector<std::complex<Floattype> > y; +<span class="lnr"> 73 </span> +<span class="lnr"> 74 </span> <span class="Statement">while</span> (x_list.back().real() < b.real() +<span class="lnr"> 75 </span> || x_list.back().imag() < b.imag()) +<span class="lnr"> 76 </span> { +<span class="lnr"> 77 </span> <span class="Comment">// Get values for this step</span> +<span class="lnr"> 78 </span> x = x_list.back(); +<span class="lnr"> 79 </span> y = y_list.back(); +<span class="lnr"> 80 </span> +<span class="lnr"> 81 </span> <span class="Comment">// Make sure we don't step past the upper boundary</span> +<span class="lnr"> 82 </span> <span class="Statement">if</span> ((x + h).real() > b.real() +<span class="lnr"> 83 </span> || (x + h).imag() > b.imag()) +<span class="lnr"> 84 </span> h = b - x; +<span class="lnr"> 85 </span> +<span class="lnr"> 86 </span> <span class="Comment">// Run Runge-Kutta mid-point stepper</span> +<span class="lnr"> 87 </span> rkstep12(x, y, y1, dy); +<span class="lnr"> 88 </span> +<span class="lnr"> 89 </span> <span class="Comment">// Determine whether the step should be accepted</span> +<span class="lnr"> 90 </span> err = cnorm(dy); <span class="Comment">// Error estimate</span> +<span class="lnr"> 91 </span> tol = tau(y, h); <span class="Comment">// Tolerance</span> +<span class="lnr"> 92 </span> <span class="Statement">if</span> (err < tol) { <span class="Comment">// Step accepted</span> +<span class="lnr"> 93 </span> x_list.push_back(x+h); +<span class="lnr"> 94 </span> y_list.push_back(y1); +<span class="lnr"> 95 </span> } +<span class="lnr"> 96 </span> +<span class="lnr"> 97 </span> <span class="Comment">// Check that we havn't hit the max. number of steps</span> +<span class="lnr"> 98 </span> <span class="Statement">if</span> (x_list.size() == n_max) { +<span class="lnr"> 99 </span> std::cerr << <span class="Constant">"Error, the max. number of steps "</span> +<span class="lnr">100 </span> << <span class="Constant">"was insufficient</span><span class="Special">\n</span><span class="Constant">"</span> +<span class="lnr">101 </span> << <span class="Constant">"Try either increasing the max. number "</span> +<span class="lnr">102 </span> << <span class="Constant">"of steps, or decreasing the precision "</span> +<span class="lnr">103 </span> << <span class="Constant">"requirements</span><span class="Special">\n</span><span class="Constant">"</span>; +<span class="lnr">104 </span> <span class="Statement">return</span>; +<span class="lnr">105 </span> } +<span class="lnr">106 </span> +<span class="lnr">107 </span> <span class="Comment">// Determine new step size</span> +<span class="lnr">108 </span> std::complex<Floattype> multiplicator (<span class="Constant">2.0f</span>, <span class="Constant">2.0f</span>); +<span class="lnr">109 </span> <span class="Statement">if</span> (err > <span class="Constant">0.0f</span>) +<span class="lnr">110 </span> h = h*pow(tol/err, power) * safety; +<span class="lnr">111 </span> <span class="Statement">else</span> +<span class="lnr">112 </span> h = multiplicator*h; +<span class="lnr">113 </span> } +<span class="lnr">114 </span>} +<span class="lnr">115 </span> +<span class="lnr">116 </span> +<span class="lnr">117 </span><span class="Comment">// Return the number of steps taken</span> +<span class="lnr">118 </span>Inttype ODE::steps() +<span class="lnr">119 </span>{ +<span class="lnr">120 </span> <span class="Statement">return</span> x_list.size(); +<span class="lnr">121 </span>} +<span class="lnr">122 </span> +<span class="lnr">123 </span><span class="Type">void</span> ODE::print() +<span class="lnr">124 </span>{ +<span class="lnr">125 </span> <span class="Statement">for</span> (Inttype i=<span class="Constant">0</span>; i<x_list.size(); ++i) { +<span class="lnr">126 </span> std::cout << x_list[i] << <span class="Special">'\t'</span>; +<span class="lnr">127 </span> <span class="Statement">for</span> (Inttype j=<span class="Constant">0</span>; j<y_list[<span class="Constant">0</span>].size(); ++j) +<span class="lnr">128 </span> std::cout << y_list[i][j] << <span class="Special">'\t'</span>; +<span class="lnr">129 </span> std::cout << <span class="Special">'\n'</span>; +<span class="lnr">130 </span> } +<span class="lnr">131 </span>} +<span class="lnr">132 </span> +<span class="lnr">133 </span><span class="Comment">// Write the x- and y-values to file</span> +<span class="lnr">134 </span><span class="Type">void</span> ODE::write(<span class="Type">const</span> <span class="Type">char</span>* filename) +<span class="lnr">135 </span>{ +<span class="lnr">136 </span> std::ofstream outstream; +<span class="lnr">137 </span> +<span class="lnr">138 </span> <span class="Comment">// Open outfile for write</span> +<span class="lnr">139 </span> outstream.open(filename); +<span class="lnr">140 </span> <span class="Statement">if</span> (!outstream) { +<span class="lnr">141 </span> std::cerr << <span class="Constant">"Error, can't open output file '"</span> +<span class="lnr">142 </span> << filename << <span class="Constant">"'.</span><span class="Special">\n</span><span class="Constant">"</span>; +<span class="lnr">143 </span> <span class="Statement">return</span>; +<span class="lnr">144 </span> } +<span class="lnr">145 </span> +<span class="lnr">146 </span> <span class="Comment">// Write output values</span> +<span class="lnr">147 </span> <span class="Statement">for</span> (Inttype i=<span class="Constant">0</span>; i<x_list.size(); ++i) { +<span class="lnr">148 </span> outstream << x_list[i].real() << <span class="Special">'\t'</span>; +<span class="lnr">149 </span> outstream << x_list[i].imag() << <span class="Special">'\t'</span>; +<span class="lnr">150 </span> <span class="Statement">for</span> (Inttype j=<span class="Constant">0</span>; j<y_list[<span class="Constant">0</span>].size(); ++j) { +<span class="lnr">151 </span> outstream << y_list[i][j].real() << <span class="Special">'\t'</span>; +<span class="lnr">152 </span> outstream << y_list[i][j].imag() << <span class="Special">'\t'</span>; +<span class="lnr">153 </span> } +<span class="lnr">154 </span> outstream << <span class="Special">'\n'</span>; +<span class="lnr">155 </span> } +<span class="lnr">156 </span> +<span class="lnr">157 </span> <span class="Comment">// Close file</span> +<span class="lnr">158 </span> outstream.close(); +<span class="lnr">159 </span> +<span class="lnr">160 </span> <span class="Statement">if</span> (verbose == <span class="Constant">true</span>) +<span class="lnr">161 </span> std::cout << <span class="Constant">"Output written in '"</span> << filename << <span class="Constant">"'.</span><span class="Special">\n</span><span class="Constant">"</span>; +<span class="lnr">162 </span>} +<span class="lnr">163 </span> +</pre> +</body> +</html> diff --git a/exam/html/ode.h.html b/exam/html/ode.h.html @@ -0,0 +1,114 @@ +<!DOCTYPE HTML PUBLIC "-//W3C//DTD HTML 4.01//EN" "http://www.w3.org/TR/html4/strict.dtd"> +<html> +<head> +<meta http-equiv="content-type" content="text/html; charset=UTF-8"> +<title>~/code/numeric/exam/ode.h.html</title> +<meta name="Generator" content="Vim/7.3"> +<meta name="plugin-version" content="vim7.3_v10"> +<meta name="syntax" content="cpp"> +<meta name="settings" content="number_lines,use_css,pre_wrap,expand_tabs"> +<style type="text/css"> +<!-- +pre { white-space: pre-wrap; font-family: monospace; color: #ffffff; background-color: #000000; } +body { font-family: monospace; color: #ffffff; background-color: #000000; } +.lnr { color: #ffff00; } +.Statement { color: #ffff00; } +.Type { color: #00ff00; } +.Constant { color: #ffff00; } +.PreProc { color: #ff40ff; } +.Comment { color: #00ffff; } +--> +</style> +</head> +<body> +<pre> +<span class="lnr"> 1 </span><span class="Comment">// Make sure header is only included once</span> +<span class="lnr"> 2 </span><span class="PreProc">#ifndef ODE_H_</span> +<span class="lnr"> 3 </span><span class="PreProc">#define ODE_H_</span> +<span class="lnr"> 4 </span> +<span class="lnr"> 5 </span><span class="PreProc">#include </span><span class="Constant"><vector></span> +<span class="lnr"> 6 </span><span class="PreProc">#include </span><span class="Constant"><complex></span> +<span class="lnr"> 7 </span><span class="PreProc">#include </span><span class="Constant">"typedefs.h"</span> +<span class="lnr"> 8 </span> +<span class="lnr"> 9 </span><span class="Comment">// ODE class</span> +<span class="lnr">10 </span><span class="Type">class</span> ODE { +<span class="lnr">11 </span> +<span class="lnr">12 </span> <span class="Comment">// Values and functions only accessible from the class internally</span> +<span class="lnr">13 </span> <span class="Statement">private</span>: +<span class="lnr">14 </span> +<span class="lnr">15 </span> <span class="Comment">// System of ordinary differential equations to solve</span> +<span class="lnr">16 </span> std::vector<std::complex<Floattype> > +<span class="lnr">17 </span> (*f)(<span class="Type">const</span> std::complex<Floattype> x, +<span class="lnr">18 </span> <span class="Type">const</span> std::vector<std::complex<Floattype> > &y); +<span class="lnr">19 </span> +<span class="lnr">20 </span> <span class="Comment">// Points to be evaluated</span> +<span class="lnr">21 </span> std::vector<std::complex<Floattype> > x_list; +<span class="lnr">22 </span> +<span class="lnr">23 </span> <span class="Comment">// Limits of range to evaluate</span> +<span class="lnr">24 </span> <span class="Type">const</span> std::complex<Floattype> a; <span class="Comment">// Lower</span> +<span class="lnr">25 </span> <span class="Type">const</span> std::complex<Floattype> b; <span class="Comment">// Upper</span> +<span class="lnr">26 </span> +<span class="lnr">27 </span> <span class="Comment">// Step size</span> +<span class="lnr">28 </span> std::complex<Floattype> h; +<span class="lnr">29 </span> +<span class="lnr">30 </span> <span class="Comment">// Results stored in 2D: vector of vectors</span> +<span class="lnr">31 </span> std::vector<std::vector<std::complex<Floattype> > > y_list; +<span class="lnr">32 </span> +<span class="lnr">33 </span> <span class="Comment">// Maximum number of steps to evaluate, defined by y size</span> +<span class="lnr">34 </span> <span class="Type">const</span> Inttype n_max; +<span class="lnr">35 </span> +<span class="lnr">36 </span> <span class="Comment">// Accuracy requirement values</span> +<span class="lnr">37 </span> <span class="Type">const</span> Floattype delta; <span class="Comment">// Absolute</span> +<span class="lnr">38 </span> <span class="Type">const</span> Floattype epsilon; <span class="Comment">// Relative</span> +<span class="lnr">39 </span> +<span class="lnr">40 </span> <span class="Comment">// Tolerance estimator</span> +<span class="lnr">41 </span> Floattype tau(<span class="Type">const</span> std::vector<std::complex<Floattype> > &y, +<span class="lnr">42 </span> <span class="Type">const</span> std::complex<Floattype> h); +<span class="lnr">43 </span> +<span class="lnr">44 </span> <span class="Comment">// Runge-Kutta mid-point stepper prototype</span> +<span class="lnr">45 </span> <span class="Type">void</span> rkstep12(<span class="Type">const</span> std::complex<Floattype> x0, +<span class="lnr">46 </span> <span class="Type">const</span> std::vector<std::complex<Floattype> > &y0, +<span class="lnr">47 </span> std::vector<std::complex<Floattype> > &y1, +<span class="lnr">48 </span> std::vector<std::complex<Floattype> > &dy); +<span class="lnr">49 </span> +<span class="lnr">50 </span> <span class="Comment">// Runge-Kutta driver function parameters</span> +<span class="lnr">51 </span> <span class="Type">const</span> Floattype power; +<span class="lnr">52 </span> <span class="Type">const</span> Floattype safety; +<span class="lnr">53 </span> +<span class="lnr">54 </span> <span class="Comment">// Runge-Kutta driver prototype</span> +<span class="lnr">55 </span> <span class="Type">void</span> rkdriver(); +<span class="lnr">56 </span> +<span class="lnr">57 </span> +<span class="lnr">58 </span> <span class="Comment">// Values and functions accessible from the outside</span> +<span class="lnr">59 </span> <span class="Statement">public</span>: +<span class="lnr">60 </span> +<span class="lnr">61 </span> <span class="Comment">// Constructor, some parameters with default values</span> +<span class="lnr">62 </span> ODE(std::vector<std::complex<Floattype> > +<span class="lnr">63 </span> (*f_in)(<span class="Type">const</span> std::complex<Floattype> x, +<span class="lnr">64 </span> <span class="Type">const</span> std::vector<std::complex<Floattype> > &y), +<span class="lnr">65 </span> <span class="Type">const</span> std::vector<std::complex<Floattype> > y_start, +<span class="lnr">66 </span> <span class="Type">const</span> std::complex<Floattype> a_in, +<span class="lnr">67 </span> <span class="Type">const</span> std::complex<Floattype> b_in, +<span class="lnr">68 </span> <span class="Type">const</span> Floattype h_start = <span class="Constant">0.01f</span>, +<span class="lnr">69 </span> <span class="Type">const</span> Inttype max_steps = <span class="Constant">1e4</span>, +<span class="lnr">70 </span> <span class="Type">const</span> Floattype delta_in = <span class="Constant">1e-3f</span>, +<span class="lnr">71 </span> <span class="Type">const</span> Floattype epsilon_in = <span class="Constant">1e-3f</span>, +<span class="lnr">72 </span> <span class="Type">const</span> Floattype power_in = <span class="Constant">0.25f</span>, +<span class="lnr">73 </span> <span class="Type">const</span> Floattype safety_in = <span class="Constant">0.95f</span> +<span class="lnr">74 </span> ); +<span class="lnr">75 </span> +<span class="lnr">76 </span> <span class="Comment">// Return the number of steps taken</span> +<span class="lnr">77 </span> Inttype steps(); +<span class="lnr">78 </span> +<span class="lnr">79 </span> <span class="Comment">// Print the x- and y-values to stdout</span> +<span class="lnr">80 </span> <span class="Type">void</span> print(); +<span class="lnr">81 </span> +<span class="lnr">82 </span> <span class="Comment">// Write the x- and y-values to file</span> +<span class="lnr">83 </span> <span class="Type">void</span> write(<span class="Type">const</span> <span class="Type">char</span>* filename); +<span class="lnr">84 </span> +<span class="lnr">85 </span>}; +<span class="lnr">86 </span> +<span class="lnr">87 </span><span class="PreProc">#endif</span> +</pre> +</body> +</html> diff --git a/exam/html/plotA.gp.html b/exam/html/plotA.gp.html @@ -0,0 +1,32 @@ +<!DOCTYPE HTML PUBLIC "-//W3C//DTD HTML 4.01//EN" "http://www.w3.org/TR/html4/strict.dtd"> +<html> +<head> +<meta http-equiv="content-type" content="text/html; charset=UTF-8"> +<title>~/code/numeric/exam/plotA.gp.html</title> +<meta name="Generator" content="Vim/7.3"> +<meta name="plugin-version" content="vim7.3_v10"> +<meta name="syntax" content="gp"> +<meta name="settings" content="number_lines,use_css,pre_wrap,expand_tabs"> +<style type="text/css"> +<!-- +pre { white-space: pre-wrap; font-family: monospace; color: #ffffff; background-color: #000000; } +body { font-family: monospace; color: #ffffff; background-color: #000000; } +.lnr { color: #ffff00; } +.Constant { color: #ffff00; } +.Statement { color: #ffff00; } +--> +</style> +</head> +<body> +<pre> +<span class="lnr">1 </span>set terminal png +<span class="lnr">2 </span>set <span class="Statement">output</span> <span class="Constant">"plotA.png"</span> +<span class="lnr">3 </span>set xlabel <span class="Constant">"x.real"</span> +<span class="lnr">4 </span>set ylabel <span class="Constant">"y.real"</span> +<span class="lnr">5 </span>set title <span class="Constant">"Integration along a real path"</span> +<span class="lnr">6 </span>set grid +<span class="lnr">7 </span>plot <span class="Constant">"funcA.dat"</span> u 1:3 t <span class="Constant">"rk12"</span> w lp +<span class="lnr">8 </span> +</pre> +</body> +</html> diff --git a/exam/html/plotB.gp.html b/exam/html/plotB.gp.html @@ -0,0 +1,32 @@ +<!DOCTYPE HTML PUBLIC "-//W3C//DTD HTML 4.01//EN" "http://www.w3.org/TR/html4/strict.dtd"> +<html> +<head> +<meta http-equiv="content-type" content="text/html; charset=UTF-8"> +<title>~/code/numeric/exam/plotB.gp.html</title> +<meta name="Generator" content="Vim/7.3"> +<meta name="plugin-version" content="vim7.3_v10"> +<meta name="syntax" content="gp"> +<meta name="settings" content="number_lines,use_css,pre_wrap,expand_tabs"> +<style type="text/css"> +<!-- +pre { white-space: pre-wrap; font-family: monospace; color: #ffffff; background-color: #000000; } +body { font-family: monospace; color: #ffffff; background-color: #000000; } +.lnr { color: #ffff00; } +.Constant { color: #ffff00; } +.Statement { color: #ffff00; } +--> +</style> +</head> +<body> +<pre> +<span class="lnr">1 </span>set terminal png +<span class="lnr">2 </span>set <span class="Statement">output</span> <span class="Constant">"plotB.png"</span> +<span class="lnr">3 </span>set xlabel <span class="Constant">"x.imag"</span> +<span class="lnr">4 </span>set ylabel <span class="Constant">"y.imag"</span> +<span class="lnr">5 </span>set title <span class="Constant">"Integration along an imaginary path"</span> +<span class="lnr">6 </span>set grid +<span class="lnr">7 </span>plot <span class="Constant">"funcB.dat"</span> u 2:4 t <span class="Constant">"rk12"</span> w lp +<span class="lnr">8 </span> +</pre> +</body> +</html> diff --git a/exam/html/plotC.gp.html b/exam/html/plotC.gp.html @@ -0,0 +1,33 @@ +<!DOCTYPE HTML PUBLIC "-//W3C//DTD HTML 4.01//EN" "http://www.w3.org/TR/html4/strict.dtd"> +<html> +<head> +<meta http-equiv="content-type" content="text/html; charset=UTF-8"> +<title>~/code/numeric/exam/plotC.gp.html</title> +<meta name="Generator" content="Vim/7.3"> +<meta name="plugin-version" content="vim7.3_v10"> +<meta name="syntax" content="gp"> +<meta name="settings" content="number_lines,use_css,pre_wrap,expand_tabs"> +<style type="text/css"> +<!-- +pre { white-space: pre-wrap; font-family: monospace; color: #ffffff; background-color: #000000; } +body { font-family: monospace; color: #ffffff; background-color: #000000; } +.lnr { color: #ffff00; } +.Constant { color: #ffff00; } +.Statement { color: #ffff00; } +--> +</style> +</head> +<body> +<pre> +<span class="lnr">1 </span>set terminal png +<span class="lnr">2 </span>set <span class="Statement">output</span> <span class="Constant">"plotC.png"</span> +<span class="lnr">3 </span>set xlabel <span class="Constant">"x.r"</span> +<span class="lnr">4 </span>set ylabel <span class="Constant">"x.i"</span> +<span class="lnr">5 </span>set zlabel <span class="Constant">"y"</span> +<span class="lnr">6 </span>set title <span class="Constant">"Integration along a complex path"</span> +<span class="lnr">7 </span>set grid +<span class="lnr">8 </span>splot <span class="Constant">"funcC.dat"</span> u 1:2:3 t <span class="Constant">"y.real"</span> w lp, <span class="Constant">""</span> u 1:2:4 t <span class="Constant">"y.imag"</span> w lp +<span class="lnr">9 </span> +</pre> +</body> +</html> diff --git a/exam/html/plotD.gp.html b/exam/html/plotD.gp.html @@ -0,0 +1,33 @@ +<!DOCTYPE HTML PUBLIC "-//W3C//DTD HTML 4.01//EN" "http://www.w3.org/TR/html4/strict.dtd"> +<html> +<head> +<meta http-equiv="content-type" content="text/html; charset=UTF-8"> +<title>~/code/numeric/exam/plotD.gp.html</title> +<meta name="Generator" content="Vim/7.3"> +<meta name="plugin-version" content="vim7.3_v10"> +<meta name="syntax" content="gp"> +<meta name="settings" content="number_lines,use_css,pre_wrap,expand_tabs"> +<style type="text/css"> +<!-- +pre { white-space: pre-wrap; font-family: monospace; color: #ffffff; background-color: #000000; } +body { font-family: monospace; color: #ffffff; background-color: #000000; } +.lnr { color: #ffff00; } +.Constant { color: #ffff00; } +.Statement { color: #ffff00; } +--> +</style> +</head> +<body> +<pre> +<span class="lnr">1 </span>set terminal png +<span class="lnr">2 </span>set <span class="Statement">output</span> <span class="Constant">"plotD.png"</span> +<span class="lnr">3 </span>set xlabel <span class="Constant">"Precision"</span> +<span class="lnr">4 </span>set ylabel <span class="Constant">"No. of steps"</span> +<span class="lnr">5 </span>set title <span class="Constant">"Number of steps required for a given absolute and relative precision"</span> +<span class="lnr">6 </span>set grid +<span class="lnr">7 </span>set <span class="Statement">log</span> xy +<span class="lnr">8 </span>plot <span class="Constant">"funcD.dat"</span> u 1:2 t <span class="Constant">"rk12"</span> w lp +<span class="lnr">9 </span> +</pre> +</body> +</html> diff --git a/exam/html/typedefs.h.html b/exam/html/typedefs.h.html @@ -0,0 +1,41 @@ +<!DOCTYPE HTML PUBLIC "-//W3C//DTD HTML 4.01//EN" "http://www.w3.org/TR/html4/strict.dtd"> +<html> +<head> +<meta http-equiv="content-type" content="text/html; charset=UTF-8"> +<title>~/code/numeric/exam/typedefs.h.html</title> +<meta name="Generator" content="Vim/7.3"> +<meta name="plugin-version" content="vim7.3_v10"> +<meta name="syntax" content="cpp"> +<meta name="settings" content="number_lines,use_css,pre_wrap,expand_tabs"> +<style type="text/css"> +<!-- +pre { white-space: pre-wrap; font-family: monospace; color: #ffffff; background-color: #000000; } +body { font-family: monospace; color: #ffffff; background-color: #000000; } +.lnr { color: #ffff00; } +.Constant { color: #ffff00; } +.Type { color: #00ff00; } +.PreProc { color: #ff40ff; } +.Comment { color: #00ffff; } +--> +</style> +</head> +<body> +<pre> +<span class="lnr"> 1 </span><span class="Comment">// Make sure header is only included once</span> +<span class="lnr"> 2 </span><span class="PreProc">#ifndef TYPEDEFS_H_</span> +<span class="lnr"> 3 </span><span class="PreProc">#define TYPEDEFS_H_</span> +<span class="lnr"> 4 </span> +<span class="lnr"> 5 </span><span class="Comment">// Define whether the program should output values of matrices</span> +<span class="lnr"> 6 </span><span class="Comment">//const bool verbose = false;</span> +<span class="lnr"> 7 </span><span class="Type">const</span> <span class="Type">bool</span> verbose = <span class="Constant">true</span>; +<span class="lnr"> 8 </span> +<span class="lnr"> 9 </span><span class="Comment">// Choose length variable type</span> +<span class="lnr">10 </span><span class="Type">typedef</span> <span class="Type">unsigned</span> <span class="Type">int</span> Inttype; +<span class="lnr">11 </span> +<span class="lnr">12 </span><span class="Comment">// Choose floating-point precision</span> +<span class="lnr">13 </span><span class="Type">typedef</span> <span class="Type">double</span> Floattype; +<span class="lnr">14 </span> +<span class="lnr">15 </span><span class="PreProc">#endif</span> +</pre> +</body> +</html> diff --git a/exam/html/vector_arithmetic.h.html b/exam/html/vector_arithmetic.h.html @@ -0,0 +1,116 @@ +<!DOCTYPE HTML PUBLIC "-//W3C//DTD HTML 4.01//EN" "http://www.w3.org/TR/html4/strict.dtd"> +<html> +<head> +<meta http-equiv="content-type" content="text/html; charset=UTF-8"> +<title>~/code/numeric/exam/vector_arithmetic.h.html</title> +<meta name="Generator" content="Vim/7.3"> +<meta name="plugin-version" content="vim7.3_v10"> +<meta name="syntax" content="cpp"> +<meta name="settings" content="number_lines,use_css,pre_wrap,expand_tabs"> +<style type="text/css"> +<!-- +pre { white-space: pre-wrap; font-family: monospace; color: #ffffff; background-color: #000000; } +body { font-family: monospace; color: #ffffff; background-color: #000000; } +.lnr { color: #ffff00; } +.Type { color: #00ff00; } +.Statement { color: #ffff00; } +.Constant { color: #ffff00; } +.PreProc { color: #ff40ff; } +.Comment { color: #00ffff; } +--> +</style> +</head> +<body> +<pre> +<span class="lnr"> 1 </span><span class="Comment">// Make sure header is only included once</span> +<span class="lnr"> 2 </span><span class="PreProc">#ifndef VECTOR_ARITHMETIC_H_</span> +<span class="lnr"> 3 </span><span class="PreProc">#define VECTOR_ARITHMETIC_H_</span> +<span class="lnr"> 4 </span> +<span class="lnr"> 5 </span><span class="PreProc">#include </span><span class="Constant"><vector></span> +<span class="lnr"> 6 </span><span class="PreProc">#include </span><span class="Constant"><cmath></span> <span class="Comment">// for sqrt</span> +<span class="lnr"> 7 </span><span class="PreProc">#include </span><span class="Constant">"typedefs.h"</span> +<span class="lnr"> 8 </span> +<span class="lnr"> 9 </span><span class="Comment">//// Overload vector methods to allow scalar-</span> +<span class="lnr">10 </span><span class="Comment">//// and element-wise arithmetic operations</span> +<span class="lnr">11 </span> +<span class="lnr">12 </span><span class="Comment">// Scalar multiplication (same scalar for real and imaginary parts)</span> +<span class="lnr">13 </span>std::vector<std::complex<Floattype> > +<span class="lnr">14 </span> <span class="Statement">operator</span>*(<span class="Type">const</span> std::vector<std::complex<Floattype> > vec, +<span class="lnr">15 </span> <span class="Type">const</span> Floattype scalar) +<span class="lnr">16 </span>{ +<span class="lnr">17 </span> std::vector<std::complex<Floattype> > result(vec.size()); +<span class="lnr">18 </span> <span class="Statement">for</span> (Inttype i=<span class="Constant">0</span>; i<vec.size(); ++i) { +<span class="lnr">19 </span> result[i].real() = real(vec[i])*scalar; +<span class="lnr">20 </span> result[i].imag() = imag(vec[i])*scalar; +<span class="lnr">21 </span> } +<span class="lnr">22 </span> <span class="Statement">return</span> result; +<span class="lnr">23 </span>} +<span class="lnr">24 </span> +<span class="lnr">25 </span><span class="Comment">// Scalar multiplication</span> +<span class="lnr">26 </span>std::vector<std::complex<Floattype> > +<span class="lnr">27 </span> <span class="Statement">operator</span>*(<span class="Type">const</span> std::vector<std::complex<Floattype> > vec, +<span class="lnr">28 </span> <span class="Type">const</span> std::complex<Floattype> scalar) +<span class="lnr">29 </span>{ +<span class="lnr">30 </span> std::vector<std::complex<Floattype> > result(vec.size()); +<span class="lnr">31 </span> <span class="Statement">for</span> (Inttype i=<span class="Constant">0</span>; i<(Inttype)vec.size(); ++i) +<span class="lnr">32 </span> result[i] = vec[i]*scalar; +<span class="lnr">33 </span> <span class="Statement">return</span> result; +<span class="lnr">34 </span>} +<span class="lnr">35 </span> +<span class="lnr">36 </span><span class="Comment">// Scalar division </span> +<span class="lnr">37 </span>std::vector<std::complex<Floattype> > +<span class="lnr">38 </span> <span class="Statement">operator</span>/(<span class="Type">const</span> std::vector<std::complex<Floattype> > vec, +<span class="lnr">39 </span> <span class="Type">const</span> std::complex<Floattype> scalar) +<span class="lnr">40 </span>{ +<span class="lnr">41 </span> std::vector<std::complex<Floattype> > result(vec.size()); +<span class="lnr">42 </span> <span class="Statement">for</span> (Inttype i=<span class="Constant">0</span>; i<(Inttype)vec.size(); ++i) +<span class="lnr">43 </span> result[i] = vec[i]/scalar; +<span class="lnr">44 </span> <span class="Statement">return</span> result; +<span class="lnr">45 </span>} +<span class="lnr">46 </span> +<span class="lnr">47 </span><span class="Comment">// Element-wise addition</span> +<span class="lnr">48 </span>std::vector<std::complex<Floattype> > +<span class="lnr">49 </span> <span class="Statement">operator</span>+(<span class="Type">const</span> std::vector<std::complex<Floattype> > vec1, +<span class="lnr">50 </span> <span class="Type">const</span> std::vector<std::complex<Floattype> > vec2) +<span class="lnr">51 </span>{ +<span class="lnr">52 </span> std::vector<std::complex<Floattype> > result(vec1.size()); +<span class="lnr">53 </span> <span class="Statement">for</span> (Inttype i=<span class="Constant">0</span>; i<(Inttype)vec1.size(); ++i) +<span class="lnr">54 </span> result[i] = vec1[i] + vec2[i]; +<span class="lnr">55 </span> <span class="Statement">return</span> result; +<span class="lnr">56 </span>} +<span class="lnr">57 </span> +<span class="lnr">58 </span><span class="Comment">// Element-wise subtraction</span> +<span class="lnr">59 </span>std::vector<std::complex<Floattype> > +<span class="lnr">60 </span> <span class="Statement">operator</span>-(<span class="Type">const</span> std::vector<std::complex<Floattype> > vec1, +<span class="lnr">61 </span> <span class="Type">const</span> std::vector<std::complex<Floattype> > vec2) +<span class="lnr">62 </span>{ +<span class="lnr">63 </span> std::vector<std::complex<Floattype> > result(vec1.size()); +<span class="lnr">64 </span> <span class="Statement">for</span> (Inttype i=<span class="Constant">0</span>; i<(Inttype)vec1.size(); ++i) +<span class="lnr">65 </span> result[i] = vec1[i] - vec2[i]; +<span class="lnr">66 </span> <span class="Statement">return</span> result; +<span class="lnr">67 </span>} +<span class="lnr">68 </span> +<span class="lnr">69 </span><span class="Comment">// Element-wise multiplication</span> +<span class="lnr">70 </span>std::vector<std::complex<Floattype> > +<span class="lnr">71 </span> <span class="Statement">operator</span>*(<span class="Type">const</span> std::vector<std::complex<Floattype> > vec1, +<span class="lnr">72 </span> <span class="Type">const</span> std::vector<std::complex<Floattype> > vec2) +<span class="lnr">73 </span>{ +<span class="lnr">74 </span> std::vector<std::complex<Floattype> > result(vec1.size()); +<span class="lnr">75 </span> <span class="Statement">for</span> (Inttype i=<span class="Constant">0</span>; i<(Inttype)vec1.size(); ++i) +<span class="lnr">76 </span> result[i] = vec1[i] * vec2[i]; +<span class="lnr">77 </span> <span class="Statement">return</span> result; +<span class="lnr">78 </span>} +<span class="lnr">79 </span> +<span class="lnr">80 </span><span class="Comment">// Normalize vector</span> +<span class="lnr">81 </span>Floattype cnorm(<span class="Type">const</span> std::vector<std::complex<Floattype> > vec) +<span class="lnr">82 </span>{ +<span class="lnr">83 </span> Floattype res = <span class="Constant">0.0f</span>; +<span class="lnr">84 </span> <span class="Statement">for</span> (Inttype i=<span class="Constant">0</span>; i<(Inttype)vec.size(); ++i) +<span class="lnr">85 </span> res += norm(vec[i])*norm(vec[i]); +<span class="lnr">86 </span> <span class="Statement">return</span> sqrt(res); +<span class="lnr">87 </span>} +<span class="lnr">88 </span> +<span class="lnr">89 </span><span class="PreProc">#endif</span> +</pre> +</body> +</html> diff --git a/exam/mainA.cpp b/exam/mainA.cpp @@ -0,0 +1,59 @@ +#include <iostream> +#include <vector> +#include <complex> +#include <cmath> +#include "typedefs.h" +#include "check.h" +#include "ode.h" +#include "functions.h" + + +int main() +{ + // Namespace declarations + using std::cout; + using std::vector; + using std::complex; + + // Calculate machine precision + Floattype eps_machine = 1.0f; + while (1.0f + eps_machine != 1.0f) + eps_machine /= 2.0f; + + const int id = 20062213; + const char n = 10; + cout << "\nMy student id is \033[1;37m" << id + << "\033[0m, resulting in exam exercise: \033[1;31m" + << id%n << "\033[0m\n"; + cout << "Examination project:\033[1;37m ODE integration " + << "with complex numbers\033[0m\n\n"; + + cout << "\033[1;33m--- Part A: Solving along a real path ---\033[0m\n"; + complex<Floattype> a(0.0f, 0.0f); // Lower limit + complex<Floattype> b(2.0f*M_PI, 0.0f); // Upper limit + cout << "Integration path: b-a = " << b-a << '\n'; + Inttype n_eqs = 2; // Number of equations in ODE system + vector<complex<Floattype> > y_start(n_eqs); + complex<Floattype> y0(0.0f, 0.0f); + complex<Floattype> y1(1.0f, 1.0f); + y_start[0] = y0; + y_start[1] = y1; + Floattype h_start = 0.01f; + ODE realode(func1, // ODE system + y_start, // Initial values + a, // Lower limit + b, // Upper limit + h_start, // Start value of step size + 10000, // Max. number of steps + eps_machine*1e12f, // Absolute precision + eps_machine*1e12f); // Relative precision + realode.write("funcA.dat"); // Write solutions to data file + + // Report to stdout + cout << "ODE system solved in " + << realode.steps() << " steps.\n\n"; + + // Return successfully + return 0; +} + diff --git a/exam/mainB.cpp b/exam/mainB.cpp @@ -0,0 +1,50 @@ +#include <iostream> +#include <vector> +#include <complex> +#include <cmath> +#include "typedefs.h" +#include "check.h" +#include "ode.h" +#include "functions.h" + +int main() +{ + // Namespace declarations + using std::cout; + using std::vector; + using std::complex; + + // Calculate machine precision + Floattype eps_machine = 1.0f; + while (1.0f + eps_machine != 1.0f) + eps_machine /= 2.0f; + + cout << "\n\033[1;33m--- Part B: Solving along an imaginary path ---\033[0m\n"; + complex<Floattype> a(0.0f, 0.0f); // Lower limit + complex<Floattype> b(0.0f, 2.0f*M_PI); // Upper limit + cout << "Integration path: b-a = " << b-a << '\n'; + Inttype n_eqs = 2; // Number of equations in ODE system + vector<complex<Floattype> > y_start(n_eqs); + complex<Floattype> y0(0.0f, 0.0f); + complex<Floattype> y1(1.0f, 1.0f); + y_start[0] = y0; + y_start[1] = y1; + Floattype h_start = 0.01f; + ODE imagode(func1, // ODE system + y_start, // Initial values + a, // Lower limit + b, // Upper limit + h_start, // Start value of step size + 10000, // Max. number of steps + eps_machine*1e12f, // Absolute precision + eps_machine*1e12f); // Relative precision + imagode.write("funcB.dat"); // Write solutions to data file + + // Report to stdout + cout << "ODE system solved in " + << imagode.steps() << " steps.\n\n"; + + // Return successfully + return 0; +} + diff --git a/exam/mainC.cpp b/exam/mainC.cpp @@ -0,0 +1,50 @@ +#include <iostream> +#include <vector> +#include <complex> +#include <cmath> +#include "typedefs.h" +#include "check.h" +#include "ode.h" +#include "functions.h" + +int main() +{ + // Namespace declarations + using std::cout; + using std::vector; + using std::complex; + + // Calculate machine precision + Floattype eps_machine = 1.0f; + while (1.0f + eps_machine != 1.0f) + eps_machine /= 2.0f; + + cout << "\n\033[1;33m--- Part C: Solving along a complex path ---\033[0m\n"; + complex<Floattype> a(0.0f, 0.0f); // Lower limit + complex<Floattype> b(2.0f*M_PI, 2.0f*M_PI); // Upper limit + cout << "Integration path: b-a = " << b-a << '\n'; + Inttype n_eqs = 2; // Number of equations in ODE system + vector<complex<Floattype> > y_start(n_eqs); + complex<Floattype> y0(0.0f, 0.0f); + complex<Floattype> y1(1.0f, 1.0f); + y_start[0] = y0; + y_start[1] = y1; + Floattype h_start = 0.01f; + ODE ode(func1, // ODE system + y_start, // Initial values + a, // Lower limit + b, // Upper limit + h_start, // Start value of step size + 10000, // Max. number of steps + eps_machine*1e12f, // Absolute precision + eps_machine*1e12f); // Relative precision + ode.write("funcC.dat"); // Write solutions to data file + + // Report to stdout + cout << "ODE system solved in " + << ode.steps() << " steps.\n\n"; + + // Return successfully + return 0; +} + diff --git a/exam/mainD.cpp b/exam/mainD.cpp @@ -0,0 +1,68 @@ +#include <iostream> +#include <vector> +#include <complex> +#include <cmath> +#include <fstream> +#include "typedefs.h" +#include "check.h" +#include "ode.h" +#include "functions.h" + + +int main() +{ + // Namespace declarations + using std::cout; + using std::vector; + using std::complex; + + // Calculate machine precision + Floattype eps_machine = 1.0f; + while (1.0f + eps_machine != 1.0f) + eps_machine /= 2.0f; + + cout << "\n\033[1;33m--- Part D: Precision analysis ---\033[0m\n"; + complex<Floattype> a(0.0f, 0.0f); // Lower limit + complex<Floattype> b(2.0f*M_PI, 2.0f*M_PI); // Upper limit + cout << "Integration path: b-a = " << b-a << '\n'; + Inttype n_eqs = 2; // Number of equations in ODE system + vector<complex<Floattype> > y_start(n_eqs); + complex<Floattype> y0(0.0f, 0.0f); + complex<Floattype> y1(1.0f, 1.0f); + y_start[0] = y0; + y_start[1] = y1; + Floattype h_start = 0.01f; + + vector<Floattype> precs; // Vector containing precision values + vector<Inttype> steps; // Vector containing number of steps required + + for (Floattype prec=eps_machine*10.0f; prec<0.1f; prec *= 10.0f) { + ODE ode(func1, // ODE system + y_start, // Initial values + a, // Lower limit + b, // Upper limit + h_start, // Start value of step size + 100000, // Max. number of steps + prec, // Absolute precision + prec); // Relative precision + precs.push_back(prec); // Save precision + steps.push_back(ode.steps()); // Save number of steps taken + } + + // Save results to text file + std::ofstream ost; // Out stream object + ost.open("funcD.dat"); // Open outfile for write + if (!ost) { + std::cerr << "Error, can't open output file!\n"; + return 1; + } + // Write output values + for (Inttype i=0; i<precs.size(); ++i) + ost << precs[i] << '\t' << steps[i] << '\n'; + + // Close file + ost.close(); + + // Return successfully + return 0; +} diff --git a/exam/ode.cpp b/exam/ode.cpp @@ -0,0 +1,163 @@ +#include <iostream> +#include <vector> +#include <cmath> // for sqrt and pow +#include <fstream> +#include "typedefs.h" +#include "ode.h" +#include "vector_arithmetic.h" + +// Constructor +ODE::ODE(std::vector<std::complex<Floattype> > + (*f_in)(const std::complex<Floattype> x, + const std::vector<std::complex<Floattype> > &y), + const std::vector<std::complex<Floattype> > y_start, + const std::complex<Floattype> a_in, + const std::complex<Floattype> b_in, + const Floattype h_start, + const Inttype max_steps, + const Floattype delta_in, + const Floattype epsilon_in, + const Floattype power_in, + const Floattype safety_in) + : f(f_in), + a(a_in), b(b_in), + h((b_in-a_in)*h_start), + n_max(max_steps), + delta(delta_in), epsilon(epsilon_in), + power(power_in), safety(safety_in) +{ + x_list.push_back(a); + y_list.push_back(y_start); + + // Perform integration + rkdriver(); +} + + +// Estimate tolerance +Floattype ODE::tau(const std::vector<std::complex<Floattype> > &y, + const std::complex<Floattype> h_i) +{ + return abs((epsilon * cnorm(y) + delta) * sqrt(h_i/(b-a))); +} + +// Runge-Kutta mid-point stepper +void ODE::rkstep12(const std::complex<Floattype> x0, + const std::vector<std::complex<Floattype> > &y0, + std::vector<std::complex<Floattype> > &y1, + std::vector<std::complex<Floattype> > &dy) +{ + // Denominator 2 used in arithmetic operations + const std::complex<Floattype> den2 (2.0f,2.0f); + + // Evaluate function at two points + (void)f(x0,y0); + const std::vector<std::complex<Floattype> > k0 = f(x0,y0); + const std::vector<std::complex<Floattype> > k12 = f(x0 + h/den2, y0 + k0*h/den2); + + // Write results to output vectors + y1 = y0 + k12*h; + dy = (k0 - k12) * h/den2; +} + + +// ODE driver with adaptive step size control +void ODE::rkdriver() +{ + std::vector<std::complex<Floattype> > dy(n_max); + std::vector<std::complex<Floattype> > y1(n_max); + + std::complex<Floattype> x; + Floattype err, tol; + std::vector<std::complex<Floattype> > y; + + while (x_list.back().real() < b.real() + || x_list.back().imag() < b.imag()) + { + // Get values for this step + x = x_list.back(); + y = y_list.back(); + + // Make sure we don't step past the upper boundary + if ((x + h).real() > b.real() + || (x + h).imag() > b.imag()) + h = b - x; + + // Run Runge-Kutta mid-point stepper + rkstep12(x, y, y1, dy); + + // Determine whether the step should be accepted + err = cnorm(dy); // Error estimate + tol = tau(y, h); // Tolerance + if (err < tol) { // Step accepted + x_list.push_back(x+h); + y_list.push_back(y1); + } + + // Check that we havn't hit the max. number of steps + if (x_list.size() == n_max) { + std::cerr << "Error, the max. number of steps " + << "was insufficient\n" + << "Try either increasing the max. number " + << "of steps, or decreasing the precision " + << "requirements\n"; + return; + } + + // Determine new step size + std::complex<Floattype> multiplicator (2.0f, 2.0f); + if (err > 0.0f) + h = h*pow(tol/err, power) * safety; + else + h = multiplicator*h; + } +} + + +// Return the number of steps taken +Inttype ODE::steps() +{ + return x_list.size(); +} + +void ODE::print() +{ + for (Inttype i=0; i<x_list.size(); ++i) { + std::cout << x_list[i] << '\t'; + for (Inttype j=0; j<y_list[0].size(); ++j) + std::cout << y_list[i][j] << '\t'; + std::cout << '\n'; + } +} + +// Write the x- and y-values to file +void ODE::write(const char* filename) +{ + std::ofstream outstream; + + // Open outfile for write + outstream.open(filename); + if (!outstream) { + std::cerr << "Error, can't open output file '" + << filename << "'.\n"; + return; + } + + // Write output values + for (Inttype i=0; i<x_list.size(); ++i) { + outstream << x_list[i].real() << '\t'; + outstream << x_list[i].imag() << '\t'; + for (Inttype j=0; j<y_list[0].size(); ++j) { + outstream << y_list[i][j].real() << '\t'; + outstream << y_list[i][j].imag() << '\t'; + } + outstream << '\n'; + } + + // Close file + outstream.close(); + + if (verbose == true) + std::cout << "Output written in '" << filename << "'.\n"; +} + diff --git a/exam/ode.h b/exam/ode.h @@ -0,0 +1,87 @@ +// Make sure header is only included once +#ifndef ODE_H_ +#define ODE_H_ + +#include <vector> +#include <complex> +#include "typedefs.h" + +// ODE class +class ODE { + + // Values and functions only accessible from the class internally + private: + + // System of ordinary differential equations to solve + std::vector<std::complex<Floattype> > + (*f)(const std::complex<Floattype> x, + const std::vector<std::complex<Floattype> > &y); + + // Points to be evaluated + std::vector<std::complex<Floattype> > x_list; + + // Limits of range to evaluate + const std::complex<Floattype> a; // Lower + const std::complex<Floattype> b; // Upper + + // Step size + std::complex<Floattype> h; + + // Results stored in 2D: vector of vectors + std::vector<std::vector<std::complex<Floattype> > > y_list; + + // Maximum number of steps to evaluate, defined by y size + const Inttype n_max; + + // Accuracy requirement values + const Floattype delta; // Absolute + const Floattype epsilon; // Relative + + // Tolerance estimator + Floattype tau(const std::vector<std::complex<Floattype> > &y, + const std::complex<Floattype> h); + + // Runge-Kutta mid-point stepper prototype + void rkstep12(const std::complex<Floattype> x0, + const std::vector<std::complex<Floattype> > &y0, + std::vector<std::complex<Floattype> > &y1, + std::vector<std::complex<Floattype> > &dy); + + // Runge-Kutta driver function parameters + const Floattype power; + const Floattype safety; + + // Runge-Kutta driver prototype + void rkdriver(); + + + // Values and functions accessible from the outside + public: + + // Constructor, some parameters with default values + ODE(std::vector<std::complex<Floattype> > + (*f_in)(const std::complex<Floattype> x, + const std::vector<std::complex<Floattype> > &y), + const std::vector<std::complex<Floattype> > y_start, + const std::complex<Floattype> a_in, + const std::complex<Floattype> b_in, + const Floattype h_start = 0.01f, + const Inttype max_steps = 1e4, + const Floattype delta_in = 1e-3f, + const Floattype epsilon_in = 1e-3f, + const Floattype power_in = 0.25f, + const Floattype safety_in = 0.95f + ); + + // Return the number of steps taken + Inttype steps(); + + // Print the x- and y-values to stdout + void print(); + + // Write the x- and y-values to file + void write(const char* filename); + +}; + +#endif diff --git a/exam/plotA.gp b/exam/plotA.gp @@ -0,0 +1,8 @@ +set terminal png +set output "plotA.png" +set xlabel "x.real" +set ylabel "y.real" +set title "Integration along a real path" +set grid +plot "funcA.dat" u 1:3 t "rk12" w lp + diff --git a/exam/plotB.gp b/exam/plotB.gp @@ -0,0 +1,8 @@ +set terminal png +set output "plotB.png" +set xlabel "x.imag" +set ylabel "y.imag" +set title "Integration along an imaginary path" +set grid +plot "funcB.dat" u 2:4 t "rk12" w lp + diff --git a/exam/plotC.gp b/exam/plotC.gp @@ -0,0 +1,9 @@ +set terminal png +set output "plotC.png" +set xlabel "x.r" +set ylabel "x.i" +set zlabel "y" +set title "Integration along a complex path" +set grid +splot "funcC.dat" u 1:2:3 t "y.real" w lp, "" u 1:2:4 t "y.imag" w lp + diff --git a/exam/plotD.gp b/exam/plotD.gp @@ -0,0 +1,9 @@ +set terminal png +set output "plotD.png" +set xlabel "Precision" +set ylabel "No. of steps" +set title "Number of steps required for a given absolute and relative precision" +set grid +set log xy +plot "funcD.dat" u 1:2 t "rk12" w lp + diff --git a/exam/typedefs.h b/exam/typedefs.h @@ -0,0 +1,15 @@ +// Make sure header is only included once +#ifndef TYPEDEFS_H_ +#define TYPEDEFS_H_ + +// Define whether the program should output values of matrices +//const bool verbose = false; +const bool verbose = true; + +// Choose length variable type +typedef unsigned int Inttype; + +// Choose floating-point precision +typedef double Floattype; + +#endif diff --git a/exam/vector_arithmetic.h b/exam/vector_arithmetic.h @@ -0,0 +1,89 @@ +// Make sure header is only included once +#ifndef VECTOR_ARITHMETIC_H_ +#define VECTOR_ARITHMETIC_H_ + +#include <vector> +#include <cmath> // for sqrt +#include "typedefs.h" + +//// Overload vector methods to allow scalar- +//// and element-wise arithmetic operations + +// Scalar multiplication (same scalar for real and imaginary parts) +std::vector<std::complex<Floattype> > + operator*(const std::vector<std::complex<Floattype> > vec, + const Floattype scalar) +{ + std::vector<std::complex<Floattype> > result(vec.size()); + for (Inttype i=0; i<vec.size(); ++i) { + result[i].real() = real(vec[i])*scalar; + result[i].imag() = imag(vec[i])*scalar; + } + return result; +} + +// Scalar multiplication +std::vector<std::complex<Floattype> > + operator*(const std::vector<std::complex<Floattype> > vec, + const std::complex<Floattype> scalar) +{ + std::vector<std::complex<Floattype> > result(vec.size()); + for (Inttype i=0; i<(Inttype)vec.size(); ++i) + result[i] = vec[i]*scalar; + return result; +} + +// Scalar division +std::vector<std::complex<Floattype> > + operator/(const std::vector<std::complex<Floattype> > vec, + const std::complex<Floattype> scalar) +{ + std::vector<std::complex<Floattype> > result(vec.size()); + for (Inttype i=0; i<(Inttype)vec.size(); ++i) + result[i] = vec[i]/scalar; + return result; +} + +// Element-wise addition +std::vector<std::complex<Floattype> > + operator+(const std::vector<std::complex<Floattype> > vec1, + const std::vector<std::complex<Floattype> > vec2) +{ + std::vector<std::complex<Floattype> > result(vec1.size()); + for (Inttype i=0; i<(Inttype)vec1.size(); ++i) + result[i] = vec1[i] + vec2[i]; + return result; +} + +// Element-wise subtraction +std::vector<std::complex<Floattype> > + operator-(const std::vector<std::complex<Floattype> > vec1, + const std::vector<std::complex<Floattype> > vec2) +{ + std::vector<std::complex<Floattype> > result(vec1.size()); + for (Inttype i=0; i<(Inttype)vec1.size(); ++i) + result[i] = vec1[i] - vec2[i]; + return result; +} + +// Element-wise multiplication +std::vector<std::complex<Floattype> > + operator*(const std::vector<std::complex<Floattype> > vec1, + const std::vector<std::complex<Floattype> > vec2) +{ + std::vector<std::complex<Floattype> > result(vec1.size()); + for (Inttype i=0; i<(Inttype)vec1.size(); ++i) + result[i] = vec1[i] * vec2[i]; + return result; +} + +// Normalize vector +Floattype cnorm(const std::vector<std::complex<Floattype> > vec) +{ + Floattype res = 0.0f; + for (Inttype i=0; i<(Inttype)vec.size(); ++i) + res += norm(vec[i])*norm(vec[i]); + return sqrt(res); +} + +#endif diff --git a/hello/Makefile b/hello/Makefile @@ -0,0 +1,71 @@ +# Define compiler +CC=g++ + +# Define compiler flags (show all warnings) +CCFLAGS=-Wall + +# Define linker +LD=g++ + +# Define linker flags +LDFLAGS= + +# Filenames of source code +ASRC=$(shell ls *.A.cpp) +BSRC=$(shell ls *.B.cpp) + +# Filenames of object files +BOBJ=$(BSRC:.cpp=.o) + +# Remove file type extension for binary filename +ABIN=$(ASRC:.cpp=) +BBIN=hello_user.B + +# The default "all" depends on A and B + +all: A B + +correct: + # Cenerating correct answer + @echo "Hello, $(USER)" > correct.txt + +A: + # Compile source code + $(CC) $(CCFLAGS) $(ASRC) -o $(ABIN) + # Execute program and redirect stdout to file + ./$(ABIN) > out.A.txt + +B: $(BOBJ) + # Link object files together + $(LD) $(LDFLAGS) $(BOBJ) -o $(BBIN) + # Execute program and redirect stdout to file + ./$(BBIN) > out.B.txt + +# Object files for B require B source code files +.B.o: $(BSRC) + $(CC) $(CCFLAGS) -c $< -o $@ + +clean-A: + # Remove binary + rm -f $(ABIN) + +clean-B: + # Remove object files + rm -f $(BOBJ) + # Remove binary + rm -f $(BBIN) + +clean: clean-A clean-B + rm -f correct.txt + +test-A: correct A + @echo "Checking A... " + @diff --brief correct.txt out.A.txt + @echo "done" + +test-B: correct B + @echo "Checking B... " + @diff --brief correct.txt out.B.txt + @echo "done" + +test: test-A test-B diff --git a/hello/hello.B.cpp b/hello/hello.B.cpp @@ -0,0 +1,6 @@ +#include <iostream> + +void hello() +{ + std::cout << "Hello, "; +} diff --git a/hello/hello_user.A.cpp b/hello/hello_user.A.cpp @@ -0,0 +1,19 @@ +#include <iostream> +#include <string> +#include <cstdlib> // For getenv() + +int main() +{ + // Namespace definitions + using std::cout; + using std::string; + + // Get USER environment variable + string username = getenv("USER"); + + // Greet user + cout << "Hello, " << username << '\n'; + + // Successfull exit + return 0; +} diff --git a/hello/main.B.cpp b/hello/main.B.cpp @@ -0,0 +1,14 @@ +#include <iostream> + +// Prototype functions +#include "main.B.h" + +int main() +{ + // Call functions from object files + hello(); + your_username(); + + // Successfull exit + return 0; +} diff --git a/hello/main.B.h b/hello/main.B.h @@ -0,0 +1,8 @@ +// Header file containing prototype functions + +// hello.B.cpp +void hello(); + +// your_username.B.cpp +void your_username(); + diff --git a/hello/your_username.B.cpp b/hello/your_username.B.cpp @@ -0,0 +1,14 @@ +#include <iostream> +#include <string> +#include <cstdlib> // For getenv() + +void your_username() +{ + + // Get USER environment variable + std::string username = getenv("USER"); + + // Print username to stdout + std::cout << username << '\n'; + +} diff --git a/integration/Makefile b/integration/Makefile @@ -0,0 +1,71 @@ +# Define compiler +CC=g++ + +# Define compiler flags (show all warnings) +CPPFLAGS=-Wall +#CPPFLAGS=-std=c++0x + +# Define linker flags +LDFLAGS= + +# Define extra libraries to be dynamically linked +#LDLIBS+=-larmadillo +#LDLIBS+=-lstdc++ + +# Compile optimized code +#CPPFLAGS+=-O2 + +# Compile debuggable code +#CPPFLAGS+=-g + +# Compile profilable code +#CPPFLAGS+=-pg +#LDFLAGS+=-pg + +# Define linker +LD=g++ + +# Filenames of source code +SRC=$(shell ls *.cpp) + +# Filenames of object files +OBJ=$(SRC:.cpp=.o) +#OBJ=downhill_simplex.o main.o +#OBJB=downhill_simplex.o main.B.o + +# Remove file type extension for binary filename +BIN=integ +#BINB=optimB + +# The default "all" depends on A and B + +#all: A B + + +A: $(BIN) + @# Run program and write amoeba positions to file + ./$(BIN) + +#B: $(BINB) +# ./$(BINB) 2> amoebaB.dat + +$(BIN): $(OBJ) + @# Link object files together + $(LD) $(LDFLAGS) $(OBJ) -o $(BIN) $(LDLIBS) + +$(BINB): $(OBJB) + @# Link object files together + $(LD) $(LDFLAGS) $(OBJB) -o $(BINB) $(LDLIBS) + +clean: + @# Remove object files + rm -f $(OBJ) $(OBJB) + @# Remove binaries + rm -f $(BIN) $(BINB) + @# Remove datafiles and plot + #rm -f *.dat *.png + @# Remove error logs + @#rm -f amoeba.dat amoebaB.dat +edit: + vim -p Makefile *.cpp *.h + diff --git a/integration/functions.h b/integration/functions.h @@ -0,0 +1,27 @@ +#ifndef FUNCTIONS_H_ +#define FUNCTIONS_H_ + +#include <cmath> +#include "header.h" + +int ncalls = 0; + +Floattype func1(const Floattype x) +{ + ++ncalls; + return x*x; +} + +Floattype func2(const Floattype x) +{ + ++ncalls; + return log(x)/sqrt(x); +} + +Floattype func3(const Floattype x) +{ + ++ncalls; + return 4.0f*sqrt(1.0f-(1.0f-x)*(1.0f-x)); +} + +#endif diff --git a/integration/header.h b/integration/header.h @@ -0,0 +1,18 @@ +// Make sure header is only included once +#ifndef HEADER_H_ +#define HEADER_H_ + +// Define whether the program should output values of matrices +const bool verbose = false; +//const bool verbose = true; + +// Prototype for checking function +void check(const bool statement); + +// Choose length variable type +typedef long int Lengthtype; + +// Choose floating-point precision +typedef double Floattype; + +#endif diff --git a/integration/integrator.cpp b/integration/integrator.cpp @@ -0,0 +1,120 @@ +#include <iostream> +#include <cmath> +#include <vector> +#include <string> +#include <typeinfo> +#include "integrator.h" +#include "header.h" + +// Constructor +Integral::Integral(Floattype (*function)(const Floattype), + const Floattype lower_limit, + const Floattype upper_limit, + const Floattype absolute_accuracy, + const Floattype relative_accuracy) + : f_in(function), acc_pres(absolute_accuracy), eps(relative_accuracy) +{ + // Initialize number of recursions to 0 + nrec = 0; + + // Test whether input interval has infinite limits + if (std::isinf(lower_limit) == 1 && std::isinf(upper_limit) == 1) { + f = Integral::infinf; + low = 0.0f; high = 1.0f; + } else if (std::isinf(lower_limit) == 0 && std::isinf(upper_limit == 1)) { + f = Integral::fininf; + low = 0.0f; high = 1.0f; + } else if (std::isinf(lower_limit) == 1 && std::isinf(upper_limit == 0)) { + f = &Integral::inffin; + low = 0.0f; high = 1.0f; + } else { + f = Integral::f_in; + low = lower_limit; + high = upper_limit; + } + + Floattype f2 = f(low + 2.0f * (high-low)/6.0f); + Floattype f3 = f(low + 4.0f * (high-low)/6.0f); + + res = adapt(low, high, acc_pres, f2, f3); +} + +// Adaptive integrator, to be called recursively +Floattype Integral::adapt(const Floattype a, + const Floattype b, + const Floattype acc, + const Floattype f2, + const Floattype f3) +{ + if (nrec > 2147483647) + return NAN; // Return NaN if recursions seem infinite + + // Function value at end points + Floattype f1 = f(a + 1.0f * (b-a)/6.0f); + Floattype f4 = f(b + 5.0f * (b-a)/6.0f); + + // Quadrature rules + Floattype Q = (2.0f*f1 + f2 + f3 + 2.0f*f4) / 6.0f * (b-a); + Floattype q = (f1 + f2 + f3 + f4) / 4.0f * (b-a); + + Floattype tolerance = acc + eps*fabs(Q); + err = fabs(Q-q); + + // Evaluate whether result is precise + // enough to fulfill criteria + if (err < tolerance) + return Q; + else { // Perform recursions in lower and upper interval + ++nrec; + Floattype q1 = adapt(a, a+(b-a)/2.0f, acc/sqrt(2), f1, f2); + ++nrec; + Floattype q2 = adapt(a+(b-a)/2.0f, b, acc/sqrt(2), f3, f4); + return q1+q2; + } +} + +// Return result +Floattype Integral::result() +{ + return res; +} + +// Return the number of recursions taken +Lengthtype Integral::recursions() +{ + return nrec; +} + +// Print results of function integration +void Integral::show(const std::string function_description) +{ + std::cout << "Integral of function '" + << function_description + << "' in interval [" + << low << ";" << high << "] is " + << res << ",\n" + << "with an absolute accuracy of " + << acc_pres + << " and a relative accuracy of " + << eps << ".\n" + << "Integral calculated in " + << nrec << " recursions, error is " + << err << ".\n"; +} + +// Functions for variable transformations when limits are infinite +Floattype Integral::infinf(const Floattype t) +{ + return (f_in((1.0f - t) / t) + f_in(-(1.0f - t) / t)) / (t*t); +} + +Floattype Integral::fininf(const Floattype t) +{ + return f_in(low + (1.0f - t) / t) / (t*t); +} + +Floattype Integral::inffin(const Floattype t) +{ + return f_in(high - (1.0f - t) / t) / (t*t); +} + diff --git a/integration/integrator.h b/integration/integrator.h @@ -0,0 +1,68 @@ +#ifndef INTEGRATOR_H_ +#define INTEGRATOR_H_ + +#include <string> +#include "header.h" + +class Integral { + + private: + + // Input function (maybe to be transformed) + Floattype (*f_in)(const Floattype); + + // Function to be evaluated (pointer to function) + Floattype (*f)(const Floattype); + + // Integral limits + Floattype low; // Lower + Floattype high; // Upper + + // Accuracy requirement values + const Floattype acc_pres; // Absolute + const Floattype eps; // Relative + + // Number of recursions + Lengthtype nrec; + + // Adaptive integrator function to be + // called recursively + Floattype adapt(const Floattype a, + const Floattype b, + const Floattype acc, + const Floattype f2, + const Floattype f3); + + // Result of integral + Floattype res; + + // Size of error + Floattype err; + + // Functions for variable transformations + // when one or both limits are infinite + Floattype infinf(const Floattype t); // a=inf + Floattype fininf(const Floattype t); // b=inf + Floattype inffin(const Floattype t); // a=inf,b=inf + + + public: + + // Constructor + Integral(Floattype (*function)(const Floattype), + const Floattype lower_limit, + const Floattype upper_limit, + const Floattype absolute_accuracy, + const Floattype relative_accuracy); + + // Return value of result + Floattype result(); + + // Return number of recursions taken + Lengthtype recursions(); + + // Print results and statistics + void show(const std::string function_description); +}; + +#endif diff --git a/integration/main.cpp b/integration/main.cpp @@ -0,0 +1,57 @@ +#include <iostream> +#include "header.h" +#include "functions.h" +#include "integrator.h" + +int main(int argc, char* argv[]) +{ + // Namespace declarations + using std::cout; + + // Calculate machine precision + Floattype eps_machine = 1.0f; + while (1.0f + eps_machine != 1.0f) + eps_machine /= 2.0f; + + Floattype a, b, acc, eps; + + // Find the integral of function 1 + a = 0.0f, b = 1.0f; + acc = 0.001f; eps = acc; + ncalls = 0; + Integral F1(func1, a, b, acc, eps); + F1.show("f(x) = x*x"); + cout << '\n'; + + // Find the integral of function 2 + a = 1e-8f, b = 1.0f; + acc = 0.001f; eps = acc; + ncalls = 0; + Integral F2(func2, a, b, acc, eps); + F2.show("f(x) = log(x)/sqrt(x)"); + cout << "Check: Integral equal to -4:\t"; + check(fabs(F2.result() + 4.0f) < 0.1f); + cout << '\n'; + + // Find the integral of function 3 + a = 0.0f, b = 1.0f; + acc = 1e-8f; eps = acc; + ncalls = 0; + Integral F3(func3, a, b, acc, eps); + F3.show("f(x) = 4*sqrt(1-(1-x)^2)"); + cout << "Check: Integral equal to pi:\t"; + check(fabs(F3.result() - M_PI) < 0.01f); + cout << '\n'; + + // Return successfully + return 0; +} + +void check(const bool statement) +{ + using std::cout; + if (statement == true) + cout << "\t\033[0;32mPassed\033[0m\n"; + else + cout << "\t\033[1;31mFail!!\033[0m\n"; +} diff --git a/leastsq/Makefile b/leastsq/Makefile @@ -0,0 +1,67 @@ +# Define compiler +CC=g++ + +# Define compiler flags (show all warnings) +CPPFLAGS=-Wall + +# Define linker flags +LDFLAGS= + +# Define extra libraries to be dynamically linked +LDLIBS+=-larmadillo + +# Compile optimized code +#CPPFLAGS+=-O2 + +# Compile debuggable code +#CPPFLAGS+=-g + +# Compile profilable code +#CPPFLAGS+=-pg +#LDFLAGS+=-pg + +# Define linker +LD=g++ + +# Filenames of source code +SRC=$(shell ls *.cpp) + +# Filenames of object files +OBJ=$(SRC:.cpp=.o) + +# Remove file type extension for binary filename +BIN=leastsq + +# The default "all" depends on A and B + +all: A B + +A: plot.A.png + +B: plot.B.png + +plot.%.png: fit.A.dat fit.B.dat data.A.txt data.B.txt + gnuplot plotall.gp + +fit.A.dat: $(BIN) + ./$(BIN) data.A.txt fit.A.dat + +fit.B.dat: $(BIN) + ./$(BIN) data.B.txt fit.B.dat + +$(BIN): $(OBJ) + @# Link object files together + $(LD) $(LDFLAGS) $(OBJ) -o $(BIN) $(LDLIBS) + @# Execute program and redirect stdout to file + @#./$(BIN) > out.txt + +clean: + @# Remove object files + rm -f $(OBJ) + @# Remove binary + rm -f $(BIN) + @# Remove datafiles and plot + rm -f *.dat *.png +edit: + vim -p Makefile *.cpp *.h *.gp + diff --git a/leastsq/data.A.txt b/leastsq/data.A.txt @@ -0,0 +1,7 @@ +20.0 10.26 2.2 +40.0 10.52 4.12 +60.0 9.22 5.2 +80.0 17.77 3.5 +100.0 38.14 5.21 +120.0 45.56 1.2 +140.0 54.29 6.2 diff --git a/leastsq/data.B.txt b/leastsq/data.B.txt @@ -0,0 +1,10 @@ +0 1 0.2 +1 2 0.2 +2 3 0.2 +3 2 0.2 +4 4 0.2 +5 3 0.2 +6 2.5 0.2 +7 3.5 0.2 +8 4 0.2 +9 5 0.2 diff --git a/leastsq/file_o.cpp b/leastsq/file_o.cpp @@ -0,0 +1,25 @@ +#include <iostream> +#include <fstream> +#include <vector> +#include "header.h" + +// Write arrays to text file, readable by plotting tools +int write_output(const std::vector<Floattype> &a, + const std::vector<Floattype> &b, + const char *outfile) +{ + std::ofstream outstream; + + // Open outfile for write + outstream.open(outfile); + + // Write output arrays in two columns + for (Lengthtype i=0; i<(Lengthtype)a.size(); ++i) + outstream << a[i] << '\t' << b[i] << '\n'; + + // Close file + outstream.close(); + + // Return successfully + return 0; +} diff --git a/leastsq/header.h b/leastsq/header.h @@ -0,0 +1,21 @@ +// Make sure header is only included once +#ifndef HEADER_H_ +#define HEADER_H_ + +#include <vector> + +// Define whether the program should output values of matrices +const bool verbose = false; + +// Choose vector length variable type +typedef int Lengthtype; + +// Choose floating-point precision +typedef double Floattype; + +// File output function prototype +int write_output(const std::vector<Floattype> &a, + const std::vector<Floattype> &b, + const char *outfile); + +#endif diff --git a/leastsq/lsfit.cpp b/leastsq/lsfit.cpp @@ -0,0 +1,61 @@ +#include <iostream> +#include <cmath> // NaN is here +#include <armadillo> +#include "header.h" +#include "lsfit.h" +#include "qrfunc.h" + +// A number of i fitting functions, evaluated at x +Floattype LSfit::fitfuncts(char i, Floattype x) +{ + // Choose function i + switch (i) { + case 0: + return 1.0f; + break; + case 1: + return x; + break; + case 2: + return x*x; + break; + default: + std::cout << "Wrong function (i = " + << i << ") specified in fitfuncts, x = " + << x << '\n'; + return NAN; + } +} + +// Constructor +LSfit::LSfit(arma::Col<Floattype> &x, + arma::Col<Floattype> &y, + arma::Col<Floattype> &delta_y) : n(x.n_rows), m(3), A(n,m), b(n), c(m) +{ + + // Initialize b values + b = y/delta_y; + + // Initialize A values + for (Lengthtype i=0; i<n; ++i) { + for (Lengthtype k=0; k<m; ++k) { + A(i,k) = fitfuncts(k,x(i)) / delta_y(i); + } + } + + // Perform QR decomposition + QR qr(A); + + // Calculate the least-squares solution + c = qr.backsub(b); +} + +// Evaluate function at x +Floattype LSfit::eval(Floattype x) +{ + // Find the linear combination of m functions to calculate F(x) + Floattype sum = 0.0f; + for (char k=0; k<m; ++k) + sum += c(k) * fitfuncts(k, x); + return sum; +} diff --git a/leastsq/lsfit.h b/leastsq/lsfit.h @@ -0,0 +1,41 @@ +// Make sure header is only included once per object +#ifndef LSFIT_H_ +#define LSFIT_H_ + +#include <armadillo> +#include "header.h" + +// lsfit structure +class LSfit { + private: + const Lengthtype n; // Input data count + const Lengthtype m; // Fitting function count + + // Covariance matrix + arma::Mat<Floattype> A; + + // Data points + arma::Col<Floattype> b; + + // Fitting coefficients + arma::Col<Floattype> c; + + // Evaluate fitting function i at x + Floattype fitfuncts(char i, Floattype x); + + public: + + // Constructor. Arguments: input data points + LSfit(arma::Col<Floattype> &x, + arma::Col<Floattype> &y, + arma::Col<Floattype> &delta_y); + + // Destructor + //~LSfit(); + + // Return value fitted at x + Floattype eval(Floattype x); + +}; + +#endif diff --git a/leastsq/main.cpp b/leastsq/main.cpp @@ -0,0 +1,89 @@ +#include <iostream> +#include <cstdio> // For fscanf +#include <armadillo> +#include <vector> +#include "header.h" +#include "lsfit.h" + +int main(int argc, char* argv[]) +{ + // Namespace declarations + using std::cout; + + // Check that a data file is given as an input argument + if (argc < 2 || argc > 3) { + cout << "Usage: " << argv[0] << " <input file> [output file]\n" + << "If 'output file' is not specified, output data " + << "will be written to stdout.\n" + << "Example: " << argv[0] << " data.A.txt\n"; + return 1; + } + + FILE *fin; // File pointer + + // Check that the input file exists and can be read + if ((fin = fopen(argv[1], "r")) == NULL) { + cout << "Error while reading " << argv[1] << '\n'; + return 1; + } + + // First, count the number of newline characters + // for preallocation purposes + Lengthtype n = 0; + int c; + while ((c = getc(fin)) != EOF) { + if (c == '\n') + ++n; + } + fclose(fin); + cout << "Input file \"" << argv[1] << "\" consists of n=" + << n << " data points.\n"; + + // Allocate input data structures + arma::Col<Floattype> x = arma::Col<Floattype> (n); + arma::Col<Floattype> y = arma::Col<Floattype> (n); + arma::Col<Floattype> delta_y = arma::Col<Floattype> (n); + + // Read data into memory + if ((fin = fopen(argv[1], "r")) == NULL) { + cout << "Error while reading " << argv[1] << '\n'; + return 1; + } + float x_tmp, y_tmp, delta_y_tmp; + for (Lengthtype i=0; i<n; ++i) { + fscanf(fin, "%f %f %f", &x_tmp, &y_tmp, &delta_y_tmp); + x(i) = x_tmp; + y(i) = y_tmp; + delta_y(i) = delta_y_tmp; + } + fclose(fin); + + // Perform least-squares fit with LSfit class + LSfit lsfit(x, y, delta_y); + + // Evaluate fit at a fine resolution + const unsigned int res = 100; + Floattype x_min = x.min(); + cout << "x_min = " << x_min; + Floattype x_max = x.max(); + cout << ", x_max = " << x_max << '\n'; + std::vector<Floattype> xo (res); + std::vector<Floattype> yo (res); + for (unsigned int i=0; i<res; ++i) { + xo[i] = x_min + (x_max - x_min) * ((Floattype)i/res); + yo[i] = lsfit.eval(xo[i]); + } + + // Write to file if specified in as command line arguments + if (argc == 3) { + write_output(xo, yo, argv[2]); + } else { + for (unsigned int i=0; i<res; ++i) { + cout << xo[i] << '\t' << yo[i] << '\n'; + } + } + + + // Return successfully + return 0; +} diff --git a/leastsq/plotall.gp b/leastsq/plotall.gp @@ -0,0 +1,17 @@ +set terminal png # Set output file format +set output "plot.A.png" # Set output filename +set key top left +set xlabel "x" +set ylabel "y" +set title "Least-squares fit A" +set grid +plot "fit.A.dat" u 1:2 w lp title "LSfit", "data.A.txt" w errorbars title "Data" + +set terminal png # Set output file format +set output "plot.B.png" # Set output filename +set key top left +set xlabel "x" +set ylabel "y" +set title "Least-squares fit B" +set grid +plot "fit.B.dat" u 1:2 w lp title "LSfit", "data.B.txt" w errorbars title "Data" diff --git a/leastsq/qrfunc.cpp b/leastsq/qrfunc.cpp @@ -0,0 +1,90 @@ +#include <iostream> +#include <armadillo> +#include "header.h" +#include "qrfunc.h" + +// QR decomposition constructor +QR::QR(arma::mat &A) + : n(A.n_cols), + A(A), + Q(A) +{ + // Initialize output structures + R = arma::zeros<arma::mat> (n,n); + + // Perform QR decomposition straight away + decomp(); +} + +// Class deconstructor (equivalent to compiler destructor) +QR::~QR() { }; + +// Return system size +Lengthtype QR::size() +{ + return n; +} + +// QR decomposition function of Armadillo matrix. +// Returns right matrix R, and modifies A into Q. +// Uses Gram-Schmidt orthogonalization +void QR::decomp() +{ + Floattype r, s; + Lengthtype j; + for (Lengthtype i=0; i<n; ++i) { + r = dot(Q.col(i), Q.col(i)); + R(i,i) = sqrt(r); + Q.col(i) /= sqrt(r); // Normalization + for (j=i+1; j<n; ++j) { + s = dot(Q.col(i), Q.col(j)); + Q.col(j) -= s*Q.col(i); // Orthogonalization + R(i,j) = s; + } + } +} + +// Solve the square system of linear equations with back +// substitution. T is an upper triangular matrix. +arma::vec QR::backsub(arma::vec &b) +{ + Floattype tmpsum; + arma::vec x = Q.t() * b; + for (Lengthtype i=n-1; i>=0; --i) { + tmpsum = 0.0f; + for (Lengthtype k=i+1; k<n; ++k) + tmpsum += R(i,k) * x(k); + + x(i) = 1.0f/R(i,i) * (x(i) - tmpsum); + } + return x; +} + +// Calculate the (absolute value of the) determinant of +// matrix A given the Q and R components. +// det(A) = det(Q) * det(R), |det(Q) = 1| +// => |det(A)| = |det(R)| +Floattype QR::det() +{ + Floattype determinant = 1.0f; + for (Lengthtype i=0; i<n; ++i) + determinant *= R(i,i); + return fabs(determinant); +} + +// Calculate the inverse of matrix A given the Q and R +// components. +arma::mat QR::inverse() +{ + arma::mat inv = arma::zeros<arma::mat> (n, n); + // In vector z, all elements are equal to 0, except z(i) = 1 + arma::vec z = arma::zeros<arma::mat> (n); + + for (Lengthtype i=0; i<n; ++i) { + z(i) = 1.0f; // Element i changed to 1 + inv.col(i) = backsub(z); + z(i) = 0.0f; // Element i changed back to 0 + } + + return inv; +} diff --git a/leastsq/qrfunc.h b/leastsq/qrfunc.h @@ -0,0 +1,49 @@ +// Make sure header is only included once +#ifndef QRFUNC_H_ +#define QRFUNC_H_ + +#include <armadillo> +#include "header.h" + +// QR structure +class QR { + private: + // System size + const Lengthtype n; + + public: + //// Data + + // Input data + arma::mat A; + + // QR decomposition matrices + arma::mat Q; + arma::mat R; + + //// Prototype functions + + // Constructor prototype + QR(arma::mat &A); + + // Destructor + ~QR(); + + // Return system size + Lengthtype size(); + + // QR decomposition of Armadillo matrix A, returning R + // and modified A (=Q) + void decomp(); + + // Backsubstitution of triangular system + arma::vec backsub(arma::vec &b); + + // Absolute value of the determinant of matrix R + Floattype det(); + + // Inverse of matrix A + arma::mat inverse(); +}; + +#endif diff --git a/matrixmul/Makefile b/matrixmul/Makefile @@ -0,0 +1,150 @@ +# Compare speed of different programming languages +# in a matrix multiplication algorithm: AB=C. +# A is uniform with the cell values '2', B should contain +# pseudorandom numbers. +# Initialization of A, B and C in the same loop is allowed, +# but the multiplication should be done in a separate loop. +# Matrix dimensions are specified as a command line argument. + +MATRIXDIMS_SLOW = 100 200 400 800 1600 +MATRIXDIMS = $(MATRIXDIMS_SLOW) 3200 +PREFIXCMD = nice \gtime -ao + +CC=gcc +CXX=g++ +CFLAGS=-Wall -O3 +CXXFLAGS=-Wall -O3 + +performance.png: plot.gp lua-arrofarrs.dat lua-linarr.dat luajit-arrofarrs.dat luajit-linarr.dat c-arrofarrs.dat c-linarr.dat c-omp-arrofarrs.dat c-omp-linarr.dat julia.dat cpp-vectorofvectors.dat cpp-linvectors.dat python-numpy.dat octave.dat + gnuplot plot.gp + +# Lua: Matrices as arrays of arrays +lua-arrofarrs.dat: lua-arrofarrs.lua + # lua-arrofars.lua + @rm -f $@ + @for dims in $(MATRIXDIMS_SLOW); do \ + $(PREFIXCMD) $@ -f "$$dims %e" lua $< $$dims; \ + echo $$dims; \ + done + +# Lua: Matrices as linear arrays +lua-linarr.dat: lua-linarr.lua + # lua-linarr.lua + @rm -f $@ + @for dims in $(MATRIXDIMS_SLOW); do \ + $(PREFIXCMD) $@ -f "$$dims %e" lua $< $$dims; \ + echo $$dims; \ + done + +# LuaJIT: Matrices as arrays of arrays +luajit-arrofarrs.dat: lua-arrofarrs.lua + # luajit-arrofars.lua + @rm -f $@ + @for dims in $(MATRIXDIMS); do \ + $(PREFIXCMD) $@ -f "$$dims %e" luajit $< $$dims; \ + echo $$dims; \ + done + +# LuaJIT: Matrices as linear arrays +luajit-linarr.dat: lua-linarr.lua + # LuaJIT lua-linarr.lua + @rm -f $@ + @for dims in $(MATRIXDIMS); do \ + $(PREFIXCMD) $@ -f "$$dims %e" luajit $< $$dims; \ + echo $$dims; \ + done + +# C: Array of arrays +c-arrofarrs.dat: c-arrofarrs + # c-arrofarrs + @rm -f $@ + @for dims in $(MATRIXDIMS); do \ + $(PREFIXCMD) $@ -f "$$dims %e" ./$< $$dims; \ + echo $$dims; \ + done + +# C: Matrices as linear arrays +c-linarr.dat: c-linarr + # c-linarr + @rm -f $@ + @for dims in $(MATRIXDIMS); do \ + $(PREFIXCMD) $@ -f "$$dims %e" ./$< $$dims; \ + echo $$dims; \ + done + +# C: Array of arrays, OpenMP +c-omp-arrofarrs.dat: c-omp-arrofarrs + # c-omp-arrofarrs + @rm -f $@ + @for dims in $(MATRIXDIMS); do \ + $(PREFIXCMD) $@ -f "$$dims %e" ./$< $$dims; \ + echo $$dims; \ + done + +# C: Matrices as linear arrays, OpenMP +c-omp-linarr.dat: c-omp-linarr + # c-omp-linarr + @rm -f $@ + @for dims in $(MATRIXDIMS); do \ + $(PREFIXCMD) $@ -f "$$dims %e" ./$< $$dims; \ + echo $$dims; \ + done + +c-omp-arrofarrs: c-arrofarrs.c + $(CC) $(CFLAGS) -fopenmp $< -o $@ + +c-omp-linarr: c-linarr.c + $(CC) $(CFLAGS) -fopenmp $< -o $@ + +# Julia, native arrays +julia.dat: julia.jl + # julia.jl + @rm -f $@ + @for dims in $(MATRIXDIMS); do \ + $(PREFIXCMD) $@ -f "$$dims %e" julia $< $$dims; \ + echo $$dims; \ + done + +# C++: Vector of vectors +cpp-vectorofvectors.dat: cpp-vectorofvectors + # cpp-vectorofvectors + @rm -f $@ + @for dims in $(MATRIXDIMS); do \ + $(PREFIXCMD) $@ -f "$$dims %e" ./$< $$dims; \ + echo $$dims; \ + done + +# C++: Linear vectors +cpp-linvectors.dat: cpp-linvectors + # cpp-linvectors + @rm -f $@ + @for dims in $(MATRIXDIMS); do \ + $(PREFIXCMD) $@ -f "$$dims %e" ./$< $$dims; \ + echo $$dims; \ + done + +# Python: Numpy module +python-numpy.dat: python-numpy.py + # python-numpy.py + @rm -f $@ + @for dims in $(MATRIXDIMS); do \ + $(PREFIXCMD) $@ -f "$$dims %e" python $< $$dims; \ + echo $$dims; \ + done + +# Octave +octave.dat: octave.m + # octave.m + @rm -f $@ + @for dims in $(MATRIXDIMS); do \ + $(PREFIXCMD) $@ -f "$$dims %e" octave -qf $< $$dims; \ + echo $$dims; \ + done + +clean: + rm -f *.o + rm -f *.dat + rm -f *.png + +edit: + vim -p Makefile plot.gp diff --git a/matrixmul/c-arrofarrs.c b/matrixmul/c-arrofarrs.c @@ -0,0 +1,72 @@ +#include <stdlib.h> +#include <stdio.h> + +void matrixMult(double** A, double** B, double** C, int N) +{ + int i, j, k; + double sum; +#pragma omp parallel for private (j,k,sum) + for (i = 0; i<N; i++) { + for (j = 0; j<N; j++) { + sum = 0.0f; + for (k = 0; k<N; k++) { + sum += A[i][k] * B[k][j]; + } + C[i][j] = sum; + } + } +} + + +int main(int argc, char* argv[]) +{ + int i, j, N; + double** A; + double** B; + double** C; + + /* Read input argument as matrix height */ + if (argc == 2) { + N = atoi(argv[1]); + } else { + printf("Sorry, I need matrix width as command line argument\n"); + return 1; + } + + /* Allocate arrays */ + A = (double**) malloc(N * sizeof(double*)); + B = (double**) malloc(N * sizeof(double*)); + C = (double**) malloc(N * sizeof(double*)); + + /* Allocate row arrays inside arrays */ + for (i = 0; i < N; i++) { + A[i] = (double*) malloc(N * sizeof(double)); + B[i] = (double*) malloc(N * sizeof(double)); + C[i] = (double*) malloc(N * sizeof(double)); + } + + /* Fill matrices A and B with values */ + for (i = 0; i < N; i++) { + for (j = 0; j < N; j++) { + A[i][j] = 2.0; + /* We want a random value between 0 and 1 */ + B[i][j] = rand()/RAND_MAX; + } + } + + /* Perform matrix multiplication */ + matrixMult(A, B, C, N); + + /* Free memory */ + for (i = 0; i < N; i++) { + free(A[i]); + free(B[i]); + free(C[i]); + } + free(A); + free(B); + free(C); + + /* Exit with success */ + return 0; +} diff --git a/matrixmul/c-linarr.c b/matrixmul/c-linarr.c @@ -0,0 +1,58 @@ +#include <stdlib.h> +#include <stdio.h> + +void matrixMult(double* A, double* B, double* C, int N) +{ + int i, j, k; + double sum; +#pragma omp parallel for private (j,k,sum) + for (i = 0; i<N; i++) { + for (j = 0; j<N; j++) { + sum = 0.0f; + for (k = 0; k<N; k++) { + sum += A[k*N+i] * B[j*N+k]; + } + C[j*N+i] = sum; + } + } +} + + +int main(int argc, char* argv[]) +{ + int i, N; + double* A; + double* B; + double* C; + + /* Read input argument as matrix height */ + if (argc == 2) { + N = atoi(argv[1]); + } else { + printf("Sorry, I need matrix width as command line argument\n"); + return 1; + } + + /* Allocate arrays */ + A = (double*) malloc(N * N * sizeof(double*)); + B = (double*) malloc(N * N * sizeof(double*)); + C = (double*) malloc(N * N * sizeof(double*)); + + /* Fill matrices A and B with values */ + for (i = 0; i < N*N; i++) { + A[i] = 2.0; + /* We want a random value between 0 and 1 */ + B[i] = rand()/RAND_MAX; + } + + /* Perform matrix multiplication */ + matrixMult(A, B, C, N); + + /* Free memory */ + free(A); + free(B); + free(C); + + /* Exit with success */ + return 0; +} diff --git a/matrixmul/cpp-linvectors.cpp b/matrixmul/cpp-linvectors.cpp @@ -0,0 +1,49 @@ +#include <iostream> +#include <vector> +#include <cstdlib> // rand is here + +int main(int argc, char* argv[]) +{ + using std::cout; + using std::vector; + + int M, N, i, j, k; + + // Read input argument as matrix height + if (argc == 2) { + N = atoi(argv[1]); + } else { + cout << "Sorry, I need matrix width as command line argument\n"; + return 1; + } + + // Width equal to height + M = N; + + vector<double> A(M*N); + vector<double> B(M*N); + vector<double> C(M*N); + + // Fill matrices A and B with values + for (i = 0; i<N; ++i) { + for (j = 0; j<M; ++j) { + A[j*M+i] = 2.0; + // We want a random value between 0 and 1 + B[j*M+i] = rand()/RAND_MAX; + } + } + + double sum; + // Perform matrix multiplication + for (i = 0; i < N; ++i) { + for (j = 0; j < M; ++j) { + sum = 0.0f; + for (k = 0; k < M; ++k) + sum += A[k*M+i] * B[j*M+k]; + C[j*M+i] = sum; + } + } + + // Exit with success + return 0; +} diff --git a/matrixmul/cpp-vectorofvectors.cpp b/matrixmul/cpp-vectorofvectors.cpp @@ -0,0 +1,49 @@ +#include <iostream> +#include <vector> +#include <cstdlib> // rand is here + +int main(int argc, char* argv[]) +{ + using std::cout; + using std::vector; + + int M, N, i, j, k; + + // Read input argument as matrix height + if (argc == 2) { + N = atoi(argv[1]); + } else { + cout << "Sorry, I need matrix width as command line argument\n"; + return 1; + } + + // Width equal to height + M = N; + + vector< vector<double> > A(M,vector<double>(M)); + vector< vector<double> > B(M,vector<double>(M)); + vector< vector<double> > C(M,vector<double>(M)); + + // Fill matrices A and B with values + for (i = 0; i<N; ++i) { + for (j = 0; j<M; ++j) { + A[i][j] = 2.0; + // We want a random value between 0 and 1 + B[i][j] = rand()/RAND_MAX; + } + } + + double sum; + // Perform matrix multiplication + for (i = 0; i < N; ++i) { + for (j = 0; j < M; ++j) { + sum = 0.0; + for (k = 0; k < M; ++k) + sum = A[k][j] * B[i][k]; + C[i][j] = sum; + } + } + + // Exit with success + return 0; +} diff --git a/matrixmul/julia.jl b/matrixmul/julia.jl @@ -0,0 +1,13 @@ +#!/usr/bin/julia + +if (length(ARGS) == 1) + N = int(ARGS[1]) +else + println("Sorry, I need the matrix width as a command line argument\n") +end + +M = N + +A = ones(M,N)*2.0 +B = rand(M,N) +C = A*B diff --git a/matrixmul/lua-arrofarrs.lua b/matrixmul/lua-arrofarrs.lua @@ -0,0 +1,30 @@ +#!/usr/bin/lua + +N = arg[1] + +-- Initialize the matrix A with two's, +-- matrix B with random numbers +A = {} +B = {} +C = {} +for i=1,N do + A[i] = {} + B[i] = {} + C[i] = {} + for j=1,N do + A[i][j] = 2 + B[i][j] = math.random() + end +end + +-- Multiply matrix A with matrix B, +-- store result in matrix C +for i=1,N do + for j=1,N do + sum = 0.0 + for k=1,N do + sum = sum + A[i][k] * B[k][j] + end + C[i][j] = sum + end +end diff --git a/matrixmul/lua-linarr.lua b/matrixmul/lua-linarr.lua @@ -0,0 +1,25 @@ +#!/usr/bin/lua + +N = arg[1] + +-- Initialize the matrix A with two's, +-- matrix B with random numbers +A = {} +B = {} +for i=1,(N*N) do + A[i] = 2.0 + B[i] = math.random() +end + +-- Multiply matrix A with matrix B, +-- store result in matrix C +C = {} +for i=1,N do + for j=1,N do + sum = 0.0 + for k=1,N do + sum = sum + A[(k-1)*N+i] * B[(j-1)*N+k] + end + C[(i-1)*N + j] = sum + end +end diff --git a/matrixmul/octave.m b/matrixmul/octave.m @@ -0,0 +1,14 @@ +#/usr/bin/octave -q -f + +if (nargin == 1) + M = str2num(argv(){1}); +else + disp("Sorry, I need the matrix size as a command line argument"); + exit(1); +end + +N = M; + +A = ones(M,N) * 2.0; +B = rand(M,N); +C = A*B; diff --git a/matrixmul/plot.gp b/matrixmul/plot.gp @@ -0,0 +1,9 @@ +set terminal png size 1200,600 +set output "performance.png" +set xlabel "Matrix side length" +set ylabel "Execution time [s]" +set title "Random number generation and matrix multiplication" +set log xy +set grid +set key outside +plot "lua-arrofarrs.dat" title "Lua: Arrays of arrays" w lp, "lua-linarr.dat" title "Lua: Linear arrays" w lp, "luajit-arrofarrs.dat" title "LuaJIT: Arrays of arrays" w lp, "luajit-linarr.dat" title "LuaJIT: Linear arrays" w lp, "c-arrofarrs.dat" title "C: Arrays of arrays" w lp, "c-linarr.dat" title "C: Linear arrays" w lp, "c-omp-arrofarrs.dat" title "C-OpenMP: Arrays of arrays" w lp, "c-omp-linarr.dat" title "C-OpenMP: Linear arrays" w lp, "julia.dat" title "Julia" w lp, "cpp-vectorofvectors.dat" title "C++: Vectors of vectors" w lp, "cpp-linvectors.dat" title "C++: Linear vectors" w lp, "python-numpy.dat" title "Python: Numpy" w lp, "octave.dat" title "Octave" w lp diff --git a/matrixmul/python-numpy.py b/matrixmul/python-numpy.py @@ -0,0 +1,16 @@ +#!/usr/bin/env python + +import sys +import numpy + +if (len(sys.argv) == 2): + M = int(sys.argv[1]) +else: + print("Sorry, I need a matrix width as input argument!") + sys.exit(1) + +N = M + +A = numpy.ones((M,N))*2.0 +B = numpy.random.random_sample((M,N)) +C = numpy.dot(A,B) # A*B is element wise on numpy arrays diff --git a/montecarlo/Makefile b/montecarlo/Makefile @@ -0,0 +1,89 @@ +# Define compiler +CC=g++ + +# Define compiler flags (show all warnings) +CPPFLAGS=-Wall +#CPPFLAGS=-std=c++0x + +# Define linker flags +LDFLAGS= + +# Define extra libraries to be dynamically linked +#LDLIBS+=-larmadillo +#LDLIBS+=-lstdc++ + +# Compile optimized code +#CPPFLAGS+=-O2 + +# Compile debuggable code +#CPPFLAGS+=-g + +# Compile profilable code +#CPPFLAGS+=-pg +#LDFLAGS+=-pg + +# Define linker +LD=g++ + +# Filenames of source code +SRC=$(shell ls *.cpp) + +# Filenames of object files +OBJ=$(SRC:.cpp=.o) +#OBJ=downhill_simplex.o main.o +#OBJB=downhill_simplex.o main.B.o + +# Remove file type extension for binary filename +BIN=mc +#BINB=optimB + +# The default "all" depends on A and B + +#all: A B + +A: performance.png error.png + +%.png: performance.dat + gnuplot plot.gp + + + +performance.dat: $(BIN) + # Removing old performance data + rm -f performance.dat + # Running program and profiling runtime and error + nice ./$(BIN) 10000 > /dev/null + nice ./$(BIN) 20000 > /dev/null + nice ./$(BIN) 50000 > /dev/null + nice ./$(BIN) 100000 > /dev/null + nice ./$(BIN) 200000 > /dev/null + nice ./$(BIN) 500000 > /dev/null + nice ./$(BIN) 1000000 > /dev/null + nice ./$(BIN) 2000000 > /dev/null + nice ./$(BIN) 5000000 > /dev/null + + +#B: $(BINB) +# ./$(BINB) 2> amoebaB.dat + +test: $(BIN) + ./$(BIN) + +$(BIN): $(OBJ) + @# Link object files together + $(LD) $(LDFLAGS) $(OBJ) -o $(BIN) $(LDLIBS) + +$(BINB): $(OBJB) + @# Link object files together + $(LD) $(LDFLAGS) $(OBJB) -o $(BINB) $(LDLIBS) + +clean: + @# Remove object files + rm -f $(OBJ) + @# Remove binaries + rm -f $(BIN) + @# Remove datafiles and plot + rm -f *.dat *.png +edit: + vim -p Makefile *.cpp *.h *.gp + diff --git a/montecarlo/functions.h b/montecarlo/functions.h @@ -0,0 +1,20 @@ +#ifndef FUNCTIONS_H_ +#define FUNCTIONS_H_ + +#include <iostream> +#include <vector> +#include <cmath> +#include "header.h" + +int ncalls = 0; + +Floattype functionA(const std::vector<Floattype> x) { + if (x.size() != 3) { + std::cout << "Error! FunctionA must be given a 3D input point!\n"; + return NAN; + } + + return 1.0f / (1.0f-cos(x[0]) * cos(x[1]) * cos(x[2])) / M_PI / M_PI / M_PI; +} + +#endif diff --git a/montecarlo/header.h b/montecarlo/header.h @@ -0,0 +1,18 @@ +// Make sure header is only included once +#ifndef HEADER_H_ +#define HEADER_H_ + +// Define whether the program should output values of matrices +const bool verbose = false; +//const bool verbose = true; + +// Prototype for checking function +void check(const bool statement); + +// Choose length variable type +typedef long int Lengthtype; + +// Choose floating-point precision +typedef double Floattype; + +#endif diff --git a/montecarlo/main.cpp b/montecarlo/main.cpp @@ -0,0 +1,62 @@ +#include <iostream> +#include <fstream> +#include <vector> +#include <cmath> +#include <cstdlib> +#include <ctime> +#include "header.h" +#include "functions.h" +#include "montecarlo.h" + +int main(int argc, char* argv[]) +{ + // Namespace declarations + using std::cout; + + // Number of sample points as input argument + Lengthtype N; + if (argc == 1) // If no args are given.. + N = 100; // 100 points are sampled + else + N = atol(argv[1]); // Else the specified number + + cout << "Sampling function at " << N << " points.\n"; + + // Calculate machine precision + Floattype eps_machine = 1.0f; + while (1.0f + eps_machine != 1.0f) + eps_machine /= 2.0f; + + // Evaluate 3D function A at 10 points in interval [a;b] + Lengthtype d = 3; // No of dimensions + std::vector<Floattype> a(d); // Lower limits + std::vector<Floattype> b(d); // Upper limits + for (Lengthtype i=0; i<d; ++i) { // Assign equidimensional limits + a[i] = 0.0f; + b[i] = M_PI; + } + Floattype tic = clock(); + MC mc(functionA, a, b, N); + mc.plain(); // Plain Monte-Carlo integration + mc.show(); // Display results + Floattype t_elapsed = (clock() - tic)/(CLOCKS_PER_SEC); + cout << "Elapsed time: " << t_elapsed << " s.\n"; + + // Append results to performance.dat + std::ofstream file; + file.open("performance.dat", std::fstream::app); // Append to file + file << N << '\t' << t_elapsed << '\t' << mc.error() << '\n'; + file.close(); + + // Return successfully + return 0; +} + +void check(const bool statement) +{ + using std::cout; + if (statement == true) + cout << "\t\033[0;32mPassed\033[0m\n"; + else + cout << "\t\033[1;31mFail!!\033[0m\n"; +} diff --git a/montecarlo/montecarlo.cpp b/montecarlo/montecarlo.cpp @@ -0,0 +1,88 @@ +#include <iostream> +#include <vector> +#include <cstdlib> // for random functionality +#include <cmath> // NaN +#include "header.h" +#include "montecarlo.h" + +// Constructor +MC::MC(Floattype (*function)(const std::vector<Floattype>), + const std::vector<Floattype> a_in, + const std::vector<Floattype> b_in, + const Lengthtype N_in) + : d(a_in.size()), f(function), a(a_in), b(b_in), x(d), N(N_in) +{ + // Check that problem is at least 1D + if (d < 1) + std::cerr << "Error! Problem must be at least 1D\n"; + + // Check that the user has specified at least 1 random sample + if (N < 1) + std::cerr << "Error! The algorithm should evaluate at least 1 point!\n"; + + // Check input data + for (Lengthtype i=0; i<d; ++i) + if (a[i] >= b[i]) + std::cerr << "Error! a >= b in dimension " << i << '\n'; + + // Initialize result and error to NaN + Q = NAN; + err = NAN; +} + +// Plain Monte Carlo integrator +void MC::plain() +{ + + // Set volume of sample space + set_volume(); + + Floattype sum = 0.0f, sumsum = 0.0f, fx; + + for (Lengthtype i=0; i<N; ++i) { + x = random_x(); // Random sample point inside intervals + fx = f(x); + sum += fx; + sumsum += fx*fx; + } + + Floattype average = sum/N; + Floattype variance = sumsum/N - average*average; + + // Write results + Q = average * V; + err = sqrt(variance/N)*V; +} + + +// Calculate volume by multiplying interval in each dimension +void MC::set_volume() +{ + V = 1.0f; + for (Lengthtype i=0; i<d; ++i) + V *= b[i] - a[i]; +} + +// Draw pseudo-random position in sample space +std::vector<Floattype> MC::random_x() +{ + std::vector<Floattype> pos(d); + for (Lengthtype i=0; i<d; ++i) { + // Random number in [a;b] interval in dimension + pos[i] = (b[i] - a[i]) * ((Floattype)rand()/RAND_MAX) + a[i]; + } + return pos; +} + + +// Print results +void MC::show() +{ + std::cout << "Integral Q = " << Q << ", error = " << err << '\n'; +} + +// Return the error +Floattype MC::error() +{ + return err; +} diff --git a/montecarlo/montecarlo.h b/montecarlo/montecarlo.h @@ -0,0 +1,58 @@ +#ifndef MONTECARLO_H_ +#define MONTECARLO_H_ + +#include <vector> +#include "header.h" + +class MC { + + private: + + // Length of vectors + const Lengthtype d; + + // Function to be evaluated (pointer to function) + Floattype (*f)(const std::vector<Floattype>); + + // Integral limits in dimension d + const std::vector<Floattype> a; + const std::vector<Floattype> b; + + // n-dimensional point + std::vector<Floattype> x; + + // Number of samples + const Lengthtype N; + + // Integration result + Floattype Q; + + // Error + Floattype err; + + // Volume + Floattype V; + void set_volume(); + + // Draw random position in sample space + std::vector<Floattype> random_x(); + + public: + + // Constructor + MC(Floattype (*function)(const std::vector<Floattype>), + const std::vector<Floattype> a_in, + const std::vector<Floattype> b_in, + const Lengthtype N_in); + + // Plain Monte Carlo integrator + void plain(); + + // Print result and error + void show(); + + // Return the error + Floattype error(); +}; + +#endif diff --git a/montecarlo/plot.gp b/montecarlo/plot.gp @@ -0,0 +1,17 @@ +set terminal png +set output "performance.png" +set xlabel "Sample points N" +set ylabel "Execution time [s]" +set title "Performane analysis" +set log xy +set grid +plot "performance.dat" u 1:2 title "Plain Monte-Carlo" w lp + +set terminal png +set output "error.png" +set xlabel "Sample points N" +set ylabel "Error" +set title "Error analysis" +set log xy +set grid +plot "performance.dat" u 1:3 title "Plain Monte-Carlo" w lp, 1/sqrt(x) w l title "O(1/sqrt(N))" diff --git a/notes.rst b/notes.rst @@ -0,0 +1,398 @@ +# vim: set tw=80:noautoindent + +For Jacobi exercise: Use atan2(). + +11-4: Introduction +##### +Deadline for exercises: ~2 weeks after assignment + +All exercises must be done on Lifa: + lifa.phys.au.dk + Ethernet: vlifa01 or vlifa02 +You need RSA keys to login from outside net. +Login from inside: NFIT credentials. + +Alternative server: molveno +Dedicated for numerical course. + +Lifa will have approx. 5 year old software, molveno is updated. + +All phys.au.dk servers have NFS shared home folders. + +Dmitri answers: +http://www.phys.au.dk/~fedorov/numeric/ + +Bash: Last command starting with e.g. 'v': !v + +Exercises are weighted 75% and the exam 25%. You do need to complete at least +51% of the exam. + +02: >51% +4: >60% +7: >80% +10: >90% +12: 100% + +16-4: Interpolation +#### +Makefiles_: +For all project, use makefiles, the project is evaluated by `make clean; make`. + +General structure of Makefile: + Definitions of variables + all: my_target + Target(s): Dependenc(y,ies) + <-tab-> Command(s) + +Strings are denoted without decoration, variables as $(VAR). +To use the CC system linker to link C++ objects, add libraries: + LDLIBS += -lstdc++ +If you use e.g. g++ as the linker, the above command is not needed. + +If an object file is required as a dependency, the Makefile will issue a CC +command, compiling the source code file with the same basename. `make -p | less` +will show all build rules, even build in. + +Interpolation_: +When you have a discrete set of datapoints, and you want a function that +describes the points. + +If a function is analytical and continuous, and an infinitely large table of +datapoints in an interval, the complete analytical function can be found from a +taylor-series of _all_ derivatives in the interval. + +k-Polynomial: k+1 unknowns (variables). High polynomials tend to oscillate. + +Cubic interpolation: The interval between points can be described by a function +of three polynomials. + +Spline interpolation: Also made of polynomials. First order spline (a0+x*a1), +simply connects the data points linearly. The splines for each inter-datapoint +interval must have the same values at the data points. The derivates are +discontinuous. + +Solving linear splines: +Datapoints: {x_i,y_i}, i = 1,...,n +Splines (n-1): S_i(x) = a_i + b_i (x - x_i) + S_i(x_i) = y_i (n equations) + S_i(x_{i+1}) = S_{i+1}(x_i) (n-2 equations (inner points)) +Unkowns: a,b: 2n-2 variables. Solution: + => a_i = y_i + S_i(x) = y_i + b_i (x-x_i) + S_i(x_{i+1}) = y_i + b_i (x_{i+1} - x_i) + => b_i = (y_{i+1} - y_i) / (x_{i+1} - x_i) + +Start of by searching which interval x is located in. Optionally, remember the +interval, as the next value will probably lie in the same interval. + +Binary search: Is x between x_1 and x_n? If not, error: Extrapolation. +Continuously divide the interval into two, and find out if the x is in the +larger or smaller part. +DO A BINARY SEARCH. + +Second order spline: + S_i(x) = a_i + b_i(x-x_i) + c_i(x-x_i)^2 + Unknowns: 3(n-1) = 3*n-3 + Equations: 3*n-4 + The derivatives are also continuous. + +Solution: + a_i = y_i + \delta x = (x-x_i) + y_i + b_i \delta x_i + c_i \delta x_i^2 = y_{i+1} + b_i + 2 c_i \delta x_i = b_{i+1} +Suppose you know b_i, you can find c_i. From that you can find b_{i+1}, and in +turn, c_{i+1}. Through recursion you find all unknowns by stepping forward. +The backwards solution can be found from the last data-point (y_n) by solving +the two equations with c_{n-1} and b_{n-1} as the two unkowns. + +Symmetry can be used as an extra condition. + + +18-4: +############## +Interpolation exercise, plotting optionally in gnuplot, or graph (from +plotutils): + graph -T png points.dat > plot.png +In makefile, example: + plot.png: points.dat + graph --display-type png $^ > $@ +Each dataset in points.dat needs a header, e.g. # m=1, S=0 + +lspline.cc: Linear spline +qspline.cc: Quadratic spline +cspline.cc: Cubic spline + +Linear spline: + S(x) = S_i(x) it x in [x_i,x_(i+1)] + S_i(x) = a_i + b_i (x-x_i) + S_i(x) = y_i + (\delta y_i)/(\delta x_i) (x-x_i) + S_i(x_i) = y_i + \delta y_i = y_(i+1) - y_i + S_i (x_(i+1)) = y_(i+1) + +In C++: + std::vector<double> x,y,p; +Maybe typecasting? Could be fun. + +Procedural programming: +-- + struct lspline { + int n; + vector<double> x,y,p; + }; + +Make a function: + struct lspline* construct_lspline(vector<double>&x, vector<double>&y) + double evaluate_lspline(struct lspline * asdf, double z) + free_lspline(); + +Object-oriented programming: +-- +If you want to take the structural approach, you keep the functions and +structures seperate. If you take a OO approach, put the functions inside the +structure (or class): + struct lspline { + int n; + vector<double> x,y,p; + lspline(...,..); // Constructor, same name as struct + double eval(double z); + }; + struct lspline ms(x,y); + ms.eval(z); + +See Dmitri's cubic spline example which uses OO. + +Functional programming (in C++11), compile with -std=c++0x: +-- +The functions can return functions: + + #include <functional> + using namespace std; + + function<double(double)> flspline (vector<double>&x, vector<double>&y); + + auto my_spline = flspline(x,y); + my_spline(5,0); + + +System of linear equations: +------- +A*x = b + +A_i,j x_j = b_i + +Solve by finding A^(-1): x = A^(-1) * b +Numerically, you calculate the inverse by solving Ax=b. +We will assume that the matrixes are not singular. + + Turn the system into a triangular form. + The main diagonal is non-zero, all lower values are 0, and upper values are + denoted T_nn. + T_nn * x_n = b_n => x_n = 1/T_nn * b_nn + T_(n-1,n-1) x_(n-1) + T_(n-1,n) x_n = b_(n-1) + T_ii x_i + sum^n_(j=i+1) T_(i,j) x_j = b_i + x_i = 1/T_ii (b_i - sum^n_(j=i+1) T_ij, x_j) + +The simplest triangulation is by Gauss elimination. Numerically, the simplest +method is LU decomposition (Lower Upper). + A = LU, where L=lower triangle, U=upper triangle. + n^2 equations. + L and U contain "(n^2 - n)/2 + n" elements. + L+U = n^2/2 + n/2 = (n(n+1))/2 + + The diagonals in L are all equal to 1: L_ii = 1. + See Dolittle algorithm in the lecture notes, which with the LU system is the + most used, and fastest method for solving a linear equation. + + Ax = b + LUx = b, Ux=y + +Another method: QR decomposition: R=Right triangle (equal to U). + A = QR + Q^T Q = 1 (orthogonal) + + Ax = b + QRx = b + Rx = Q^T b + + Gram-Schmidt (spelling?) orthogonalization: + Consider the columns of your matrix A. Normalize them. Orthogonalize all other + columns to the first column. + + Normalizing the column: ||a_1|| = sqrt(a_1 dot a_i) + Orthoginalize columns: a_2/||a_1|| -> q_1 + + Numerically: + for j=2 to m: + a_j - dot(dot(q_1,a_j),q_1) -> a_j + a_2 -> a_2/||a_2|| = q_2 + +Find inverse matrix: + A A^-1 = diag(1) + + + +30-4: Diagonalization +######################### + +Runtime comparison: Do a number of comparisons with different matrix sizes +etc.numeric-2012 + +Easiest way to diagonalize a matrix: Orthogonal transformation + A -> Q^T A Q = D +Q matrix can be built with e.g. QR decomposition. Rotation: Computers to cyclic +sweep, which is faster than the classic rotation. +Cyclic: Zero all elements above the diagonal, and do your rotations until the +matrix is diagonal. The matrix is converged if none of the eigenvalues have +changed more than machine precision. You will destroy the upper half plus the +diagonal. If you store the diagonal in another vector, you can preserve the +matrix values. +In the beginning, you look for the largest element in each row, and create an +index which is a column that keeps the numbers of the largest elements. +With an index, you can perform the operation fast. You have to update the index +after each operation. + +For very large matrixes, > system memory: +Krylov subspace methods: You use only a subspace of the whole space, and +diagonalize the matrix in the subspace. + + +02-5: Ordinary Differential Equations (ODEs) +############################################ + + dy/dx = f(x,y) + +Usually coupled equations: + dy_1/dx = f_1(x,y_1,y_2) + dy_2/dx = f_1(x,y_1,y_2) + [y_1, ..., y_n] = y + [f_1(x,y), ..., f_n(x,y) = f(x,y) + y' = f(x,y) +x is usually time. If f has the above form, it is a ODE. +Sine function: + y'' = -y +In this form, it is not a ODE, because it has a second derivative. +You can redifine it to be an ODE: + => y = u_1 -> u'_1 = u_2 + y' = u_2 -> u'_2 = -u_1 + +Solving ODEs with a starting condition: + y' = f(x,y) + y(x_0) = y_0 + a->b (Shape of the function in the interval) +You can calculate the derivative at (x_0,y_0). You then make a step towards +x_0+h, and apply Euler's method: + y' \approx (\delta y)/(\delta x) = (y(x_0+h) - y_0)/h + y(x_0+h) \approx y(x_0) + h*f(x_0,y_0) + y(x) = y_0 + y'_0(x-x0) + 1/(2!) y''_0(x_0)(x-x_0)^2 + ... (Taylor exp.) +You can find the higher-order terms by sampling points around (x_0,y_0). +Estimate the new derivative half a step towards h, which is used for the first +point. You sample your function, and fit a polynomial. Higher order polynomiums +tend to oscillate. +When solving ODE's the many sampled points create a tabulated function. If you +want to do something further with the function, you interpolate the points by +e.g. splines. +Only three steps are needed if the equation is a parabola. Other functions are +called multi-step methods. You do not want to go to high in the +Taylor-expansion, as there lies a danger with the higher order terms +oscillating. +If the function is changing fast inside an interval, the step size h needs to be +small. This is done by a driver. +The Runge-kutta is single step, and the most widely used. + +Absolute accuracy (delta) vs. relative accuracy (epsilon) can behave +significantly differently when the machine precision is reached. +The user usually specifies both, and the accuracy is satisfied, when one of the +conditions are met (the tolerance, tau). + tau = delta + epsilon*||y|| + +Stepper: Must return y(x+h), e (error estimate). +Driver: x->x+h. If e is smaller than tau, it accepts the step. +The driver finds the size of the next h: + h_next = h(tau/e)^power * safety + power = 0.25 + safety = 0.95 + The driver thus increases the step size if the error was low relative to the + tolerance, and vice-versa. + + Runge principle for determining the error: + e = C*h^(p+1) + You first do a full h step, then two h/2 steps. + 2*c*(h/2)^(p+1) = c * (h^(p+1))/(2^p) + For the error, you calculate y_1 from full step, and then y_1 from two half steps: + e = ||y_1(full step) - y_1(2 half steps)|| + You can also instead use RK1 and RK2, and evaluate the difference between the + two for the same step. + +The effectiveness of the driver and stepper is determined by how many times you +call the right-hand side of the ODEs. + +Save the x and y values in a C++ vector, and add dynamically using the +push_back() function. + + +07-5: Nonlinear equations and optimization +##### +Pipe stdin to program: + cat input.txt | ./my_prog +Redirect stdout to file: + ./my_prog > out.txt +Redirect stderr to file: + ./my_prog 2> err.txt + +Jacobian matrix: Filled with partial derivatives. Used for linearizing the +problem. f(x+\delta x) \approx f + J \delta x + +Machine precision: double: 2e-16. It is the largest number where 1 + eps = 1. + +Quasi-Newton methods: Might be exam exercise. + + +14-5: Numerical integration +##### +Examination problem: FFT or Multiprocessing + +Functions for calculating: + \int_a^b f(x) dx + +Often not possible analytically, thus numerical aproach. Most differential +operations are possible analytically. + +Riemann integral: Numerical approximation. +You divide the interval into subintervals (t_1, ..., t_n). +Riemann sum: + S_n = \sum_{i=1}^n f(x_i) \Delta x_i +Approaching the integral as the interval approaches 0.0. The geometric +interpretation corresponds to rectangles with height f(x_i) and width \Delta +x_i. The function passes trough the middle of the top side of the rectangles. +The integral exists also for discontinuous functions. + +Rectangle method: Stable, does no assumptions + \int_a^b f(x) dx \approx \sum_{i=1}^n f(x_i) \Delta x_i +Trapezum method: + Instead of one point inside the interval, two points are calculated, and the + average is used. + \int_a^b f(x) dx \approx \sum_{i=1}^n (f(x_{i+1} )+ f(x_i))/2 * \Delta x_i +Adaptive methods: + The \Delta x_i is adjusted to how flat the function is in an interval. + As few points as possible are used. +Polynomial functions: + Analytical solutions exist. Taylor series expansion. w_i: Weight. + \sum_{i=1}^n w_i f(x_i) => \int_a^b f(x) dx + `---- quadrature -----` + n functions: {\phi_1(x),...,\phi_n(x)} + Often the polynomials are chosen: {1,x,x^2,...,x^{n-1}} + \sum_{i=1}^n w_i \phi_k(x_i) = I_k + I_k = \int_a^b \phi_k (x) dx +Gauss methods: + Also use the quadratures as tuning parameters, as the polynomial method. + w_i,x_i as tuning parameters: 2*n tuning parameters. + \int_a^b f(x) w(x) dx + Functions like 1/sqrt(x) or sqrt(x) are handled through the weight function. + See the lecture notes for details. Both nodes and weights are adjusted. + +The method to use is dependent on the problem at hand. + + + diff --git a/notes/book.pdf b/notes/book.pdf Binary files differ. diff --git a/notes/eigen.pdf b/notes/eigen.pdf Binary files differ. diff --git a/notes/fft.pdf b/notes/fft.pdf Binary files differ. diff --git a/notes/integration.pdf b/notes/integration.pdf Binary files differ. diff --git a/notes/interp.pdf b/notes/interp.pdf Binary files differ. diff --git a/notes/krylov.pdf b/notes/krylov.pdf Binary files differ. diff --git a/notes/leastsq.pdf b/notes/leastsq.pdf Binary files differ. diff --git a/notes/lineq.pdf b/notes/lineq.pdf Binary files differ. diff --git a/notes/montecarlo.pdf b/notes/montecarlo.pdf Binary files differ. diff --git a/notes/odes.pdf b/notes/odes.pdf Binary files differ. diff --git a/notes/roots.pdf b/notes/roots.pdf Binary files differ. diff --git a/notes/sfun.pdf b/notes/sfun.pdf diff --git a/optimization/Makefile b/optimization/Makefile @@ -0,0 +1,69 @@ +# Define compiler +CC=g++ + +# Define compiler flags (show all warnings) +CPPFLAGS=-Wall -std=c++0x + +# Define linker flags +LDFLAGS= + +# Define extra libraries to be dynamically linked +LDLIBS+=-larmadillo -lstdc++ + +# Compile optimized code +#CPPFLAGS+=-O2 + +# Compile debuggable code +#CPPFLAGS+=-g + +# Compile profilable code +#CPPFLAGS+=-pg +#LDFLAGS+=-pg + +# Define linker +LD=g++ + +# Filenames of source code +SRC=$(shell ls *.cpp) + +# Filenames of object files +#OBJ=$(SRC:.cpp=.o) +OBJ=downhill_simplex.o main.o +OBJB=downhill_simplex.o main.B.o + +# Remove file type extension for binary filename +BIN=optim +BINB=optimB + +# The default "all" depends on A and B + +all: A B + + +A: $(BIN) + @# Run program and write amoeba positions to file + ./$(BIN) 2> amoeba.dat + +B: $(BINB) + ./$(BINB) 2> amoebaB.dat + +$(BIN): $(OBJ) + @# Link object files together + $(LD) $(LDFLAGS) $(OBJ) -o $(BIN) $(LDLIBS) + +$(BINB): $(OBJB) + @# Link object files together + $(LD) $(LDFLAGS) $(OBJB) -o $(BINB) $(LDLIBS) + +clean: + @# Remove object files + rm -f $(OBJ) $(OBJB) + @# Remove binaries + rm -f $(BIN) $(BINB) + @# Remove datafiles and plot + #rm -f *.dat *.png + @# Remove error log + rm -f amoeba.dat amoebaB.dat +edit: + vim -p Makefile *.cpp *.h + diff --git a/optimization/downhill_simplex.cpp b/optimization/downhill_simplex.cpp @@ -0,0 +1,173 @@ +#include <iostream> +#include <armadillo> +#include <vector> +#include <functional> +#include "downhill_simplex.h" + +using namespace std; +using namespace arma; + +// Amoeba constructor +amoeba::amoeba(function<double(vec)> fun, vector<vec> simplex) + : d(simplex.size()-1), f(fun), value(zeros<vec>(d)), p_ce(zeros<vec>(d)) +{ + p_ce_o = p_ce; + + for (int i=0; i<d+1; ++i) + p.push_back(simplex[i]); + for (int i=0; i<d+1; ++i) + value[i] = f(p[i]); + + update(); + +} + +// Update amoeba parameters +void amoeba::update() +{ + + p_ce_o = p_ce; + + // Find highest point + hi = 0; + for (int i=1; i<d+1; ++i) + if (value[i] > value[hi]) + hi = i; + + // Find lowest point + lo = 0; + for (int i=1; i<d+1; ++i) + if (value[i] < value[lo]) + lo = i; + + // Find centroid + p_ce = zeros<vec>(d); + for (int i=0; i<d+1; ++i) + if (i != hi) + p_ce += p[i]; + p_ce /= d; + + // Write amoeba position to stderr + pos(); +} + +// Find size of simplex +double amoeba::size() +{ + double s = 0; + for (int i=0; i<d+1; ++i) + if (i != lo) { + double n = norm(p[i] - p[lo], 2.0f); + if (n>s) + s=n; + } + return s; +} + +void amoeba::downhill(double simplex_size_goal) +{ + // Find operation to update position + while (size() > simplex_size_goal) { + + vec p_re = p_ce + (p_ce - p[hi]); // Try reflection + double f_re = f(p_re); + if (f_re < value[lo]) { + vec p_ex = p_ce + 1.0f * (p_ce - p[hi]); // Try expansion + double f_ex = f(p_ex); + if (f_ex < f_re) { + value[hi] = f_ex; + p[hi] = p_ex; // Accept expansion + update(); continue; // Start loop over + } + } + + if (f_re < value[hi]) { + value[hi] = f_re; + p[hi] = p_re; // Accept reflection + update(); continue; // Start loop over + } + + vec p_co = p_ce + 0.5f * (p[hi] - p_ce); // Try contraction + double f_co = f(p_co); + if (f_co < value[hi]) { + value[hi] = f_co; + p[hi] = p_co; // Accept contraction + update(); continue; // Start loop over + } + + for (int i=0; i<d+1; ++i) + if (i != lo) { + p[i] = 0.5f * (p[i] + p[lo]); // Do reduction + value[i] = f(p[i]); + } + update(); continue; + } +} + +void amoeba::downhill_mod(double simplex_size_goal) +{ + // Find operation to update position + while (size() > simplex_size_goal) { + + vec p_re = p_ce + (p_ce - p[hi]); // Try reflection + + double f_re = f(p_re); + if (f_re < value[lo]) { + vec p_ex_n = p_ce + 1.0f * (p_ce - p[hi]); // Try expansion + double f_ex_n = f(p_ex_n); + vec p_ex_o = p_ce_o + 1.0f * (p_ce_o - p[hi]); // Try expansion, old val + double f_ex_o = f(p_ex_o); + + // Find out which expansion gave the lowest value + vec p_ex; double f_ex; + if (f_ex_n < f_ex_o) { // New val + p_ex = p_ex_n; + f_ex = f_ex_n; + } else { // Old val + p_ex = p_ex_o; + f_ex = f_ex_o; + } + + if (f_ex < f_re) { + value[hi] = f_ex; + p[hi] = p_ex; // Accept expansion + update(); continue; // Start loop over + } + } + + if (f_re < value[hi]) { + value[hi] = f_re; + p[hi] = p_re; // Accept reflection + update(); continue; // Start loop over + } + + vec p_co = p_ce + 0.5f * (p[hi] - p_ce); // Try contraction + double f_co = f(p_co); + if (f_co < value[hi]) { + value[hi] = f_co; + p[hi] = p_co; // Accept contraction + update(); continue; // Start loop over + } + + for (int i=0; i<d+1; ++i) + if (i != lo) { + p[i] = 0.5f * (p[i] + p[lo]); // Do reduction + value[i] = f(p[i]); + } + update(); continue; + } +} + +// Write simplex points to stderr +void amoeba::pos() +{ + for (int i=0; i<d+1; ++i) + std::cerr << p[i][0] << '\t' + << p[i][1] << '\n'; +} + +// Return lowest point of simplex +vec amoeba::low() +{ + return p[lo]; +} diff --git a/optimization/downhill_simplex.h b/optimization/downhill_simplex.h @@ -0,0 +1,61 @@ +#ifndef DOWNHILL_SIMPLEX_H_ +#define DOWNHILL_SIMPLEX_H_ + +#include <functional> +#include <vector> +#include <armadillo> +using namespace std; +using namespace arma; + +class amoeba { + + public: + + // Constructor + amoeba(function<double(vec)> fun, vector<vec> simplex); + + // Return lowest simplex position + vec low(); + + // Move amoeba downhill with relevant method + void downhill(double simplex_size_goal); + + // Move amoeba downhill with relevant method, modified method + void downhill_mod(double simplex_size_goal); + + private: + + // Number of vertexes + int d; + + // Maxima points of function + int hi, lo; + + // Vertex points + vector<vec> p; + + // Functions to be evaluated + std::function<double(vec)> f; + + // Function values at points + vec value; + + // Centroid + vec p_ce; + + // Centroid from previous step + vec p_ce_o; + + // Filename of output file + + // Private class functions + void update(); + + // Returns size of the amoeba + double size(); + + // Write amoeba position vertexes to stderr + void pos(); +}; + +#endif diff --git a/optimization/functions.h b/optimization/functions.h @@ -0,0 +1,24 @@ +#include <functional> +#include <armadillo> +using namespace arma; +using namespace std; + +int ncalls = 0; + + +// Rosenbrock's function +function<double(vec)> rosenbrock = [&ncalls] (vec p) { + ++ncalls; + double x = p[0], y = p[1]; + return (1.0f - x) * (1.0f - x) + + 100.f * (y - x*x) * (y - x*x); +}; + +// Himmelblau's function +function<double(vec)> himmelblau = [&ncalls] (vec p) { + ++ncalls; + double x = p[0], y = p[1]; + return (x*x + y - 11.0f) * (x*x + y - 11.0f) + + (x + y*y - 7.0f) * (x + y*y - 7.0f); +}; + diff --git a/optimization/header.h b/optimization/header.h @@ -0,0 +1,12 @@ +// Make sure header is only included once +#ifndef HEADER_H_ +#define HEADER_H_ + +// Define whether the program should output values of matrices +const bool verbose = false; +//const bool verbose = true; + +// Prototype for checking function +void check(const bool statement); + +#endif diff --git a/optimization/main.B.cpp b/optimization/main.B.cpp @@ -0,0 +1,93 @@ +#include <iostream> +#include <armadillo> +#include <functional> +#include "header.h" +#include "functions.h" +#include "downhill_simplex.h" +using namespace arma; +using namespace std; + +int main(int argc, char* argv[]) +{ + // Namespace declarations + using std::cout; + + // Calculate machine precision + double eps = 1.0f; + while (1.0f + eps != 1.0f) + eps /= 2.0f; + + cout << "\n\033[1;36m## Minimization with downhill-simplex, part B ##\033[0m\n"; + + // Try amoeba on Rosenbrock's valley function + cout << "\n\033[1;33m--- Rosenbrock's valley function ---\033[0m\n"; + int d = 2; + vec p(2); p[0]=5; p[1]=6; + if (verbose == true) + p.print("Initial simplex is chosen around the point:"); + vector<vec> simplex; + vector<vec> simplex_mod; + for(int i=0; i<d+1; ++i) { + simplex.push_back(p); + simplex_mod.push_back(p); + } + double dx = 1; + for(int i=0; i<d; ++i) { + simplex[i][i] += dx; + simplex_mod[i][i] += dx; + } + amoeba A(rosenbrock, simplex); + ncalls = 0; + A.downhill(10.0f*eps); + if (verbose == true) + A.low().print("Lowest point:"); + cout << "Amoeba converged after " << ncalls << " calls\n"; + amoeba A_mod(rosenbrock, simplex_mod); + ncalls = 0; + A_mod.downhill_mod(10.0f*eps); + if (verbose == true) + A_mod.low().print("Lowest point (modified downhill)"); + cout << "Amoeba converged after " << ncalls << " calls with modified method\n"; + + + // Try amoeba on Himmelblau's function + cout << "\n\033[1;33m--- Himmelblau's function ---\033[0m\n"; + vec p2(2); p2[0]=5; p2[1]=6; + if (verbose == true) + p2.print("Initial simplex is chosen around the point:"); + vector<vec> simplex2; + vector<vec> simplex2_mod; + for(int i=0; i<d+1; ++i) { + simplex2.push_back(p2); + simplex2_mod.push_back(p2); + } + double dx2 = 1; + for(int i=0; i<d; ++i) { + simplex2[i][i] += dx2; + simplex2_mod[i][i] += dx2; + } + amoeba A2(himmelblau, simplex2); + ncalls = 0; + A2.downhill(10.0f*eps); + if (verbose == true) + A2.low().print("Lowest point:"); + cout << "Amoeba converged after " << ncalls << " calls\n"; + amoeba A2_mod(himmelblau, simplex2_mod); + ncalls = 0; + A2_mod.downhill_mod(10.0f*eps); + if (verbose == true) + A2_mod.low().print("Lowest point (modified downhill)"); + cout << "Amoeba converged after " << ncalls << " calls with modified method\n"; + + // Return successfully + return 0; +} + +void check(const bool statement) +{ + using std::cout; + if (statement == true) + cout << "\t\033[0;32mPassed\033[0m\n"; + else + cout << "\t\033[1;31mFail!!\033[0m\n"; +} diff --git a/optimization/main.cpp b/optimization/main.cpp @@ -0,0 +1,68 @@ +#include <iostream> +#include <armadillo> +#include <functional> +#include "header.h" +#include "functions.h" +#include "downhill_simplex.h" +using namespace arma; +using namespace std; + +int main(int argc, char* argv[]) +{ + // Namespace declarations + using std::cout; + + // Calculate machine precision + double eps = 1.0f; + while (1.0f + eps != 1.0f) + eps /= 2.0f; + + cout << "\n\033[1;36m## Minimization with downhill-simplex, part A ##\033[0m\n"; + + // Try amoeba on Rosenbrock's valley function + cout << "\n\033[1;33m--- Rosenbrock's valley function ---\033[0m\n"; + int d = 2; + vec p(2); p[0]=5; p[1]=6; + p.print("Initial simplex is chosen around the point:"); + vector<vec> simplex; + for(int i=0; i<d+1; ++i) + simplex.push_back(p); + double dx = 1; + for(int i=0; i<d; ++i) + simplex[i][i] += dx; + amoeba A(rosenbrock, simplex); + ncalls = 0; + A.downhill(10.0f*eps); + A.low().print("Lowest point:"); + cout << "Amoeba converged after " << ncalls << " calls\n" + << "(Newton's root-finding method did this in 4181 calls)\n"; + + // Try amoeba on Himmelblau's function + cout << "\n\033[1;33m--- Himmelblau's function ---\033[0m\n"; + vec p2(2); p2[0]=5; p2[1]=6; + p2.print("Initial simplex is chosen around the point:"); + vector<vec> simplex2; + for(int i=0; i<d+1; ++i) + simplex2.push_back(p2); + double dx2 = 1; + for(int i=0; i<d; ++i) + simplex2[i][i] += dx2; + amoeba A2(himmelblau, simplex2); + ncalls = 0; + A2.downhill(10.0f*eps); + A2.low().print("Lowest point:"); + cout << "Amoeba converged after " << ncalls << " calls\n" + << "(Newton's root-finding method did this in 33 calls)\n"; + + // Return successfully + return 0; +} + +void check(const bool statement) +{ + using std::cout; + if (statement == true) + cout << "\t\033[0;32mPassed\033[0m\n"; + else + cout << "\t\033[1;31mFail!!\033[0m\n"; +} diff --git a/roots/Makefile b/roots/Makefile @@ -0,0 +1,74 @@ +# Define compiler +CC=g++ + +# Define compiler flags (show all warnings) +CPPFLAGS=-Wall -std=c++0x +CXXFLAGS=-Wall -std=c++0x + +# Define linker flags +LDFLAGS= + +# Define extra libraries to be dynamically linked +LDLIBS+=-larmadillo -lstdc++ + +# Compile optimized code +#CPPFLAGS+=-O2 + +# Compile debuggable code +#CPPFLAGS+=-g + +# Compile profilable code +#CPPFLAGS+=-pg +#LDFLAGS+=-pg + +# Define linker +LD=g++ + +# Filenames of source code +SRC=$(shell ls *.cpp) + +# Filenames of object files +OBJ=$(SRC:.cpp=.o) + +# Remove file type extension for binary filename +BIN=roots + +# The default "all" depends on A and B + +#all: A B + +#A: plot.A.png + +#B: plot.B.png + +#plot.%.png: fit.A.dat fit.B.dat data.A.txt data.B.txt +# gnuplot plotall.gp + +#fit.A.dat: $(BIN) +# ./$(BIN) data.A.txt fit.A.dat + +#fit.B.dat: $(BIN) +# ./$(BIN) data.B.txt fit.B.dat +# + +run: $(BIN) + ./$(BIN) 2> error.log + +$(BIN): $(OBJ) + @# Link object files together + $(LD) $(LDFLAGS) $(OBJ) -o $(BIN) $(LDLIBS) + @# Execute program and redirect stdout to file + @#./$(BIN) > out.txt + +clean: + @# Remove object files + rm -f $(OBJ) + @# Remove binary + rm -f $(BIN) + @# Remove datafiles and plot + #rm -f *.dat *.png + @# Remove error log + rm -f error.log +edit: + vim -p Makefile *.cpp *.h + diff --git a/roots/functions.h b/roots/functions.h @@ -0,0 +1,84 @@ +#include <functional> +#include <armadillo> +using namespace arma; +using namespace std; + +int ncalls = 0; + +// System of equations (first task) +function<vec(vec)> sys_2_eqs = [&ncalls] (vec p) { + ++ncalls; + double A = 10000.0f; + vec f(2); + double x = p[0], y = p[1]; + f[0] = A * x * y - 1.0f; + f[1] = exp(-x) + exp(-y) - 1.0f - 1.0f/A; + return f; +}; + +// Rosenbrock's function +function<vec(vec)> rosenbrock = [&ncalls] (vec p) { + ++ncalls; + vec f(1); + double x = p[0], y = p[1]; + f[0] = (1.0f - x) * (1.0f - x) + + 100.f * (y - x*x) * (y - x*x); + return f; +}; + +// Gradient of Rosenbrock's function +function<vec(vec)> rosenbrockGrad = [&ncalls] (vec p) { + ++ncalls; + vec f(2); + double x = p[0], y = p[1]; + f[0] = 2.0f * (1.0f - x) * (-1.0f) + + 100.0f * 2.0f * (y - x*x) * (-1.0f) * 2.0f * x; + f[1] = 100.0f * 2.0f * (y - x*x); + return f; +}; + +// Return derivatives of Rosenbrock's function (Jacobian matrix) +mat rosenbrockJacobian(vec p) { + mat J(2,2); + double x = p[0], y = p[1]; + J(0,0) = 2.0f + 1200.0f*x*x - 400.0f*y; + J(0,1) = -400.0f*x; + J(1,0) = -400.0f*x; + J(1,1) = 200.0f; + return J; +}; + + +// Himmelblau's function +function<vec(vec)> himmelblau = [&ncalls] (vec p) { + ++ncalls; + vec f(1); + double x = p[0], y = p[1]; + f[0] = (x*x + y - 11.0f) * (x*x + y - 11.0f) + + (x + y*y - 7.0f) * (x + y*y - 7.0f); + return f; +}; + +// Gradient of Himmelblau's function +function<vec(vec)> himmelblauGrad = [&ncalls] (vec p) { + ++ncalls; + vec f(2); + double x = p[0], y = p[1]; + f[0] = 2.0f * (x*x + y - 11.0f) * 2.0f * x + 2.0f*(x + y*y - 7.0f); + f[1] = 2.0f * (x*x + y - 11.0f) + 2.0f * (x + y*y - 7.0f) * 2.0f * y; + return f; +}; + +// Return derivatives of Himmelblau's function (Jacobian matrix) +mat himmelblauJacobian(vec p) { + mat J(2,2); + double x = p[0], y = p[1]; + J(0,0) = 4.0f * x*2.0f * x + + 4.0f * (x*x + y - 11.0f) + 2.0f; + J(0,1) = 4.0f * x+4.0f * y; + J(1,0) = 4.0f * x+4.0f * y; + J(1,1) = 4.0f * y * 2.0f * y + + 4.0f * (y*y + x - 7.0f) + 2.0f; + + return J; +} diff --git a/roots/header.h b/roots/header.h @@ -0,0 +1,20 @@ +// Make sure header is only included once +#ifndef HEADER_H_ +#define HEADER_H_ + +// Define whether the program should output values of matrices +const bool verbose = false; +//const bool verbose = true; + +// Choose vector length variable type +typedef int Lengthtype; + +// Choose floating-point precision +//typedef float Floattype; +typedef double Floattype; +//typedef long double Floattype; + +// Prototype for checking function +void check(const bool statement); + +#endif diff --git a/roots/main.cpp b/roots/main.cpp @@ -0,0 +1,79 @@ +#include <iostream> +#include <armadillo> +#include <functional> +#include "header.h" +#include "functions.h" +using namespace arma; +using namespace std; + +vec newton(function<vec(vec)> f, vec x_0, vec dx, Floattype eps); +vec newtonJac(function<vec(vec)> f, vec x_0, vec dx, Floattype eps, + mat (*J)(vec)); + +int main(int argc, char* argv[]) +{ + // Namespace declarations + using std::cout; + + // Calculate machine precision + Floattype eps = 1.0f; + while (1.0f + eps != 1.0f) + eps /= 2.0f; + + cout << "\nFinding the solution to the two-equation linear system:\n"; + vec x1(2); x1[0] = 2.0f; x1[1] = 2.1f; + vec dx1(2); dx1[0] = 1e-6f; dx1[1] = 1e-6f; + ncalls = 0; + vec root1 = newton(sys_2_eqs, x1, dx1, eps*10.f); + root1.print("Solution x:"); + sys_2_eqs(root1).print("f(x):"); + cout << "It took " << ncalls << " function calls\n"; + + cout << "\nFinding the minumum of the Rosenbrock's valley function:\n"; + vec x2(2); x2[0] = 5.0f; x2[1] = 6.0f; + vec dx2(2); dx2[0] = 1e-6f; dx2[1] = 1e-6f; + ncalls = 0; + vec root2 = newton(rosenbrockGrad, x2, dx2, eps*10.f); + root2.print("Solution x:"); + rosenbrock(root2).print("Rosenbrock at x:"); + cout << "It took " << ncalls << " function calls\n"; + + cout << "\nFinding the minumum of the Rosenbrock's valley function, Jacobian matrix predefined:\n"; + vec x2J(2); x2J[0] = 5.0f; x2J[1] = 6.0f; + vec dx2J(2); dx2J[0] = 1e-6f; dx2J[1] = 1e-6f; + ncalls = 0; + vec root2J = newtonJac(rosenbrockGrad, x2J, dx2J, eps*10.f, rosenbrockJacobian); + root2J.print("Solution x, Jacobian:"); + rosenbrock(root2J).print("Rosenbrock at x, Jacobian:"); + cout << "It took " << ncalls << " function calls\n"; + + cout << "\nFinding the minumum of the Himmelblau's function:\n"; + vec x3(2); x3[0] = 5.0f; x3[1] = 6.0f; + vec dx3(2); dx3[0] = 1e-6f; dx3[1] = 1e-6f; + ncalls = 0; + vec root3 = newton(himmelblauGrad, x3, dx3, eps*10.f); + root3.print("Solution x:"); + himmelblau(root3).print("Himmelblau at x:"); + cout << "It took " << ncalls << " function calls\n"; + + cout << "\nFinding the minumum of the Himmelblau's function, Jacobian matrix predefined:\n"; + vec x3J(2); x3J[0] = 5.0f; x3J[1] = 6.0f; + vec dx3J(2); dx3J[0] = 1e-6f; dx3J[1] = 1e-6f; + ncalls = 0; + vec root3J = newtonJac(himmelblauGrad, x3, dx3, eps*10.f, himmelblauJacobian); + root3J.print("Solution x:"); + himmelblau(root3J).print("Himmelblau at x, Jacobian:"); + cout << "It took " << ncalls << " function calls\n"; + + // Return successfully + return 0; +} + +void check(const bool statement) +{ + using std::cout; + if (statement == true) + cout << "\t\033[0;32mPassed\033[0m\n"; + else + cout << "\t\033[1;31mFail!!\033[0m\n"; +} diff --git a/roots/newton.cpp b/roots/newton.cpp @@ -0,0 +1,89 @@ +#include <functional> +#include <armadillo> +#include "qrfunc.h" + +using namespace arma; + +// Newton's method +vec newton(std::function<vec(vec)> f, vec x_0, vec dx, double eps) +{ + vec x = x_0; + int n = x.size(); + mat A(n,n); + vec y(n); + vec fx(n); + vec fy(n); + vec df(n); + vec Dx(n); + + do { + fx = f(x); + for (int j=0; j<n; ++j) { + x[j] += dx[j]; + df = f(x) - fx; + + for (int i=0; i<n; ++i) + A(i,j) = df[i]/dx[j]; + + x[j] -= dx[j]; + } + + // Perform QR decomposition + QR qr(A); + vec fx_neg = -1.0f*fx; + Dx = qr.backsub(fx_neg); + + double lambda = 2.0f; + do { + lambda /= 2.0f; + y = x + Dx * lambda; + fy = f(y); + } while (norm(fy,2.0f) > (1.0f-lambda/2.0f)*norm(fx,2.0f) && lambda > 0.02f); + + x = y; + fx = fy; + std::cerr << "Newton: f(x).norm() = " << norm(fx, 2.0f) << '\n'; + } while (norm(Dx,2.0f) > norm(dx,2.0f) && norm(fx,2.0f) > eps); + + // Return solution + return x; +} + +// Newton's method - the user supplies his own derivatives +vec newtonJac(std::function<vec(vec)> f, vec x_0, vec dx, double eps, + mat (*J)(vec)) +{ + vec x = x_0; + int n = x.size(); + vec y(n); + vec fx(n); + vec fy(n); + vec Dx(n); + + fx = f(x); + + do { + + // Get Jacobian matrix + mat Jx = J(x_0); + + // Perform QR decomposition + QR qr(Jx); + vec fx_neg = -1.0f*fx; + Dx = qr.backsub(fx_neg); + + double lambda = 2.0f; + do { + lambda /= 2.0f; + y = x + Dx * lambda; + fy = f(y); + } while (norm(fy,2.0f) > (1.0f-lambda/2.0f)*norm(fx,2.0f) && lambda > 0.02f); + + x = y; + fx = fy; + std::cerr << "Newton: f(x).norm() = " << norm(fx, 2.0f) << '\n'; + } while (norm(Dx,2.0f) > norm(dx,2.0f) && norm(fx,2.0f) > eps); + + // Return solution + return x; +} diff --git a/roots/qrfunc.cpp b/roots/qrfunc.cpp @@ -0,0 +1,90 @@ +#include <iostream> +#include <armadillo> +#include "header.h" +#include "qrfunc.h" + +// QR decomposition constructor +QR::QR(arma::mat &A) + : n(A.n_cols), + A(A), + Q(A) +{ + // Initialize output structures + R = arma::zeros<arma::mat> (n,n); + + // Perform QR decomposition straight away + decomp(); +} + +// Class deconstructor (equivalent to compiler destructor) +QR::~QR() { }; + +// Return system size +Lengthtype QR::size() +{ + return n; +} + +// QR decomposition function of Armadillo matrix. +// Returns right matrix R, and modifies A into Q. +// Uses Gram-Schmidt orthogonalization +void QR::decomp() +{ + Floattype r, s; + Lengthtype j; + for (Lengthtype i=0; i<n; ++i) { + r = dot(Q.col(i), Q.col(i)); + R(i,i) = sqrt(r); + Q.col(i) /= sqrt(r); // Normalization + for (j=i+1; j<n; ++j) { + s = dot(Q.col(i), Q.col(j)); + Q.col(j) -= s*Q.col(i); // Orthogonalization + R(i,j) = s; + } + } +} + +// Solve the square system of linear equations with back +// substitution. T is an upper triangular matrix. +arma::vec QR::backsub(arma::vec &b) +{ + Floattype tmpsum; + arma::vec x = Q.t() * b; + for (Lengthtype i=n-1; i>=0; --i) { + tmpsum = 0.0f; + for (Lengthtype k=i+1; k<n; ++k) + tmpsum += R(i,k) * x(k); + + x(i) = 1.0f/R(i,i) * (x(i) - tmpsum); + } + return x; +} + +// Calculate the (absolute value of the) determinant of +// matrix A given the Q and R components. +// det(A) = det(Q) * det(R), |det(Q) = 1| +// => |det(A)| = |det(R)| +Floattype QR::det() +{ + Floattype determinant = 1.0f; + for (Lengthtype i=0; i<n; ++i) + determinant *= R(i,i); + return fabs(determinant); +} + +// Calculate the inverse of matrix A given the Q and R +// components. +arma::mat QR::inverse() +{ + arma::mat inv = arma::zeros<arma::mat> (n, n); + // In vector z, all elements are equal to 0, except z(i) = 1 + arma::vec z = arma::zeros<arma::mat> (n); + + for (Lengthtype i=0; i<n; ++i) { + z(i) = 1.0f; // Element i changed to 1 + inv.col(i) = backsub(z); + z(i) = 0.0f; // Element i changed back to 0 + } + + return inv; +} diff --git a/roots/qrfunc.h b/roots/qrfunc.h @@ -0,0 +1,49 @@ +// Make sure header is only included once +#ifndef qrfunc_h_ +#define qrfunc_h_ + +#include <armadillo> +#include "header.h" + +// QR structure +class QR { + private: + // System size + const Lengthtype n; + + public: + //// Data + + // Input data + arma::mat A; + + // QR decomposition matrices + arma::mat Q; + arma::mat R; + + //// Prototype functions + + // Constructor prototype + QR(arma::mat &A); + + // Destructor + ~QR(); + + // Return system size + Lengthtype size(); + + // QR decomposition of Armadillo matrix A, returning R + // and modified A (=Q) + void decomp(); + + // Backsubstitution of triangular system + arma::vec backsub(arma::vec &b); + + // Absolute value of the determinant of matrix R + Floattype det(); + + // Inverse of matrix A + arma::mat inverse(); +}; + +#endif