Polska wersja

Scientific Computing 2019-20

Tuesday lecture 1015am in room 5870, lab 1215pm room 2041(lab)
The solutions of homeworks and projects and reports (instead of an oral exam if one wishes..) should be uploaded on moodle - details below.
Attention please Please upload the homework and project solutions as octave scripts to moodle .
CHAT I started a chat channel on chat.mimuw.edu.pl: #on-sciecomp-19-20-tue-10-14; here there is a link. I will try on Tue 1015-1345 reply your questions (and during my office hours as well).

Exam

Project (25%) and lab short problems (solved in the lab) 50% - oral exam (or a theoritical report) 25%. During the exam one should show (now due to pandemia upload to moodle) a short report (2-3 pages details in project section) and the code (working) and answer a few questions concerning the material presented in the lectures.

ORAL EXAM or a theoretical report

There two options for an oral exam : classical - we set a date and time and connect on google meet or a short report on a given subject extending knowledge from our course - e.g. describe how to solve in LAPACK eignvalue problem what functions are available how to use them a short code using them etc or sparse matrices what formats are available with some details... I will propose a few subjects, see below.

September session

The Rules are the same. Upload everything on moodle and let me know by an email.
If one wants to pass an oral exam on google meet, please contact me by an email.

    Questions/problems for the oral exam. Three questions. I will ask one question and the examinee may pick the other two.
  1. Sparse matrices formats in octave and in general.
  2. Octave: how to solve in octave the basic computational problems of numerical linear algebra such as solving linear systems, LLSP, eigenvalue problems, SVD etc
  3. Octave: how to compute in octave integrals, double or triple integrals?
  4. Octave: how to solve nonlinear systems of equations? How to compute the inverse of an univariate function? How to compute a value of an implicitely defined function?
  5. Octave: how to solve in octave ODEs?
  6. Numerical solvers for ODEs : simple examples, what is stiffness, adaptive step refinement, multistep schems, Runge's schemes etc (no theory - just ideas as it was presented in the lectures)
  7. Numerical libraries: how to call a fortran function in C? Matrices in fortran: how are they stored in the memory.
  8. Numerical libraries: BLAS what type of functions are in this lib?
  9. Numerical libraries: LAPACK : what does this lib contain? i.e. what computational tasks can we solve using LAPACK?
  10. Numerical libraries: other numerical libs - which library one can use for solving/computing: integrals, nonlinear equations, ODEs. Examples of general numerical libs (contaning solvers almost all problems)
  11. Finite Difference Method - the main idea described on the example of BVP: -u''=f on [a,b] u(a)=al;u(b)=beta. How to treat a Neumann boundary condition, e.g., u'(b)=beta?
  12. Finite Element Method (the linear element - the simplest one) - the idea described on the example of BVP: -u''=f on [a,b] u(a)=al;u(b)=beta (optional)
  13. Numerical Codes Optimization Techniques: memory hierarchy, Amdahl's law, basic techniques like loop unrolling etc, optimization compiler options. (optional)
  14. make, makefile - what is this? What are the basic rules? etc How to write a simple makefile? (optional)
Optional - means that one can refuse to reply to this question/problem - then I will ask another.
    Proposed subjects of reports: (if one do not want to pass the oral exam)
  1. nonlinear solvers: Picard and fix point iterations - describtion and properties (that's a bit outside of the main scope of the course - however nonlinear solvers are very important in scie. computing)
  2. implicit Runge-Kutta IVP solvers - idea - a few examples
  3. quadratic FEM in 1D - in the 14th lecture there was described linear FEM - describe the quadratic FEM for the same problem more or less in a similar manner in particular write a resulting linear system for equidistant mesh - assuming that the rhs is computed by the Simpson composite rule (on the same mesh naturally)
  4. Sparse matrice formats - other than the 3 ones,i.e. COO,CSR,CSC, discussed on our course.
  5. Sparse direct solvers - an idea (advanced - difficult)
  6. UMF pack - what type of problems can be solved, discussion of functions - examples of usage
  7. eigenproblem solvers in LAPACK - what type of eigenproblems can be solved, discussion of functions (drivers) - examples of usage
  8. SVD in LAPACK - what is SVD, type of SVD type problems can be solved, discussion of functions (drivers) - examples of usage
  9. KINSOL - library - what type of problems can be solved, discussion of functions - examples of usage
  10. CVODE - library - what type of problems can be solved, discussion of functions - examples of usage
  11. make - makefile - advance usage... automatic variables, pattern rules etc
  12. how to measure a real computation time in C/C++? (that's should be in our course but there was not enough time)

Lecture I will give here a brief program:
  1. Lect no 4 and 5:
    • numerical linear algebra reminder and howtos in octave - LLSP and SVD,
    • sparse matrices (howtos in octave and some basic formats);
    • nonlinear solvers in octave and how to compute a value of the onverse function or implicit function.
    • An introduction to the lab problem from Chapter 3 of Lecture notes (I will write in the lab section the main things in English)
  2. Lect 6:
  3. Lect 7 (Apr 7, 2020)
      Numerical integration in octave. Read documentation in particular of function quad() quadl()etc
    • how to compute a double integral or triple integrals using specific octave function Or just one of quad family function (which one will work? Answer in octave docs)
    • read docs of trapz() function... And how composite trapezoid and Simpson rules work..
    • which function is better for computing an integral of a function given as solution of ode quad()or trapz(). I.e. E.g. F(s)=y(1;s) where y solves IVP y'=f(t,y) with y(0)=s. Compute an integral of F on [0,1]. To use quad() one has to define function which computes F(s) and we pass its handle to quad(). To use trapz just compute the values of IVP for no of times and use them in trapz. Hint : how many times IVP is solved in the both cases...
  4. Lect 8. (Apr 21, 2020)
  5. Lect 9 (Apr 28, 2020) C programming language basics repeat Today lecture is on C language - we assume that on some level of your previous education learned some basics of C programming language - so it is just a reminder or repeat. In scientific computing we need in general only the basic structures and things - we will focus on computations not on the interfaces etc Below you have a list of basics which we may need later in a few future labs.
    1. basic variable types:
      • integer: int, char, signed char, unsigned char;
      • real types: float - single precision, double - double precision
      • arrays:

        int ip[3], float x[6];
        ip[2]=-12;x[5]=12.34;

      • structures: e.g.

        struct complex {
        float re;
        float im;};
        struct complex z;
        z.re=0.0;z.im=1.0;

      • pointers ... how to define one, what are basic usages etc : e.g.

        int *k,i=0;
        ..
        k=&i;*k=23;
        printf("%d\n",i);

      • static, global variables...
      • file type...
    2. loops:

      for(k=0;k<10;k++){
      ...
      }
      while(..){
      ....
      } ;
      do ... while(...) etc

    3. Conditionals:

      if(...) { ..}
      else {...}

    4. boolean operators e.g.

      ==, !=, ! as NOT etc

    5. functions - how to define one, what types can be returned by a function, what types can be parameters of a function... what a prototype of a functions is and why it is important to write ones even if compiler does not require them. Recurrence - howto? Define a factorial function recurentially.
    6. casting: e.g.

      int k=(int)3.4;

      Here

      int k=3.4;

      would give the same result but the 1st version makes the code clearer.
    7. dynamical allocation of memory: malloc(), calloc(), free()
    8. how to represent a real/complex matrix or a vector...
    9. precompiler commands: #include, #define, #if etc what are macros and what is a difference between a macro and a function: e.g.

      double square(double x){ return(x*x);}

      is a function and

      #define square1(x) (x)*(x)

      is a macro computing the square of x... what is the difference? Hint: can we apply both the function and the macro to any numerical type e.g. float or integer one?
    10. standard header files e.g. #include, #include etc and the content of other standard libraries... in particular repeat basic output/input functions like: printf(), scanf() etc
    11. how to compile a simple one file source code, or a bit more complicated projects comprising many files
    12. file operations how to read a text file, or create and wrtte to a text file some numerical data e.g. a vector or a matrix, how to save/load a memory block (bytewise):

      fopen(),fclose(), fscanf(),fprintf(), fgetchar(), fread(),fwrite(), feof() etc

    and that's it - this list contains basic which you should know - you should already learned all that in the previous years if not, please catch up with that - there are thousands of tutorials on line, wiki pages etc Anyway I will provide solutions for some of the lab problems later after the lab date naturally...
  6. Lect 10 (May 5, 2020) How to use a fortran function (library function) in C.
    1. Numerical types in fortran and their counterparts in C languae.
    2. Vectors, matrices - how are represented in Fortran and how to represent them in C in order to use a fortran numerical functions
    3. How to write a prototype of a fortran function knwoing only the source code of this function in fortran.
    4. How to call such a fortran function.
    5. How to compile (link) our C code with a fortran function or library
    6. A simple example namely DNRM2() fortran function from BLAS LEVEL 1 numerical lib - how t owrite a prototype and a simple usage examle.
    7. A more complex example namely how to use a BLAS level 2 function (a subroutine what means a function not returning any value) sgemv() function - how to write a prototype and how to use it - again fro small size data - this function computed y=alpha*Ax+beta*y or y=alpha*A^Tx+beta*y for a given vectors x,y, scalars alpha,beta and a matrix A. The vector y is overwritten.
    More details are in pdf files with the notes, see below Section 5.1 of the notes.
  7. Lect 11 (May 12,2020) Numerical libraries
    1. BLAS basic linear algebra subprograms - please read thru the docs of tbis libeary, it us a free lub, docs are readily available. There are 3 levels of BLAS functions, naming is similar as in Lapack.
    2. LAPACK Linear Algebra Package. The best free linear algebra lib for dense matrices. Read about naming, types of functions, numerical algebra problems that may be solved. In general lapack needs optimized blas lib, then it is fast.
    3. eigenvalue problem for sparse matrices: ARPACK
    4. ode solving libraries: e.g. lsode lib, cvode, odeint etc
    5. nonlinear solvers e.g. kinsol, minpack, petsc etc
    6. sparse numerical linear algebra packages; sparse direct solver e.g. umfpack, sparse,
    7. iterative linear solvers e.g. Petsc (actually just a part of PETSC contains linear iterative solvers and preconditioners), hypre etc
    8. PDEs - PETSC - Portable, Extensible Toolkit for Scientific Computation (PETSc), is a suite of data structures and routines for the scalable (parallel) solution of scientific applications modeled by partial differential equations.
    9. general numerical libs NAG (many languages in particular C,Fortran,C++) , GNU scientific library, Harwell Subroutine Library (fortran) etc
    10. FFTW- (FFT in the West) fast fourier and related transforms ©
    11. QUADPACK- numerical integration in one dimension (octave use it in the function quad())
    12. Intel MKL, Intel Math Kernel Library (in C), a library of optimized math routines for science, engineering, and financial applications, written in C/C++ and Fortran. Core math functions include BLAS, LAPACK, ScaLAPACK, sparse solvers, fast Fourier transforms, and vector math.
    More details are in pdf files with the notes, see below Section 5.1 of the notes.
  8. Lect 12 (May 19,2020) Boundary value problems in 1D i.e. the 2nd order ODE with two boundary values at he ends of an interval, how to solve them...
  9. Lect 13 (May 26,2020)
  10. Lect 14 (June 2,2020) and the right Neumann one:

    u(b)=be

    We then consider a new space Vl with functions being zero at the left end (Dirichlet bnd condition - in general in FEM- Dirichlet conditions go into the definition of the space) and the weak formulation (again integrating by parts):

    u \in Vl \int_a^b u'v' +c uvdx = \int_a^b f v dx +u'(b)v(b) for all v in Vl

    New linear FEM space (functions that are zero at the left end only):

    Vlh=\{u in C: u(a)=0,u linear on [x_k,x_{k+1}]\}

    and variational problem:

    u \in Vlh \int_a^b uh'v' +c uhvdx = \int_a^b f v dx +be*v(b) for all v in Vlh

    Nodal basis (gk) k=1,...,N (one more extra nodal function relatedto the end b!) and we get the system

    Al * U =Fl

    Al has an extra row and column (but is still symmetric and positive definite and tridiagonal - problem) Fl has also an extra entry i.e Fl=[F;Fl(N)] the last one namely: Fl(N)= \int_a^b f g_N dx + be*g_N(b) which can be approximated (integral by trapezoidal rule on [x_{N-1},x_N]) as 0.5f(b)(x_N-x_{N-1})+be in case of equaidistant mesh: 0.5*f(b)h +be - we get that this equation (the whole one not just the rhs..) is equivalent to the FDM approximation equation of Neumann cond by the backward difference + correction by utilizing the ODE equation at b...
  11. Other classical FEM spaces: we can get a higher order approximations of both problems by just exchanging discrete spaces, e.g. for zero D. bnd. conditions we can take a space of continuous piecewise quadratic functions which are zero at the ends getting classical quadratic FEM space and then method... The unknowns are usually taken as the values at the nodes and at the midpoints. The order of convergence is higher than the one of the linear FEM in an apropriate norm if the solution is smooth enough. Naturally we have to use a quadrature rule which utilizes all respective nodal points, i.e. the nodes x_k and midpoints 0.5(x_k+x_{k+1}). The most natural seems the composite Simpson rule.
  12. Lect 15 (June 9,2020) Some miscalanous things connected with numerical code writing:
  13. Some notes concerning the lecture: SC19_20.pdf

    Office hours:

    Plan. (room 5010 - IV floor). - during classes only, otherwise please contact me by e-mail

    Homeworks - lab problems

    Each lab problem will have a due date. Lab short problems which are late (showed me after the due date), will obtain the graded no of points multiplied by 0.8 (if late up to 2 labs; or 0.5 if later). I will naturally postpone all dedlines of the home problems due to the extraordinary situation.

    Next lab

    Octave scripts or source C files

    Projects

    References

    1. Piotr Krzy�anowski, Obliczenia in�ynierskie i naukowe. Skuteczne, szybkie, efektowne, Wydawnictwo Naukowe PWN, 2011
    2. Piotr Krzy�anowski, Obliczenia naukowe. Html lecture notes (in Polish). Pdf version is also available there. University of Warsaw. 2010.

    Lab