Introduction to C++ and Object-Oriented Programming for Scientific Computing

1. Introduction to C++ and Object-Oriented Programming for Scientific Computing Shermane Austin,

2. Programming Exercise • Compute sin(x) using Taylor series • Note: convert degrees to radians radians = degrees * pi/180 • How many terms are needed? How might the question be answered? • Implement as a function, mysin(x) and compute sine for –pi to + pi

3. Boolean Data Type • Type used to generate true/false evaluation • Data type: bool • Examples: bool done = false; bool errordetected = false; while (!done && !errordetected) { } int a, b; … bool b1 = a==b; bool greater (int a, int b) {return a > b;}

4. Bool conversions • Bools are converted to ints in logical and arithmetic operations: bool a=true; bool b = true; bool c = a + b; //a+b = 2, and is non-zero , c=true When a result is converted back to bool type, Non-zero value = true and zero = false

5. File I/O • using the fstream library • File read: ifstream • File output: ofstream

6. File Input Example #include <fstream> #include <string> using namespace std; int main() { ifstream fin(“myfile.txt”); int ValueCount = 0; float Sum = 0, Value; while (fin >> Value) { Sum += Value; ++ValueCount; }

7. File1.cc //File I/O Example 1 - reading file input and writing file output #include <fstream> #include <iostream> #include <string> using namespace std; int main() { ifstream fin("mydata.txt"); int count = 0; float sum=0, value; while(fin >> value) { sum += value; ++count; } if(count > 0) { ofstream fout("output.txt"); float average = sum / count; fout << "Average: " << average << endl; } else { cerr << "No list to average" << endl; } return 0; }

8. //File I/O Example 2 - prompt user for filename #include <fstream> #include <iostream> #include <string> using namespace std; int main() { string FileName; cout << "Enter name of file to process: "; cin >> FileName; //fstream lib uses member of string, c_str() to send as argument ifstream fin(FileName.c_str()); int count = 0; float sum=0, value; while(fin >> value) { sum += value; ++count; } if(count > 0) { ofstream fout("output.txt"); float average = sum / count; fout << "Average: " << average <<endl; } else { cerr << "No list to average" << endl; } return 0; }

9. File I/O - checks if (! fin) { cerr << "Cannot open " << FileName << endl; exit(1); }

10. Standard Input, Output, Error • Standard Input – keyboard • Use cin • May be redirected • Standard Output – terminal • Use cout • May be redirected • Standard Error – terminal • Use cerr

11. Programming Exercise Data transformation Store the following xy-data in a file 0.1 1.1 0.2 1.8 0.3 2.2 0.4 1.8 Create a simple function to transform the y data (your choice) Write a C++ program to read the file, transform the y data using your function and write the new xy data to a new file.

12. Scope • Statement block – statements in {} • Statement blocks are legal anywhere a statement can be placed • Objects are used as per scope

13. Global Scope • Example: • int i • int main() • ... • } • void f() • i++ • char i

14. Function parameters • Default: pass by value • & - using to reference • constant parameter – function cannot change the object, e.g. • void myfunction (const int a, int b, int c) • b = a + 3; //legal • a = c + 5; //illegal assignment

15. Examining function arguments • Modify the simple routine to compute the area of a circle by changing the value of the radius in the function. • In your main program, print out the value of the radius after calling the function. • What happens?

16. More on functions -void as return type signifies subroutines similar to Fortran -return values are typed: double f(double x) { return sin(x)*pow(x,3.2); } -Default transfer of arguments: “call by value”, e.g. x1 = 3.2; q = f(x1); f takes a copy x of x1 and the value of x1 in the calling block (main or another function) cannot be changed.

17. Call by reference -Changing values of arguments requires call by reference -C uses pointers, int n=8; … somefunc(&n); /* &n is a pointer to n */ … void somefunc(int *i) { *i = 10; // n is changed to 10 } -Pointers also work in C++ but more commonly, references are used, e.g. int n=8; somefunc(n); void somefunc(int& i) // reference to i { i = 10; // n is changed to 10 }

18. Reference Use Always uses references for large objects This function implies a copy of x: void somefunc(MyArray<double> x) //example of a template { ... } Copying is inefficient if the object is large Instead use a reference in the function: void somefunc(MyArray<double>& x) { // danger: elements in the array can be changed x(5) = 10; } Manipulation of the array can be avoided using the const keyword: void somefunc(const MyArray<double>& x) { // can NOT manipulate the entries in x x(5) = 10; // illegal r = x(1); // ok }

19. Function Overloading • Functions with the same name that behave differently (parameters are used to determine which function to use) • E.g. • void Swap( int &x, int &y) • void Swap(float &x, float &y) • STL – standard template library – includes functions

20. Memory Allocation - C++ • Dynamic memory allocation – conserves memory, e.g. • int myarray[maxsize] is a static memory allocation, for duration of program. • C uses malloc and free – requires size of data type and use of pointers • C++ handles this very simply with functions new() and delete()

21. Dynamic memory allocation in C++ • C++ new and delete • double* x = new double[n]; //creates an array of n elements// • delete [] x; // deallocates memory for the array • // allocate a single variable: • double* p = new double; • delete p; • double* x = new double[n]; • .

22. Object-Oriented Programming Paradigm • Encapsulation • Objects are analogous to data structures • Data hiding or protection • Only member functions can access data • Inheritance • Derive new objects from existing objects • SparseMatrix from Matrix object • Polymorphism

24. C++ Classes A class is a collection of data structures and operations on them An object is an instance of a class The MatDense object is a good example: data: matrix size + array entries operations: creating a matrix, accessing matrix entries, matrix-vector products,.. A class is essentially a new data type A class can contain other objects. Complex types can be created that are more easily manipulated.

25. Class construct • Class – programmer-defined data type • Terminology • Data or attributes • Member functions • Define: • attributes or data • access: public, private (also protected and friend)

26. Generic Types of Member Functions • Constructors – create new objects • Inspectors – view attributes of objects • Mutators – change attributes of objects • Facilitators – operations on objects • Destructors – destroy objects

27. Class Rational Public Interface: Add(), Subtract(), Multiply(), Divide(), Equal(), LessThan(), Insert(), Extract() Data Members: NumeratorValue, DenominatorValue Other members: Numerator(), Denominator(), SetNumerator(), SetDenominator()

28. Client Program Example #include <iostream> #include <string> #include “rational.h” using namespace std; int main() { Rational r(3/4), s(2/3); Rational t(r); Rational Sum = r + s; Rational Product = r * s; cout << r << “+” << s << “=“ << Sum << endl; return(0); }

29. Rational Class Interface Description //rational.h #include <iostream> #include <string> Using namespace std; //Class description Class Rational { public: // member functions //constructors Rational(); Rational( int numer, int denom=1); //some arithmetic and stream facilitators …

30. Example: Rational Class • Want to represent rational numbers • Think of as a new data structure • Member data – numerator and denominator • Member functions would include: • Constructors • Arithmetic functions a/b + c/d = (ad+bc)/bd a/b – c/d = (ad-bc)/bd a/b * c/d = ac/bd a/b / c/d = ad/bc • Copy • Compare • Display • Extract value • Inspect values • Set values

31. auxilary functions and operators • operator overloading

32. An example Initialize A and x in a single function: void init (MyArray<double>& A, MyList<double>& x) { const int n = x.size(); int i,j; for (j=1; j<=n; j++) { x(j) = j/2.0; // or double(j)/2.0 for (i=1; i<=n; i++) { A(i,j) = 2.0 + double(i)/double(j); } } }

33. Complex Arithmetic Develop a class for complex numbers Note: C++ already has a class complex in its standard template library (STL), use that instead of the complex example #include <complex> using namespace std; … complex<double> z(5.3,2.1); complex<double> y(0.3); cout << z*y + 3;

34. Using the new Complex Class #include “complex.h" void main () { Complex a(0,1); // imaginary unit Complex b(2), c(3,-1); Complex q = b; cout << "q=" << q << ", a=" << a << ", b=" << b << "\n"; q = a*c + b/a; cout << "Re(q)=" << q.Re() << ", Im(q)=" << q.Im() << "\n"; }

35. Basic contents of Complex class Data members: real and imaginary part (of the complex number) Functions: construct or initialize complex numbers Complex a(0,1); // imaginary unit Complex b(2), c(3,-1); Write out complex numbers: cout << "a=" << a << ", b=" << b << "\n"; Perform arithmetic operations: q = a*c + b/a;

36. Declaration of Complex Class class Complex { private: // this is the default double re, im; // real and imaginary part public: Complex (); // Complex c; Complex (double re, double im = 0.0); // Complex a(4,3); Complex (const Complex& c); // Complex q(a); ~Complex () {} Complex& operator= (const Complex& c); // a = b; double Re () const; // double real_part = a.Re(); double Im () const; // double imag_part = a.Im(); double abs () const; // double m = a.abs(); // modulus friend Complex operator+ (const Complex& a, const Complex& b); friend Complex operator- (const Complex& a, const Complex& b); friend Complex operator* (const Complex& a, const Complex& b); friend Complex operator/ (const Complex& a, const Complex& b); }; //friend means that stand-alone functions can work on private parts (re, im)

37. Simple member functions Extract the real and imaginary part (recall: these are private, i.e., invisible for users of the class; here we get a copy of them for reading) double Complex:: Re () const { return re; } double Complex:: Im () const { return im; } Computing the modulus: double Complex:: abs () const { return sqrt(re*re + im*im); }

38. Using const const variables cannot be changed: const double p = 3; p = 4; // ILLEGAL!! compiler error... const arguments (in functions) cannot be changed: double myabs (const Complex& c) { return sqrt(c.re*c.re + c.im*c.im); } c cannot be changed, e.g. c.re=0 is illegal const Complex arguments can only call const member functions: double myabs (const Complex& c) { return c.abs(); } // ok because c.abs() is a const function Without const, i.e., double Complex:: abs () { return sqrt(re*re + im*im); } the compiler would not allow the c.abs() call in myabs

39. Overloaded operators C++ allows us to define + - * / for arbitrary objects The meaning of + for Complex objects is defined in the function Complex operator+ (const Complex& a, const Complex& b); // a+b The compiler translates c = a + b; into c = operator+ (a, b); i.e., the overhead of a function call It would be more efficient to have the function body inlined, i.e., copied into the calling code This is enabled by the inline keyword: inline Complex operator+ (const Complex& a, const Complex& b) { return Complex (a.re + b.re, a.im + b.im); }

40. Using inline If operator+, operator= and the constructor Complex(r,i) all are inline functions, c = a + b; that is, c.operator= (operator+ (a,b)); is transformed to c.re = a.re + b.re; c.im = a.im + b.im; by the compiler, i.e., no function calls This is particularly crucial for loops, e.g., Complex s; // have also Complex a and b for (i = 1; i <= huge_n; i++) { s = s + a; a = a/3.0; } Without inlining s=s+a we introduce two function calls inside a loop, which prevent aggressive optimization by the compiler

41. Constructors Constructors have the same name as the class The declaration statement Complex q; calls the member function Complex() A possible implementation is Complex:: Complex () { re = im = 0.0; } meaning that declaring a complex number means making the number (0,0) Alternative: Complex:: Complex () {} Downside: no initialization of re and im

42. Constructors with arguments The declaration statement Complex q(-3, 1.4); calls the member function Complex(double, double) A possible implementation is Complex:: Complex (double re_, double im_) { re = re_; im = im_; }

43. The Assignment operator a = b implies a call to a.operator= (b) – this is the definition of assignment Implement operator= as a part of the class: Complex& Complex:: operator= (const Complex& c) { re = c.re; im = c.im; return *this; } If you forget to implement operator=, C++ will make one (just memberwise copy – dangerous if you perform dynamic memory management!)

44. Copy constructor The statements Complex q = b; Complex q(b); makes a new object q, which becomes a copy of b Simple implementation in terms of the assignment: Complex:: Complex (const Complex& c) { *this = c; } this is a pointer to “this object”, *this is the present object, so *this = c means setting the present object equal to c, i.e., this->operator= (c)

45. Output function Output format of a complex number: (re,im), i.e., (1.4,-1) Desired user syntax: std::cout << c; any_ostream_object << c; The effect of << for a Complex object is defined in ostream& operator<< (ostream& o, const Complex& c) { o << "(" << c.Re() << "," << c.Im() << ") "; return o;} The input operator (operator>>) is more complicated (need to recognize parenthesis, comma, real numbers)

46. Header file We divide the code of class Complex into a header file Complex.h (extension .h) and a file Complex.cpp with the body of the functions (extension .cpp)

47. Example: class MyVector Class MyVector: a vector Data: plain C array Functions: subscripting, change length, assignment to another vector, inner product with another vector, ... This examples demonstrates many aspects of C++ programming Note: this is mainly an educational example; for professional use one should use a ready-made vector class (std::valarray for instance)

48. MyVector functionality I Create vectors of a specified length: MyVector v(n); Create a vector with zero length: MyVector v; Redimension a vector to length n: v.redim(n); Create a vector as a copy of another vector w: MyVector v(w); Extract the length of the vector: const int n = v.size();

49. MyVector functionality II Extract an entry: double e = v(i); Assign a number to an entry: v(j) = e; Set two vectors equal to each other: w = v; Take the inner product of two vectors: double a = w.inner(v); or alternatively a = inner(w,v);

50. MyVector functionality III Write a vector to the screen: v.print(std::cout); Arithmetic operations with vectors: // MyVector u, y, x; double a u = a*x + y; // ’DAXPY’ operation The proposed syntax is defined through functions in class MyVector Class MyVector holds both the data in the vector, the length of the vector, as well as a set of functions for operating on the vector data Users can only operate on the vector data through the offered functions MyVector objects can be sent to Fortran and C functions as follows: // v is MyVector call_my_F77_function (v.getPtr(), v.size(), ...) // array length