Home

FYSA120 C++ Numerical Programming Vesa Apaja November 11

1. using my_cvec lib1 part3 section4 vector lt complex lt double gt gt e A fast algorithm in a slow computer computes faster than a slow algorithm in a fast computer Those are my principles and if you don t like them well I have others Groucho Mars 187
2. vector lt double gt v 1 4 1 6 0 2 1 8 0 1 1 5 or do many push back s cout lt lt original vector lt lt endl vector_out v auto it min_element v begin v end vector lt double gt iterator it cout lt lt minimum element lt lt xit lt lt endl it max_element v begin v end cout lt lt maximum element lt lt xit lt lt endl 63 find element with some value it find v begin v end 0 2 dips end cout lt lt failed to find value lt lt endl else cout lt lt found value lt lt xit lt lt endl reverse some elements cout lt lt reverse starting from lt lt xit lt lt endl reverse it v end vector_out v j cout lt lt sorting lt lt endl sort v begin v end vector_out v return 0 64 Iterators may look a bit messy but they are easily hidden from view Especially if all you want is the number where the iterator points at Example 17 numerics easyusemin cpp C Standard Library algorithms another version of minimum search Finding min element trying to simplify usage include lt iostream gt include lt vector gt include lt algorithm gt include lt iterator gt namespace my create my own namespace double min element const std vector lt double gt amp v no explicit iterator We need only iterator return std min element v begin v end j j calling routine is clean and sim
3. S M l 1520 169 Obviously there are 2 possible spin combinations for a chain N spins long These can be economically coded in binary L4 4 1 00 00 0 dad 4 7 00 01 1 dad st 4 00 10 2 dad st f 00 11 3 ARNET A 20 1 Next write the Hamiltonian in the basis of these states as a matrix The matrix elements between states a and b are Has a H b a b 0 2 1 We impose periodic boundary conditions and tie the chain to a loop so that if i N then j i 1 1 This can be done neatly using the modulo operation j mod i 1 N or in C j i 1 ZN 170 The Hamiltonian matrix has plenty of zeros The non zero elements Hap are only the following 1 Term a 57 57 b states must be the same a b Thus Hab Haa and this happens in four cases 1 1 1 case 1 Haa 0 0 S755 0 0 ESE 4 1 1 1 case 2 Aaa bsec ios S785 mcus 5 5 E 1 1 1 case 3 Hag 1i 05 97 SF 14 05 Gl Se dey jo 1 1 case 4 Haa beluslie Si 5j elis gu Toa 4 2 Term a 3 S7 5 b The only non zero element from this term is 1 z case 5 Hab 1 05 15 57 55 Me UL N 3 Term aj S7 Sj b The only non zero element from this term is lc case 6 Hab 0i tjd 5 87 lei Ojos Now we are ready to compute 171 Algoritm Go through all states a Cases 1 and 4 Set Haa
4. so that the function main gives out the integer 0 out where Hmm to the operating system Safe to leave out entirely 15 2 1 Meaning of include lt iostream gt The line include lt iostream gt tells the compiler that the code has something defined in the header iostream which is part of the C standard library Two ways to include headers Hinclude lt iostream gt part of standard library include myheader h a home made header The difference is the header search path the latter uses the current path first after that the include path To see what headers were included try g M code cpp and to include a path g Imypath code cpp You get pretty far with these headers e include lt iostream gt input and output cout cin cerr clog here c means console which means screen include lt iomanip gt formatted output if you like printf O forget this include lt cmath gt C or include math h C or C mathematical functions sin cos log sqrt pow z is coded pow x 2 include complex complex numbers include lt fstream gt file handling also ofstream ifstream include lt string gt character strings include random pseudo random numbers include functional function objects Take a peek at what headers exist and what s in them www cplusplus com reference 16 Function styles double f double x
5. tal double f double x 1 j if it s short double f double x I ll be using double precision real numbers a lot as they are accurate enough for most needs Here the function f returns a number of type double the argument x is also a double 17 Example 2 basic fun_before cpp Function defined before main include lt iostream gt include lt cmath gt using namespace std double f double x return sin x cos xx x j int main 1 double x y gem y f x cout lt lt f lt lt x lt lt lt lt y lt lt endl return 0 f 2 333 0 482627 The compiler knows all it needs to know about the function f when it s met in main 18 Example 3 basic fun_after cpp Function defined after main include lt iostream gt include lt cmath gt using namespace std double f double int main double x y 2 owas y f x cout lt lt f lt lt x lt lt lt lt y lt lt endl return 0 l double f double x j return sin x cos xsx 5 When the compiler meets f inside main it wants to know what it takes as arguments and what it returns Therefore you need to provide it with a sneak peak double f double called prototype or function declaration You have to give the prototype also if the function is defined in another file Failing to declare a function results the error message fun after cpp 9 error
6. esa Ma E 1 10 7399900 H a a 0 25 only Boost else H a a 0 25 bbit abit ff ora gt e ard s x 43 lt lt Ws bbit flip i bbit flip j b static_cast lt int gt bbit to_ulong Hiab 0 53 23The or comments show how to do it using bitwise operators funny effective and unreadable for non experts 173 23 1 Add numbers to file names Sometimes it s convenient to have file names contain numerical parameters One option is to use stringstream objects see basic stringstream_ex cpp It s not convenient because one has to dig a C string out of the stringstream object s using s srt c_str P4 An easier way is shown below Example 61 basic string_to_filename cpp Open files res_ suffix is a real number Principle how to get a computed number to a file name include lt iostream gt include lt fstream gt include lt math h gt include lt string gt using namespace std int main const string str res_ for unsigned i 1 i lt 5 i ofstream out str to_string cos i out close opened files res_0 540302 res_ 0 416147 res_ 0 989992 res_ 0 653644 24A valid file name has to end with null character 0 like a character string in C 174 24 Lambda functions expressions I have already mentioned lambda functions or expressions in a few examples They are nice They are fast A lambda function is com
7. if b a and the bits i and j in states a and b are equal Cases 2 and 3 Set Haa i if the bits 7 and j are not equal in states a and b Cases 5 and 6 Set Hab gt if flipping spins 7 and j in state a gives state b Pseudocode iN Code A Sandvik Boston University H 0 set all elements of H to zero do a 0 2 N 1 N body states a do i 0 N spin i j i 1 N j i 1 periodic boundary conditions if i N j 1 if alil a jl H a a 1 4 else H a a 1 4 b flip a i j flip bits i ja j in a to get state b H a b 1 2 end if end do end do The Hamiltonian matrix is large full diagonalization may be slow but a few lowest eigenvalues can be found effectively using e g the Lanczos algorithm for example IETL has one 172 The C Standard Library has a convenient container for binary data std bitset Next choose a C matrix representation such as Boost Hermitian matrix typedef ublas hermitian_matrix lt cmplx ublas lower gt Matrix void FillHermitianMatrix Matrix amp H std bitset lt std numeric_limits lt int gt digits gt abit std bitset lt std numeric_limits lt int gt digits gt bbit dus ab 5 J Gto loo int NN H size int N log2 NN H clear zero H only Boost for a 0 a lt NN a abit std bitset lt std numeric limits lt int gt digits gt a for i 0 i lt N i j i4 1 N periodic boundary conditions Tf abit bil or at Ca
8. 0 0 afat 1 0 0 return GSL SUCCESS 154 GSL has driver functions in the header gsl_odeiv2 h which are much easier to use than the basic functions in the header gsl_odeiv h Whichever some preparatory steps need to be done e Choose the integration algorithm that is how the next point x t dt y t dt is computed Plenty of choice for this stepper here rk2 rk4 rkck rk8pd rk2imp rk4imp bsimp rklimp msadams and msbdf e Choose the control criterion that is when are shorter steps needed For example absolute error 1078 relative error 0 e GSL wants all information about the problem at hand in one data structure of type gsl odeiv2 system gsl_odeiv2_system sys func jac 2 amp mu This collects information about the functions their derivatives Jacobi matrix order and the parameter s Any number of parameters are passed as a pointer void params This naked pointer doesn t let the compiler check your data properly so better be careful A GSL driver function is set up auto solver g Sl odeiv2 driver alloc y new amp sys gsl odeiv2 step rk8pd 1le 6 1e 6 0 0 and it s called later like this int status gsl odeiv2 driver apply solver amp t ti y Passed values are the algorithm solver start point t end point ti and initial values y The solution at time ti is in the return value of y l here are four drivers to choose fr
9. If the code can have the situation where for example the numbers p 2 and q 1 happen to be exactly the same element you have aliasing Two pointers point to the same data In other words you had two different names for the same thing If you are absolutely sure this can never happen during the lifetime of pointers p and q let the compiler know I solemnly swear not to write a line of code where a pointer points to the same object as p int restrict p This is something the compiler can t deduce by itself it has to trust your word for it If you lie the code may run but the results are arbitrary In C use the restrict keyword often but with caution Further information demystifying the restrict keyword 49 9 C Standard Library A closer look 9 1 std vector container Example 11 basic vector cpp Create an empty container and push a few elements to it Create empty vector container and push data in include lt iostream gt include lt vector gt include lt cmath gt using namespace std int main vector lt double gt x for unsigned i 0 i lt 6 i x push back pow 1 2 5 77 push 1 2 to vector x cout lt lt i lt lt i lt lt x lt lt x i cout lt lt size of x lt lt x size lt lt endl cout lt lt final size lt lt x size lt lt n cout lt final capacity lt lt x capacity lt lt n return 0 final size
10. f was not declared in this scope Here we encounter the concept of scope C 11 standard made an addition the return type of a function can be declared like this auto f double gt double trailing return type Pretty useless in this example If you however think of how the return types of class methods become clear only after going through all the types of the arguments you may appreciate how this helps the compiler There is another use of foo bar means foo bar 19 2 2 Scope The scope tells the range of visibility Which parts of the program is an object defined and also who may know the thing even existd The basic ideas are global objects visible everywhere and local objects visible in a limited scope The scope limits are curly brackets boundaries of a block double z I m visible to anyone mess me up at will int main i double x I m visible between E double function f double k double x I m not the same x as in main z 10 but I m the one and only z A local loop variable exist only in the loop unsigned j 5 for unsigned i 0 i lt 10 i fif i j break cout lt lt i j when i lt lt i lt lt endl WON T WORK When the program leaves the loop the variable i ceases to exist The last line is not allowed Sounds modern invented in the 50 s 20 2 3 Simple file operations C idea cr
11. A possible GSL FFT header ifndef GSL_FFT_HPP define GSL FFT HPP include gsl gsl errno h include lt gsl gsl_fft_complex h gt namespace my_GSL_FFT 1 template typename T gt int fft T amp data int direction int status const size t stride l Size t n ifdef ARMA Armadillo data has no method size use n_elem n data n elem else n data size endif doublex pdata reinterpret_cast lt doublex gt amp data 0 if direction gt 0 status gsl_fft_complex_radix2_forward pdata stride n j else status gsl fft complex radix2 backward pdata stride n i if status GSL SUCCESS return 1 return 0 j j Hendif 144 22 2 1 Passing a pointer to complex data C 11 complains if you try to take the address of a complex number like this doublex pdata amp data 0 real may not compile Here data 0 real is an rvalue a temporary whose life is about to end see chapter 8 3 It does not have an identity that lives long enough to be used The first complex element does have an identity an lvalue so I tell the compiler to use that and pretend it s pointing to a double using reinterpret_cast The type change done by reinterpret cast is in your responsibility it basically tells the compiler trust me I know what I m doing doublex pdata reinterpret cast doublex amp data 0 Armadillo vector does not conform with a standard contain
12. calls r s amp result amp error gsl monte plain free s display_results method result error method MISER auto s gsl_monte_miser_alloc dim gsl_monte_miser_integrate amp G a b dim calls r s amp result amp error gsl monte miser free s display_results method result error method VEGAS warmup auto s gsl_monte_vegas_alloc dim warmup gsl_monte_vegas_integrate amp G a b dim 10000 r s amp result amp error display_results method result error 166 method VEGAS do gsl_monte_vegas_integrate amp G a b dim calls 5 r s amp result amp error display_results method result error method VEGAS continue while fabs gs1_monte_vegas_chisq s 1 0 gt 0 5 gsl_monte_vegas_free s return 0 double func double x size_t dim void params Static const double A 1 0 M_PI 4 M_PI x M_PI return A 1 0 cos x 0 cos x 1 cos x 2 void display_results string title double result double error const double exact 1 3932039296856768591842462603255 cout lt lt setiosflags ios fixed cout lt lt setfill lt lt setw 50 lt lt lt lt endl cout lt lt Method lt lt title lt lt endl cout lt lt setfill lt lt setw 50 lt lt lt lt end1 cout precision 8 cout lt lt result lt lt result lt lt lt lt error lt lt endl cout lt l
13. interpolate using natural cubic my spline x y xx yy esplxzme output xx yy interpolate using natural Akima my spline x y xx yy Akima output xx yy return 0 a out gt tt gnuplot M Peis 2 w p pe 7 po 1 8 E mor polares PEMPE takima you could use std generate XX spline spline x size x i 1 x O ixdx 2 3 3 wy lis Ton e eee 5 79 ab 24 y l9 161 20 15 10 ALS XX SSSSLA RTE LINO SA known points cubic akima 22 5 GSL Monte Carlo integration Let s compute an N dimensional integral over a hypercube by ba bn dai dry f din F 21 09 00 s a 2 an d a GSL offers three ways see Wikipedia Monte Carlo sampling 1 Plain sample random points inside thye hyper cube a very crude method 2 MISER Press amp Farrar stratified sampling algorithm sample points from areas that give largest error 3 VEGAS Lepage combines both stratified sampling and importance sampling sample points from areas that affect the result the most The example code evaluates the integral fa fa fa 1 39320392 z 1 TE 7 Jo 0 H o 1 cos z cos y cos z Prerequisites e Random number generator e Integrand Data Type gsl monte function This data type defines a general function with parameters for Monte Carlo integration double f double x size t dim void params this function sh
14. lt 5x vector Wu for auto elem x doubleout elem cout lt lt endl return 0 This is just as good as the std for_each example and compiles faster 89 This is where for_each shines e Access only part of the elements for each x begin x begin 2 doubleout output two values from the beginning e Access all elements or some elements and use anything a class can contain Example 32 numerics for_each limited sum cpp std for_each can do things range for loops can t include lt iostream gt include lt vector gt include lt algorithm gt struct LimitedSum vold operator ant a if 4 gt 1 sm 1 int sum 0 js int main stadir vector intix 122 324 LimitedSum lim std for_each x begin x end LimitedSum std cout lt lt Limited sum lt lt lim sum lt lt n Notice how we access the member of the class LimitedSum after a call to std for_each The return value of for_each is the very same class object that the function object in the argument is Function object are discussed more in chapter Without storing the return value we had no access to the data member sum See the details on the next page 90 12 3 std for each in detail One way std for_each can be implemented is this Example 33 advanced foreach_template cpp std for_each works essentially like this template lt class Iter class Func gt Func for_each
15. lt lt lt lt cos x lt lt lt lt tan x lt lt std endl j y int main TrigFuns comp comp 000 return 0 Here comp is an object but used as if it were a function hence it s a function object 99 Why use a function object Why not an ordinary function 1 Function objects are not passed as pointers so the compiler can easily inline it insert to its place As you saw in exerc 1 pdf one way to send a function to a function is as a function pointer void apply double g double double x std cout lt lt g x lt lt std end1 but it s also simple to wrap the function in a class and make it a function object include lt iostream gt include lt cmath gt class Func public double operator double x return sin x E void apply Func g double x std cout lt lt g x lt lt std endl int main Func f apply f 10 0 return 0 2 A class can contain data members such as counters accessible only to class methods A function object can easily perform complicated tasks For example a function object can compute the cosine of all input and simultaneously compute the sum of these cosines this is done in numerics forarch functor2 cpp 100 Example 37 numerics foreach_functor cpp A function object used with std for each include lt iostream gt include lt vector gt include lt algorithm gt include lt
16. no point writing result swap a b nonsense T c std move a create class T object c and move the contents of a to object c std move make sure you use the correct move function from the std namespace std move a fix a to a movable object one whose contents can be stolen What can be swapped with this template C 11 Type T shall be move constructible and move assignable or have swap defined for it 8C 11 replaced copy T c a for move T c std move a search the net for C move semantics 66 Many C Standard Library algorithms are compact and efficient Take for example this one this page has only utility functions Example 19 numerics algo_permutations cpp C Standard Library permutation algorithm Finding permutations using C Standard Library algorithm next_permutation Ti include lt iostream gt include lt vector gt include lt algorithm gt include lt iterator gt using namespace std IGAC to prime Ore a vector template lt typename T gt void vector_out const vector lt T gt v for anbto ziy contro 1 cout lt lt endl int factorial const int n ime racy ls if m lt 1 return l else fact n factorial n 1 recursion return fact 67 int main constant em vector lt int gt state for int i 0 i lt N i state push_back i sort state begin state end make sure next permutation starts ok
17. subtract average from all elements deprecated in C 11 transform x begin x end xx begin bind2nd minus lt double gt average C 11 version using a lambda function transform x begin x end xx begin average const double amp elem return elem average sigma sqrt inner product xx begin xx end xx begin 0 0 N 1 75 Example 24 numerics mystatistics ct 1l cpp get_stats using plain C 11 C 11 version Computes average and standard deviation for data stored in std vector include lt vector gt include lt cmath gt namespace mystat_C11 void get_stats const std vector lt double gt amp x double amp average double amp sigma average sum i 1 N x i N fi Ml sum i 1 N x i x 2 I sigma sgrt LC N 1 int N x size average 0 sigma 0 for auto x_i x average x_i average N for auto x_i x sigma pow x_i average 2 sigma sqrt sigma N 1 j I used the namespace mystatC11 to tell from the previous srd version The header in file mystatistics_c 11 hpp has the same namespace You may disagree but to me this version is much easier to grasp that the C Standard Library version The code and the math formulas are more easily compared The range for loop is neat and safe For completeness I present the main program to test both variants 76 Example 25 ma
18. 0 a 0 B 0 c 0 but a compiler can be so stupid that it adds up the whole million element vector A B C puts it to D and only then looks for D 0 MBlitz still works behind scenes Last time I looked Scipy Scientific Python module used parts of the Blitz library IBvValarray tries to avoid intermediates temporaries using proxy objects Most libraries prefer the expression template technique due to its generality 108 Continue reading if you the previous discussion didn t drop out of your cache Think how templates can help to avoid temporaries If a compiler meets a line of code it can t immediately recognize it starts looking for a suitable template model Finding one it instantiates the template i e brings it alive This template can itself instantiate other templates and so on Apart from instantiating another template a template can instantiate another copy of itself recursive template This recursion is in the heart of many expression templates For identification templates have template parameters which tell exactly what kind of model to instantiate template lt parameters gt blaablaa Letting templates instantiate new templates lets one express the addition of objects as a tree A B C The vertices could represent operators Let A B C D be class Array objects just some class data could be in a std vector Then the tree could be Array A B C D DE ANC
19. 11 loops 73 ivector The principle of get_stats that uses std algorithms The algorithm std accumulate computes the sum of elements unless told to do something else The container std vector has has method size which gives the number of elements For standard deviation we compute how much elements differ from the average using the algorithm std transform There is also std bind2nd which is left for the reader It s there because std transform can do only operations for one or two elements so the problem of feeding in the average as well is solved using a lambda function old days one used std bind2nd but it s deprecated If you at this point feel exhausted don t be alarmed I feel the same 74 Example 23 numerics mystatistics cpp get_stats using std algorithms Computes average and standard deviation for data stored in std vector include lt vector gt include lt algorithm gt include lt numeric gt include lt functional gt include lt cmath gt namespace mystat void get_stats const std vector lt double gt amp x double amp average double amp sigma average sum i 1 N x i N Wi V sum i 1 N x_i lt x gt 72 Se SE ox m IGE eigna gopi com ere im E ete cites i dd N 1 N 1 using namespace std int N x size s average accumulate x begin x end 0 0 N sum up values vector lt double gt xx x xx x lt x gt
20. 6 size of the container increased automatically final capacity 8 more space reserved than needed just in case you extend the vector 50 9 1 1 Iterators Iterators are things pointing to elements of a container bit like custom made pointers Example 12 basic vector2 cpp Print all elements of a std vector using iterator print all but 2 elements of a vector container using iterator it include lt iostream gt include lt vector gt include lt cmath gt using namespace std int main vector lt double gt x puse 22 for unsigned i 0 i lt 6 i x push_back pow i 2 use iterator to print all but the two last elements for auto it x begin it x end 2 it cout lt lt it lt lt j cout lt lt endl To print all elements use the range for loop for auto elem x cout lt lt elem lt lt cout lt lt endl return 0 7 output fi Q 1 AE ffi 1 AO IG 25 51 The keyword auto it is a message to the compiler figure out yourself the type of it you can I recommend C 11 auto type let the compiler figure out the type In the example auto spared you from typing the ugly type vector lt double gt iterator it type of an iterator One way to fill a std vector is std fi11 it s very fast if you need to reset an existing std vector std vector lt double gt y 5 std fill y begin y end piod stds
21. A class Bird may have an object crow which tells the compiler that a crow is a bird The method comp energy is conveniently a class method a member function Rough outlook class ClassName here private things pub iici here public things and member functions used to manipulate and show contents of private date this line may already contain objects of this class One possible class describing the state of the system just a draft improve it later class SystemState double E puc void comp_energy compute energy double energy const return E read energy statel state2 This defines the class and immediately constructs two objects state1 and state2 You can do this later as well SystemState statel 22 4 Member function with const or noexcept after it C 11 has several compiler directives that facilitate code optimization Class MyClass pub iici void memberi const does not change any data member of the class except mutable void member2 noexcept does not throw any exceptions void memder3 const noexcept does neither throw nor change any data member unless mutable throw refers to a mechanism used to make exceptions errors or warnings More about this in chapter for throw catch noexcept can potentially lead to faster code see chapter In short if a function cannot fail use noexcept In numerics you seldom want
22. DUp EE GB Be Pc Ee edd qme ewe Reg mee hc ee dra 44 Lob aw Pu EUR EOS aoe GE Foy Boh edo anG 45 8 5 The restrict keyword in C read on spare time esop dec des 48 9 C Standard Library A closer look 50 9 1 std vector contaimer 2 s gt s skarr ee a a a Y We vo aa 50 O I lt rat rS x6 2 9o ae o9 w a a e PR EE Re D e 51 pa Gh Sek p a Gre OG a a i NON aoe A eae 53 9 2 Moving not COpying s e 4x pia 3c Rok OY Ew ee EBRD 4 EUR ROS R3 XU Pon SUR SOR RU 55 9 3 std valarray class x ouo a b ke OE Ge EUR Ok Bee ee we Ae EUR E OR d Ree Rn 58 MEM MM 59 9 5 Stream iterators read on spare ANN 61 Panpa PH A be AE A a 62 9 6 1 std min_element max element find sort reverse 63 9 6 2 std swapis a template ee 66 TT 72 9 8 std complex complex numbers and arithmetics eA 78 79 11 Operator overloading all for readability 82 88 D2 Stade stor each era oy Soy See zu Ry Eee e emp M3 Ee REGE ble ee RUE ae ee 88 12 2 When to use std for each 7 2 less sS ss oss 89 12 9 std for each in detalll ua 0 ose ko Rowe BOY X a Fox rie REOR B JR e e RO Re Ph s 91 i44 Webokg b X SRORSECRCEO REOR xg SORS EP we eee g Gig Rec AUR HU Rede 93 LES duo UPPER TE UL 94 12 6 C Standard Library algorithms stateful objects and std ref o 95 13 A few things that can potentially speed up your code 13 1 noexcept no throw quarantee ee 13 2 constexpr Compile time constant expressions 13
23. a 5 b 10 integ integ a b return 0 no args to integ integrating from O to 1 two args to integ integrating from 5 to 10 80 Example 28 basic function_overload2 cpp Functions second argument is by default 1 0 Function fun has a default value 1 0 for the second parameter b include lt iostream gt using namespace std double fun double a double b 1 0 IMPORTANT LINE int main double a 5 b 10 fun a fun a b return 0 j double fun double a double b if b 1 cout lt lt a lt lt a lt lt default case b 1 lt lt endl j else cout lt lt a lt lt a lt lt not default case b lt lt b lt lt endl j return 1 0 Just test a 5 default case b 1 a 5 not default case b 10 Default value is given only in the function declaration This can make the default value hard to find 10 Technically it s possible to set default values in function definition but I strongly advice you not to You code will not be 81 11 Operator overloading all for readability In section 9 8 we learned that the complex number multiplication is done correctly by the operator The way this was achieved is operator overloading an operator can be told to do a slightly different operation depending on the data type Operator overloading can greatly improve code readability Operator overloading is something you don t have to learn to do but you will appre
24. a container Recent compilers should produce as fast code with std valarray as with std vector because the use the very efficient expression template technique The closest relative to std valarray is in the Boost library std valarray boost ublas vector To be honest std valarray is not popular among C programmers and some recommend to stay away from 1 58 9 4 typedef and using give a good name to your data type IMHO the lack of built in numerical data structures such as vectors and matrices prevents C to be a perfect language for numerics I face a dilemma How to code vectors and matrices A few options 1 Should I stick to std valarray or std vector and live with the deficiences In a small project that needs to be more portable than numerically less ambitious Yes In a serious numerical project No go to option 3 2 Should I write a class for a vector myself Never write your own classes for vectors or matrices You may see C numerics cources in the web and books telling you how to make a matrix class Don t you ll make mistakes waste your time and get a slow program Sure these are deceivingly simple tasks at a first glance but once you hit the speed bumps you can hear your tailpipe clank to the floor Go to option 3 3 Should I use a vector or matrix from an external external to C standard library YES But not directly utilize typedef or using This has the probl
25. are used for optimization Genetic Algorithm Library e Bartlomiej Filipek gives examples how to use the standard algorithms instead of raw loops Top Beautiful Cplusplus std Algorithms Examples 186 28 Farewell words for C numerical programmers e Prefer C Standard Library containers over dynamic arrays Don t transform one to another e For math matrices and vectors use Armadillo MTL4 Blaze don t code them yourself e Use C Standard Library algoritms not loops e Forget type safety it s the least of C numerics problems e Move don t copy For random number generators use std ref in std generate e Pass functions as function objects lambda functions or use std function Avoid pointers to functions e Keep the daily stuff the interface code clean Hide ugly details library calls to headers fft data l n overloaded fft library call in a header gsl fft complex radix2 forward data l n Sac JO Wes as winy 2 clo sea e Use namespace protection for you own functions my get_data e Templates are you friends e Use operator overloading to make cleaner code Test all aspects outfile lt lt myobject overloaded lt lt more readable than a function call A B M c overloaded and B M complex matrices c real A multiply_complexmat_realvec sum_complexmat_complexmat B M c C Oh no e Make data types short and descriptive with typedef or using
26. be able to use them you have to e learn how to call fortran subroutines from C tricky business e find a third party wrapper C code that calls fortran or the C versions of the libraries Calling old fortran 77 FORTRAN library routines is pretty safe but calling for example new fortran 95 parts of the MKL library has spawned a warning Caution Avoid calling BLAS 95 LAPACK 95 from C C Such calls require skills in manipulating the descriptor of a deferred shape array which is the Fortran 90 type Moreover BLAS95 LAPACK95 routines contain links to a Fortran RTL 19The netlib versio of CLAPACK was f2c d translated from fortran to C I find this approach and buying with photocopied money more than a bit troublesome Why not use the real thing 127 A C programmer would prefer a package written in C There was once LAPACK 4 but the project faded away around year 2000 It was superseced by TNT Template Numerical Toolkit which too looks pale Boost is very much alive but mostly not numerics stuff The Boost collection boost ublas is outdated limited and the interface is too far from math notation only my personal opinion Here is an incomplete list of projects providing well coded C algorithms e OpenBlas and GotoBlas2 original code by Kazushige Goto is optimized to processor architectures Ne halem Sandy Bridge Haswell AMD Bulldozer Piledriver etc e Armadillo NICTA Australia 20
27. call swap a b Ineffective code does copying void swap int amp a int amp b is 5 Cm o GE C style swap handles pointers function call swap amp a amp b void swap int a intx b xsv eg c a a b b c One benefit of passing by value is that the original value is safe only a copy is sent out In C you can give the reference const attribute to protect it from changes Caveat some libraries seem to ignore this attribute You can never be absolutely sure A malicious software package can also pretend it s using input by value only but snatch the reference instead and change your input data The function calls look the same Not that I know of any such packages 41 8 2 Can a reference be unsafe Alas a reference is not completely safe it s still possible to shoot yourself in the leg A dangling reference is one way to mis use an address of an object after the object has ceased to exist Example 9 basic dangling_reference cpp A dangling reference Unsafe code example ff Twin ia Weillesedincl e Ely QUe include lt iostream gt int amp get number ims mil s amg e sp ml mil return ae mal s int main int number get number std cout lt lt number is lt lt number lt lt n Think where the reference r_n1 is referring to The line inside get_number back in main C r nl int amp nl
28. choose 0 if type compare akima 0 type compare Akima 0 choose 1 switch choose case 0 spline gsl_spline_alloc gsl_interp_cspline x size break case 1 spline gsl_spline_alloc gsl_interp_akima x size break gsl_spline_init spline amp x 0 amp y 0 x size 159 auto posy yy begin for auto posx xx begin posx xx end posx posyt gsl_spline_eval spline posx acc gsl_spline_free spline gsl interp accel free acc j j endif Example 59 numerics gsl spline cpp Test code for GSL spline header AES ISO sb dca SS pico pS onto Si include lt iostream gt include lt iomanip gt include lt vector gt include lt algorithm gt home made call to GSL spline include gsl_spline hpp using namespace std typedef vector lt double gt vec void output const vec amp x const vec y cout lt lt fixed lt lt setprecision 8 for auto i 0 i lt x size i cout lt lt x i lt lt lt lt y i lt lt endl cout lt lt n n n int main veo x V0 et 7 known points x y 160 double t b D05 for unsigned i 0 i lt x size i xa te yli 10 0 tx t vt t 1 0 y 4 10 0 1ift one point up to demonstrate locality nonlocality output x y interpolated points xx vec xx 100 yy 100 double dx x x size 1 x 0 x for unsigned i 0 i lt xx size i
29. cout to cout lt lt endl Use the keyboard to give characters and press ctrl d when you are done Here std copy is an algorithm and I m going to tell more about them next Tistream stream in ostream stream out 61 9 6 C Standard Library algorithms and utilities The C Standard Library has many useful methods search algorithms sorting algorithms etc As soon as you know what you need visit the pages www cplusplus com reference algorithm www cplusplus com reference utility Examples If v and w are std vector containers and it is an iterator of that container type then it std max element v begin v end returns the iterator it pointing to the largest element the largest element is it std sort v begin v end sorts the container std swap v w swaps the contents of containers v and w C 11 has this in the header lt utility gt 62 9 6 1 std min element max element find sort reverse Example 16 numerics algo_minmax cpp C Standard Library algorithms minimum maximum and search of an element Finding elements from a vector container include lt iostream gt include lt vector gt include lt algorithm gt include lt iterator gt using namespace std utility to print out a vector template lt typename T gt void vector out const vector lt T gt v for auto xv cobren cout lt lt endl int main
30. do in C In C the common practise is to put to headers only function declarations 12 What if get_stats is part of a huge pile of code where another get_stats happens to exist with the same type of arguments You get an ambiquity error also called name collision Let s use namespace encapsulation Define our own namespace where get_stats lives so that we can be sure which of the many get_stats should be invoked Example 22 numerics mystatistics hpp A better get_stats header ifndef MYSTATISTICS HPP define MYSTATISTICS HPP include lt vector gt namespace mystat j Hendif void get stats const std vector lt double gt amp double amp double amp A remark about style and good habits refrain from taking the whole std namespace unnecessarily using namespace std don t do this namespace mystat This is impolite because if you ever give this header to a college to be used in his her code the poor fellow gets the whole std namespace wanted or not This easily leads to to name collisions if the college was not carefully protecting his her cute and short function names such as get Namespace encapsulation saves the day mystat get is always recognizable as a mystat function and myclasses is your vector not std vector The next task is to code the function itself I present here two styles the first utilizes std algorithms the other plain C
31. end1 eigvec j i lt lt trans x lt lt trans x lt lt endl Armadillo calls the LAPACK routine dsyev double symmetric eigenvalue so link e g with llapack 130 One possible output 1 6804 0 5919 1 2605 0 5919 0 6704 1 3971 1 2605 1 3971 0 7296 1 7146 0 9135 1 2307 0 9279 0 7969 1 0895 Oth eigenvalue 1 30315 x 0 3020 check Ax lambda 0 3020 ith eigenvalue 0 803214 x 0 1275 check Ax lambda 0 1275 2th eigenvalue 0 281196 x 0 3701 check Ax lambda 0 3701 1 7146 0 9135 1 2307 0 2832 1 4111 0 1113 0 1113 0 5699 0 5699 0 3639 0 3639 0 9279 0 7969 1 0895 1 4111 0 3134 0 0491 0 0491 0 7894 0 7894 0 2307 0 2307 131 8003 8003 1882 1882 3048 3048 5035 5035 0197 0197 7645 7645 18 2 Blaze example Example 47 numerics blazel cpp Blaze Simple vector and matrix operations source Blaze web page https code google com p blaze lib wiki Getting Started eh An Example Involving Matrices added main include lt blaze Math h gt using namespace blaze int main Instantiating a dynamic 3D column vector DynamicVector lt int gt x 3UL lO 4 S e ae Instantiating a dynamic 2x3 row major matrix preinitialized with 0 Via the function ealik operator three values of the matrix are explicitly set to get the matrix M Caoa M Oey DynamicMatrix
32. failing functions 5 Code rotting and deprecated functions You want to replace an old function with a better one but don t want to upset your boss by making his legacy code fail to compile Use C 14 deprecated C 14 deprecated void foo int You can add a message too deprecated Sorry Boss newfoo replaces foo A call to foo causes g to emit warning void foo int is deprecated Sorry Boss newfoo replaces foo Wdeprecated declarations and a C 11 compiler emits the do not quite understand message warning deprecated attribute directive ignored Wattributes 23 To work the class SystemState still needs the definition of the method comp_energy If it s short it can be defined in situ on the spot in side the class as we did for energy but you are free to define it later member function of class SystemState void SystemState comp_energy E 3 12515e 10 just anything for testing j The variable E used byt the method is exactly the same E as in the object not just any double E Now the energy of the state can be computed and read int main statel comp_energy compute statel s energy cout statei energy end1 print state1 s energy return 0 This introduced one benefit of a class data encapsulation which means that data can be protected from inadvertent changes or hackers Just make it private and give
33. function does not return at all so no need to have a return in a function that only throws In the example the message is output to stream std cerr similar to std cout but specialized for error outputs This makes it possible to separate error output from normal output In bash typing a out 2 gt program err causes normal std cout output to screen and std cerr output to go to file program err There is also std clog for log outputs 139 22 Gnu Scientific Library GSL GSL in wikipedia GSL is free and written in C For C users the linkage is made easy The library header files automatically define functions to have extern C linkage when included in C programs This allows the functions to be called directly from C Example 51 numerics gsl_bessel cpp Bessel function JO Compilation g gsl_bessel cpp lgsl lgslcblas or g gsl bessel cpp gsl config libs try gsl config libs in linux include lt iostream gt include lt iomanip gt include lt gsl gsl_sf_bessel h gt using namespace std int main void double x 5 0 double y gsl_sf_bessel_JO x pontes ID We Y or using iomanip cout lt lt fixed lt lt setprecision 18 std COMES SOS ESAS lt lt endi return 0 JO 5 0 177596771314338264 140 22 1 GSL statistics double gsl stats mean const double data size_t stride size t n Function This function returns the arithmetic mea
34. gener define generator void initrng void Suele BOOL Ide UNE y afi first auto seed static cast uint fast32 t std time 0 cout lt lt seed seed endl gener seed seed first false j j double unirand void Stawie Bo wiles vine 115 static function lt double gt rnd dd irse initrng uniform_real_distribution lt double gt unif_dist 0 1 rnd bind unif_dist gener Wish falser return rnd double gaussrand void Er A MO static function lt double gt rnd if first 1 initrng normal_distribution lt double gt norm_dist 0 1 rnd bind norm_dist gener first false return rnd double gaussrand2 void warning does not initialize generator static normal_distribution lt double gt norm_dist 0 1 return norm dist gener double exprand void statici oo M sib CS Static function lt double gt rnd df forest initrng exponential_distribution lt double gt expo 116 rnd bind expo gener first false return rnd 117 16 Boost library Boost library has served as a test bench for ideas that may be included in a future C standard C 11 took many of its new features from Boost It offers many extensions to the C standard library algorithms special functions differential equation solvers and many more The documentation is http www boost org doc lib
35. is op C first should translate as expression X lt Array plus Array gt C second should translate as expression X lt X lt Array plus Array gt plus Array gt Here X is a vertex for objects left and right the operator plus in between The code advanced recursive template cpp shows how such a tree structure X X Array plus Array plus Array is created The example advanced expression template cpp computes the sum of std vector s using expression templates l6Excuse me for mixing template arguments and template parameters I honestly don t know their difference 109 What do we gain Expression templates can delay the evaluation until the sign is reached lazy evaluation The compiler goes through the whole tree at compilation time instantiating vertices X as many as needed Reaching the end it knows exactly what is needed such as only D 0 The idea is that operations return expressions not results Expressions are something that can be further manipulated during compilation Compiling the code with a recursive template essentially generates changes to the code or the code changes itself This is the essence of template metaprogramming templates can generate new code during program compilation Alas compilation takes longer ee a o arguments exist IBM ven ae on their linux ee pages publ ldi f em aptly I Should say A template argument for a template template parameter
36. its name only just like in C e The compiled fortran subroutine is in object code file dgemm o or inside a library with an unerscore dgemm is shown as dgemm e fortran is always pass by reference so all arguments have to be pointers Here double a points to the start of a double array 134 As if this were not enough there is one more problem matrix storage in fortran and in Matlab is column major but in C C row major Below is a figure about the two storage habits 123 4 math notation e 6 7 5 fortran Matlab row major 15263748 C C column major 12345678 Both are natural but moving between languages you have to change indexing before and after a call to dgemm For square matrices this means transpose If you have defined a static array in C style like this const int n 100 double a n n the fortran call must include const int n Finally by default indices begin from 0 C C or from 1 fortran fortran 1 N meaning V 1 V 2 V N C 0 N 1 meaning V 0 V 1 V N 1 135 20 Fixed size arrays in C plain array and std array Two different things that can be called arrays e Plain C style array A static array is made like this double array1 10 memory allocation for 10 elements int array_int 3 6 12 memory allocation and fill with values Dynamic arrays are created with keywords new and deleted with delete I m not te
37. long int constexpr factorial int n i int main return n gt 0 n x factorial n 1 1 one return statement is allowed constexpr long int fact13 factorial 13 13 is computed in compile time If you look at the assembly code after compilation g S code cpp and look at code s you see that indeed fact13 is set equal to 6227020800 constexpr means that everything the expression may be known at compile time If all is known everything is done in compile time If not the function is treated as a normal function evaluated in run time It a win win situation Although compilation takes longer Numerical constants are by definition already constant expressions so no need to replace const s with constexpr s there Extra information const means that once initialized the value cannot be changed constexpr is more than const it s a constant expression Functions declared constexpr can compute the allocator template parameter of e g std list lt gt see stackoverflow difference between constexpr and const If you know C metaprogramming and Haskell you might be interested to read What Does Haskell Have to Do with C 97 For curiosity the same task with template metaprogramming Template metaprogram include lt iostream gt just for testing template lt int N gt struct Factorial 1 static const long int value N Factorial lt N 1 gt value E tem
38. main const int n 100000000 ofstream fileout timing double time reftime for int nc 1 nc 5 nc omp_set_num_threads nc override OMP_NUM_THREADS double sum 0 auto t0 high resolution clock now pragma omp parallel for reduction sum for int i 0 i lt n i sum 1 cos i thread sums put together in the end auto ti high_resolution_clock now auto d duration_cast lt milliseconds gt t1 t0 time d count 1000 0 cout lt lt sum lt lt sum lt lt cores lt lt nc lt lt timing lt lt time lt lt s lt lt endl if nc 1 reftime time fileout lt lt nc lt lt lt lt time lt lt lt lt reftime time lt lt endl fileout close 182 numerics openmp reduction cpp Timings on AMD FX 8350 Bulldozer eight core processor 8 0 T0 c timing ideal linear speedup 6 0 5 0 Speedup 4 0 3 0 2 0 1 0 Number of cores cores time seconds speed up factor 1 4 336 1 2 2 135 2 0309 3 1 451 2 9882 4 1 085 3 9963 5 0 908 4 7753 6 0 760 5 7052 7 0 665 6 5203 8 0 583 7 4373 The apparent super optimal speed up for two cores is due to inaccurate one core timing used as a reference OpenMP has some overhead don t expect speed ups for short tasks 183 26 Tips and tricks Here means some lines of code The triple dot has a C meaning it s an ellipsis e How to pause exe
39. oe Sees CTETUR ae up PMCID MONETE oras dorado RE E NE dE dE 23 Physics example Spin 1 2 Heisenberg chain 23 1 Add numbers to file names a a a 24 Lambda functions expressions 5 openmp parallel programming 6 Tips and tricks 7 Some more C in the net N N N 28 Farewell words for C numerical programmers 140 141 143 145 150 158 163 169 174 175 180 184 186 187 1 Course itinerary 8 Lectures lecturer Vesa Apaja YN212 Nanoscience Center 2nd floor 4 Demos and 2 groups Group 1 starts on Nov 5 Group 2 starts on Nov 3 held in Physics Dep Computer class FL349 Demos are supervised by Sampsa Kiiskinen warm up demo 1 solved during exercise Requirements At least half of the exercises of the demos 2 3 4 5 and 6 plus one programming task No exam You may use a machine of your choice you need a fairly recent C compiler GSL and armadillo libraries the Physics Department server calc phys jyu fi Log in using your JYUNET login name and password In calc window type module add gcc This loads a recent version of gcc along with the C compiler g The system compiler is outdated thanks to the concervative dino age is safer package policy of RHEL 1 1 A Brief History of C e 70 s and 80 s FORTRAN 77 annoyed users with its fixed size tables Pascal became popular Borland introduced a fast compiler Dennis Ritchie creates C language with dynamic
40. of amp amp is widely considered unfortunate but live with it Scott Meyers calls T amp amp a universal reference some call it forwarding reference To make sense out of a thing like amp amp amp an rvalue reference to a reference one has a rule Reference Collapsing rule amp always wins Tis T T amp is T amp T amp amp 8 is T amp T amp amp is T amp T amp amp amp amp is TRE These rules seem odd but they work Return to the wrapper task and see how the argument types are cast wrapper 1 0 int amp y int amp amp z AN GAMA cU f std forward lt int amp amp gt 1 0 std forward lt int amp amp amp y std forward lt int amp amp amp amp z will reference collapse to f int amp amp int amp y int amp amp z rvalue reference reference rvalue reference This is the way we wanted f s arguments to be Task accomplished 47 8 5 The restrict keyword in C read on spare time You probably encounter this subjet in C coded scientific programs Wikipedia In C or C as mandated by the strict aliasing rule pointer arguments in a function are assumed to not alias if they point to fundamentally different types except for char and void which may alias to any other type Some compilers allow the strict aliasing rule to be turned off so that any pointer argument may alias any other pointer arguments In this case the compiler must assume that any accesses th
41. safe printf discussion Boost has worked on the issue quiet a while see on choices to be made in print I suggest that as soon as you define a class think of how the objects should be printed Write a member function to do that or overload lt lt to do the job The former way looks like this class Thing mytype data some data type can be self defined as here public print some code how to print data E int main Thing blob input data create and initialize a Thing blob print output the Thing You can improve this example easily to use arbitrary stream for output and hide the formatting to the member function Same for reading as blob read By the way printf is an example of a function that can have a variable number of arguments this is known as variadic function C 11 introduced a template that can have any number of template arguments variadic template O ao ss Here the is literally in the code it s known as the ellipsis This can be quite handy and I m just learning to appreciate it It s also template meta programming and that makes us all feel groovy Note The C Standard Library has already a nice variadic template std tuple 125 Example 44 advanced variadic_template cpp Variadic template to sum squares of numbers Any number of arguments to sum of squares include lt iostream gt include lt cmath gt using std c
42. tables e From 1979 on Bjarne Stroustrup creates C as a sort of improved C but is it The jury is still out The name C was coined in 1983 e Characteristic to C strong typing of variables supposedly reduces programming mistakes Object oriented programming possible though not forced data belonging together is collected together to an object that can be moved around as one big chunk Compiled language for speed e Early 90 s No standard compiler manufacturers make their own decisions e 1998 Horray C Standard An extensive standard library e 2011 C 11 Standard working name was C 0x Easier ways to do very common tasks C manuals written before about 2010 are better left to collect dust C is a multi purpose language and was never intended to specialize on numerics You may also come across the fact that old languages have the burden of legacy code Backwards compatibility and committee decision making tends to make old languages contain everything but the kitchen sink Maybe we the numerical people would actually need that sink 1 2 Future of C This timeline is given by the ISO C committee TS stands for Technical Specification TR for Technical Report C 98 C 03 C 11 C 14 C 17 major TC bug fixes only major minor major l a l 98 99 00 01 02 03 04 05 06 07 08 09 10 11 HH de File System TS J ne Library TR aka Ts b Fundamenta
43. them is a big thing in C 11 Moving is a Robin Hood operation poor std move rich pilfers rich from it s resources and hands them to poor Way more effective than copying money or gold The std move rich does not actually move rich it makes it movable Many times the compiler cannot tell whether an object is movable so with std move you can tell that an object has resources that can be pilfered In other words std move is for turning lvalues to rvalues see section 8 3 so that you can call the move constructor The move constructor is the Robin Hood code A move constructor and a move assignment look like this X X X amp amp other 01 10 move constructor X amp X operator X amp amp other C 11 move assignment operator You can be even more dramatic If objects of a class should never be copied you can forbid copying by deleting the copy constructor and the copy assignment only in C 11 class TM T const T amp delete forbid copying T amp operator const T amp delete The next lengthy example tries to elucidate situations when an object is copied and when moved For transparency the copy and move constructors as well as the copy and move assignments print out a message so you can see which is invoked 95 Example 14 advanced_move_constructor cpp Copying or moving include lt iostream gt include lt vector gt class T PUDI
44. till y begin jl y end 2 2 0 math y 1 0 2 E es Jf menos GO s 0 OO il 2 CU To create a std vector and set values the cleanest ways are Std vector lt douple gt yo 1005 2 7 maths y 1 0 1 05 13051205170 std vector lt double gt y 1 0 1 07 1 0 1 0 1 0 C 11 universal initialization C 11 lets you initialize almost anything with the universal initialization C 11 has also a range for or range based for loop std vector lt double gt y 100 for auto amp elem y elem 1 0 notice the amp use reference to elements The loop goes goes through all elements of y without you worrying about how many there are Be careful with the ambersand amp for auto elem y elem 1 0 this does nothing at all for auto elem y cout lt lt elem lt lt this works but doesn t try to change y 6 Can you do algebra with the length of a std vector x resize 0 sets the length to zero so x size is 0 Ok but x sizeO 1 is 18446744073709551615 not 1 Lengths are unsigned integers and can t store negative numbers 52 9 1 2 Storing objects into std vector If you know C well check out numerics better_vector_of_class_objects cpp Example 13 numerics vector of class objects cpp std vector storing objects g std c 11 vector of class objects cpp include lt iostream gt include lt vector gt include l
45. 00 0 0 01 5 cout lt lt final point lt lt x lt lt endl use method do_step to follow the solution step by step ofstream out ode dat const double dt 0 01 for double t 0 0 5 t 100 0 t dt stepper do_step harmonic_oscillator x t dt coupespecc Ol end oupsecp eE lt lt Mende out close 120 Example 43 numerics boost_ode_adaptive cpp Adaptive step size solution using boost numeric odeint ES SEACE po C OM retro CO Adaptive integration using Boost numeric odeint compile g std c 11 boost ode simple cpp include lt iostream gt include lt iomanip gt include lt vector gt include boost numeric odeint hpp gt using state type std vector lt double gt using stepper_type boost numeric odeint runge kutta cash karp54 state type namespace my using namespace std void system const state type amp y state type amp dydt const double t j void output double t double y double exact cout lt lt fixed lt lt setprecision 16 static bool first true ii rirse 1 cout lt lt setw 20 lt lt t lt lt setw 20 lt lt Boost solution cout lt lt setw 20 lt lt exact solution lt lt setw 20 lt lt error n first false dydt 0 txy 0 i cout lt lt setw 20 lt lt t lt lt setw 20 lt lt y lt lt setw 20 lt lt exact lt lt setw 20 lt lt y exact lt lt endl 121 int main 1
46. 14 10 30 Version 4 500 Singapore Sling can utilize MKL and OpenBlas to mention a few e MTLA Simunova a small c software company in Dresden speed tests e I T4 4 Sweden version 4 2 0 came out 2010 long break version 4 3 1 came out 6 7 2013 e Blaze 20 6 2014 version 2 1 A project initiated by Klaus Iglberger main developers in Erlangen Germany A link to a fine article on Smart Expression Templates Sandia A huge collection of methods actively maintained version 11 12 1 Oct 24th 2014 Eigen Eigenvalue problem specialist version 3 2 2 came out 4 8 2014 e ETL The Iterative Eigensolver Template Library is another eigenvalue specialist Uses Boost LAPACK BLAS ATLAS and possibly Blitz The web page was updated 2004 so I wonder if there is any development Armadillo MTL4 and Blaze have a simple interface close to Matlab Octave They use generic programming and boldly apply the expression template technique to optimize e g matrix and vector operations In this group Blaze is a newcomer and blazingly fast It works parallel in shared memory machines such as your laptop Header only libraries only need to be copied to your machine Include the library path in compilation is all you need to do 128 18 1 Armadillo examples Example 45 numerics arma_matrix_multi cpp Armadillo Matrix product source Armadillo web page include lt iostream gt include lt armadillo gt
47. 3 Function objects functors 9 scs seda km a RO RO RO eS Se CR RUE ex OWOX AAA 13 4 Four ways to pass a function to a function oo a a 13 5 Cache dalal 4 sacres acy em ak Be SOEUR S Y a aed GR E BS oe de UE edm XD REOR XR Re Ge deos 13 6 Use emplace back instead of push backK a 13 7 If available use the methods of containers rather than algorithms 13 8 Expression templates 2 22 ee 14 The best random number generator 15 C 11 generation of pseudo random numbers 15 1 std bind to simplify usage ec poo e om UR o ES FE oe e FOR gu X EE eed o Rt EUROS 16 Boost library 16 1 Boost ordinary differential equations ODE 17 Formatted output 17 1 Formatted output using include iomanip e 17 2 printf type safety and variadic functions and templates 020000002 eee 18 Linear algebra which library to use 18 1 Armadillo examples a 18 2 Blaze example 4x 9 ma koe bU e RR mex o3 door e Bee ee d ee ne baal es Res X UR UH eed 19 Calling C or fortran from C 20 Fixed size arrays in C plain array and std array 21 Exception handling with throw and catch 111 112 114 118 119 123 123 125 127 129 132 134 136 138 22 Gnu Scientific Library GSL 22 1 E 43 bei bales 3 9 SASS R AS ES eee amp She SOGCRG OR ee SEEPS wee Med 22 2 GSL Fast Fourier Transform FET uu oue ox RU o UR ee ee Od eee SE Se
48. B2 a o erm als tells r n1 is a reference to an integer which is n1 Nothing wrong with that but n1 is a local variable so it s no longer alive after you exit the function get number After returning to main the reference r_n1 points to where ni used to be What makes dangling references perilous is that sometimes the reference happens to dangle over the right spot 42 Example 10 basic dangling_reference3 cpp A dangling reference that the g compiler will notice Unsafe code example g warns at complile time include lt iostream gt int amp get number 1 ime amd 100 return n1 j int main int number get number std cout lt lt number is lt lt number lt lt n You don t necessarily need a function call to get a dangling reference Example basic dangling reference2 cpp shows what happens if you store the reference to a vector element and then resize the vector There is no guarantee the stored reference is pointing to a valid location any more valgrind a out often detects such memory problems My own valgrind user manual Any message from valgrind q a out means your code has a bug 43 8 3 lvalue and rvalue and rvalue references C 11 has actually five kinds of references and I introduce two them so that you can read C 11 code e lvalue reference Marked with a single ambersand amp for example int amp a All references I ve b
49. FYSA120 C Numerical Programming Vesa Apaja November 23 2015 Contents 1 Course itinerary 6 a BMG Eth Geek Ghia ees Gees ye Gee ects ee ee ee ee 7 Peat ee ene ee Ge eee mee Be eg Gr es Beet ng ees Won Ge ae He ee ace es 8 1 3 How transferable is C 0 0 0 aaa a ll s sss 9 1 4 About this document o a aaa a sss 10 1 5 Why C why not Java fortran Lisp Haskell French 7 2 2 llle 11 LA bf todd ao domo 4o e ae eh oe Ga Fox e eee eee d S4 E RR Nodo vun s 13 IT Oe im the Net do Pe ea ad we e Swe Ge Oo qum A wee do da a 14 1 8 C in Matlab or Octava 14 1 9 Cor C in Python A 14 2 A really Brief introduction to C 15 2 1 Meaning of include lt iostream gt 2 222 2 lll 4 e s s s s 16 2 9 SCOPE sarka ug sos ae ede eee nds RR ade keene ee had eRe REOR a a 20 2 3 Simple file operations ee 21 3 Class 22 4 Member function with const or noexcept after it 23 5 Code rotting and deprecated functions 23 5 1 Delegating a Constructor ead Dn Spare HIE sese a ARA oe ee RORCR GR 30 5 2 Own data structures pairs of numbers 2 les 31 6 Templates generic instructions and algorithms 33 7 C libraries 35 7 1 C Template libraries just alist o e 35 7 2 C Standard Library ee 36 8 C reference variables 38 8 1 Why would a reference be safer than a pointer 222A 39 8 2 Can a reference be unsafe oaa a a els sss os ss 42
50. Iter first Iter last Func f while first last f first irst return f M Oslo return moy The third argument f can be a class as long as f first is defined This hints that the name of a class can sometimes be used the same way as a function they are called function objects or functors for_each argument third is class Func f and the return value is the same f The pass by value feature Notice also how the third argument is passed by value as class Func This means the argument object Func has a one way ticket it does not return as an argument But the function does return it as a return value So in as argument out as a return value Example If you send in a function object and change something in that object you have to use the return value object The example numerics foreach_functor2 cpp shows the principle it computes the cosines of elements done by the function object and collect their sum as a data member in the object The gcc version 4 9 1 has the file include c 4 9 1 parallel for_each h Here the parallel gives away that for_each can be parallelized All elements are pushed separately though the same function so why not do it in parallel 91 The page Draft of Technical Specifaction tells just how draft the parallel part is still in 2015 Working Draft Technical Specification for C Extensions for Parallelism Note this is an early draft It s known to be incom
51. T2 amp y TOV EZ wrapper const T1 amp x T2 amp y T3 amp z wrapper Tix const 128 y T3 amp z wrapper Tig 12 y const T3 amp z wrapper const Ti amp x const T2 amp y T3 amp z wrapper const T1 amp x T2 y const T3 amp z wrapper Tig x const 12 y const 13 2 wrapper const Ti amp x const T2 y const T3 amp z This gets impossible with increasing number of arguments for n arguments you would need 2 overloads And it only gets worse once you add the possibility of rvalue references We could want f int x int amp y int amp amp z so to perfectly forward the int amp amp you would need more overloads C 11 solves this task with the type cast std forward solution template lt typename Ti typename T2 typename T3 gt void wrapper Ti amp amp x T2 amp amp y T3 amp amp z f std forward lt T1 gt x std forward T2 wv std forward lt T2 gt z 46 As Scott Myers says std move moves nothing and std forward forwards nothing They do a type cast std move is a single minded Robin Hood cast poor std move rich rich becomes an rvalue movable property ready to be stolen std move does this always poor std move another poor really always a single minded cast How did std forward do the wrapper task First T amp amp is not always an rvalue reference There are cases where one instead does type deduction This recycling
52. TE int count std vector lt double gt x TO default s constructor T int count std vector lt double gt x_ count count_ x x_ constructor T noexcept default destructor copy assignment auto operator T amp rhs amp gt T amp x rhs x count rhs count std cout lt lt copy assignment called n copy constructor const Me rha x rhs x Count rns count Sidi sc out cau copy cOn site tor lle di mue move assignment auto operator T amp amp rhs amp noexcept gt T amp x std move rhs x count std move rhs count ls count Us rhs x resizes x to 0 length 56 std cout lt lt move assignment called n move constructor T T amp amp rhs noexcept x std move rhs x count std move rhs count stadir Cout lt move Constructor calles M mE We int main Toxt137 0172 2115 Std scouts x2 NS D IUE il 5 258 SE ee xs std Cout a TE SIT a maxx Dr 6 7 same as To std cout lt lt T x4 std move x1 An als x4 std move x1 std cout lt lt T x5 std move x2 n Tesi Sstdh mower x2 i 57 9 3 std valarray class As you may have noticed std vector is far from ideal for numerics you can hardly think of anything more horrendous std vector is not a mathematical vector it s a list that can expand and shrink This led one to introduce to C the class std valarray not
53. ach the cord and tie it to anything else A C reference cannot be detached from its variable int amp c 39 Amuse yourself with these almost relevant examples Pass by value C or C Write the numbers 5 and 10 on a piece of paper show it to you college and tell him her to copy the numbers to his her own piece of paper in reverse order Swap failed your paper still holds 5 and 10 not 10 and 5 Pass by reference C Write the numbers 5 and 10 on a piece of paper hand it to you college without letting go of the paper Tell him her to swap the numbers on your paper then pull the piece of paper back Swap achieved Robustness Your college immediately sees if your hand is empty or there are no numbers on the paper Pass by pointer C tai C Write the numbers 5 and 10 on a piece of paper and tell your college the paper is in a specific drawer Tell him her to swap the numbers written on the paper and to put it back to the drawer Swap achieved Robustness your college picks a wrong piece of paper from the drawer and wonders why you want to swap two telephone numbers 40 Remarks 1 If you have many colleges and or 10000 numbers copying them around may take a lot of time and paper Prefer passing references it s fast and economical 2 For the swap to function properly you can t pass numbers 5 and 10 by value that is copies of the values Two working ways e C style swap handles references function
54. al distribution like this the name normal dist is my own std normal_distribution lt double gt normal dist 0 1 112 4 Give the generator a seed only once suffix u means unsigned gener seed 4835267u same sequence every time you run the code or pick the seed from system clock gener seed static cast uint fast32 t std time 0 at least 32 bits different sequence every time you run the code if time 0 changed or rely on std ramdom device gener seed std random_device This is not without problems see a discussion in cpps random_device html 5 EITHER i bind the generator and the distribution together auto normal_random std bind normal_dist gener do this once here auto is actually static std function lt double void gt and use the combination like this double random normal_random OR ii get the random number directly from a call to distribution generator double random normal_dist gener Here the random number is random 113 15 1 std bind to simplify usage Before continuing with random numbers lets look what auto normal_random std bind normal_dist gener did It created a function normal_random that means call normal_dist with argument gener Example 40 basic bind example cpp std bind is versatile How to use bind to a change the order of arg
55. ate integrate with typename T for MyFunc and replaces the function f O with MyFunc operator as an inline function This makes the code very fast 102 3 C 11 Pass the function as a function object of the class template std function Recommended method double f double x return 1 0 exp 10 0 x sample integrand double integrate const std function lt double double gt amp f double a double b integrate function f x from a to b use like integrate f 1 0 2 0 integrate f x If the function has been declared earlier you can let the compiler deduce the types now double double double integrate const std function lt decltype f gt amp f double a double b 4 Pass the function pointer as a template parameter double f double x return 1 0 exp 10 0 x sample integrand template lt double Tfunc double gt double integrate double a double b use function Tfunc as usual j use like integrate lt f gt 1 0 2 0 integrate function f Examples of each style is in the file numerics function to function speed test cpp The speed depends on whether the compiler inlines a function pointer recent ones do IMHO methods 2 and 3 looks best and conform with the C general goal to avoid naked pointers 103 13 5 Cache data Don t recompute the same data over and over again Use cached values 13 6 Use emplace_back instead of push_back Recom
56. ciate it if someone has done a good job overloading operators As an example without overloading adding two complex numbers c1 and c2 would read something like this cs vadd cd cas With overloaded operation it reads es Ci 23 The compiler does change the to function call but it s out of sight The code is readable and just like math If you overload an operator make sure it works as expected This is related to The law of least astonishment The program should behave in a way that least astonishes the user Steve Oualline How Not To Program in C Here is a story of a code that didn t obey the law of least astonishment I was surfing the net one day for examples on how operators can be overloaded and came across with one that overloaded the operator for complex numbers Upon testing I found that the sum c3 c1 c2 really is computed correctly I thought the implementation was ok and used it in my code Then I was astonished to get wrong results No not because c3 was computed wrong no hou It was because I didn t come to think that the operation c3 c1 c2 was changing also the value of c1 It did and in the code that followed c1 had a wrong value A classic piece of bad code is this attempt to use a macro for squaring BAD CODE Never use define this is just for educational purposes define SQUARE x x x x iiu me int b SQUARE a you would think thi
57. cmath gt using namespace std class TakeCos public void operator double amp x x cos x function object Js int main vectorcedouble lios cout vector x US for auto ex cout lt lt e lt NS cout endl for each x begin x end TakeCos COMES yector colr m for auto e x cout lt e lt lt cout lt lt endl return 0 If the task is this simple or it s supposedly used only here it s more convenient to use a lambda function introduced later in chapter 24 One example of a lambda function was already in numerics mystatistics cpp 101 13 4 Four ways to pass a function to a function 1 One way is to pass a function pointer double f double x return 1 0 exp 10 0 x sample integrand double integrate double f double double a double b integrate function f x from a to b i use like integrate f 1 0 2 0 integrate f x 2 Wrap the function in a class and make it a function object Write integrate as a template class MyFunc jo ILILE S double operator double x return 1 0 exp 10 0 x define a function object y template lt class T gt double integrate T f double a double b use function f as usual use like integrate Myfunc 1 0 2 0 integrate MyFunc class function object Notice how this code does not have any explicit MyFunc class objects The compiler instantiates the templ
58. cution cout lt lt PRESS ENTER TO CONTINUENn cin ignore e Infinite loops can be typed as while true or more to the point as Hdefine ever for ever 1 break to get out SE I heartily recommend stackoverflow com discussions This is THE tip As what comes to having fun with C this one is absolutely priceless what is the worst real world macros pre processor abuse youve ever com e Linux Long compiler messages can be piped to more or less Notice the amp after the pipe character g std c 11 program cpp amp more e Linux It s tedious to type g std c 11 all the time Write a makefile use an IDE or set an alias alias gtt g std ct 11 bash put this line to SHOMEJ bashrc and type SHOME bashrc alias g gt std c 11 csh or tcsh put this line to HOME cshrc and type source HOME cshrc and you only need g program cpp 184 e std numeric limits from lt limits gt tells what traits a type has hence they are type_traits For example std numeric_limits lt int gt max gives 2147483647 If you have in mind a class that can use an optimized algorithm and another related one that cannot you can give the objects a trait Based on that trait the code can automatically decide if the optimized or the default algoritm is applied See advanced traits cpp 185 27 Some more C in the net e Genetic algorithms
59. de lt map gt using namespace std void vector_out vector lt int gt v 1 copy v begin v end from ostream iterabor inb cout j 77 to cout lt lt endl int main 4 const int N 6 single fermion states const int Nfermions 3 number of fermions cout lt lt All lt lt Nfermions lt lt fermion states for N lt lt N lt lt endl vector lt int gt state for int i 0 i lt N i if i lt Nfermions 69 state push_back 1 else state push back 0 j sort state begin state end int is l do cout lt lt state lt lt i lt lt is vector_out state aL while next_permutation state begin state end coute found lt lt i l_ lt fermion states VH return 0 70 All 3 fermion states for N 6 state 4 1 is 000111 state 4 2 is 001011 state 4 3 is 001101 state 18 is 110010 state 4 19 is 110100 state 4 20is 111000 found 20 fermion states What would for example 110010 stand for Electrons can have spin up 1 or down so the states could be coded as occupations of up down pairs 110010 would mean T 00 7 0 meaning the 1st one body state has up and down electrons the 2nd is empty and the 3rd has a spin up electron Each state has only zeroes and ones so it would be very economical to encode them in binary 71 9 7 Header guards and namespace encapsulation This section demonstrates how to write a function for hom
60. dl using std vector wector cdouble y 2 8151 51675 138 3 IO mean gsl stats mean amp v 0 1 5 variance gsl stats variance amp v 0 1 5 cout lt lt mean lt lt mean lt lt endl cout lt lt variance lt lt variance lt lt endl Passing a std vector as a pointer amp v 0 is not pretty but does the job All that the GSL function needs is the start address and the length of the data elements of a std vector are stored in consecutive memory slots 142 22 2 GSL Fast Fourier Transform FFT I want to transform complex data with length 2 N Zso like this fft data direction direction 1 forward 1 backward Here direction is 1 for a Fourier transform and 1 for an inverse transform Notice that I want to keep this syntax no matter what library does the FFT The decision to use GSL s FFT is made in a include which loads the header given below The variable type of data is flexible and the template works at least for std vector std array fixed size not plain array std valarray armadillo vector and Blaze vector Armadillo vectors needs some special care The compiler option DARMA chooses the first option in ifdef ARMA Armadillo data has no method size use n_elem n data n_elem gt else n data size endif This solves the problem of missing size method but it s a hacky solution 143 Example 53 gsl fft hpp
61. e made statistics with C Standard Library algorithms and without This is just for demonstration there are much better statistical library routines than this Here we push numbers to a std vector and compute the statistical mean and standard deviation of the data in the function get_stats Let s start with a header Headers are for the compiler with some information to you Java doesn t need them Example 21 numerics mystatistics1 hpp First attempt as a header for get_stats ifndef MYSTATISTICS HPP define MYSTATISTICS HPP include lt vector gt void get_stats const std vector lt double gt amp double amp double amp Hendif The preprosessor directives ifndef MYSTATISTICS HPP define MYSTATISTICS HPP endif make up a header guard They make sure this piece of code is not processed more than once To use this header stored in the file numerics mystatistics hpp put in the beginning of the program the line include mystatistics hpp This line may well be in many program units hence the header guard ifndef stands for if not defined If MYSTATISTICS_HPP is not defined define MYSTATISTICS_HPP and process the rest The next time the compiler tries to include mystatistics hpp it already has the code processed and does nothing The file suffix hpp is one way to tell that this is a header file not to be processed unless include d In C one has the suffix h and it will also
62. e of the std vector element you just tell how it should be made Example 39 numerics emplace_parameter_arguments cpp emplace_back can construct with parameters ii emplace_back with parameter arguments if include lt iostream gt include lt cmath gt include lt vector gt using namespace std struct MyObj double value MyObj const double amp par double angle noexcept value parxsin angle E int main const int N 1000000 const double par 3 866 vector lt My0bj gt v for auto i 0 i lt N i 1 v emplace back par i M PI N not emplacing an object but forwarding parameters to constructor cout lt lt first 10 values n for auto it v begin it lt v begin 10 it cout lt lt it value lt lt COMA 106 13 7 If available use the methods of containers rather than algorithms The general idea is that a container specific algorithm can take advantage of the container properties and gain speed Sometimes generic algorithms are not available at all the list container has no random access iterator so std sort can t work Instead there is a sort method in list std list lt int gt mylist Hi sea ES mylist sort calls the sort member function 107 13 8 Expression templates This chapter gives you an idea how complicated it is to master C and why I never can What I hope you to learn from all of this is use good libraries What do expre
63. eate a stream either for writing ofstream or for reading ifstream o as out i as in and f as file To a stream we stuff things using the operator lt lt and read with the operator gt gt The console output cout console input cin work like streams Example 4 basic fileio cpp write to a file include lt iostream gt include lt fstream gt using namespace std int main ofstream myfile output myfile lt lt Text to file output lt lt endl myfile close return 0 This declared and opened a stream in the one line with ofstream ofstream myfile output The file output contains now the text Text to file output and myfile close myfile s close there may be many other close functions closes the file The object myfile of the class ofstream has a method member function close 21 3 Class Say you want an object to be the state of the system The state holds as variables the number of microstates N their occupations n1 n2 ny and the energy of the state E and things you didn t come to think of at this stage In addition the value of E should be set only by the method comp_energy A collection like this can be stored as a class A class is a collection of data and methods to manipulate the data A class is just an abstract type a bit like int only tells what an integer is You need an object a manifestation of the class
64. ecified after all a class must have some constructor to be useful at all The only thing a default constructor does is to reserve some space for an object nothing else C language specifies that If you give a constructor for a class the default constructor vanishes This is exactly what happened in the example a colourless Car construction was not there Objects can be copied Car buick cadillac 27 and buick picks the colour from the cadillac Wait a moment what function did the copying It was the copy constructor provided by the compiler automati cally 28 As you guessed there is a method called destructor which deletes an object It has the strange tilde sign Car cout lt lt Car demolished n Compilers create this automatically then it s called a default destructor You can also provide one yourself and include a clear message when an object is destroyed To summarize a class has a set of building methods add the prefix default if you let the compiler create a method constructor how an object is created when or if created copy constructor how an object is copied resulting two similar objects with different name copy assignment makes a copy from right to left move constructor tells how the information of an object is transfered to another the object essentially pilfers steals the resources of the other object leaving it in default statp move assignmen
65. een using so far e rvalue reference Marked with a double ambersand amp amp for example T amp amp But wait not all amp amp s are rvalue references read the next chapter There is also xvalue glvalue and prvalue I m cutting corners so for deeper understanding read for example c11 tutorial explaining the ever elusive lvalues and rvalues by Danny Kalev Two things you need to consider can the thing be moved and does it have an identity Expressions program statements that have a value have long or short lifetimes The long living expressions have an identity so that you can call them the short living expressions die soon and have no identity An lvalue is a non movable expression or one with identity It lives long or cannot be moved An rvalue is a movable expression or one without identity It may be about to die such as a temporary Evaluating x a c b the compiler can create a temporary to store the result of a c that s an rvalue Simplified rule lvalue is a variable rvalue is a temp With an rvalue reference you can prolong the rvalue s lifetime For example int amp amp temporarily_rich return 5000 number 5000 is an rvalue no identity temporarily poor amp amp std move temporarily_rich std cout lt lt poor lt lt std endl 44 8 4 The strange T amp amp and the perfect forwarding problem This line of thinking is not my idea it s been frequently used
66. em of being not part of the standard Ha so the whole library may cease to conform with the standard in the future The developers retire or more probably someone writes a better and faster library and your vocal colleges push you to using that instead There is a way to overcome these difficulties It s not a cure but a workaround With typedef or using you can give a name to your own data structure Write the code so that the decision about a data type is in one place in case you change your mind typedef existing type new type name Example typedef double v data call double by the name v data typedef std vector v data stl vector vector of v data is called stl vector typedef std vector stl vector gt stl matrix vector of vectors is a matrix If you change your mind and want to store integers to st1 vector and stl matrix just change the first line to typedef int v data Another way to do the same using new type name existing type using v_data double call double by the name v data using stl_vector std vector lt v_data gt vector of v data is called stl vector using stl matrix std vector stl vector vector of vectors is a matrix Example you find that using library smack fictitious that has a nice smack vector you could build upon this alias using my_vector smack vector If you decide to replace smack vector with a close relatives
67. equations details in 22 3 Then use your great mathematical intuition and tell if the solution has kinks or is otherwise badly behaving Choose the integration algorithm the stepper that solves the next point Fixed step size routines Simple and fast accuracy of the solution is your responsibility Adaptive step size routines Try to reach an accuracy goal absolute and relative error limit Boost has a seasoned library for solving ODE S as does GSL Boost calls the current values the state which contains x t x t Hence the type state type Example 42 numerics boost_ode_simple cpp Fixed step size solution using boost numeric odeint include lt iostream gt include boost numeric odeint hpp gt include lt fstream gt using namespace std using state type std vector lt double gt using stepper type boost numeric odeint runge_kutta4 lt state_type gt ES OMS 369 G5 cc ECC ME amc Cs SoLL GO MEC OTITNIK GI 2 GE Ws y E x t gt gam yCt notation E ol Y A 7 Ce cese OT 1 Ge elche lal void harmonic_oscillator const state_type amp x state_type amp dxdt const double t const double gam 0 15 ba O x dxdt 1 x 0 gam x 1 119 int main stepper_type stepper state type x 1 0 2 04 7 amitial values x 0 1 x7 0 2 integrate all the way to final time Eimer consti MwA y memomile_ oscullero E 0 0 1
68. er specifically it has no size method The least I want is to add a third argument to the clean call fft data direction to tell the data size I have chosen to use the preprocessor and define ARMA if Dm using Armadillo vectors Test code numerics gsl fft arma main cpp is shown below one using std vector is in numerics gsl fft main cpp 145 Example 54 numerics gsl_fft_arma_main cpp FFT test using Armadillo modified from GSL manual C version Compile the below is printed wrong in latex if ss seed gel fx anma maria ay md esusig lilos tinclude lt iostream gt include lt cmath gt include lt complex gt include lt fstream gt define ARMA tinclude gsl_fft hpp include lt iomanip gt include lt armadillo gt using namespace std using namespace arma using my_GSL_FFT fft short name for data type using vectype Col lt complex lt double gt gt overload lt lt to print complex numbers in vectype as I like them here std ostream amp operator lt lt std ostream amp os const vectype amp v using namespace std os lt lt scientific lt lt setprecision 8 some I O manipulation int Sie for auto x v os lt lt setw 5 lt lt i lt lt lt setw 20 lt lt x real lt lt lt lt setw 20 lt lt x imag endi return os 146 int main void const int n 128 const int ni 10 Keep less than n vectype data n ofstream myfi
69. ermine the sequences are not random Ok why not run one mill with a mother seed to give you seeds for many rng s That too has not been tested to give sufficiently random results Algorithmic generation of pseudo random numbers is a tricky thing 94 12 6 C Standard Library algorithms stateful objects and std ref Algorithms std for_each and std generate are not useless in context of stateful objects There is a simple remedy to the copy problem For stateful objects use a reference wrapper std vector v 100 std generate v begin v end std ref gen gen is stateful object std ref is a helper function to generate a std reference wrapper meaning see std ref or std reference_wrapper Std ref 10 is std reference_wrapper lt int gt Example 35 numerics generate_id_vector cpp Fill a std vector with unique id numbers from an id generator include lt iostream gt include lt algorithm gt include lt functional gt include lt vector gt using namespace std class IdGen int id 7 4 objects stare joxtilelstc 3 IdGen id 0 cout lt lt constructed an IdGen n constructor int operator void return id functor lidgen int main vector lt int gt vi 20 20 generate v begin v end generate w begin w end std cout lt lt v n std ref idgen try these without std std ref idgen tor tanto rsy conte
70. es not use references so I may talk about pointers as well Just because you may one day find yourselves on the drivers seat of a car with manual transmission GSL is written in C so that is at least one point where pointers are needed Usually I put my faith in libraries acknowledging they are also written by humans Specifically I hope this course i gives you self confidence to write your own C program that solves a numerical task using a chosen library ii helps you read C programs and sometimes even understand how they function on a good day even C programs may look familiar It s easier to make a working code faster than to make a fast code to work A code optimized for speed can be a very unfriendly beast to debug Most codes have portions that need not to be fast at all The English in these lecture notes is not perfect but I remind you how Yoda the Lucasfilm trademark character the master of The Force after over 874 years says Your father he is but defeat him you must This is BTW called OSV or Object Subject Verb word order C you learn 10 1 5 Why C why not Java fortran Lisp Haskell French Depends on what you want to do Maybe French is the best language if you want to read old mathematical manuscripts or write poems The learning curve in C is steep and it was designed to do all but numerics Fortran was designed for numerics making it more compact We would probably never
71. from library crunge you edit these lines and not touching code using my vector Be cautious crunge vector and smack vector may be incompatible and support different operations Rather than inventing a name my_vector it s better to use your own namespace my and hide the definitions there Then my vector is the thing to use We ll come to this later 60 9 5 Stream iterators read on spare time Stream is an object representing an I O input output channel Streams are considered one of the best qualities of C by C programmers and one of the worst by C programmers A stream iterator is something that goes through a stream So a stream iterator is a boat Think of the following task Write a program that reads an arbitrary string of characters from the keyboard Sounds simple but in many programming languages this is very lengthy to do safely Remember user may hit any key and and keep on hitting for a month Below is a C suggestion Example 15 advanced streamiterl cpp Read character to a std vector include lt iostream gt include lt vector gt include lt algorithm gt include lt iterator gt using namespace std int main vector lt string gt s input from standard input cin copy istream_iterator lt string gt cin from istream_iterator lt string gt end back_inserter s to output to standard output cout copy s begin s end from ostream iterator string
72. have heard about Linus Torwalds had he written the linux operating system in fortran Relax if you know C you can more easily learn Java and fortran I m not talking about the ancient FORTRAN 77 phew but fortran 2008 fortran 2015 coming They all have similar structures like classes This excludes functional languages such as Haskell because their way of thinking is quite different from the imperative languages C or C To mention a few differences Java has no operator overloading appreciated by numerical programmers Speed is crucial C programs can challenge and sometimes beat fortran I suspect that a lousy fortran programmer can write faster code than a lousy C programmer Speed comparisions age quickly but you get hints of the situation P If you plan to run the program a few times pay attention to time spent on programming If you plan to run the program hundreds or millions of times pay attention to time spent by the program Numerics is the latter It s number crunching squeezing the last drops of the computer juice http javarevisited blogspot fi 2011 08 why java does not support operator html engineering linkedin com 49 linkedin coding competitions graph coloring haskell c and java 11 12 1 6 Easy tasks Yilun Zhang Learning data science kB Zhiraz Bet unnar Oel ike Merchant Hello world in Python print Hello world Java l public class HelloWorld public static void main S
73. hows how a std cout object eats output manipulations width of the field 10 the output form fixed and the number of decimals 6 cout lt lt fixed lt lt setprecision 6 6 decimals cout lt lt setw 10 lt lt x lt lt setw 10 lt lt y lt lt setw 10 lt lt z lt lt end1 field width is 10 cout lt lt setw 10 lt lt x lt lt setw 10 lt lt y lt lt setw 10 lt lt z lt lt enal This is awful but at least the same formatting works also with writing to a file it s a stream just as cout ofstream output data out output lt lt fixed lt lt setprecision 6 output lt lt setw 10 lt lt x lt lt setw 10 lt lt y lt lt setw 10 lt lt z lt lt n data out 1 542234 12 423400 0 121300 13 000000 4 234000 1 000000 123 See numerics output_formatting cpp for more examples Beware that after formatting some setting are still on fixed while some are immediately forgotten setw C has the famous printf function which is quite readable but not considered a good C practise I suggest that if you are already good with printf use it Pros of a stream object is type safety the compiler can always tell if a data type cannot be sensibly dealt with and flexibility it works the same way on screen and on file 124 17 2 printf type safety and variadic functions and templates There are several attempts to code a type safe printf in C for example by Andrei Alexandrescu see type
74. in end vector lt gt v v begin points to the begin of the container v v end points to one step past the end of the container v v begin v end Why one step past the end It was chosen like that because of cleaner loops for auto it v begin it lt v end d e algorithms are frequently needed ways to process data sort find min_element max_element reverse 36 Example 8 basic std swap cpp Swap two numbers using std swap include lt iostream gt include lt utility gt int main using namespace std double a 5 double b 10 cout lt lt before lt lt a lt lt lt lt b lt lt endl swap a b cout lt lt after lt lt a lt lt lt lt b lt lt endl return 0 before 5 10 after 10 5 Someone familiar with C might be concerned about what values a and b have after the call to swap How are the arguments passed to the function and what comes back 37 8 C reference variables Three ways to pass data to functions e Pass by value void dothis int function prototype c 15 dothis c work with number 15 not with c cannot change contents of c e C Pass by reference void dothat int amp function prototype c 15 dothat c work with c can change contents of c e Pass a pointer void doodd int function prototype GS line doodd amp c works w
75. in teaching C 11 by the smart ones like Scott Meyers and Eli Bendersky I give here a try The task Write a wrapper function wrapper x y z that perfectly forwards here was the keyword all arguments to function f x y z Simple ha It s just a simple template try one template lt typename T1 typename T1 typename T3 gt void wrapper Ti x T2 y 13 y Ee ec ee D Then you realize that this will always call by value so the arguments x y z are copied over This is no good if you can have declared f int amp x int amp y int z Clever you you rewrite this version try two template lt typename T1 typename T1 typename T3 gt void wrapper Ti amp x T2 amp y T3 amp y Elx yz Then you realize that this won t do if you can call f with rvalues things that are movable or have no identity like f 1 0 3 0 The compiler complains about an invalid initialization of non const reference from an rvalue It cannot make a reference required by T1 amp x that points to number 1 0 because the number has no identity to refer to You could cure the problem by adding a const to an argument To any argument because f could have argument types 1 0 int amp y int amp z or f int amp x 3 0 int amp z or 45 At this stage you smell something burning You would have to write overloaded wrapper functions for all possible combinations in this case wrapper TRE Se
76. in_mystatistics cpp Main program to test mystatistics cpp and mystatistics_c 11 cpp compile g std c 11 main_mystatistics cpp mystatistics_c 11 cpp f MY SUEMGILS SIPICIDIDMSUI SIS jaye include lt iostream gt include lt vector gt include mystatistics hpp zxmiesliudeaemwsiterbhisibsies ester ede p include mirala pps int main i 1 hi using namespace std vectorcedouble x 5 595 3 56 Ser ODA double avei ave2 sigmai sigma2 cout rie myutil vector_out x cout lt lt std algoritm version vs C 11 version lt lt endl mystat get stats x avei sigmal std version mystat Cii get stats x ave2 sigma2 C 11 version cout lt lt average lt lt avel lt lt vs lt lt ave2 lt lt endl cout lt lt standard deviation lt lt sigmal lt lt vs lt lt sigma2 lt lt endl std algoritm version vs C 11 version average 7 53333 vs 7 53333 standard deviation 2 02452 vs 2 02452 TT 9 8 std complex complex numbers and arithmetics All common complex operations multiplication division addition real part real etc are there already Example 26 numerics complex_ex cpp Basic operations with complex numbers include lt iostream gt include lt complex gt int main 1 using namespace std complex lt double gt c c2 ci complex lt double gt 1 5 2 2 c2 complex lt double gt 1 0 3 3 cout lt lt c1
77. include lt iostream gt include lt vector gt include lt algorithm gt include lt cmath gt include lt complex gt using namespace std int main typedef complex lt double gt cmplx vector cmpl yv enpl Tilt cmpixi 220271 empire for aute xv cout lt lt x e cout lt lt original data lt lt endl take complex conjugate of elements for each begin end cmplx amp c c conj c for aubo xov Contas we cout lt lt complex conjugate lt lt endl Tulostus 11 12 22 21 33 32 11 12 22 21 93 92 Note The elements are passed by reference to complex conjugation 177 Example 65 advanced lambda4 cpp Calculate the product of std vector elements include lt iostream gt include lt vector gt include lt algorithm gt int main using namespace std vector lt double gt v 1 0 2 0 3 0 double prod 1 0 for_each v begin w eud amp prod double x prod x cout lt lt product lt lt prod lt lt endl product 6 Note prod is captured by reference 178 Example 66 advanced lambda5 cpp Chop all elements below some limit to zero include lt iostream gt include lt vector gt include lt algorithm gt using namespace std typedef vector lt double gt vec void rand_vector vec amp v generate v begin v end return rand double RAND_MAX int vector_chop vec amp v con
78. into 3b eed do cout lt lt state lt lt i lt lt is vector_out state srt E while next_permutation state begin state end cout lt lt N lt lt factorial N lt lt endl return 0 state 4 1 is 01234 state 2is 0124 3 state 3 is 01324 There are N permutations computed for checking in the recursive function factorial Pl Now next permutation state begin state end does all the work all else is supporting code It generates the next permutation of elements and returns true until it no longer finds a new permutation and returns false The loop do while next_permutation state begin state end goes on until the test while true changes to while false A recursive function calls itself The implementation is poor there is no test that the result fits an int 68 Physics example Generate all many body states with a fixed number of fermions Each spin state can hold 0 or 1 fermions Pauli rule If you have 4 spin states and 2 fermions you can have states occupied as 0011 0101 0110 1001 1010 1100 in total 6 different many body states Example 20 numerics algo fermion states cpp Fermion states by permutation Finding fermion basis states using C Standard Library algorithm next_permutation include lt iostream gt include lt vector gt include lt algorithm gt include lt iterator gt inclu
79. is example used the initialization list Sum imt 1 aut 9 ala a2 J Les which means that the value of i is given to al and the value of j is given to a2 30 5 2 Own data structures pairs of numbers Warning The example below is BAD no need to re invent the wheel Example 6 numerics class_pairs_of numbers cpp A class for pairs of numbers include lt iostream gt include lt cmath gt using namespace std Class for pairs of numbers class TwoDouble publica double x y y double distance const TwoDouble amp const TwoDouble amp int main TwoDouble pointi point2 pointi x 1 0 Class data member x can be accessed directly It was public pointi ya pointa z 20 point2 y 1 0 cout lt lt distance lt lt distance point1 point2 lt lt end1 return 0 j distance between two points x y double distance const TwoDouble amp ci const TwoDouble amp c2 return sqrt pow ci x c2 x 2 pow ci y c2 y 2 j Notice how the name chosen for the function distance is conveniently descriptive but potentially dangerous The namespace std is used in it s entirety and who knows if there already is a function called distance doing a completely wrong thing Actually C 98 does have std distance for iterators Please check out the example class pairs of numbers2 cpp to see how to use a common name like distance safely We ll come to name conflicts la
80. is the name of a class template 110 14 The best random number generator If only the best is good enough FRESH QUANTUM RANDOM NUMBER GENERATOR The fastest quantum random number generator to date Credit ICFO Delivers a very random number in every nanosecond This was needed in a quantum entanglements experiment at Delft More in einstein god dice 111 15 C 11 generation of pseudo random numbers A random number generator is a program that given a seed say 23525176471263 produces a sequence of numbers that appears random It s like a number mill In goes the seed and out comes flour of random numbers Always the same output Stages to invoke C 11 random number generation 1 Include the headers 2 include lt random gt include functional if you use std function Choose the type of random number algorithm std mt19937 gener Mersenne twister Here gener is my name for the generator only this appears in the code and changing this single line I can change to another generator linear congruential engine subtract with carry engine std linear congruential engine gener Choose from available distributions uniform real distribution parameters start and end normal distribution parameters meand value and variance Uniform distribution unif dist U 0 1 is created like this std uniform real distribution cdouble unif dist 0 1 and a norm
81. ith the pointer to c can change c at will Importrant by looking at just the function call there is no way to tell if the argument is passed by value or passed by reference Only the function prototype declaration reveals which it is The problem is that the prototype can be on line 3467 in the header file myheader for this job hpp In practise you end up guessing whether a function can change the contents of variable c dothis c doesn t but dothat c does 38 8 1 Why would a reference be safer than a pointer Imagine that the memory of a computer is a stack of drawers A pointer int c means roughty the integer c is in the upmost drawer or actually just upmost drawer take c Such pointers may be doing funny mistakes from there e The upmost drawer may have a wrong number or not an integer but a mouse trap real number e The upmost drawer can be bottomless contain a pointer to the drawer below e The freaking drawer doesn t even exist pointer to outside computer memory space A programmer is let to do all this without the compiler ever noticing any bad deeds Running the program gives odd results I once wrote a program that was supposed to do a simple calculation result was something like Full moon tonite C reference to variable c int amp c is like a cord or a thin rope tied to the integer c There is no place for mistake the end of the cord always has c because there is no way to det
82. le result YN 74 symmetric pulse abs a aba ona oO yat a E al aba n 11 ones 10 ones data ones for int i n1 1 i lt n n1 i data i 0 0 cout lt lt before fft n cout lt lt data myfile lt lt data myfile lt lt endl ff Base gt 3 forward FFT int status fft daba ll je jubar eS oa if status 0 COLES wesiile wn s return le cout lt lt after forward fft divided by sqrt n n data sqrt n cout lt lt data myfile lt lt data 147 148 Result 15 original pulse e Fourier transformed pulse e 0 5 0 20 40 60 80 100 120 The data is represented as usual in FFT in wrap around order You can easily write utilities to do the wrapping from natural order and unwrapping to natural order move data or use an index table 149 22 3 GSL differential equations Let s solve a simple equation one that we can solve analytically y t ty y 0 1 exact solution y t 2e 472 Example 55 numerics gsl_ode_simple cpp ESSE SN COM dE SAO 0 2 g std c 11 gsl ode simple cpp gsl config libs tinclude lt iostream gt include lt iomanip gt include lt cmath gt include lt vector gt tinclude lt gsl gsl_matrix h gt include lt gsl gsl_odeiv2 h gt using namespace std int func double t const double y double f void xparams 0 t y 0 yO y f 0 y dy dt t y re
83. linear algebra Eigen Library linear algebra Matrix Template Library Loki C Native Template Library PoCo C Libraries CGAL Blitz Fastflow View Template Library STXXL Standard Template Library for Extra Large Data Sets 35 7 2 C Standard Library C Standard Library comes with all C compiler suites What then is the Standard Template Library STL many people talk about I rise my hands here I ve figured out this much The STL existed before the Stan dard Library and parts of it seem to be now in the Standard Library People talk about STL and mean the Standard Library and also talk about the Standard Library when they mean STL As an end user of the C products I couldn t care less and have taken the pragmatic approach to accept it s a library that comes with the compiler suite I call thme both as the C Standard Library Period On late hours I simply think that STL means Standard Template Library official STL means STandard Library inofficial my own idea to save my boiling brain The Standard Library contains also the C Standard Library so you are able to use some C features as well But not all of C notably C programmers avoid naked pointers to the extreme and they are never really needed in C Some parts are templates multi purpose models such as e containers are collections of objects vector stack deque list map e iterators are for going through elements of containers beg
84. lling you how to use plain arrays because they don t conform with C containers and their memory management is manual labour Deleting a two dimensional dynamic array array2 is not simply delete array2 Below are two examples Example 49 basic array to function cpp Passing an C array to function in C style include lt iostream gt using namespace std youd f int d conet ant sized C style empty tells compiler d is an array d is passed by reference You are dealing with the original array not with it s copy for int i 0 i lt sized i cout lt lt i lt lt lt lt d i lt lt endl int main int 3254755 j 3 return 0 136 Example 50 basic array_to_function cpp Passing an C array to function in C style include lt iostream gt using namespace std void f int d const int sized C style think of d as a pointer for int i 0 i lt sized i cout lt lt i lt lt lt lt d i lt lt endl int main int 3 512 4 5 j 3 return Q e There is one more data type suitable for numerical data in C std array It s a container for start From cppreference std array can be initialized like this std array lt int 3 gt al 1 2 3 double braces required Std array int d 22 122 3 except after sbdosarray std string 02 asi 1 std stringa but Notice how std array type and dimension is defined as
85. ls TS Performance TR Parallelism TS Array TS Concepts TS Networking TS Tx Memory TS Concurrency TS 1 3 How transferable is C We d better first settle what is meant with transferability Do we transfer source code or binary Java is well known for its tranferability thanks to a clever combination of a programming language and the Java Virtual Machine JVM C is just a programming language The questions one should ask is how much work do you have to do to make a program 10 years from now It s fairly common that people in science write the same task already a third time Example The first version was in fortran then in C and now in Python planning to write it in Haskell Getting any faster No but the structure is nice and the user interface is cute More modular and easier to expand You hope Why does one reprogram It has to do with sociology some programming languages are more fashionable 1 4 About this document The author is not a C guru and never will become one I see programming languages as tools that help me understand Nature just like mathematics Programming and mathematics deserve to be disciplines for their own sake but I have no desire of becoming a car manufacturer or repairman I just want to drive I spend some time in introducing you some C 11 features it s a whole lot different to what old C used to be As s general goal I try to avoid pointers at least naked pointers and use references C do
86. lt return OSs continues 86 int main 1 using namespace std MyClass testcllass 252 3014 150 2035 e una versal initralization cout testclass endl uses overloaded ofstream out MyClass out out testclass endl same output to file j on screen and in file MyClass out a 2 20000000 b 3 10000000 wv 1 00000000 2 00000000 3 00000000 The overloaded operator lt lt can also write to a file Without overloading cout lt lt testclass cannot work As s reward we have a clean output of objects of a self made class without a visible call to a function This may sound a small achievement but in numerics you learn to appreciate clean formulas Assume you have a matrix A vectors a b and c and hte job is to compute c Ab a With overloading this will at best look likes c Axb a which is easier to decipher than a nested function call c vec_add matrix_multiply A b a Now that you asked the overloaded operator lt lt is compiled with this logic Operation cout lt lt testclass means find from the class of the object cout that ll be ostream the method lt lt and call it with arguments cout testclass This translates to something not completely unlike ostream cout testclass Next the compiler searches the methods ostream lt lt to find a function whose arguments are ostream amp MyClass amp or without the ambersand a
87. lt f IO Ofi t yi t ya t brO f anc at d Oy t nr NEC RUN mU Au Je do ur 1 i on n x t anta ya t dfdt o y t df tyi t y2 t dfdt 1 153 Implementation e The function func 4 top rows of the table computes x t y t from known values x t y t The results are stored to table f elements f 0 f 1 e The function jac O 4 bottom rows of the table hold the Jacobi matrix This information is used only in higher order solvers the more knowledge the less function evaluations A twist to bookkeeping jac is supposed to fill a 1 dimensional table dfdy where the Jacobi matrix is stored as m i j dfdy i dimensio j Hence the odd calls to gsl_matrix_ If you find this hard to follow fill dfdy as you please A plain for loop is not a bad idea Example 56 gsl func jac h int func double t const double y double f void xparams double mu double x params t 0 y 1 ft 1 y 0 mu y 1 y 0 y 0 1 i return GSL SUCCESS int jac double t const double y double dfdy double dfdt void params double mu x double x params g sl matrix view dfdy mat gsl_matrix_view_array dfdy 2 2 gsl_matrix m amp dfdy_mat matrix gol Matrixsses m 0 0 0 0 gsl_matrix_set m 0 1 1 0 gsl matrix set m 1 0 2 0 muxy O y 1 1 0 gsl matrix set m 1 1 mux y 0 y 0 1 0 dfdt 0
88. lt lt c1 lt lt endl cout lt lt c2 lt lt c2 lt lt endl Hija Ge cule ea e cea cout lt lt real c2 lt lt real c2 lt lt endl cout lt lt imag c2 lt lt imag c2 lt lt end l GO lt lt cop ccu ens cout cix c2 clx c2 endl cout lt lt conj c1 lt lt conj c1 lt lt endl cout lt lt c1 c2 lt lt c1 c2 lt lt endl return 0 7 In math the multiplication of complex numbers is x a ib y c id zy ac bd i ad bc So how does the compiler know that if x and y are type std complex then x y means this operation It s called operator overloading but let s first take a look at the simpler function overloading 78 10 Function overloading optional arguments and default arguments Function overloading in C means you can assign different but related tasks under one function name This is nothing new in math the exponent of a real number z is exp x and the exponent of a complex number c is exp c Even exp M of matrix M is under the same name exp C has no function overloading only math functions have been overloaded Designers of C have apparently a different opinion of good practise The compiler has to distinguish which task you want the function to perform It does this by means of the types and number of arguments a bit like in math Moreover function overloading makes two things possible e Optional arguments are ones you don t always have t
89. lt int gt A 2UL 3UL 0 ACTO E AU 4 A e Os Performing a dense matrix dense vector multiplication DynamicVector lt int gt y A x 132 Printing the resulting vector Sbd comes y Sn lt lt Nao Instantiating a static column major matrix The matrix is directly initialized as LEX Ed cue M C0 2 I KA Q StaticMatrix lt int 3UL 2UL columnMajor gt B 3 0 1 1 2 0 Performing a dense matrix dense matrix multiplication DynamicMatrix lt int gt C A B Printing the resulting matrix Seles com lt lt YC Ama lt lt lt Wal g std c 11 Ipath to blaze blazel cpp calc module add gcc g std c 11 I usr local blaze blaze1 cpp 133 19 Calling C or fortran from C Here are some principles about mixing fortran C and C Example 48 Fortran subroutine dgemm is part of BLAS it computes the matrix matrix product SUBROUTINE DGEMM TRANSA TRANSB M N K ALPHA A LDA B LDB BETA C LDC DOUBLE PRECISION ALPHA BETA INTEGER K LDA LDB LDC M N CHARACTER TRANSA TRANSB DOUBLE PRECISION A LDA B LDB C LDC Declaration in C extern C i void dgemm_ char transa char transb intx m int n int k doublex alpha doublex a int lda double b int ldb doublex beta doublex c intx ldc e extern C means compile C style It tells the C compiler to forget about function overloading and identify the function by
90. ma omp parallel num_threads 3 parallel over 3 threads pragma omp for following for loop is parallel pragma omp for private k each thread has it s own copy of k pragma omp for shared m all threads use the same m pragma omp for ordered schedule dynamic parts of for must happen in due order pragma omp ordered following statement must be done in ordered fashion A loop variable such as i in the previous example is automatically private meaning each thread get its own value of i from the openmp scheduler The scheduler decides who computes what In linux the number of thread is set to four usign environment variables setenv OMP NUM THREADS 4 csh or tcsh shell export OMP NUM THREADS 4 bash shell These are for run time only you can change these as needed without re compiling the code export OMP NUM THREADS 2 a out export OMP NUM THREADS 16 a out runs a out first using two threads cores and then using 16 cores The function omp set num threads 16 does the same overriding OMP NUM THREADS See the next example OMP NUM THREADS and omp set num threads are only suggested maximum number of threads 181 Example 68 numerics openmp reduction cpp for loop reduction with timing include lt iostream gt include lt vector gt include lt cmath gt include lt fstream gt include lt omp h gt include lt chrono gt using namespace std using namespace std chrono int
91. mendation Use always if available emplace back will contruct the value at the end of the std vector while push back will contruct it someplace else and move it to the vector The latter simply has to be slower Example 38 numerics emplace_vector cpp How emplace_back beats push back in speed include lt iostream gt include lt cmath gt include lt vector gt include lt chrono gt using namespace std using namespace std chrono struct MyThing 1 int idat std vector lt double gt x MyThing int idat_ vector lt double gt x noexcept idat idat_ x x_ E 104 int main 1 const int N 1000000 vector lt MyThing gt v auto tO high resolution clock now for auto i 0 i lt N i v push_back MyThing i vector lt double gt 1 0 2 0 3 0 4 0 5 0 6 0 auto ti high_resolution_clock now auto d duration cast milliseconds gt t1 t0 cout vector push back took lt lt d count lt lt ms n v clear remove reservations in the container tO high_resolution_clock now for auto i 0 i lt N i v emplace_back i vector lt double gt 1 0 2 0 3 0 4 0 5 0 6 0 j ti high resolution clock nov d duration cast milliseconds gt t1 t0 cout lt lt vector emplace back took lt lt d count lt lt ms n 105 emplace_back can also call the constructor with parameters needed in contructing You don t even consruct the valu
92. mp It can find one as a method with the friend attribute the one we just wrote 87 12 C Standard Library More algorithms 12 1 std for each The algorithm std for_each performs a given operation to all elements As arguments you give the beginning the end and what to do Example 30 numerics vector print foreach cpp std for_each and printing some elements of a std vector include lt iostream gt include lt vector gt include lt algorithm gt using namespace std void doubleout double y cout lt lt lt lt y int main vector lt double x IMPIMPI 3y cout lt lt x vector NNU for_each x begin x end doubleout for_each x begin x begin 2 doubleout cout lt lt end l return 0 This applies doubleout to all elements In general the applied function can be anything as long it s declared anything double x i e it must eat doubles 88 12 2 When to use std for_each Now you may wonder what more std for_each has to offer than the range for loops in C 11 After all you can print all elements of a container more neatly with a range for loop Example 31 numerics vector print range for cpp Range for loop and printing all elements of a std vector include lt iostream gt include lt vector gt using namespace std void doubleout double y cout lt lt lt lt y int main vector lt double xl comm lt
93. mplate cpp and basic static assert cpp The latter shows how to check types in compile time to catch attempts to use a wrong template by mistake Skip the rest unless you are familiar with templates or curious Take a peek at the series of articles An Idiots Guide to Cplusplus Templates Still curious See what strange things can be done at l ype Rich Style for Cplusplus The original source is Bjarne Stroustrup s Youtube lecture It tells the story about C 11 code like this int main Speed spd 1m 1 2s Acceleration acc 8m 3 3s2 i Own types Speed and Acceleration are defined simply with typedef or using see 9 4 This we know But the units Stroustrup claims that the template implementation of the unit system is so clever that the compiler can do type checking and find programming mistakes that show as unit mismatch Moreover the compiled code has only the double s but no units around after compilation No speed penalty Anyone who after lengthy calculations on paper has got the final result energy 234 4 kg m values greatly automated unit checking 34 7 C libraries 7 1 C Template libraries just a list Standard Template Library Boost C Libraries STLSoft C Libraries Electronic Arts Standard Template Library Adobe Source Libraries Intel Threading Building Blocks C Templated Image Processing Library Database Template Library Windows Template Library Armadillo C Library
94. n of data a dataset of length n with stride stride The arithmetic mean or sample mean is denoted by ji and defined as M yn where z are the elements of the dataset data For samples drawn from a gaussian distribution the variance of ji is o N double gsl stats variance const double data size_t stride Function size t n This function returns the estimated or sample variance of data a dataset of length n with stride stride The estimated variance is denoted by 6 and is defined by 42 1 212 0 Won LA where 2 are the elements of the dataset data Note that the normalization factor of 1 N 1 results from the derivation of 6 as an unbiased estimator of the population variance 0 For samples drawn from a gaussian distribution the variance of 6 itself is 20 N This function computes the mean via a call to gs1 stats mean If you have already computed the mean then you can pass it directly to gs1 stats variance m 141 Example 52 numerics gsl_statistics cpp Arithmetic mean and standard deviation include lt iostream gt include lt vector gt tinclude lt gsl gsl_statistics h gt using namespace std int main void using plain array double data 5 17 2 185 16 5 18 3 12 6 double mean variance mean gsl_stats_mean data 1 5 variance gsl_stats_variance data 1 5 cout lt lt mean lt lt mean lt lt endl cout lt lt variance lt lt variance lt lt en
95. n to two coupled 1st order equations Define a new variable y x t y t y t x t uy t x t 1 An n th order differential equation is solved as n coupled 1st order differential equations The code solves two unknown functions x t y t point by point starting from initial values at t 0 Once a point x t y t has been solved the derivatives in that point can be computed and the next point x t dt y t dt can be solved and so on The Van Der Pol oscillator position x t can have very sharp turns for some values of u Near sharp turns we apply adaptive stepsize In order to maintain numerical accuracy the algorithm takes shorter steps in t 1 f this is neglected the solution goes astray after every steep turn It will follow a solution curve but not the one we started with until after the next turn it picks yet another new curve Since the numerical accuracy is limited anyhow you can be sure this happens sooner or later in t 152 The biggest challenge is bookkeeping We have i the math on paper ii 1st order differential equation notation in GSL and iii the notation in GSL solver The GSL manual uses this general notation dy t HO lO an y E Lam Don t spend too much time studying this but tabulated the three notations are related like this A B C x t yi t y 0 y t y t yl v t y t an fim FO t y t v t EIOS 1 28 fat walt vo
96. nction multMatrix lhs rhs result res and use it as D A x B x C It s also simple to create you own matrix type and overload fortran type matrix double pointer data end type Matrix interface operator x module procedure multMatrix end interface operator x function multMatrix lhs rhs result res and write D A B C but then you have to dig the data from the objects in fortran it s done A data 85 Example 29 advanced myclass_overload cpp operator lt lt overloaded to print self made class objects How to teach lt lt to print a self made class object include include include include include lt iostream gt lt vector gt lt iterator gt lt iomanip gt lt fstream gt class MyClass double a b std vector lt double gt v public Mes italia abd on MyClass double a_ double b_ std vector lt double gt v_ overload lt lt and tell it how to print MyClass objects define lt lt as friend to let it access private data a a_ b b_ v v_ friend std ostream amp operator lt lt std ostream amp os const MyClass amp obj ie std ostream amp operator lt lt std ostream amp os const MyClass amp obj using namespace std os lt lt fixed lt lt setprecision 8 some I O manipulation os lt lt a lt lt obj a lt lt b lt lt obj b lt lt endl os lt lt v for auto ele obj v os lt lt ele lt
97. nstructor of the class Generator std generate may copy the third argument std for_each may copy of the third argument The random number generator is copied too Why not copy a random number generator rng A rng is just another program Given a seed it can produce a nearly random number sequence The sequence depends only on the seed the algorithm is deterministic and the same seed gives exactly the same number sequence That s why it s often called a pseudo random number generator If you copy the rng you have two identical number mills If you compare the two number sequences they produce you find that within each sequence the numbers are almost random but the numbers in the two sequences are badly correlated Turning the crank of two similar mills in a simulation code can give you exciting but wrong results Every generator that is not giving a constant output must have a state In other words a generator that can give non constant output has to remember where it is A clock that cannot keep track of time is a stopped clock Copying a stateful generator has to be done with care 1 Ah why not give the generator a new seed and get a new random sequence different from the other That way you would have several number mills with different seeds Bad idea The rng algorithms have been tested to produce a number sequence random only in relation to numbers within the same sequence The basic problem is that the seeds that det
98. o give For example estimate x may do a slightly different calculation than estimate x a No need to call them estimatei x and estimate2 x a e Default arguments unless given the argument has its default value Imagine the boredom and messy program if you always have to call the function myfun x y alpha beta gamma with all five arguments even though you know you in most cases have alpha 1 beta 5 gamma 14 5 In C you can set these values as defaults and call the function simply myfun x y Important difference Optional arguments are used in the function only if they are present Their presence or absence causes different code to be executed Default arguments are always used in the function and they must have some values default or given Don t overdo function overloading There s no real benefit in calling two functions by the same name if their tasks have nothing in common 79 Example 27 basic function_overload cpp Function arguments a and b are optional Function integ can do two different things depending on arguments include lt iostream gt using namespace std double integ void cout lt lt no args to integ integrating from 0 to i lt lt end1 return 120 77 just test double integ double amp a double amp b cout lt lt two args to integ integrating from lt lt a lt lt to lt lt b lt lt endl return 2 0 just test int main double
99. o write ones yourself because most C libraries use templates to achieve excellent performance and multi purpose applicability A template is an abstract model of what to do std sort is a template to sort arguments many kinds of arguments not just numbers std pair is a template to combine two things many kinds of things std swap is in a template to swap two things many kinds of things For example std vector lt int gt x tells the compiler find a model for a vector in the std library and implement instantiate it using int as data type and call it x The angular bracket is for template arguments Basic templates are simple Suppose you have a function that is written for double s void f double x cout lt lt argument is lt lt x lt lt n do something more but you know it should work the same way also for int s and string s then options are to overload f more on overloading functions in chapter and manually write the three functions double f double x int f int x and string f string x much typing or make a template template lt typename T gt void T s cout lt lt argument is lt lt x lt lt n do something more i 33 Even this simple template is not omnipotent there are lots of type T objects that cannot be printed by std cout and those give a long compiler error message For a more serious example see basic simple te
100. olor on luokan Car j sen se ei ole public joten se on private get_color on ainoa toimiva menetelma lukea auton vari cadillac is black_and_white ferrari is red buick is black_and_white cout lt lt buick color lt lt endl Car demolished Car demolished tuottaa k nt j n virheilmoituksen Car demolished car cpp 5 10 error std string Car color is private 26 Notice how the class Car has a function with the very same name as the class itself Car is a class constructor used to create an object a Car and to give it some properties This constructor is used if and only if the color is given in parenthesis The compiler looks how we try to create a Car and find a match from the list of possible ways to constructor a Car Object is a manifestation of a class a thing that has all the properties defined in the class Car lada black 7 Car as a class lada is am object Meaning lada is a Car and it s black When the compiler comes to this line it fetches a constructor in the class Car answering the needs a function that looks like this Car string Such a function exists so the compiler will create a black lada It may become as a surprise that constructor of a colourless Car fails Car volvo fails because Car has no constructor Car By the way this would be the so called default constructor one made by the compiler if no other constructor is sp
101. om 155 Set initial values and find the solution between t 0 100 Example 57 numerics gsl_ode_part cpp include lt stdio h gt tinclude gsl gsl errno h tinclude lt gsl gsl_matrix h gt tinclude lt gsl gsl_odeiv2 h gt tinclude esl mae Jae lal int main void double mu 10 gsl_odeiv2_system sys func jac 2 amp mu auto solver gsl_odeiv2_driver_alloc_y_new amp sys gsl_odeiv2_step_rk8pd le 6 1e 6 0 0 abla al double t 0 0 t1 100 0 double y 2 1 0 9 0 r5 for i 1 i lt 100 i double ta i ti 97 00 05 int status gsl_odeiv2_driver_apply solver amp t ti y status GS SUCCESS NT printf error return value d n status break C printf is here a lot shorter than C iomanip would be printt f 5e V 50 heat t VIO yl lis j gsl odeiv2 driver free solver return 0 156 2 5 2 NNNM 0 5 bl 0 5 A 1 5 PA A 2 ao af o Y 0 10 20 30 40 50 60 70 80 90 100 Van Der Pol oscillator with y 10 The line is just a linear interpolation between solved points 157 22 4 GSL interpolation It s convenient to hide the GSL specific parts to a header I chose to use std vector for data The example recognizes two spline types not very ambitious e cspline natural cubic spline Changes in one point causes non local changes to the spline uns
102. only carefully chosen methods the permit to touch the data Secondly upon improving the code you can improve the methods without disrupting the workings of the rest of the code A class has properties that can be inherited meaning another class can inherit the properties and methods of a parent class Smart base classes can be inherited and you can recycle code In the C language itself a lot of things are inherited for example ofstream header file inherits this chain of mre general purpose headers ios base ios ostream ofstream On the bottom there is a general purpose class ios_base The next example sorry that this has nothing to do with numerics makes a class Car and an object cadillac 24 Example 5 basic car cpp Car class How to create a Car class with member function get_color include lt iostream gt using namespace std class Car string color private member public string pet_color const return color Gar string col color col CONS tubo Car cout lt lt Car demolished n destructor UE int main Car cadillac black_and_white create a Car with a color Car ferrari red coute Caisliililee ae s cout lt lt cadililac get_colior lt lt endl get cadiliac s color Cout ur errari cout ferrari get color endl Car buck cadillac copy cadillac to buick all properties COULD cout lt lt buick get_color lt lt endl 25 c
103. ould return the value f x params for the argument x and parameters params where x is an array of size dim giving 163 the coordinates of the point where the function is to be evaluated size_t dim the number of dimensions for x void params a pointer to the parameters of the function The integrand is passed as a function pointer ROOT CERN C package includes among plenty of other things Monte Carlo integration and a wrapper to pass a function object to GSL 164 Example 60 numerics gsl_monte carlo cpp see http www gnu org software gs1 manual html_node Monte Carlo Examples htm1 ADS SN gol nont carlo CM SIC One t te Ol oI include lt iostream gt include lt iomanip gt include lt string gt tinclude lt gsl gsl_math h gt tinclude lt gsl gsl_monte h gt tinclude lt gsl gsl_monte_plain h gt tinclude lt gsl gsl_monte_miser h gt tinclude lt gsl gsl_monte_vegas h gt using namespace std double func double x size_t dim void params void display_results string title double result double error int main 1 coust size t dim j double result error string method gsl_rng r gsl_monte_function G double a dim 0 0 double b dim M PI size_t calls 500000 random number generator gsl_rng_env_setup r gsl_rng_alloc gsl_rng_default 165 method plain auto s gsl_monte_plain_alloc dim gsl_monte_plain_integrate amp G a b dim
104. out using std pow double sum_of_squares return 0 end of recursion template lt typename T typename Ts gt double sum_of_squares T first Ts rest double output to be sure pow takes a double or float as base int base will be implicitly converted return pow first 2 sum_of_squares rest recursive i int main cout lt lt sum_of_squares 1 2 lt lt n cout lt lt sum_of_squares 1 2 3 4 5 lt lt n cout lt lt sum_of_squares 1 1 2 2 3 3 lt lt n cout lt lt sum_of_squares 152 153 1 ntc In the coming standard C 1y g std c 1y you could write Mi EXPE sunmor scueares C sales Ws ESE For more examples see http en wikipedia org wiki Variadic template which contains also an example of how one could improve printf in C 11 Don t worry about runtime speed C variadic templates are expanded at compile time 18In C variadic functions are resolved at run time so there is a speed penalty 126 18 Linear algebra which library to use THE answer is LAPACK and BLAS written in fortran and available in Processor manufacturers have their optimized versions ACML MKL and an independent freely distributed automatically optimized ATLAS C versions CLAPACK CBLAS are heavily used in csf GSL the Gnu Scientific Library offers an extensive package of C numerics There are optimized LAPACK and BLAS libraries by Intel MKL and AMD ACML To
105. piled to a function object and easily inlined Boost has it s own lambda s parts of the ideas came to be in C 11 Usage Define a temporary function on the spot exactly where it s used It can remain unnamed Bonus No class needed to get a function object The lambda remains local and you can t misuse it elsewhere The brackets in a lambda are for capturing things from outside Here is a list of what is captured and how capture nothing x amp y capture x by value y by reference amp capture all external variables by reference capture all external variables by value amp x capture x by value all else by reference amp x capture x by reference all else by value Example 62 advanced autolambda cpp A locally applied function func g std c 11 autolambda cpp include lt iostream gt using namespace std int main auto func PO cout lt lt Hello world n gt funcio 175 Example 63 advanced lambdal cpp Output a container and the cubes of the elements include lt iostream gt include lt vector gt include lt algorithm gt include lt cmath gt using namespace std int main vector lt int gt v 11 22 33 for_each v begin v end Tinten cout lt lt n lt lt lt lt pow n 3 lt lt endl s 11 1331 22 10648 33 35937 Note int n is passed by value 176 Example 64 advanced lambda2 cpp
106. plate lt gt Struct Factorial lt l gt static const long int value 1l J int main i const long int facti3 Factorial lt 13 gt value std cout lt lt facti3 lt lt std endl The number 13 is fed to Factorial as a template parameter Factorial lt 13 gt The compiler instantiates the template Factorial lt 13 gt which tells to instantiate the template Factorial lt 12 gt and so on until Factorial lt 1 gt is instantiated At this points the compiler notices that there is a template specialization corresponding to argument lt 1 gt and the recursive instantiation ends Compilers can t handle very deep recursive instantiation and builtin integer data type can t hold large factorials anyhow 98 13 3 Function objects functors Recommendation Use frequently in numerical code Function objects or functors for short are rather popular in numerics You can define a function in a class that is not a method member function This function is called when the class name is used as a function name Don t be put off this is actually very simple Example 36 numerics functor cpp A function object to compute sin x cos x and tan x Function object functor include lt iostream gt include lt cmath gt class TrigFuns jonuiodLstce 8 double operator double x std cout lt lt x lt lt x lt lt std endl std cout lt lt sin x cose Bau Ge in 2 std cout lt lt sin x
107. ple int main 1 std vector lt double gt v1 1 153 2 8 dd std cout lt lt minimum element my min element v std endl return 0 minimum element 1 3 Be cautious when coding specialized versions of C Standard Library algorithms Here I used the namespace my to protect the home made function min element This is necessary because somebody else may use that simple descriptive name and once you incorporate his her program to your you end up with a name conflict 65 9 6 2 std swap is a template Previously given swap using reference variables is for swapping one type of things only What if we want to swap two real numbers or two other type of variables Terribly boring to write each variable type it s own version of swap Write a template Actually it s been done already std swap a b looks like this Like many standard library codes this too were refurnished in C 11 The user sees only the improved performance Example 18 advanced swap_template cpp std swap template template lt class T gt void swap T amp a T amp b j template class T size t N gt void swap T amp a N T amp b N T c std move a a std move b b std move c for size_t i 0 i lt N i swap a i b i C 98 version was TUE I T cla a b b c Vir template this is a template class T means T is a class no need to know what class void swap means swap returns nothing
108. plet and incorrekt and it has lots of bad fomatting This kind of put me off and I had to go for coffee 92 12 4 std generate algorithm One way to generate values to a container Example 34 numerics generate_random_vector cpp Fill a std vector with random numbers include lt iostream gt include lt algorithm gt include lt vector gt using namespace std double double_random poor random numbers 0 1 return rand x 1 0 RAND_MAX avoid int int int main 4 vector lt double gt v 20 fill vector with random numbers generate v begin v end double_random output for unsigned i 0 i lt v size i cout lt lt i lt lt lt lt v i lt lt endl j j The next chapter takes a look at why generate may be dangerous for random number generation 93 12 5 C Standard Library algorithms take care of copies As the previous example showed std generate is a nice way to fill a container with random numbers There is a potential risk however The functioning of std generate is equivalent to this template class Forwardlterator class Generator gt void generate ForwardIterator first ForwardIterator last Generator gen 1 j while first last x first gen The third argument is Generator gen and Generator is a class The generator is passed by value so a copy of gen is used l he compiler creates a copy of gen using the copy co
109. rough these pointers can alias This can prevent some optimizations from being made In C99 the restrict keyword was added which specifies that a pointer argument does not alias any other pointer argument In Fortran procedure arguments and other variables may not alias each other unless they are pointers or have the target attribute and the compiler assumes they do not This enables excellent optimization and is one major reason for Fortran s reputation as a fast language Note that aliasing may still occur within a Fortran function For instance if A is an array and i and j are indices which happen to have the same value then Afi and A j are two different names for the same memory location Fortunately since the base array must have the same name index analysis can be done to determine cases where A i and A j cannot alias In pure functional languages function arguments may usually alias each other but all pointers are read only T hus no alias analysis needs to be done In short If you use pointers and you are sure two pointers never ever point to the same data you can let compiler generate much faster C code by using the restrict keyword 48 Example pointers p and q point to integers int p int q and will be pointing to elements of a table of integers like this q 3 q 2 q 1 q p p D p 2 pe3 q 1 p 4 q 2 p 5 EM
110. s and a very interesting math part is in The examples are thorough e g how to use your own vector type in a equation is shown here As long as your data structure conforms with a standard or a Boost container you can just throw it into Boost Especially the library odeint for solving ordinary differential equations is comprehensive Boost is in most parts a header only library In linux this means you just copy it to say directory home myusername boost_1_59_0 The example program chaotic system cpp resides in home myusername boost_1_59_0 libs numeric odeint examples and it can be translated like this bash mkdir tmp cd tmp export BOOSTDIR home myusername boost_1_59_0 In s BOOSTDIR libs numeric odeint examples chaotic_system cpp In s S BOOSTDIR libs numeric odeint examples gram schmidt hpp g std c 11 O3 ISBOOSTDIR chaotic system cpp Of course installing Boost to your include path makes life simple if you have enough admin rights Boost examples in calc phys jyu fi module add gcc g I usr local boost 1 59 0 usr local boost 1 59 0 libs numeric odeint examples solar_system cpp a out gt solar system dat and in gnuplot p solar_system dat u 2 4 w 1 u 5 7 wl u 8 10 wl u 11 13 wl u 14 16 w 1 u 17 19 wl 118 16 1 Boost ordinary differential equations ODE If your equation is nth order second order or higher you first separate it to n coupled first order differential
111. s U combas ny Sib deen Olt lt lt lt em E for ato oo conte t 1 compas Nn e ref 95 13 A few things that can potentially speed up your code You may be interested to check out Wiki C Performance improving features If you are a good programmer find out what Move Semantics and Perfect Forwarding mean in C and dive into the pool of Smart Pointers Read books and postings by Scott Meyers 13 1 noexcept no throw quarantee Recommendation Use frequently in numerical code If your function never throws an exception the keyword noexcept may let the compiler to optimize your code more If it does throw your code will terminate Ha you lied Still noexcept is one of the very latest features of C so don t believe just any blog posts just give it a try Usage void myfunction noexcept Mba 13 2 constexpr Compile time constant expressions Recommendation Use frequently in numerical code Not quaranteed to give any speedup but an interesting concept Computation of factorials recursively can be traditionally done like this Common way to compute a factorial int factorial comest ime no abes Face Il 32 if once 12 return 1 13We return to throw catch in chapter 21 96 else fact nxfactorial n 1 recursion return fact i C 11 is able to perform tasks during compilation coded as template metaprogramming or with constexpr
112. s home page After all he created the language 1 8 C in Matlab or Octave You can call a compiled C or C programs from Matlab using MEX files Matlab executable files Matlab can access C C library routines and may work through a numerical bottleneck faster See instructions in mex files html Octave is a free Matlab lookalike with an almost identical syntax but different internals Octave has interface to several languages including C The compiled files used in Octave are called oct files see instruction in 1 9 Cor C in Python 3 Extending Python 3 with C ctypes 2 A really Brief introduction to C A list of some basic elements is in C sheet pdf Example 1 basic mainl cpp Main program include lt iostream gt here cout is defined using namespace std all names like cout come from this list int main 1 ALL BEGINS AH Lg integers i and j de set values to i and j j 20 int k i j define k on the fly cout lt lt k lt lt k lt lt endl output to screen ALL ENDS Story include lt iostream gt Take Standard Input Output Streams Library Header that defines the standard input output stream objects using namespace std unless told otherwise the keyword is defined in std namespace main parenthesis indicate that the main program is a function and can have arguments int main it returns an integer The function main can end with return 0 or return 0
113. s squares 2 and then increments it by one No The result is 2 3 6 Astonishing 83 There are some rules and limitations to operator overloading e Think of operators as functions with one or two arguments They are called unary and binary operators If your operation needs two arguments take one existing binary operator and overload that e You can t invent new operators my_clever new operator is no good Only some of the existing ones can be overloaded FF 4 amp 4 lt gt s k f Y a f le b gt lt c Soa ls lt gt ke gt new delete e The order of execution prevails is executed before You may find tempting to overload an operator to compute powers Especially since the math form x and the function call pow x y look so very different Resist the temptation and use pow or you ll be astonished 84 Curiosity In case you are interested fortran has operator overloading too Funnily in C you can t invent your own operators in fortran you can overload existing ones if you apply it to your own data type but if you apply it to an existing data type you must invent a new operator So working on 2D tables double A B where the data type exists you must invent an operator such as x the strange colours come from C settings fortran interface operator x module procedure multMatrix end interface operator x fu
114. ssion templates do faster than a for loop If you ever find out please tell me If you have no desire to become a C guru skip this chapter and do something useful All this said here goes The basic idea behind expression templates is to use operator overloading to build parse trees Expression templates were used by Todd Veldhuizen in his matrix library Blitz in mid 90 s ever since applied in practically all numerical C libraries see Todd Veldhuizen Techniques for Scientific C Early expression templates were only able to cure one bad C side effect namely that one should not create temporaries in every corner But this only taught C behave reasonably Since then we have learned quite a bit and realized that avoiding temporaries is not the only essence in speed Where did the avoid temporaries goal come about C has this wonderful thing called operator overloading The problem is that if you apply is straightforwardly the compiler easily creates temporaries Consider how the addition operator can be translated to function calls D A B C means add A B C so set temporary M add A B finally set D add M C If A B C and D fill the fast cache memory then the intermediate result M drops something out of cache to slow RAM 10x slower or more Another example is adding three vectors and using only one component D O A B4 C 0 You can easily expand this the most effective way D
115. st double lim int count 0 for each v begin v end amp count lim double amp x if x lt lim x 0 0 count return count int main vec v 7 rand_vector v for auto xcv Cout x AN const double limit 0 4 cout lt lt chopped lt lt vector_chop v limit lt lt values below lt lt limit lt lt endl for auto zir comas lt Nns 179 25 openmp parallel programming Your computer most probably has multiple cores or threads openmp offers an easy way to parallelize loops see for example http bisqwit iki fi story howto openmp Example 67 numerics openmp ex cpp Openmp parallel loop too small a task to get any speed up though g std c il pipe 03 march native fopenmp mfpmath sse msse2 openmp ex cpp include lt iostream gt include lt cmath gt include lt vector gt int main 1 using namespace std vector lt double gt tab 200 const int N tab size pragma omp parallel for for int i 0 i N i threads have their own private i 1 tab sin 2xM_PI i N cout lt lt i lt lt endl will produce mess but shows parallelism j cout end1 for auto t tab contestes gt cout lt lt endl The preprosessor directive pragma omp parallel for tells what to parallelize 180 You need to read more from a better source but here are some other pragmas pragma omp parallel Somewhere later comes a parallel computation prag
116. t exact lt lt exact lt lt endl cout lt lt diff lt lt fabs result exact lt lt endl 167 Method plain result 1 41220870 0 01343586 exact 1 39320393 ldiff 0 01900477 result 1 39132158 0 00346056 exact 1 39320393 ldiff 0 00188235 result 1 39267259 0 00341041 exact 1 39320393 ldiff 0 00053134 result 1 39328139 0 00036248 exact 1 39320393 diff 0 00007746 168 23 Physics example Spin 1 2 Heisenberg chain The Hamiltonian of a spin 1 2 chain has spin spin interactions of two nearest neighbour spins N N A JY Si Sia DO ISP SE SUSE S753 gel i 1 N 1 E J ICAA S754 SS i 1 Here operators S and S the operators measuring spin x and y components respectively were combined to rising and lowering a k a ladder operators S S iS The basis states used here are s m which are the eigenstates of S spin squared S s m s s 1 s m and the eigenstates of S S s m m s m Units are chosen so that A 1 Rising and lowerin spin z component works like this S s m v s s 1 m m aie 1 s m T 1 S s m m s m The Heisenberg chain has spin 1 2 and the states are 11 s st 1 1 PEPE The only difference is in the spin z component m so we keep track on that only The system state is a chain of spin up and down values The ladder operators change the state S t 0 S 1
117. t cmath gt using namespace std class WaveFunction pubik double energy vector lt double gt density E int main vector lt WaveFunction gt basis a vector of WaveFunctions WaveFunction wf for int i 0 i lt 10 i make a 10 wavefunction basis wf energy i i for int j 0 j lt 5 j wf density push back sqrt j i basis push_back wf wf density clear REMEMBER THIS or wf density keeps growing Op use LO Conwy for auto wf basis wf goes through elements of basis cout lt lt energy lt lt wf energy lt lt endl cout lt lt density for auto den wf density cout lt lt den lt lt den goes through a density in wf cout lt lt endl 53 D As you see you can fill a std vector with almost any type of data Here the container contains object of the self made type WaveFunction blue boxes and a number energy and a yet a std vector density basis energy energy energy energy density density LL 1E YE JE OUL density density LJLILILUILU L JE AE IE IE These are exactly the same thing class WaveFunction pPublktice class all is private by default double energy vector lt double gt density is struct WaveFunction double energy struct all is public by default vector lt double gt density We 54 9 2 Moving not copying You may have noticed that moving objects instead of copying
118. t pilfer the resources on the right and gives them to the left destructor tells how an object should be destroyed C 11 has a rule of five saying that if you need to write your own destructor chances are your class is so wierd that you need to provide all five destructor copy constructor move constructor copy assignment operator move assignment operator because none of the default version generated by the compiler meet the needs of your class 5The example used a clumsy constructor a better way would have been to use an initialization list 20 5 1 Delegating a Constructor read on spare time C 11 introduced a way to delegate a constructor The idea is to let for example a constructor without any arguments to delegate the construction to another constructor with with a single argument by telling what a default argument is The points is to avoid copying the same piece of code over and over again in all constructors Here is a bit artificial example this doesn t save much code Lf C 98 EMO ala class Sum class Sum PUDE publici Sum a1 0 a2 0 s ai a2 Sum Sum 0 Sum int E al i a1 0 fs al a2 Sum int i Sum i 0 f Sum int i int j al i a2 j s al a2 Sum int i int j a1 i a2 j s al a2 private private int al a2 s imit El m2 E E In the C 11 code the object mysum to be created without an argument Sum mysum is delegated to Sum 0 who delegates it to Sum 0 0 Th
119. table Natural refers to setting end point derivatives to zero e akima or Akima natural Akima spline Changes in one point causes only local changes to the spline stable Steps 1 Reserve space for accelerator speeds up searches auto accel gsl_interp_accel_alloc 2 Reserve space for interpolation gsl_spline spline gsl_spline_alloc gsl_interp_cspline 3 Initialize the interpolation with known points x y now in a std vector gsl spline init spline amp x 0 amp y 0 x size 4 Compute the interpolated y values to yy iterator posy at points xx iterator posx posy gsl spline eval spline posx acc If needed the 1st and 2nd derivative could be computed gsl spline eval deriv gsl spline eval deriv2 158 Example 58 numerics gsl_spline hpp wrapper to GSL spline spline type cspline is natural cubic spline spline type akima or Akima is natural Akima spline ifndef GSL SPLINE HPP define GSL SPLINE HPP include lt iostream gt include lt cstring gt include lt vector gt include lt gsl gsl_errno h gt tinclude lt gsl gsl_spline h gt namespace my using namespace std using vec vector lt double gt void spline const vec amp x const vec amp y vec amp xx vec amp yy const string amp type auto acc gsl interp accel alloc gsl_splinex spline int
120. template parameters Since it s a container it has many built in methods just like std vector has Fixed size makes compile time checks possible since C 11 also template arguments can be checked For that C has std static assert see basic static assert cpp std static assert can greatly benefit code development by letting you know if you by mistake contradict your own ideas 137 21 Exception handling with throw and catch The C Standard Library has a class dedicated to exception handling std exception One part of it is runtime_error One way to make your own error handling process is to inherit the class class MyException public std exception e and add a new property In numerics exceptions are often simple I m using this handling Coy my_function catch char constx e cerr lt lt e lt lt endl return l ja funktiossa on rivi void my_function void if test throw test failed Error happens if test is true and the exception with message test failed is thrown and we leave the function Later the exception is caught with catch 20The function my_function throws and exception and wishes some exception handling routine will take care of it properly 138 If throw is executed the function execution terminates it s still more gentle than halting the program execution Any code after throw is not executed and essentially the
121. ter in this course 3l The previous example about pairs of numbers obviously works but it is a very common structure And common structure are found in C Standard Library If you ever need a pair of well practically anything use the standard structure std pair lt typel type2 gt Example 7 numerics pairs of anything cpp std pair and std make pair using std pair to make pairs of almost anything include lt iostream gt include lt vector gt int main std pair lt int int gt iipair iipair std make_pair 10 12 std cout lt lt iipair first lt lt lt lt iipair second lt lt std endl std pair lt int double gt idpair 1 120 324 std cout lt lt idpair first lt lt lt lt idpair second lt lt std endl H waero 0E pallets EA A 19 std vector lt std pair lt int std vector lt double gt gt gt ivpairs std vector lt double gt v 10 1 20 2 30 3 40 4 testing vector for int i 0 i lt 10 i ivpairs push_back std make_pair i v for auto icivpairs std cout lt lt i lt lt i first lt lt vector for auto j i second std coute std cout lt lt std endl 32 6 Templates generic instructions and algorithms A template is a generic model how things fed to the model should be handled Templates are advertised as one of the best things C has to offer as always some disagree You ll benefit from them even though you have no idea how t
122. tring args 3 System out println Hello world 4 C 1 include lt iostream gt 2 int main 3 std cout lt lt Hello world n C include lt stdio h gt 2 int main void 3 printf Hello world n 4 JS 1 myfile js 2 console log Hello World 4 command line node myfile js MySQL DELIMITER SSCREATE FUNCTION hello world RETURNS TEXT 2 DELIMITER SELECT hello world COMMENT D 13 Python Numpy Scipy is often an unbeat able combination Eigenvalues and eigenvectors of a random 5x5 matrix import numpy as np import scipy linalg as la A np random randint 0 10 25 reshape 5 5 print matrix A n A la eig A print eigenvalues e_vals e_vecs Le vals how can scipy be fast if it is written in an interpreted language like python 1 7 C in the Net e www ohjelmointiputka net oppaat php in Finnish mureakuha is out of action competition killed it e www cplusplus com e www greenteapress com thinkcpp Book 2012 Don t read old books Anything written before 2011 should be taken with a grain of salt It s futile to force any good outlook for C programs just be systematic with your style When to use capital letters and when to use underscores is up to you Think carefully what to put in header files ll speak up my opinion later About style and good C practices I recommend Bjarne Stroustrup
123. turn GSL SUCCESS j void output double t double y double exact cout lt lt fixed lt lt setprecision 16 static bool first true eii sel cout lt lt setw 20 lt lt t lt lt setw 20 lt lt gsl solution cout lt lt setw 20 lt lt exact solution lt lt setw 20 lt lt error n tarse falso j cout lt lt setw 20 lt lt t lt lt setw 20 lt lt y lt lt setw 20 lt lt exact lt lt setw 20 lt lt y exact lt lt endl 150 int main double exact const int n 50 of points doubles tO 0 0 tl 10 065 time start and end double dt ti t0 n 1 time step double y 1 2 0 initial value table with one entry gsl_odeiv2_system sys func nullptr 1 nullptr not using a Jacobian auto driver gsl_odeiv2_driver_alloc_y_new amp sys gsl_odeiv2_step_rk8pd le 13 1 e 13 0 0 output t0 y 0 y 0 for int i 0 i lt m i auto t t0 ix dt begin of t interval if gsl_odeiv2_driver_apply driver amp t t dt y GSL_SUCCESS cout lt lt FAILED near lt lt t lt lt endl return 1 exact 2 0xexp 0 5 txt output t y 0 exact gsl odeiv2 driver free driver return 0 j 151 The next one is more involved it s directly from the GSL manual The 2nd order nonlinear Van Der Pol oscillator equation is a t uz t z t 1 a t 0 The numerical solver is for 1st order equations so we split this 2nd order equatio
124. uments b turn a 3 argument function to a 2 argument function include lt iostream gt include lt functional gt void f const double amp x const double amp y const double amp z int main std cout lt lt called f with arguments lt lt x lt lt lt lt y lt lt lt lt z lt lt std endl using namespace std placeholders for _1 _2 pd le ee es ae a change the order of arguments with bind auto nvi ctd bind f 23 7 97 1 bind return type is std function my ta d n 9c ge b create a two argument function from f Auto p std bindi 1 52 393993 9 7 7 bind 2rd argument to fixed value 333 3 AZ 114 The next examples are a set of easy to use helper functions to initialize a generator and get uniform distribution unirand normal distribution gaussrand and gaussrand2 or exponentially distributed exprand random numbers The goal here was to hide all inconveniences to the header numerics random hpp include random hpp double random gaussrand clean and simple The header is by no means perfect especially the initialization is clumsy esl Example 41 numerics random cpp C 11 random number generation std mt19937 random number generator ASS SS ta ua cion an ds cle band include lt iostream gt include lt random gt include lt ctime gt include lt fstream gt include lt functional gt using namespace std std mt19937
125. using namespace std const double e abs le 15 e rel le 15 absolute and relative error goal auto stepper boost numeric odeint make_controlled lt stepper_type gt e_abs e rel state type y 2 0 condition y 0 2 gives y t 2 0 exp 0 5 t t const int n 50 solve 50 points const double ti 0 0 t2 10 0 const double dt t2 t1 n 1 for int i 0 ien i auto t ti ixdt auto exact 2 0xexp 0 5 txt my output t y 0 exact auto steps integrate_adaptive stepper my system y t t dt 0 01 cout lt lt steps lt lt steps lt lt endl how many sub division steps were needed Such a simple my system begs to be replaced with a lambda expression auto steps integrate adaptive stepper const state type amp y state type amp dydt const double t dydt 0 t y 0 y t tdt 0 01 This same ODE will be solved using GSL in chapter 22 3 122 17 Formatted output It s important to keep numerical output readable The following example is mixing columns pretty badly XxX yz 1 542234 12 4234 0 1213 13 0 4 234 1 00 and this is easily the default outcome Another thing to remember is that the method may limit the accuracy of the result not the number of stored and filed decimals It s up to you not to publish 6 decimals if the method is reliable only up to 2 decimals 17 1 Formatted output using include lt iomanip gt The next s
126. using namespace std using namespace arma int main mat A randu lt mat gt 4 5 mat is double 4 5 mat B randu lt mat gt 4 5 B B cout A trans B endl cout lt lt Axtrans B lt lt endl return 0 This calls internally a CBLAS routine since the attempt to compile as g arma matrix multi cpp gives the error message among a plethora of other lines undefined reference to cblas dgemm Here dgemm stands for LAPACK BLAS subroutine name meaning double general matrix matrix You get CBLAS as part of GSL or OpenBlas or ATLAS 11 arma matrix multi cpp gsl config libs std c 11 arma matrix multi cpp lopenblas g std c 11 arma matrix multi cpp L usr lib64 atlas lcblas 129 Example 46 numerics arma_eigenvalues cpp Armadillo Eigenvalues of a symmetric matrix include lt iostream gt include lt iomanip gt include lt armadillo gt using namespace std using namespace arma int main mat A randu lt mat gt 5 5 vec eigval mat eigvec A A trans A a way to make a symmetric matrix cout lt lt A An lt lt A lt lt endl eig_sym eigval eigvec A this does all work vec x eigval for unsigned i 0 i lt A n_rows i cout lt lt i lt lt th eigenvalue x eigvec col i x j cout lt lt X x Axx eigval i cout lt lt check Ax lambda return 0 lt lt eigval i lt lt

FYSA120 C++ Numerical Programming Vesa Apaja November 11

Contents

Download Pdf Manuals

Related Search

Related Contents