#+SETUPFILE: ../packages/org-themes/src/readtheorg_inline/readtheorg_inline.theme

#+STARTUP: org-pretty-entities entitiespretty
#+PROPERTY: header-args :comments no
#+OPTIONS: timestamp:nil
#+OPTIONS: h:3 num:t toc:3
#+TITLE: The pugs user manual
#+OPTIONS: author:nil date:nil
#+OPTIONS: tex:t

#+LANGUAGE: en
#+ATTR_LATEX: :width 4cm
#+HTML_HEAD_EXTRA: <style> pre.src-pugs:before { content: 'pugs'; } </style>
#+HTML_HEAD_EXTRA: <style> pre.src-pugs-error:before { content: 'invalid pugs'; } </style>
#+HTML_HEAD_EXTRA: <style> .footref{ color: #2980b9; font-size: 100%; }  </style>
#+LATEX_CLASS_OPTIONS: [10pt]
#+LATEX_HEADER: \usepackage[hmargin=2.5cm,vmargin=1.5cm]{geometry}
#+LATEX_HEADER: \usepackage{ae,lmodern}
#+LATEX_HEADER: \usepackage[OT1]{fontenc}
#+LATEX_HEADER: \usepackage{booktabs}
#+LATEX_HEADER: \usepackage{siunitx}
#+LATEX_HEADER: \usepackage[table]{xcolor}
#+LATEX_HEADER: \usepackage{colortbl}
#+LATEX_COMPILER: pdflatex --shell-escape

#+LATEX_HEADER_EXTRA: \usepackage{amsmath}
#+LATEX_HEADER_EXTRA: \usepackage{amsthm}
#+LATEX_HEADER_EXTRA: \usepackage{amssymb}
#+LATEX_HEADER_EXTRA: \usepackage{xcolor}
#+LATEX_HEADER_EXTRA: \usepackage{mdframed}

#+LATEX_HEADER_EXTRA: \BeforeBeginEnvironment{tabular}{\rowcolors[]{2}{orange!5}{orange!10}}
#+LATEX_HEADER_EXTRA: \BeforeBeginEnvironment{minted}{\begin{mdframed}[linecolor=blue,backgroundcolor=blue!5]}
#+LATEX_HEADER_EXTRA: \AfterEndEnvironment{minted}{\end{mdframed}}
#+LATEX_HEADER_EXTRA: \BeforeBeginEnvironment{verbatim}{\begin{mdframed}[linecolor=gray,backgroundcolor=green!5]}
#+LATEX_HEADER_EXTRA: \AfterEndEnvironment{verbatim}{\end{mdframed}}
#+LATEX_HEADER_EXTRA: \newtheorem{note}{Note}
#+LATEX_HEADER_EXTRA: \BeforeBeginEnvironment{note}{\begin{mdframed}[linecolor=orange,backgroundcolor=orange!5]}
#+LATEX_HEADER_EXTRA: \AfterEndEnvironment{note}{\end{mdframed}}
#+LATEX_HEADER_EXTRA: \newtheorem{warning}{Warning}
#+LATEX_HEADER_EXTRA: \BeforeBeginEnvironment{warning}{\begin{mdframed}[linecolor=red,backgroundcolor=red!5]}
#+LATEX_HEADER_EXTRA: \AfterEndEnvironment{warning}{\end{mdframed}}

#+LATEX_HEADER_EXTRA: \usepackage{mathpazo}

* Introduction

~pugs~[fn:pugs-def] is a general purpose solver collection built to
approximate solutions of partial differential equations. It is mainly
(but not only) designed to deal with finite-volume methods.

~pugs~ is a parallel software that uses ~MPI~[fn:MPI-def] and multi-threading
for parallelism. Multi-threading is achieved through an encapsulation
of some [[https:github.com/kokkos/kokkos][Kokkos]] mechanisms.

The philosophy of ~pugs~ is to provide "simple" numerical tools that are
assembled together through a high-level language (a DSL[fn:DSL-def]
close to the mathematics) to build more complex solvers. This approach
is inspired by the success of [[http://freefem.org][FreeFEM]], which use a similar approach.

Before detailing the leading concepts and choices that we have made to
develop ~pugs~, we give a simple example.

** An example
For instance, the following code builds a uniform Cartesian grid of
$]-1,1[^2$ made of $20\times20$ cells and transforms it by displacing the
mesh nodes according to a user-defined vector field $T: \mathbb{R}^2
\to \mathbb{R}^2$.
#+NAME: introduction-example
#+BEGIN_SRC pugs :exports both :results output
  import mesh;
  import math;
  import writer;

  let pi:R, pi = acos(-1);
  let theta:R^2 -> R, x -> 0.5*pi*(x[0]*x[0]-1)*(x[1]*x[1]-1);
  let M:R^2 -> R^2x2, x -> [[cos(theta(x)), -sin(theta(x))],
                            [sin(theta(x)),  cos(theta(x))]];
  let T: R^2 -> R^2, x -> x + M(x)*x;
  let m:mesh, m = cartesian2dMesh([-1,-1], [1,1], (20,20));

  m = transform(m, T);

  write_mesh(gnuplot_writer("transformed", 0), m);
#+END_SRC
The example is quite easy to read.
- First, some *modules* are loaded: the ~mesh~ module, which contains some
  mesh manipulation functions. The ~math~ module provides a set of
  classical mathematical functions ($\sin$, $\cos$, ...). The ~writer~
  module is used to save meshes or discrete functions to files using
  various formats.
- The second block of data defines variables of different types
  - ~pi~ is a real value that is initialized by an approximation of $\pi$.
  - ~theta~ is the real value function $\theta$ defined by
    \begin{equation*}
      \theta: \mathbb{R}^2 \to \mathbb{R},\quad\mathbf{x} \mapsto
      \frac{\pi}{2} (x_0^2-1)(x_1^2-1)
    \end{equation*}
  - ~M~ is the $\mathbb{R}^{2\times2}$  matrix field $M$ defined by
    \begin{equation*}
      M: \mathbb{R}^2 \to \mathbb{R}^{2\times2},\quad\mathbf{x} \mapsto
      \begin{pmatrix}
        \cos(\theta(\mathbf{x})) & -\sin(\theta(\mathbf{x}))  \\
        \sin(\theta(\mathbf{x}))  & \cos(\theta(\mathbf{x}))
      \end{pmatrix}
    \end{equation*}
  - ~T~ is the vector field $T$ defined by
    \begin{equation*}
      T: \mathbb{R}^2 \to \mathbb{R}^2, \quad\mathbf{x} \mapsto (I+M(\mathbf{x}))\mathbf{x}
    \end{equation*}
  - Finally ~m~ is defined as the uniform Cartesian mesh grid of
    $]-1,1[^2$. The last argument: ~(20,20)~ sets the number of cells in
    each direction: $20$.
- The third block is the application of the transformation $T$ to the
  mesh. Observe that if the resulting mesh is stored in the same
  variable ~m~, the old one was not modified in the process. It is
  important to already have in mind that the ~pugs~ language *does not
  allow* the modifications of values of variables of *non-basic*
  types. This is discussed in the section [[high-level-types]].
- Finally, the last block consists in saving the obtained mesh in a
  ~gnuplot~ file. The result is shown on Figure [[fig:intro-example]].

#+NAME: intro-transform-mesh-img
#+BEGIN_SRC gnuplot :exports results :file (substitute-in-file-name "${PUGS_SOURCE_DIR}/doc/intro-transform-mesh.png")
  reset
  unset grid
  unset border
  unset key
  unset xtics
  unset ytics
  set terminal png truecolor enhanced
  set size square
  plot '<(sed "" $PUGS_SOURCE_DIR/doc/transformed.0000.gnu)' w l
#+END_SRC

#+CAPTION: Obtained transformed mesh
#+NAME: fig:intro-example
#+ATTR_LATEX: :width 0.38\textwidth
#+ATTR_HTML: :width 300px;
#+RESULTS: intro-transform-mesh-img

Even if this first example is very simple, some key aspects can
already be discussed.
- There is no predefined constant in ~pugs~. Here a value is provided
  for ~pi~.
- There are two kinds of variable in ~pugs~: variables of basic types
  and variable of high-level types. This two kinds of variable behave
  almost the same but one must know their differences to understand
  better the underlying mechanisms and choices that we made. See
  section [[basic-types]] and [[high-level-types]] for details.
- Also, there are two types of functions: *user-defined* functions and
  *builtin* functions. In this example, ~theta~, ~M~ and ~T~ are user-defined
  functions. All other functions (~cos~, ~cartesian2dMesh~,...) are
  builtin functions and are generally defined when importing a
  module. These functions behave similarly, one should refer to
  [[functions]] for details.

** Concepts and design

As it was already stated, ~pugs~ can be viewed as a collection of
numerical methods or utilities that can be assembled together, using a
user friendly language, to build simulation scenarios.

Utilities are tools that are often used in numerical
simulations. Examples of such utilities are mesh manipulations,
definition of initial conditions, post-processing, error
calculations,...

*** A C++ toolbox driven by a user friendly language

Numerical simulation packages are software of a particular
kind. Generally, in order to run a calculation, one has to define a
set of data and parameters. This can simply be definition of a
discretization parameter such as the mesh size. One can also specify
boundary conditions, equations of state, source terms for a specific
model. Choosing a numerical method or even more, setting the model
itself, is common in large codes.

In ~pugs~, all these "parameters" are set through a
DSL[fn:DSL-def]. Thus, when ~pugs~ is launched, it actually executes a
provided script. A ~C++~ function is associated to each instruction of
the script. The ~C++~ components of ~pugs~ are completely unaware one of
the others. ~pugs~ interpreter is responsible of data flow between the
components: it manages the data transfer between those ~C++~ components
and ensures that the workflow is properly defined.

**** Why?

In this section we motivate the choice of a language and not of a more
standard approach.

***** Data files are evil

There are lots of reasons not to use data files. By data file, we
refer to a set of options that describe physical models, numerical
methods or their settings.

- Data files are not flexible. This implies in then one hand that
  application scenarios must be known somehow precisely to reflect
  possible option combinations and in the other hand even defining a
  specific initial data may require the creation of a new option and
  the associated code (in ~C++~ for instance). \\
  Usually, the last point is addressed by adding a local interpreter
  to evaluate user functions.
- Data files may contain irrelevant information. Actually, it is quite
  common to allow to define options that are valid but irrelevant to a
  given scenario. This policy can be changed but it is generally not
  an easy task and requires more work from the user (which can be a
  good thing).
- Generally data files become rapidly obsolete. An option was not the
  right one, or its type changed to allow other contexts... This puts
  pressure on the user.
- Even worst, options meaning can depend on other
  options. Unfortunately, this happens commonly. For instance, a
  global option can change implicitly the treatment associated to
  another one. This is dangerous since writing or reading the data
  file requires an important and up to date knowledge of the code's
  internals.
- Another major difficulty when dealing with data files is to check
  the compatibility of provided options.

***** Embedded "data files" are not a solution

Using directly the general purpose language (~C~, ~C++~, ~Fortran~,...) used
to write the code can be tempting. It has the advantage that no
particular treatment is necessary to build a parser (to read data
files or a script), but it presents several drawbacks.

- The first one is probably that it allows to much freedom. While
  defining the model and numerical options, the user has generally
  access to the whole code and can change almost everything, even
  things that should not be changed.
- Again, one can easily have access to irrelevant options and it
  requires a great knowledge of the code to find important ones.
- With that regard, defining a simulation properly can be a difficult
  task. For instance, in the early developments of ~pugs~ (when it was
  just a raw ~C++~ code) it was tricky to change boundary conditions for
  coupled physics.
- Another difficulty is related to the fact that code's internal API
  is likely to change permanently in a research code. Thus valid
  constructions or setting may become rapidly obsolete.  In other
  words keeping up to date embedded "data file" might be difficult.
- Finally it requires recompilation of pieces of code (which can be
  large in some cases) even if one is just changing a simple
  parameter.

***** Benefits of a DSL

Actually, an honest analysis cannot conclude that a DSL is the
solution to all problems. However, it offers some advantages.

- First, it allows a fine control on what the user can or cannot
  perform. In some sense, it offers a chosen level of flexibility.
- It allows to structure the code in the sense that new developments
  have to be designed not only focusing on the result but also on the
  way it should be used (its interactions with the scripting language).
- In the same idea, it provides a framework that drives the desired
  principle of "do simple things and do them well".
- There are no hidden dependencies between numerical options: the DSL
  code is easier to read (than data files) and is less likely to
  become obsolete (this is not that true in early developments since
  the language itself and some concepts are still likely to change).
- The simulation scenario is *defined* by the script, it is the
  responsibility of the user to check its meaning (not to the charge
  of the code).

***** ~pugs~ language purpose

~pugs~ language is used to assemble ~C++~ tools which are well
tested. These tools should ideally be small pieces of ~C++~ code that
do one single thing and do it well.

Another purpose of the language is to allow to perform high-level
calculations. In other words, the language defines a data flow and
checks that each ~C++~ piece of code is used correctly. Since each piece
of code acts as a pure function (arguments are unchanged by calls),
the calling context is quite easy to check.

Finally it aims at simplifying the definition of new methods since
common utilities are available directly in scripts.

***** The framework: divide and conquer

~pugs~ is a research oriented software, thus generally the user is also
a developer. If this paragraph is indeed more dedicated to the
developer, by reading it, the user will have a better understanding of
the development choices and the underlying policy that the code
follows.

Actually the development framework imposed by the DSL tends to guide
writing of new methods.

- In the process of writing a *new numerical methods*, one must create
  *new functions in the language*. Moreover, if a new method is close to
  an existing one, it is generally *better* to use completely new
  underlying ~C++~ code than to patch existing methods. Starting from a
  *copy* of the existing code ~C++~ is *encouraged* for developments. This
  may sound weird since classical development guidelines encourage
  inheritance or early redesign. Actually, this policy is the result
  of the following discussion.
  - Using this approach, one *separates* clearly the *design* of a
    numerical method and the *design* of the code.
  - From the computer science point of view, early design for new
    numerical methods is generally wrong: usually one cannot
    anticipate precisely enough eventual problems or method
    corrections.
  - It is much more difficult to introduce bugs in existing methods,
    since previously validated methods are unchanged!
  - For the same reason, existing methods performances are naturally
    unchanged by new developments.
  - Also, when comparing methods, it is better to compare to the
    original existing code.
  - If the new method is abandoned or simply not kept, the existing
    code is not polluted.
  - Finally when a method is validated and ready to integrate the
    mainline sources of the code, it is easier to see differences with
    existing ones and *at this time* one can redesign the ~C++~ code
    checking that results are unchanged and performances not
    deteriorated. At this time, it is likely that the numerical method
    design is finished, thus (re)designing the source code makes more
    sense.
- Another consequence is that utilities are not be developed again and
  again.
  - This implies an /a priori/ safer code: utilities are well tested and
    validated.
  - It saves development time obviously.
  - The code of numerical methods is not polluted by environment
    instructions (data initialization, error calculation,
    post-processing,...)
  - The counterpart is somehow classical. In the one hand, the
    knowledge of existing utilities is required, this document tries
    to address a part of it. In the other hand, if the developer
    requires a new utility, a good practice is to discuss with the
    other ones to check if it could benefit to them. Then one can
    determine if it should integrate rapidly or not the main
    development branch.

***** Why not python or any other scripting language?

As it was already pointed out above, general purpose languages offer
to much freedom: it is not easy to protect data. For instance in the
~pugs~ DSL, non basic variables are constant (see paragraph
[[high-level-types]]). It is important since it prevents the user from
modifying data in inconsistent ways. Also, one must keep in mind that
constraining the expressiveness is actually a strength. As said
before, one can warranty coherence of the data, perform calculations
without paying attention to the parallelism aspects,... Observe that
it is not a limitation: if the DSL's field of application needs to be
extended, it is always possible. But these extensions should never
break the rule that a DSL must not become a general purpose language.

Providing a language close to the application reduces the gap between
the application and its expression (code): Domain Specific Languages
are made for that.

Finally, python is ugly.

*** A high-level language

Following the previous discussion, the reader should now understand
the motivations that drove the design choices that conduct to build
~pugs~ as a ~C++~ toolbox driven by a user friendly language.

#+begin_verse
Keep it simple and relevant!
#+end_verse

In the development process of an application, the easy step is
generally the implementation itself. It is even more true when the
developer feels that changes to the code are natural and that the
modifications themselves look easy. Obviously, any experienced
programmer knows that writing the code is only the first step of a
much longer process which requires rigorous tests and validations (the
enrichment of a non-regression database). Generally one also desires
to ensure that the development is used within the correct bounds which
requires to implement data verification. In a perfect world, an
up-to-date documentation of the functionality and its domain of
validity should be provided.

This is even more true when defining a language (or a DSL). Enriching
a language syntax (or grammar) is not something that must be done to
answer a specific need. It must not be done /because it is possible to
do it/!

#+begin_verse
When designing a language, the difficulty is not to offer new functionalities,\\
it is generally to decide not to offer them.\\
--- Bjarne Stroustrup, C++ conference 2021.
#+end_verse

If the grammar of ~pugs~ is still likely to be extended, it should *never*
integrate low-level instructions. Low-level instructions give too much
freedom and thus are a source of errors. Several things are already
done to forbid this kind of evolution. The constness of high-level
data is a good illustration. For instance, meshes or discrete
functions *cannot* be modified. This is not only a security to protect
the user from doing "dangerous" manipulations, but it also permits to
define high-level optimizations.
- Since meshes are constant objects, one can for instance compute
  geometric data on demand. These data are kept into memory as long as
  the mesh lives.
- Forbidding the modification of values of a discrete function ensures
  that parallel communication instructions should never appear in a
  ~pugs~ script.

Another benefit of not providing low-level instructions is that the
scripts are more easy to write and read, and it is more difficult to
write errors.

* Language

** Variables

In order to simplify the presentation, before going further, we
introduce the syntax that allows to print data to the terminal. It
follows ~C++~ streams construction for convenience. For instance
#+NAME: cout-preamble-example
#+BEGIN_SRC pugs :exports both :results output
  cout << "2+3 = " << 2+3 << "\n";
#+END_SRC
produces the following output
#+results: cout-preamble-example
The code is quite obvious for ~C++~ users, note that ~"\n"~ is the
linefeed string (there is no character type in ~pugs~, just strings).

Actually, ~cout~ is itself a variable, we will come to this later.

~pugs~ is a strongly typed language. It means that a variable *cannot*
change of type in its lifetime.


*** Declaration and affectation syntax

**** Declaration of simple variables

To declare a variable ~v~ of a given type ~V~, one writes
#+BEGIN_SRC pugs :exports source
  let v:V;
#+END_SRC

This instruction is read as
#+begin_verse
Let $v\in V$.
#+end_verse

Actually,
- ~let~ is the declaration keyword,
- ~v~ is the variable name,
- ~:~ is a separation token (it can be read as "/in/" in this context),
- ~V~ is the identifier of the type, and
- ~;~ marks the end of the instruction.

For instance to declare a real variable ~x~, one writes
#+NAME: declare-R
#+BEGIN_SRC pugs :exports both :results none
  let x:R;
#+END_SRC
The instructions that follow the declaration of a variable can use it
while defined in the same scope, but it cannot be used before. Also,
after its declaration, one cannot declare another variable with the
same name in the same scope (see [[blocks-and-life-time]]).
#+NAME: redeclare-variable
#+BEGIN_SRC pugs-error :exports both :results output
  let x:R; // first declaration
  let x:R; // second declaration
#+END_SRC
produces the following error
#+results: redeclare-variable

**** Affectation of simple variables

To affect the value of an expression ~expression~ to variable ~v~ one
simply uses the ~=~ operator. Thus assuming that a variable ~v~ has
already been /declared/, one writes simply
#+BEGIN_SRC pugs :exports source
  v = expression;
#+END_SRC
There is not to much to comment, reading is quite natural
- ~v~ is the variable name,
- ~=~ is the affectation operator,
- ~expression~ is some code that provides a value of the same type as ~v~
  (it can be another variable, an arithmetic expression, the result of
  a function,...), and
- ~;~ marks the end of the instruction.

For instance, one can write
#+NAME: simple-affectations-example
#+BEGIN_SRC pugs :exports both :results output
  import math; // to load the sin function

  let a:N;
  let b:N;
  let x:R;
  a=3;
  b=2+a;
  x=sin(b)+a;

  cout << "a = " << a << " b = " << b << " x = " << x << "\n";
#+END_SRC
In this example, we import the ~math~ module which provides the ~sin~
function and we use the ~N~ data type of natural integers
($\mathbb{N}\equiv\mathbb{Z}_{\ge0}$).

Running the example gives the following result.
#+results: simple-affectations-example

#+BEGIN_warning
Actually, *separating* the *declaration* from the *initialization* of a
variable is quite *dangerous*. This is prone to errors and can lead to
the use of undefined values. During the compilation of scripts, ~pugs~
detects uninitialized values. For instance,
#+NAME: uninitialized-variable
#+BEGIN_SRC pugs-error :exports both :results output
  let x:R;
  let y:R;
  y=x+1;
#+END_SRC
produces the following compilation error
#+results: uninitialized-variable

For more complex constructions, it can be very difficult to detect it
at /compile time/, this is why it is *encouraged* to use variable
definition (see [[definition-simple-variables]]).

Observe nonetheless that ~pugs~ checks at /run time/ that used variables
are correctly defined. If not, a /runtime/ error is produced.
#+END_warning

**** Definition of simple variables<<definition-simple-variables>>

The best way to define variables is the following.  To define a
variable ~v~ of a given type ~V~, from an expression one writes
#+BEGIN_SRC pugs :exports source
  let v:V, v = expression;
#+END_SRC
- ~let~ is the declaration keyword,
- ~v~ is the variable name,
- ~:~ is a separation token (can be read as "/in/" in this context),
- ~V~ is the identifier of the type,
- ~,~ is the separator that indicates the beginning of the affectation,
- ~expression~ is some code that provides a value of type ~V~ (it can be
  another variable, an expression, the result of a function,...), and
- ~;~ marks the end of the instruction (here the definition of ~v~).

A practical example is
#+NAME: simple-definition-example
#+BEGIN_SRC pugs :exports both :results output
  import math; // to load the sin function

  let a:N, a=3;
  let b:N, b=2+a;
  let x:R, x=sin(b)+a;

  cout << "a = " << a << " b = " << b << " x = " << x << "\n";
#+END_SRC
which produces the result
#+results: simple-definition-example

*** Blocks and lifetime of variables<<blocks-and-life-time>>

In pugs scripts, variables have a precise lifetime. They are defined
within scopes. The main scope is implicitly defined and sub-scopes are
enclosed between curly bracket pairs: ~{~ and ~}~. Following ~C++~, a
variable exists (can be used) as soon as it has been declared and
until the end of the scope (where it has been declared). At this point
we give a few examples.

**** A variable cannot be used before its declaration
#+NAME: undeclare-variable
#+BEGIN_SRC pugs-error :exports both :results output
  n = 3;
  let n:N;
#+END_SRC
#+results: undeclare-variable

**** A variable cannot be used after its definition scope
#+NAME: out-of-scope-variable-use
#+BEGIN_SRC pugs-error :exports both :results output
  {
    let n:N, n = 3;
    n = n+2;
  }
  cout << n << "\n";
#+END_SRC
#+results: out-of-scope-variable-use

**** Variable name can be reused in an enclosed scope
#+NAME: nested-scope-variable-example
#+BEGIN_SRC pugs :exports both :results output
  let n:N, n = 0; // global variable
  {
    cout <<  "global scope n = " << n << "\n";
    let n:N, n = 1;  // scope level 1 variable
    {
      cout <<  "scope level 1 n = " << n << "\n";
      let n:N, n = 2;  // scope level 2 variable
      cout <<  "scope level 2 n = " << n << "\n";
    }
    {
      cout <<  "scope level 1 n = " << n << "\n";
      let n:N, n = 4;  // scope level 2.2 variable
      cout <<  "scope level 2.2 n = " << n << "\n";
    }
    cout <<  "scope level 1 n = " << n << "\n";
  }
  cout << "global scope n = " << n << "\n";
#+END_SRC
#+results: nested-scope-variable-example
This example is self explanatory. Obviously such constructions are
generally *bad ideas*. This kind of constructions can appear in loops
where the variables defined in blocks follow the same lifetime rules.

*** Basic types<<basic-types>>

Basic types in ~pugs~ are boolean (~B~), natural integers (~N~), integers (~Z~),
real (~R~), small vectors (~R^1~, ~R^2~ and ~R^3~), small matrices (~R^1x1~, ~R^2x2~
and ~R^3x3~) and strings (~string~).

#+BEGIN_note
Observe that while mathematically, obviously $\mathbb{R} = \mathbb{R}^1
= \mathbb{R}^{1\times1}$, the data types ~R~, ~R^1~ and ~R^1x1~ are different
in ~pugs~ and are *not implicitly* convertible from one to the other!

This may sound strange but there are few reasons for that.
- First, these are the reflect of internal ~pugs~ ~C++~-data types that
  are used to write algorithms. In its core design pugs aim at writing
  numerical methods generically with regard to the dimension. One of
  the ingredients to achieve this purpose is to use dimension $1$
  vectors and matrices when some algorithms reduce to dimension $1$
  instead of ~double~ values. To avoid ambiguity that may arise in some
  situations (this can lead to very tricky code), we decided to forbid
  automatic conversions of these types with ~double~. When designing the
  language, we adopted the same rule to avoid ambiguity.
- A second reason is connected to the first one. Since ~pugs~ aims at
  providing numerical methods for problems in dimension $1$, $2$ or
  $3$, this allows to distinguish the nature of the underlying objects.
  - It is natural to consider that the coordinates of the vertices
    defining a mesh in dimension $d$ are elements of $\mathbb{R}^d$,
  - or that a velocity or a displacement are also defined as
    $\mathbb{R}^d$ values.
  Thus using ~R^1~ in dimension $1$ for this kind data precise their
  nature in some sense .
#+END_note

**** Expression types

The language associates statically some types to several special
expressions.
- Special boolean (type ~B~) expressions. Two *keywords* allow to define
  boolean values: ~true~ and ~false~.
- Integers (type ~Z~) are defined as a contiguous list of digits.
  - For instance, the code ~0123~ is interpreted as the integer $123$.
  - However the sequence ~123 456~ is *not interpreted* as $123456$ but as
    the two integers $123$ and $456$.
- Real values (type ~R~) use the same syntax as in ~C++~. For instance,
  the following expressions are accepted to define the number $1.23$.
#+BEGIN_SRC pugs :exports both
  1.23;
  0.123E1;
  .123e+1;
  123e-2;
  12.3E-1;
#+END_SRC
- Small vectors values (types ~R^1~, ~R^2~ and ~R^3~) are written through
  square brackets. Each component of vector values must be a scalar
  expression.
#+BEGIN_SRC pugs :exports both
  [1];     // R^1 value
  [1,2];   // R^2 value
  [1,2,3]; // R^3 value
#+END_SRC
- Small matrices values (types ~R^1x1~, ~R^2x2~ and ~R^3x3~) are written by
  lines through square brackets. Each line is enclosed between square
  brackets.
#+BEGIN_SRC pugs :exports both
  [[1]];                     // R^1x1 value
  [[1,2],[3,4]];             // R^2x2 value
  [[1,2,3],[4,5,6],[7,8,9]]; // R^3x3 value
#+END_SRC
- ~string~ values are defined as the set of characters enclosed between
  two double quotes ( ~"~ ). The string /Hello world!/ would be simply
  written as ~"Hello world!"~. Strings support the following escape
  sequences (similarly to ~C++~):
  #+LATEX: \definecolor{contiYellow}{RGB}{255,165,0}
  #+LATEX: \rowcolors[]{2}{contiYellow!5}{contiYellow!20}
  | escape | meaning         |
  |--------+-----------------|
  | ~\'~     | single quote    |
  | ~\"~     | double quote    |
  | ~\?~     | question mark   |
  | ~\\~     | backslash       |
  | ~\a~     | audible bell    |
  | ~\b~     | backspace       |
  | ~\f~     | new page        |
  | ~\n~     | new line        |
  | ~\r~     | carriage return |
  | ~\t~     | horizontal tab  |
  | ~\v~     | vertical tab    |
  These special characters are not interpreted by ~pugs~ itself but
  interpreted by the system when an output is created. They are just
  allowed in the definition of a ~string~.

**** Variables of basic types are stored by value

This means that each variable of basic type allocates its own memory
to store its data. This is the natural behavior for variables. It is
illustrated in the following example
#+NAME: basic-type-value-storage
#+BEGIN_SRC pugs :exports both :results output
  let a:N, a = 3;
  let b:N, b = a;

  a = 1;
  cout << "a = " << a << " b = " << b << "\n";
#+END_SRC
which produces the expected result
#+results: basic-type-value-storage
When defining ~b~, the *value* contained in ~a~ is copied to set the value
of ~b~. Thus changing ~a~'s value does not impact the variable ~b~.

**** Variables of basic types are mutable

In ~pugs~ the only variables that are mutable (their value can be
*modified*) are of basic types. Executing the following code
#+NAME: basic-type-mutable-value
#+BEGIN_SRC pugs :exports both :results output
  let a:N, a = 2;
  a += 3;
  cout << "a = " << a << "\n";
#+END_SRC
gives
#+results: basic-type-mutable-value
which is not a surprise. However, the use of the ~+=~ operator results
in the modification of the stored value. There is no copy.

Actually, this is not really important from the user point of
view. One just have to keep in mind that, as it will be depicted
after, high-level variables *are not mutable*: their values can be
*replaced* by a new ones but *cannot be modified*.

*** Implicit type conversions<<implicit-conversion>>

In order to avoid ambiguities, in ~pugs~, there is *no implicit*
conversion in general.

#+BEGIN_note
Actually, there are only three situations for which implicit type
conversion can occur: when the value is
- given as a parameter of a function,
- used as the returned value of a function, or
- used to define a tuple
This will be illustrated in section [[functions]] and [[tuples]]
#+END_note

This means that all affectations, unary operators and binary operators
are defined explicitly for supported types.

Here is a table of implicit type conversions *when allowed*.

| expected type | convertible types                              |
|---------------+------------------------------------------------|
| ~N~             | ~B~, ~Z~                                           |
| ~Z~             | ~B~, ~N~                                           |
| ~R~             | ~B~, ~N~, ~Z~                                        |
| ~R^1~           | ~B~, ~N~, ~Z~, ~R~                                     |
| ~R^2~           | ~0~ (special value)                              |
| ~R^3~           | ~0~ (special value)                              |
| ~R^1x1~         | ~B~, ~N~, ~Z~, ~R~                                     |
| ~R^2x2~         | ~0~ (special value)                              |
| ~R^3x3~         | ~0~ (special value)                              |
| ~string~        | ~B~, ~N~, ~Z~, ~R~, ~R^1~, ~R^2~, ~R^3~, ~R^1x1~, ~R^2x2~, ~R^3x3~ |

*** Operators

**** Affectation operators

In the ~pugs~ language, the affectation operators are the following.
| operator | description          |
|----------+----------------------|
| ~=~        | affectation operator |
|----------+----------------------|
| ~+=~       | increment operator   |
| ~-=~       | decrement operator   |
| ~*=~       | multiply operator    |
| ~/=~       | divide operator      |

#+BEGIN_note
It is important to note that in ~pugs~ language, affectation operators
have *no* return value. This is a notable difference with ~C~ or
~C++~. Again, this is done to avoid common mistakes that can be
difficult to address. For instance, the following ~C++~ code is valid
but does not produce the expected result.
#+BEGIN_SRC C++ :exports source
  bool b = true;
  // do things ...
  if (b=false) {
    // do other things ...
  }
#+END_SRC
Obviously the mistake is that the test should have been ~(b==false)~,
since otherwise, the conditional block is never executed (in ~C++~, the
result of an affectation is the value affected to the variable, which
is always false in that case.

This cannot happen with ~pugs~. A similar example
#+NAME: no-affectation-result
#+BEGIN_SRC pugs-error :exports both :results output
  let b:B, b=true;
  // do things
  if (b=false) {
    // do other things
  }
#+END_SRC
produces the following error
#+results: no-affectation-result

Actually, affectations are /expressions/ in ~C++~. In ~pugs~, affectations
are /instructions/.

#+END_note

***** List of defined operator ~=~ for basic types.

As already mentioned, operator ~=~ is defined for *all* types in ~pugs~ if
the right hand side expression has the *same* type as the left hand side
variable. This is true for basic types as well as for high-level
types.

We now give the complete list of supported ~=~ affectations. The lists
are sorted by type of left hand side variable.

- ~B~: boolean left hand side variable. One is only allowed to affect boolean
  values.
  | ~B =~ allowed expression type |
  |-----------------------------|
  | ~B~                           |

- ~N~: natural integer ($\mathbb{N}$ or $\mathbb{Z}_{\ge0}$) left hand side variable.
  | ~N =~ allowed expression type |
  |-----------------------------|
  | ~B~                           |
  | ~N~                           |
  | ~Z~  (for convenience)        |

- ~Z~: integer ($\mathbb{Z}$) left hand side variable.
  | ~Z =~ allowed expression type |
  |-----------------------------|
  | ~B~                           |
  | ~N~                           |
  | ~Z~                           |

- ~R~: real ($\mathbb{R}$) left hand side variable.
  | ~R =~ allowed expression type |
  |-----------------------------|
  | ~B~                           |
  | ~N~                           |
  | ~Z~                           |
  | ~R~                           |

- ~R^1~: vector of dimension 1 ($\mathbb{R}^1$) left hand side variable.
  | ~R^1 =~ allowed expression type |
  |-------------------------------|
  | ~B~                             |
  | ~N~                             |
  | ~Z~                             |
  | ~R~                             |
  | ~R^1~                           |

- ~R^2~: vector of dimension 2 ($\mathbb{R}^2$) left hand side variable.
  | ~R^2 =~ allowed expression type               |
  |---------------------------------------------|
  | ~R^2~                                         |
  | ~0~  (special value)                          |
  | list of 2 scalar (~B~, ~N~, ~Z~ or ~R~) expressions |
  An example of initialization using an $\mathbb{R}^2$ value or the special value ~0~ is
  #+NAME: affectation-to-R2-from-list
  #+BEGIN_SRC pugs :exports both :results output
    let u:R^2, u = [-3, 2.5];
    let z:R^2, z = 0;
    cout << "u = " << u << "\n";
    cout << "z = " << z << "\n";
  #+END_SRC
  which produces
  #+RESULTS: affectation-to-R2-from-list
  Observe that ~0~ is a special value and *treated as a keyword*. There is no conversion from
  integer values. For instance:
  #+NAME: R2-invalid-integer-affectation
  #+BEGIN_SRC pugs-error :exports both :results output
    let z:R^2, z = 1-1;
  #+END_SRC
  produces the compilation error
  #+results: R2-invalid-integer-affectation

- ~R^3~: vector of dimension 3 ($\mathbb{R}^3$) left hand side variable.
  | ~R^3 =~ allowed expression type               |
  |---------------------------------------------|
  | ~R^3~                                         |
  | ~0~  (special value)                          |
  | list of 3 scalar (~B~, ~N~, ~Z~ or ~R~) expressions |
  An example of initialization is
  #+NAME: affectation-to-R3-from-list
  #+BEGIN_SRC pugs :exports both :results output
    let u:R^3, u = [-3, 2.5, 1E-2];
    let z:R^3, z = 0;
    cout << "u = " << u << "\n";
    cout << "z = " << z << "\n";
  #+END_SRC
  the output is
  #+RESULTS: affectation-to-R3-from-list

- ~R^1x1~: matrix of dimensions $1\times1$ ($\mathbb{R}^{1\times1}$) left hand side variable.
  | ~R^1x1 =~ allowed expression type |
  |---------------------------------|
  | ~B~                               |
  | ~N~                               |
  | ~Z~                               |
  | ~R~                               |
  | ~R^1x1~                           |

- ~R^2x2~: matrix of dimension $2\times2$ ($\mathbb{R}^{2\times2}$) left hand side variable.
  | ~R^2x2 =~ allowed expression type             |
  |---------------------------------------------|
  | ~R^2x2~                                       |
  | ~0~  (special value)                          |
  | list of 4 scalar (~B~, ~N~, ~Z~ or ~R~) expressions |
  An example of initialization using an $\mathbb{R}^{2\times2}$ value or
  the special value ~0~ is
  #+NAME: affectation-to-R2x2-from-list
  #+BEGIN_SRC pugs :exports both :results output
    let u:R^2x2, u = [[-3, 2.5],
                      [ 4, 1.2]];
    let z:R^2x2, z = 0;
    cout << "u = " << u << "\n";
    cout << "z = " << z << "\n";
  #+END_SRC
  which produces
  #+RESULTS: affectation-to-R2x2-from-list

- ~R^3x3~: matrix of dimension $3\times3$ ($\mathbb{R}^{3\times3}$) left hand side variable.
  | ~R^3x3 =~ allowed expression type             |
  |---------------------------------------------|
  | ~R^3x3~                                       |
  | ~0~  (special value)                          |
  | list of 9 scalar (~B~, ~N~, ~Z~ or ~R~) expressions |
  An example of initialization is
  #+NAME: affectation-to-R3x3-from-list
  #+BEGIN_SRC pugs :exports both :results output
    let u:R^3x3, u = [[ -3, 2.5, 1E-2],
                      [  2, 1.7,   -2],
                      [1.2,   4,  2.3]];
    let z:R^3x3, z = 0;
    cout << "u = " << u << "\n";
    cout << "z = " << z << "\n";
  #+END_SRC
  the output is
  #+RESULTS: affectation-to-R3x3-from-list

- ~string~ left hand side variable. Expressions of any basic type can be
  used as the right hand side.
  | ~string =~ allowed expression type |
  |----------------------------------|
  | ~B~                                |
  | ~N~                                |
  | ~Z~                                |
  | ~R~                                |
  | ~R^1~                              |
  | ~R^2~                              |
  | ~R^3~                              |
  | ~R^1x1~                            |
  | ~R^2x2~                            |
  | ~R^3x3~                            |
  | ~string~                           |

***** List of defined operator ~+=~ for basic types.

- ~B~: the ~+=~ operator is not defined for left hand side boolean
  variables.

- ~N~: natural integer ($\mathbb{N}$ or $\mathbb{Z}_{\ge0}$) left hand side variable.
  | ~N +=~ allowed expression type |
  |------------------------------|
  | ~B~                            |
  | ~N~                            |
  | ~Z~  (for convenience)         |

- ~Z~: integer ($\mathbb{Z}$) left hand side variable.
  | ~Z +=~ allowed expression type |
  |------------------------------|
  | ~B~                            |
  | ~N~                            |
  | ~Z~                            |

- ~R~: real ($\mathbb{R}$) left hand side variable.
  | ~R +=~ allowed expression type |
  |------------------------------|
  | ~B~                            |
  | ~N~                            |
  | ~Z~                            |
  | ~R~                            |

- ~R^1~: vector of dimension 1 ($\mathbb{R}^1$) left hand side variable.
  | ~R^1 +=~ allowed expression type |
  |--------------------------------|
  | ~R^1~                            |

- ~R^2~: vector of dimension 2 ($\mathbb{R}^2$) left hand side variable.
  | ~R^2 +=~ allowed expression type |
  |--------------------------------|
  | ~R^2~                            |

- ~R^3~: vector of dimension 3 ($\mathbb{R}^3$) left hand side variable.
  | ~R^3 +=~ allowed expression type |
  |--------------------------------|
  | ~R^3~                            |

- ~R^1x1~: matrix of dimensions $1\times1$ ($\mathbb{R}^{1\times1}$) left hand side variable.
  | ~R^1x1 +=~ allowed expression type |
  |----------------------------------|
  | ~R^1x1~                            |

- ~R^2x2~: matrix of dimension $2\times2$ ($\mathbb{R}^{2\times2}$) left hand side variable.
  | ~R^2x2 +=~ allowed expression type |
  |----------------------------------|
  | ~R^2x2~                            |

- ~R^3x3~: matrix of dimension $3\times3$ ($\mathbb{R}^{3\times3}$) left hand side variable.
  | ~R^3x3 +=~ allowed expression type |
  |----------------------------------|
  | ~R^3x3~                            |

- ~string~ left hand side variable. Expressions of any basic type can be
  used as the right hand side.
  | ~string +=~ allowed expression type |
  |-----------------------------------|
  | ~B~                                 |
  | ~N~                                 |
  | ~Z~                                 |
  | ~R~                                 |
  | ~R^1~                               |
  | ~R^2~                               |
  | ~R^3~                               |
  | ~R^1x1~                             |
  | ~R^2x2~                             |
  | ~R^3x3~                             |
  | ~string~                            |

***** List of defined operator ~-=~ for basic types.

- ~B~: the ~-=~ operator is not defined for left hand side boolean variables.

- ~N~: natural integer ($\mathbb{N}$ or $\mathbb{Z}_{\ge0}$) left hand side variable.
  | ~N -=~ allowed expression type |
  |------------------------------|
  | ~B~                            |
  | ~N~                            |
  | ~Z~  (for convenience)         |

- ~Z~: integer ($\mathbb{Z}$) left hand side variable.
  | ~Z -=~ allowed expression type |
  |------------------------------|
  | ~B~                            |
  | ~N~                            |
  | ~Z~                            |

- ~R~: real ($\mathbb{R}$) left hand side variable.
  | ~R -=~ allowed expression type |
  |------------------------------|
  | ~B~                            |
  | ~N~                            |
  | ~Z~                            |
  | ~R~                            |

- ~R^1~: vector of dimension 1 ($\mathbb{R}^1$) left hand side variable.
  | ~R^1 -=~ allowed expression type |
  |--------------------------------|
  | ~R^1~                            |

- ~R^2~: vector of dimension 2 ($\mathbb{R}^2$) left hand side variable.
  | ~R^2 -=~ allowed expression type |
  |--------------------------------|
  | ~R^2~                            |

- ~R^3~: vector of dimension 3 ($\mathbb{R}^3$) left hand side variable.
  | ~R^3 -=~ allowed expression type |
  |--------------------------------|
  | ~R^3~                            |

- ~R^1x1~: matrix of dimensions $1\times1$ ($\mathbb{R}^{1\times1}$) left hand side variable.
  | ~R^1x1 -=~ allowed expression type |
  |----------------------------------|
  | ~R^1x1~                            |

- ~R^2x2~: matrix of dimension $2\times2$ ($\mathbb{R}^{2\times2}$) left hand side variable.
  | ~R^2x2 -=~ allowed expression type |
  |----------------------------------|
  | ~R^2x2~                            |

- ~R^3x3~: matrix of dimension $3\times3$ ($\mathbb{R}^{3\times3}$) left hand side variable.
  | ~R^3x3 -=~ allowed expression type |
  |----------------------------------|
  | ~R^3x3~                            |

- ~string~: the ~-=~ operator is not defined for left hand side string variables.

***** List of defined operator ~*=~ for basic types.

- ~B~: the ~*=~ operator is not defined for left hand side boolean variables.

- ~N~: natural integer ($\mathbb{N}$ or $\mathbb{Z}_{\ge0}$) left hand side variable.
  | ~N *=~ allowed expression type |
  |------------------------------|
  | ~B~                            |
  | ~N~                            |
  | ~Z~  (for convenience)         |

- ~Z~: integer ($\mathbb{Z}$) left hand side variable.
  | ~Z *=~ allowed expression type |
  |------------------------------|
  | ~B~                            |
  | ~N~                            |
  | ~Z~                            |

- ~R~: real ($\mathbb{R}$) left hand side variable.
  | ~R *=~ allowed expression type |
  |------------------------------|
  | ~B~                            |
  | ~N~                            |
  | ~Z~                            |
  | ~R~                            |

- ~R^1~: vector of dimension 1 ($\mathbb{R}^1$) left hand side variable.
  | ~R^1 *=~ allowed expression type |
  |--------------------------------|
  | ~B~                              |
  | ~N~                              |
  | ~Z~                              |
  | ~R~                              |

- ~R^2~: vector of dimension 2 ($\mathbb{R}^2$) left hand side variable.
  | ~R^2 *=~ allowed expression type |
  |--------------------------------|
  | ~B~                              |
  | ~N~                              |
  | ~Z~                              |
  | ~R~                              |

- ~R^3~: vector of dimension 3 ($\mathbb{R}^3$) left hand side variable.
  | ~R^3 *=~ allowed expression type |
  |--------------------------------|
  | ~B~                              |
  | ~N~                              |
  | ~Z~                              |
  | ~R~                              |

- ~R^1x1~: matrix of dimensions $1\times1$ ($\mathbb{R}^{1\times1}$) left hand side variable.
  | ~R^1x1 *=~ allowed expression type |
  |----------------------------------|
  | ~B~                                |
  | ~N~                                |
  | ~Z~                                |
  | ~R~                                |

- ~R^2x2~: matrix of dimension $2\times2$ ($\mathbb{R}^{2\times2}$) left hand side variable.
  | ~R^2x2 *=~ allowed expression type |
  |----------------------------------|
  | ~B~                                |
  | ~N~                                |
  | ~Z~                                |
  | ~R~                                |

- ~R^3x3~: matrix of dimension $3\times3$ ($\mathbb{R}^{3\times3}$) left hand side variable.
  | ~R^3x3 *=~ allowed expression type |
  |----------------------------------|
  | ~B~                                |
  | ~N~                                |
  | ~Z~                                |
  | ~R~                                |

#+BEGIN_note
Observe that for these small matrix types ($\mathbb{R}^{d\times d}$) the
construction ~A *= B;~ where ~B~ is a matrix of the same type as ~A~ is not
allowed. The main reason for that is that for $d>1$ this operation has
no interests since it requires a temporary. One will see bellow that
it is possible to write ~A = A*B;~ if needed.
#+END_note

- ~string~: the ~*=~ operator is not defined for left hand side string variables.

***** List of defined operator ~/=~ for basic types.

- ~B~: the ~/=~ operator is not defined for left hand side boolean variables.

- ~N~: natural integer ($\mathbb{N}$ or $\mathbb{Z}_{\ge0}$) left hand side variable.
  | ~N /=~ allowed expression type |
  |------------------------------|
  | ~N~                            |
  | ~Z~  (for convenience)         |

- ~Z~: integer ($\mathbb{Z}$) left hand side variable.
  | ~Z /=~ allowed expression type |
  |------------------------------|
  | ~N~                            |
  | ~Z~                            |

- ~R~: real ($\mathbb{R}$) left hand side variable.
  | ~R /=~ allowed expression type |
  |------------------------------|
  | ~N~                            |
  | ~Z~                            |
  | ~R~                            |

- ~R^d~: the ~/=~ operator is not defined for left hand side vector (of
  dimension $d\in\{1,2,3\}$) variables.

- ~R^dxd~: the ~/=~ operator is not defined for left hand side matrix (of
  dimension $d\times d$ with $d\in\{1,2,3\}$) variables.

**** Unary operators

The ~pugs~ language allows the following tokens as unary operators
| operator | description          |
|----------+----------------------|
| ~not~      | not operator         |
| ~+~        | plus unary operator  |
| ~-~        | minus unary operator |
|----------+----------------------|
| ~++~       | increment operator   |
| ~--~       | decrement operator   |
|----------+----------------------|
| ~[]~       | access operator      |

The ~not~, ~+~ and ~-~ operators apply to the *expression* on their right. ~++~
and ~--~ operators apply only to a *variable* that can be positioned
before (pre increment/decrement) or after the token (post
increment/decrement). These operators are also inspired from their ~C++~
counterparts for commodity.

The ~+~ unary operator is a convenient operator that is *elided* when
parsing the script.

For basic types, when the operators ~not~, ~+~, ~-~, ~++~ or ~--~ are defined,
they return a value of the same type as the argument (except for the
operator ~-~ if the argument is a ~B~ or a ~N~, then the result is a
~Z~). These operators can be defined for high-level types.
- The ~not~ operator is only defined for boolean values (~B~).
- The ~-~ unary operator is defined for numeric basic types: ~B~,
  ~N~, ~Z~, ~R~, ~R^1~, ~R^2~, ~R^3~, ~R^1x1~, ~R^2x2~ and ~R^3x3~. It is not defined
  for ~string~ variables.
- Pre and post increment operators, ~--~ and ~++~, are defined for integer
  types: ~N~ and ~Z~. They are not defined for ~B~, ~R~, ~R^1~, ~R^2~, ~R^3~, ~R^1x1~,
  ~R^2x2~, ~R^3x3~ and ~string~ variables.

Note that the pre increment/decrement operators behave slightly
differently than their ~C++~ counterparts since they are not allowed to
be chained. In ~C++~, the following code is allowed
#+BEGIN_SRC C++ :exports source
  int i = 0;
  int j = ++ ++i;
#+END_SRC
In ~pugs~, it is forbidden:
#+NAME: double-pre-incr-result
#+BEGIN_SRC pugs-error :exports both :results output
  let i:N, i=0;
  let j:N, j = ++ ++i;
#+END_SRC
produces the compilation error
#+results: double-pre-incr-result
Again, this is done to simplify the syntax and to avoid weird
constructions.

- Access operators are only defined for small vectors ~R^d~ and small
  matrices ~R^dxd~. To avoid use of uninitialized variables (or
  partially uninitialized variables), these are ~read-only~ access
  operators. Their syntax is the following.
#+NAME: Rd-Rdxd-access-operator
#+BEGIN_SRC pugs :exports both :results output
  let x:R^2,   x = [1,2];
  let A:R^3x3, A = [[1,2,3],[4,5,6],[7,8,9]];

  cout << "x[0] = " << x[0] << "\nx[1] = " << x[1] << "\n";
  cout << "A[0,0] = " << A[0,0] << "\nA[2,1] = " << A[2,1] << "\n";
#+END_SRC
This code produces
#+results: Rd-Rdxd-access-operator

**** Binary operators

Syntax for binary operators follows again a classical structure: if
~exp1~ and ~exp2~ denotes two expressions and if ~op~ denotes a binary
operator, one simply writes ~exp1 op exp2~.

Here is the list of binary operators
| keyword | operator              |
|---------+-----------------------|
| ~and~     | logic and             |
| ~or~      | logic or              |
| ~xor~     | logic exclusive or    |
|---------+-----------------------|
| ~==~      | equality              |
| ~!=~      | non-equality          |
| ~<~       | lower than            |
| ~<=~      | lower or equal than   |
| ~>~       | greater than          |
| ~>=~      | greater or equal than |
|---------+-----------------------|
| ~<<~      | shift left            |
| ~>>~      | shift right           |
|---------+-----------------------|
| ~+~       | sum                   |
| ~-~       | difference            |
| ~*~       | product               |
| ~/~       | division              |

Binary operators can be defined for high-level types. For basic types,
they follow a few rules.

- Logical operators ~and~, ~or~ and ~xor~ are defined for boolean operands
  (type is ~B~) only. The result of the expression is a boolean.
  #+begin_src latex :results drawer :exports results
    \begin{equation*}
      \left|
        \begin{array}{rl}
          \mathtt{and}:&\quad \mathbb{B} \times \mathbb{B} \to \mathbb{B}\\
          \mathtt{or}:& \quad\mathbb{B} \times \mathbb{B} \to \mathbb{B}\\
          \mathtt{xor}:& \quad \mathbb{B} \times \mathbb{B} \to \mathbb{B}
        \end{array}
      \right.
    \end{equation*}
  #+end_src
- Comparison operators ~==~, ~!=~, ~<~, ~<=~, ~>~ and ~>=~ are defined for all
  basic scalar type and return a boolean value.
  #+begin_src latex :results drawer :exports results
    \begin{equation*}
      \forall \mathbb{S}_1, \mathbb{S}_2 \in \{\mathbb{B},\mathbb{N},\mathbb{Z},\mathbb{R}\},
      \quad
      \left|
        \begin{array}{rl}
          \mathtt{==}:& \mathbb{S}_1 \times \mathbb{S}_2 \to \mathbb{B}\\
          \mathtt{!=}:& \mathbb{S}_1 \times \mathbb{S}_2 \to \mathbb{B}\\
          \mathtt{<}: & \mathbb{S}_1 \times \mathbb{S}_2 \to \mathbb{B}\\
          \mathtt{<=}:& \mathbb{S}_1 \times \mathbb{S}_2 \to \mathbb{B}\\
          \mathtt{>}: & \mathbb{S}_1 \times \mathbb{S}_2 \to \mathbb{B}\\
          \mathtt{>=}:& \mathbb{S}_1 \times \mathbb{S}_2 \to \mathbb{B}
        \end{array}
      \right.
    \end{equation*}
  #+end_src
  When comparing a boolean value (type ~B~) with another scalar value
  type (~N~, ~Z~ or ~R~), the value ~true~ is interpreted as $1$ and the value
  ~false~ as $0$.
\\
  For vector and matrix basic types, the only allowed operators are ~==~
  and ~!=~.
  #+begin_src latex :results drawer :exports results
    \begin{equation*}
      \forall d \in \{1,2,3\},\quad
      \left|
        \begin{array}{rl}
          \mathtt{==}:& \mathbb{R}^d \times \mathbb{R}^d \to \mathbb{B}\\
          \mathtt{==}:& \mathbb{R}^{d \times d} \times \mathbb{R}^{d \times d} \to \mathbb{B}\\
          \mathtt{!=}:& \mathbb{R}^d \times \mathbb{R}^d \to \mathbb{B}\\
          \mathtt{!=}:& \mathbb{R}^{d \times d} \times \mathbb{R}^{d \times d} \to \mathbb{B}
        \end{array}
      \right.
    \end{equation*}
  #+end_src
\\
  This is also the case for ~string~ values: only allowed operators are
  ~==~ and ~!=~.
  #+begin_src latex :results drawer :exports results
    \begin{equation*}
      \left|
        \begin{array}{rl}
          \mathtt{==}:& \mathtt{string} \times \mathtt{string} \to \mathbb{B}\\
          \mathtt{!=}:& \mathtt{string} \times \mathtt{string} \to \mathbb{B}
        \end{array}
      \right.
    \end{equation*}
  #+end_src

- Shift operators ~<<~ and ~>>~ are not used to define binary operators
  between two basic types.

- Arithmetic operators (~+~, ~-~, ~*~ and ~/~) are defined for a set of
  combinations of basic types. We classify them by their returned
  types.

  - Their is no arithmetic operation that returns a boolean ~B~.

  - Operators that return a natural integer ~N~.
    #+begin_src latex :results drawer :exports results
      \begin{equation*}
        \forall \mathbb{S}_1, \mathbb{S}_2 \in \{\mathbb{B},\mathbb{N}\},
        \quad
        \left|
          \begin{array}{rl}
            \mathtt{+}:& \mathbb{S}_1 \times \mathbb{S}_2 \to \mathbb{N}\\
            \mathtt{*}:& \mathbb{S}_1 \times \mathbb{S}_2 \to \mathbb{N}
          \end{array}
        \right.
      \end{equation*}
    #+end_src
    Boolean values (type ~B~) are not allowed as the right operand of a
    division.
    #+begin_src latex :results drawer :exports results
      \begin{equation*}
        \forall \mathbb{S} \in \{\mathbb{B},\mathbb{N}\},
        \mathtt{/}: \mathbb{S} \times \mathbb{N} \to \mathbb{N}
      \end{equation*}
    #+end_src

    Observe that ~-~ is *not* in the list.

  - Operators that return an integer ~Z~.
    #+begin_src latex :results drawer :exports results
      \begin{equation*}
        \forall \mathbb{S}_1, \mathbb{S}_2 \in \{\mathbb{B},\mathbb{N},\mathbb{Z}\},
        \mbox{ such that }\mathbb{S}_1 = \mathbb{Z}\mbox{ or }\mathbb{S}_2 = \mathbb{Z},
        \quad
        \left|
          \begin{array}{rl}
            \mathtt{+}:& \mathbb{S}_1 \times \mathbb{S}_2 \to \mathbb{Z}\\
            \mathtt{-}:& \mathbb{S}_1 \times \mathbb{S}_2 \to \mathbb{Z}\\
            \mathtt{*}:& \mathbb{S}_1 \times \mathbb{S}_2 \to \mathbb{Z}\\
            \mathtt{/}:& \mathbb{S}_1 \times \mathbb{S}_2 \to \mathbb{Z}
          \end{array}
        \right.
      \end{equation*}
    #+end_src
    The divide operator does not allow boolean values as a right
    operand.
    #+begin_src latex :results drawer :exports results
      \begin{align*}
        \forall \mathbb{S}_1 \in \{\mathbb{B},\mathbb{N},\mathbb{Z}\}&\mbox{ and } \mathbb{S}_2 \in \{\mathbb{N},\mathbb{Z}\},
        \mbox{ such that }\mathbb{S}_1 = \mathbb{Z}\mbox{ or }\mathbb{S}_2 = \mathbb{Z},\\
        &\mathtt{/}: \mathbb{S}_1 \times \mathbb{S}_2 \to \mathbb{Z}.
      \end{align*}
    #+end_src
    Finally the following  operator is also defined
    #+begin_src latex :results drawer :exports results
      \begin{equation*}
        \mathtt{-}: \mathbb{N} \times \mathbb{N} \to \mathbb{Z}
      \end{equation*}
    #+end_src
    The difference of two natural integers returns a /signed/ integer.

  - Operators that return a real ~R~.
    #+begin_src latex :results drawer :exports results
      \begin{equation*}
        \forall \mathbb{S}_1, \mathbb{S}_2 \in \{\mathbb{B},\mathbb{N},\mathbb{Z},\mathbb{R}\},
        \mbox{ such that }\mathbb{S}_1 = \mathbb{R}\mbox{ or }\mathbb{S}_2 = \mathbb{R},
        \quad
        \left|
          \begin{array}{rl}
            \mathtt{+}:& \mathbb{S}_1 \times \mathbb{S}_2 \to \mathbb{R}\\
            \mathtt{-}:& \mathbb{S}_1 \times \mathbb{S}_2 \to \mathbb{R}\\
            \mathtt{*}:& \mathbb{S}_1 \times \mathbb{S}_2 \to \mathbb{R}
          \end{array}
        \right.
      \end{equation*}
    #+end_src
    The divide operator does not allow boolean values as a right
    operand.
    #+begin_src latex :results drawer :exports results
      \begin{align*}
        \forall \mathbb{S}_1 \in \{\mathbb{B},\mathbb{N},\mathbb{Z},\mathbb{R}\}&\mbox{ and } \mathbb{S}_2 \in \{\mathbb{N},\mathbb{Z},\mathbb{R}\},
        \mbox{ such that }\mathbb{S}_1 = \mathbb{R}\mbox{ or }\mathbb{S}_2 = \mathbb{R},\\
        &\mathtt{/}: \mathbb{S}_1 \times \mathbb{S}_2 \to \mathbb{Z}.
      \end{align*}
    #+end_src

  - Operators that return a small vector ~R^d~.
    #+begin_src latex :results drawer :exports results
      \begin{equation*}
        \forall d \in\{1,2,3\},
        \quad
        \left|
          \begin{array}{rl}
            \mathtt{+}:& \mathbb{R}^d \times \mathbb{R}^d \to \mathbb{R}^d\\
            \mathtt{-}:& \mathbb{R}^d \times \mathbb{R}^d \to \mathbb{R}^d\\
            \mathtt{*}:& \mathbb{B} \times \mathbb{R}^d \to \mathbb{R}^d\\
            \mathtt{*}:& \mathbb{N} \times \mathbb{R}^d \to \mathbb{R}^d\\
            \mathtt{*}:& \mathbb{Z} \times \mathbb{R}^d \to \mathbb{R}^d\\
            \mathtt{*}:& \mathbb{R} \times \mathbb{R}^d \to \mathbb{R}^d\\
            \mathtt{*}:& \mathbb{R}^{d\times d} \times \mathbb{R}^d \to \mathbb{R}^d
          \end{array}
        \right.
      \end{equation*}
    #+end_src

  - Operators that return a small matrix ~R^dxd~.
    #+begin_src latex :results drawer :exports results
      \begin{equation*}
        \forall d \in\{1,2,3\},
        \quad
        \left|
          \begin{array}{rl}
            \mathtt{+}:& \mathbb{R}^{d\times d} \times \mathbb{R}^{d\times d} \to \mathbb{R}^{d\times d}\\
            \mathtt{-}:& \mathbb{R}^{d\times d} \times \mathbb{R}^{d\times d} \to \mathbb{R}^{d\times d}\\
            \mathtt{*}:& \mathbb{B} \times \mathbb{R}^{d\times d} \to \mathbb{R}^{d\times d}\\
            \mathtt{*}:& \mathbb{N} \times \mathbb{R}^{d\times d} \to \mathbb{R}^{d\times d}\\
            \mathtt{*}:& \mathbb{Z} \times \mathbb{R}^{d\times d} \to \mathbb{R}^{d\times d}\\
            \mathtt{*}:& \mathbb{R} \times \mathbb{R}^{d\times d} \to \mathbb{R}^{d\times d}\\
            \mathtt{*}:& \mathbb{R}^{d\times d} \times \mathbb{R}^{d\times d} \to \mathbb{R}^{d\times d}
          \end{array}
        \right.
      \end{equation*}
    #+end_src

  - Operators that return a ~string~.
    #+begin_src latex :results drawer :exports results
      \begin{equation*}
        \forall \mathbb{S} \in \{\mathbb{B},\mathbb{N},\mathbb{Z},\mathbb{R}, \mathbb{R}^1,  \mathbb{R}^2,  \mathbb{R}^3, \mathbb{R}^{1\times1}, \mathbb{R}^{2\times2}, \mathbb{R}^{3\times3}\},
        \quad
        \left|
          \begin{array}{rl}
            \mathtt{+}:& \mathbb{S} \times \mathtt{string} \to \mathtt{string}\\
            \mathtt{+}:& \mathtt{string} \times \mathbb{S} \to \mathtt{string}\\
            \mathtt{+}:& \mathtt{string} \times \mathtt{string}  \to \mathtt{string}
          \end{array}
        \right.
      \end{equation*}
    #+end_src

**** Operator precedence and associativity

To avoid confusions, the operators precedence in ~pugs~ language follows
the same rules as in ~C++~.

This is summarized in the following table, where ~a~ and ~b~ denotes two
expressions.
| Precedence | Operator |
|------------+----------|
|          1 | ~a++~      |
|            | ~a--~      |
|            | ~a[]~      |
|------------+----------|
|          2 | ~++a~      |
|            | ~--a~      |
|            | ~not a~    |
|------------+----------|
|          3 | ~a*b~      |
|            | ~a/b~      |
|------------+----------|
|          4 | ~a+b~      |
|            | ~a-b~      |
|------------+----------|
|          5 | ~a<<b~     |
|            | ~a>>b~     |
|------------+----------|
|          6 | ~a<b~      |
|            | ~a<=b~     |
|            | ~a>b~      |
|            | ~a>=b~     |
|------------+----------|
|          7 | ~a==b~     |
|            | ~a!=b~     |
|------------+----------|
|          8 | ~a=b~      |
|            | ~a+=b~     |
|            | ~a-=b~     |
|            | ~a*=b~     |
|            | ~a/=b~     |

As already said, we forbid some constructions to try to avoid spurious
constructions. By construction associativity of many operators makes
no sense in ~pugs~ language (affectations and increment/decrement
operators for instance).

For all other operators, the associativity rule is *left to right*. Thus
the following code
#+BEGIN_SRC pugs :exports source
  - 1 + 3 - 4 + 2;
#+END_SRC
is equivalent to
#+BEGIN_SRC pugs :exports source
  (((- 1) + 3) - 4) + 2;
#+END_SRC

Obviously ~pugs~ allows the use of parenthesis to write
expressions. It is enough to give a simple example.
#+NAME: parenthesis-precedence
#+BEGIN_SRC pugs :exports both :results output
  cout << " 2 + 3  * 4 = " <<  2 + 3  * 4 << "\n";
  cout << "(2 + 3) * 4 = " << (2 + 3) * 4 << "\n";
#+END_SRC
the output is
#+RESULTS: parenthesis-precedence

*** High-level types<<high-level-types>>

Aside from the basic types described in the previous section, ~pugs~
also deals with "high-level" types. This term is more to understand as
"non-basic types". The ~pugs~ language is not object oriented to keep it
simple.

To fix ideas, let us give examples of a few high-level data
types. These can be meshes (the ~mesh~ type), output streams (the
~ostream~ type), boundary descriptors (~boundary~), quadrature formula
(~quadrature~), discrete functions (~Vh~),...

One can already see that the complexity of these types may vary a lot.

The main difference between these types and basic types is that,
high-level types are not available in directly the language but they
are loaded on demand (using the ~import~ keyword in the preamble of the
script).

The second difference is that data of these types are *constant*. More
precisely, the content of high-level variables can be replaced by a
new one *but* it cannot be modified. For this reason, the following
operators can never be applied to variables of this kind
| forbidden operators | description                 |
|---------------------+-----------------------------|
| ~++~                  | pre/post increment operator |
| ~--~                  | pre/post decrement operator |
| ~+=~                  | assignment by sum           |
| ~-=~                  | assignment by difference    |
| ~*=~                  | assignment by product       |
| ~/=~                  | assignment by quotient      |

We conclude by stating that if access operator ~[]~ can eventually be
overloaded for high-level types, it should be done with care. It is
not recommended.

Finally, the last difference lies in the fact that high-level types
use shallow copies and not value copies as it is the case for basic
types. This is transparent to the user and provides the intuitive
(similar) behavior since data of high-level variables are constant. To
illustrate this, let us consider the following example.

#+BEGIN_SRC pugs :exports both
  import mesh;

  let m1:mesh, m1 = cartesian2dMesh(0, [1,1], (10,10));
  let m2:mesh, m2 = m1;
  let m3:mesh, m3 = cartesian2dMesh(0, [1,1], (10,10));

  m3 = m1;
#+END_SRC

In this example, we are dealing with 3 ~mesh~ variables.
- First, ~m1~ is defined as a uniform cartesian mesh in dimension 2 of
  $]0,1[^2$ made of $10\times10$ identical square cells. Let us call this
  resident mesh $\mathcal{M}_a$ for clarity.
- Next, the variable ~m2~ is defined as a copy of the variable ~m1~. It is
  a copy of the variable, but its content is *the same*. The variable ~m2~
  is also referring to the mesh $\mathcal{M}_a$: there is no
  duplication of the mesh in memory.
- Then, the variable ~m3~ is defined to refer to a *new* mesh (called
  $\mathcal{M}_b$. It is /similar/ to $\mathcal{M}_a$, but it is *not* the
  same one! When ~m3~ is defined two meshes $\mathcal{M}_a$ and
  $\mathcal{M}_b$ (and two distinct connectivities) are resident in
  memory.
- Finally, the last instruction (~m3 = m1;~) sets ~m3~ to also refer
  $\mathcal{M}_a$. Since no other variable refers to $\mathcal{M}_b$,
  this mesh is destroyed (memory is freed). At the end of the program,
  all the variables ~m1~, ~m2~ and ~m3~ are referring to $\mathcal{M}_a$
  which is the only mesh that resides in memory.

*** Compound types

The ~pugs~ language allows to deal with compound types. The idea is to
define a list of variables as a member of a product space (each
variable belongs to one of them, each variable has a simple type:
basic or high-level).

**** Compound declaration

Let us provide an example to fix ideas.
#+NAME: compound-declaration
#+BEGIN_SRC pugs :exports both :results output
  let (A,x,n): R^2x2*R^3*N;
  A = [[1,2],[3,4]];
  x = [2,4,6];
  n = 2;

  cout << "A = " << A << "\nx = " << x << "\nn = " << n << "\n";
#+END_SRC
the output is
#+RESULTS: compound-declaration

This is completely equivalent to declaring the variables one after the
other.

**** Compound definition

One can also use the following definition instruction
#+NAME: compound-definition
#+BEGIN_SRC pugs :exports both :results output
  let (A,x,n): R^2x2*R^3*N, (A,x,n) = ([[1,2],[3,4]], [2,4,6], 2);

  cout << "A = " << A << "\nx = " << x << "\nn = " << n << "\n";
#+END_SRC
which gives the same result
#+RESULTS: compound-definition

A potential mistake with this construction, is that variables
must appear in the same order in the declaration part and in the
affectation part of a definition. For instance, the following code is
invalid.

#+NAME: invalid-compound-definition
#+BEGIN_SRC pugs-error :exports both :results output
  let (x,y):R*R, (y,x) = (0,1);
#+END_SRC
It produces the following error
#+results: invalid-compound-definition
which is easy to fix.

Another potential mistake is that all variables of the list are marked
as defined at the same time. Thus one cannot use construction like
this:
#+NAME: undeclared-compound-definition
#+BEGIN_SRC pugs-error :exports both :results output
  let (x,y):R*R, (x,y) = (0,2+x);
#+END_SRC
It produces the following error
#+results: undeclared-compound-definition
While the variable ~x~ is defined *before* ~y~, this kind of construction is
forbidden. From a technical point of view, this behavior would be easy
to change (allow to use the fresh value of ~x~ in the definition of ~y~),
but this make the code unclear and this is not the purpose of compound
types.

#+BEGIN_note
Observe that there is no implicit conversion when dealing with
compound types. The ~=~ operators are used sequentially to set the
different variables.
#+END_note

**** Compound affectation

The last way to use compound types is illustrated by the following
example.
#+NAME: compound-affectation
#+BEGIN_SRC pugs :exports both :results output
  let A: R^2x2;
  let x: R^3;
  let n: N;

  (A,x,n) = ([[1,2],[3,4]], [2,4,6], 2);

  cout << "A = " << A << "\nx = " << x << "\nn = " << n << "\n";
#+END_SRC
It produces again
#+results: compound-affectation

#+BEGIN_note
Observe that the only operator allowed for this kind of construction
is the operator ~=~.
#+END_note

**** Use of compound types

Actually using compound types the way it is presented in this
paragraph is not recommend. The purpose of compounds types and
compound affectations is related to functions. As one will see in
section [[functions]], functions can return compound values, thus compound
affectations (or definitions) are needed to get returned values in
that case.

*** Tuple types<<tuples>>

The last special type construction is the ability to deal with tuples
in the ~pugs~ language. The tuples we are describing here are lists of
data of a *unique* and *simple* type (one cannot use compound types to
define tuples). Tuples in the ~pugs~ language have *no relationship* with
~C++~'s ~std::tuple~ these are two completely different notions.

The list of values given to the tuple must be *implicitly* convertible
to the type of tuple elements (see the conversion table in section
[[implicit-conversion]]). There is no ambiguity, the implicit conversions
follow the rules of operator ~=~.

**** Tuple declaration and affectation

For instance, one may write
#+NAME: tuple-declaration-affectation
#+BEGIN_SRC pugs :exports both :results output
  let x:(R);
  x = (1,2,3.4,2);
  cout << x << "\n";
#+END_SRC
Executing this code, one gets
#+results: tuple-declaration-affectation

**** Tuple definition

The definition syntax is also possible.
#+NAME: tuple-definition
#+BEGIN_SRC pugs :exports both :results output
  let x:(R^2), x = ([1,2],[3.4,2], 0, [2,3]);
  cout << x << "\n";
#+END_SRC
This code gives
#+results: tuple-definition

Observe that the special value ~0~ is used there.

**** Tuple purpose

Tuples variables are just lists of data of the same type. In the ~pugs~
language, one cannot access to a specific value of the list nor alter
one of them. This is not something that should ever change. Tuples are
not arrays! The ~pugs~ language is not meant to allow low-level
instructions.

The use case of tuples is to provide lists of data to the ~C++~
underlying methods. A classical example is to provide a set of
boundary conditions to a method.

** Statements

The ~pugs~ language supports classical statements to control the data
flow. For simplicity, these statements syntax follow their ~C++~
counterpart. The only statement that is not implemented in ~pugs~ is the
~switch...case~. This may change but in the one hand, up to now it was
never necessary (up to now, we did not encountered the need to chain
~if...else~ statements), and on the other hand, keeping the language as
simple as possible remains the policy in ~pugs~ development.

*** ~if...else~ statement

The simplest statement is the conditional statement ~if...else~.  If a
condition is satisfied, the ~truestatement~ (a single instruction or
block of instructions) is executed. In the other case, if the /optional/
~else~ keyword is used, the ~falsestatement~ (a single instruction or
block of instructions) is executed.
#+BEGIN_SRC pugs :exports code
  if (condition) truestatement
  if (condition) truestatement else falsestatement
#+END_SRC

The condition itself *must be* a boolean value (of type ~B~).

#+NAME: if-instruction
#+BEGIN_SRC pugs :exports both :results output
  if (2>1) cout << "yes: 2 > 1\n";
  if (not (2>1)) cout << "hmm... 2 <= 1 ?\n";
#+END_SRC
generates the following expected output. The second output is
hopefully not executed.
#+results: if-instruction
It is probably better to use the ~else~ keyword to obtain the same
result.

An example of the use of the ~else~ keyword follows.
#+NAME: if-else-instruction
#+BEGIN_SRC pugs :exports both :results output
  if (2<=1)
    cout << "hmm... 2 <= 1 ?\n";
  else
    cout << "ok: 2 > 1\n";
#+END_SRC
#+results: if-else-instruction

However, it is recommended even in the case of single instructions in
the conditional statements to use instruction blocks:
#+NAME: if-else-block-of-one
#+BEGIN_SRC pugs :exports both :results output
  import math;
  let sin100_positive:B;

  if (sin(100) > 0){
     sin100_positive = true;
  } else {
     sin100_positive = false;
  }

  cout << "sin(100)>0: " << sin100_positive << "\n";
#+END_SRC
#+results: if-else-block-of-one

We give a final illustration
#+NAME: if-block
#+BEGIN_SRC pugs :exports both :results output
  let m:R, m = 2.5;
  let M:R, M = 2;
  let has_swapped:B, has_swapped = false;

  if (m>M){
     let tmp:R, tmp = m;
     m = M;
     M = tmp;
     has_swapped = true;
  }

  cout << "min = " << m << " max = " << M << " has_swapped = " << has_swapped << "\n";
#+END_SRC
#+results: if-block

*** ~for~ loops

~pugs~ allows to write ~for~ loops. It follows the ~C++~ syntax
#+BEGIN_SRC pugs :exports code
  for (declarationinstruction ; condition ; postinstruction) statement
#+END_SRC
Following ~C++~, ~declarationinstruction~, ~condition~ and ~postinstruction~
are optional. The ~condition~ argument, if it is present, *must* be a
boolean value (type ~B~). If it is absent, the default value ~true~ is
used.

- The ~declarationinstruction~ is execute only *once* /before/ the beginning
  of the loop. The lifetime of the variable declared here is defined
  by the ~for~ instruction itself.

- The ~condition~ is evaluate /before/ each loop.

- The ~postinstruction~ is executed /after/ each loop.

- The ~statement~ is either a single instruction or a block of
  instructions. The ~statement~ is executed if the ~condition~ has the
  value ~true~.

For instance, one can write
#+NAME: for-block
#+BEGIN_SRC pugs :exports both :results output
  let n:N, n = 10;
  let sum:N, sum = 0;
  for (let i:N, i=1; i<=n; ++i) {
      sum += i;
  }
  cout << "sum = " << sum << "\n";
#+END_SRC
which gives as expected
#+results: for-block

The lifetime of the declared variable (in the ~declarationinstruction~
statement) is illustrated by the following example
#+NAME: for-scope-error
#+BEGIN_SRC pugs-error :exports both :results output
  for (let i:N, i=0; i<2; ++i) {
      cout << "i = " << i << "\n";
  }
  cout << "i = " << i << "\n";
#+END_SRC
Running this example produces the following error
#+results: for-scope-error

To fix the previous code, one can write
#+NAME: for-no-decl
#+BEGIN_SRC pugs :exports both :results output
  let i:N, i=0;
  for (; i<2; ++i) {
      cout << "i = " << i << "\n";
  }
  cout << "i = " << i << "\n";
#+END_SRC
One gets
#+results: for-no-decl

*** ~do...while~ loops

The second kind of loops that is supported is the ~do...while~
construction which executes /at least/ one time the enclosed ~statement~.
#+BEGIN_SRC pugs :exports code
  do statement while (condition);
#+END_SRC

The ~statement~ is either a single instruction or a block of
instructions. The ~condition~ is an expression of boolean value (type
~B~).
#+NAME: do-while-block
#+BEGIN_SRC pugs :exports both :results output
  let sum:N, sum = 0;
  let i:N, i = 0;
  do {
    sum += i;
    ++i;
  } while (i<=10);
  cout << "sum = " << sum << "\n";
#+END_SRC
It gives
#+results: do-while-block

*** ~while~ loops

The last kind of loops that is allowed in ~pugs~ language is the ~while~
loop.
#+BEGIN_SRC pugs :exports code
  while (condition) statement
#+END_SRC
The ~statement~ is either a single instruction or a block of
instructions. The ~condition~ is an expression of boolean value (type
~B~).

This time is the ~condition~ is not satisfied (~false~ when reaching the
~while~ instruction), the ~statment~ is never executed.

An example of the ~while~ loop is the following.
#+NAME: while-block
#+BEGIN_SRC pugs :exports both :results output
  let sum:N, sum = 0;
  let i:N, i = 1;
  while (i<=10) {
    sum += i;
    ++i;
  }
  cout << "sum = " << sum << "\n";
#+END_SRC
The result is
#+results: while-block


*** ~break~ and ~continue~ keywords.

These *keywords* are used to control loops behavior from enclosed loop
statements. They follow their ~C++~ counterparts.

- The ~continue~ keyword is used to skip the instructions that follow in
  the loop statement and continue the loop.

- The ~break~ keyword leaves the current loop.

An example of use of the ~continue~ keyword is
#+NAME: nested-continue
#+BEGIN_SRC pugs :exports both :results output
  for (let i:N, i = 0; i < 5; ++i) {
    cout << i << ": ";
    for (let j:N, j=0; j < 5; ++j) {
      if (j<i) continue;
      cout << j << " ";
    }
    cout << "\n";
  }
#+END_SRC
The result is
#+results: nested-continue

Replacing the ~continue~ keyword by a ~break~ (and changing the test)
#+NAME: nested-break
#+BEGIN_SRC pugs :exports both :results output
  for (let i:N, i = 0; i < 5; ++i) {
    cout << i << ": ";
    for (let j:N, j=0; j < 5; ++j) {
      if (j>i) break;
      cout << j << " ";
    }
    cout << "\n";
  }
#+END_SRC
changes the output to
#+results: nested-break

Obviously the behavior is the same using ~do...while~ or ~while~ loops.

** Functions<<functions>>

The ~pugs~ language allows the manipulation and definition of
functions. These are *mathematical functions* and are not functions as
functions in ~C++~: functions in the ~pugs~ language are *not subroutines*.

To be more precise, a function $f$ follows the following structure
\begin{align*}
    \mbox{let }f:& X_1\times \cdots \times X_n \to Y_1\times \cdots \times Y_m,\\
                 & (x_1,\ldots,x_n)\mapsto (y_1,\ldots,y_m)=f(x_1,\ldots,x_n),
\end{align*}
where $n,m\in\mathbb{N}_{\ge1}$, and where ${(X_i)}_{1\le i\le n}$ and ${(Y_i)}_{1\le
i\le m}$ are /simple/ types. Actually $X_1\times \cdots \times X_n$ and
$Y_1\times \cdots \times Y_m$ are /compound/ types.

Thus assuming that the function $f$ is defined in the language as ~f~,
one can use the following syntax to evaluate it. This is a
pseudo-code, real examples will follow. Assuming that ~(x1...,xn)~ has
been properly defined in ~X1*...*Xn~ one can write
#+BEGIN_SRC pugs :exports code
  let (y1,...,ym):Y1*...*Ym, (y1,...,ym) = f(x1,...,xn);
#+END_SRC
or if ~(y1,...,ym)~ has already been declared in  ~Y1*...*Ym~
#+BEGIN_SRC pugs :exports code
  (y1,...,ym) = f(x1,...,xn);
#+END_SRC

The ~pugs~ language handles two kinds of functions. User-defined
functions (see [[user-defined-functions]]) and builtin functions (see
[[builtin-functions]]). They behave essentially the same and are both
constant. In this context, it means that one cannot just declare or
modify a function.

*** Pure functions

In the ~pugs~ language, functions are *pure functions* in the sense that
arguments given to a function are *never* modified by the function. They
act as operators.

#+BEGIN_note
Actually these functions are not strictly /pure functions/ in the
computer science sense. The reason is that they can eventually have
side effects. As an example, it is possible to modify the random seed
used by the code. In that case, the modified value is not a variable
of the language itself but the internal random seed.
#+END_note

*** Implicit type conversion for parameters and returned values

As already said, in general, there is no implicit type conversion in
the ~pugs~ language. The only exceptions are when initializing tuples
(see [[tuples]]), passing arguments to functions and dealing with returned
values. See the conversion table in section [[implicit-conversion]].

*** User-defined functions<<user-defined-functions>>

To define user functions, the syntax mimics mathematics. The
definition of the function itself is a *single* expression. This means
that one cannot use ~if...else~, loops and there is not such keyword as
~return~ in the ~C++~.

The syntax of the definition of functions is
#+BEGIN_SRC pugs :exports code
  let f1 : X -> Y, x -> e;
  let f2 : X -> Y1*...*Ym, x -> (e1,...,em);
  let f3 : X1*...*Xn -> Y, (x1,...,xn) -> e;
  let f4 : X1*...*Xn -> Y1*...*Ym, (x1,...,xn) -> (e1,...,em);
#+END_SRC
~X~, ~Y~, ~X1~, ..., ~Xn~ and ~Y1~, ..., ~Ym~ are /simple/ data types (not
compound). The type of ~x~ is ~X~ and for lists, each argument ~xi~ belongs
to ~Xi~. The function themselves are defined by the expressions ~e~ (or
~e1~, ..., ~em~) which types are set by ~Y~ (or ~Y1~, ..., ~Ym~).


Let us give a few examples.
#+NAME: R-to-R-function
#+BEGIN_SRC pugs :exports both :results output
  let f: R -> R, x -> -x;

  cout << "f(2.3) = " << f(2.3)
       << "\nf(1) = " << f(1)
       << "\nf(true) = " << f(true)
       << "\n";
#+END_SRC
One observes the implicit conversion of the values passed by argument.
#+results: R-to-R-function

The following example illustrates the implicit conversion to the
returned type
#+NAME: R-to-R1-function
#+BEGIN_SRC pugs :exports both :results output
  let f: R -> R^1, x -> 2*x;

  cout << "f(3.2) = " << f(3.2) << "\n";
#+END_SRC
#+results: R-to-R1-function

Using compound types as input and output, one can write
#+NAME: R22-R-string-to-R-string-function
#+BEGIN_SRC pugs :exports both :results output
  let f : R^2x2*R*string -> R*string*R^3x3,
          (A,x,s) -> (x*A[0,0]*A[1,1]-A[1,0],
                      A+","+x+","+s,
                      [[A[0,0], A[0,1], 0],
                       [A[1,0], A[1,1], 0],
                       [     0,      0, x]]);
  let x : R, x=0;
  let s : string, s= "";
  let A : R^3x3, A = 0;
  (x,s,A) = f([[1,2],[3,4]], 2.3, "foo");
  cout << "x = " << x
       << "\ns = " << s
       << "\nA = " << A << "\n";
  let (y,t,A2):R*string*R^3x3, (y,t,A2) = f([[3,1],[4,2]], -5.2, "bar");
  cout << "y = " << y
       << "\nt = " << t
       << "\nA2 = " << A2 << "\n";
#+END_SRC
This meaningless example produces the following result.
#+results: R22-R-string-to-R-string-function

**** Lifetime of functions arguments

The arguments used to define a function are *local* variables that
exists only during the evaluation of the function.

Let us consider the following example
#+NAME: lifetime-of-function-args
#+BEGIN_SRC pugs :exports both :results output
  let a:R, a = 1.4;

  let plus: R -> R, a -> (a>0)*a;

  cout << "plus(1)    = " << plus(1) << "\n";
  cout << "plus(-3.2) = " << plus(-3.2) << "\n";
  cout << "a          = " << a << "\n";
#+END_SRC
This gives the expected result: the value of the variable ~a~ is
unchanged.
#+results: lifetime-of-function-args

**** Non-arguments variables in function expressions

Here we discuss rapidly of using variables (which are not arguments)
in function expressions.

#+NAME: non-arg-variables-in-functions
#+BEGIN_SRC pugs :exports both :results output
  let a:R, a = 1.4;

  let plus: R -> R, x -> (x>0)*a;
  cout << "a = " << a << "\n";
  cout << "plus(2.2)  = " << plus(2.2) << "\n";
  cout << "plus(-3.2) = " << plus(-3.2) << "\n";

  a = 2;
  cout << "a = " << a << "\n";
  cout << "plus(2.2)  = " << plus(2.2) << "\n";
  cout << "plus(-3.2) = " << plus(-3.2) << "\n";
#+END_SRC
Running the example, one gets
#+results: non-arg-variables-in-functions
While the function itself is a constant object, one sees that since
the value of ~a~ is changed, the value function is implicitly
modified. /This is a dangerous feature and should be avoid!/

Since functions themselves are variables one can use functions in
function expressions.
#+NAME: functions-in-functions
#+BEGIN_SRC pugs :exports both :results output
  let plus: R -> R, x -> (x>0)*x;
  let minus: R -> R, x -> -(x<0)*x;

  let pm : R -> R*R, x -> (plus(x), minus(x));
  let toR2: R*R -> R^2, (x,y) -> [x,y];

  cout << "pm(2) = " << toR2(pm(2)) << " pm(-3) = " << toR2(pm(-3)) << "\n";
#+END_SRC
One observes the utility function ~toR2~ that is used to perform the
output since ~cout~ does not handle compound types output. One gets
#+results: functions-in-functions

**** Lifetime of user-defined functions

Since functions are somehow variables, the lifetime of functions
follows the similar rules.

Let us give an example
#+NAME: functions-lifetime
#+BEGIN_SRC pugs :exports both :results output
  let f: R -> R, x -> 2.3*x+2;
  {
    let f: R -> R, x -> 1.25*x+3;
    cout << "block 1:  f(2) = " << f(2) << "\n";
    {
      let f: R -> R, x -> 3.25*x-0.3;
      cout << "block 2:  f(2) = " << f(2) << "\n";
    }
  }
  for (let i:N, i = 1; i<=2; ++i) {
    let f: R -> R, x -> i*x+2;
    cout << "for(i=" << i << "): f(2) = " << f(2) << "\n";
  }
  cout << "global:   f(2) = " << f(2) << "\n";
#+END_SRC
Running this example produces
#+results: functions-lifetime

**** Recursion

Since functions cannot be simply declared but must be defined, double
recursion is not possible. Thus, there is no equivalent construction
in ~pugs~ of the following ~C++~ code
#+BEGIN_SRC C++ :exports source
  int f(int i);
  int g(int i)
  {
    if (i < 0)
      return 0;
    else
      return f(i-1);
  }

  int f(int i) { return g(i); }
#+END_SRC
Moreover simple recursion is forbidden for similar reasons, since
functions are made of a single expression (one cannot use statements),
there would be no way to end the recursion. Thus, the code
#+NAME: no-recursion
#+BEGIN_SRC pugs-error :exports both :results output
  let f: N->N, n -> f(2*n);
#+END_SRC
produces the following compilation time error
#+results: no-recursion

*** Builtin functions<<builtin-functions>>

In ~pugs~ language, builtin functions are ~C++~ pieces of code that can be
called in scripts. There usage is very similar to user-defined
functions. They differ from user-defined functions in three points.
- Builtin functions can have no parameter or no returned value.
- Builtin functions are polymorphic. More precisely, this means that
  the signature of a builtin function is also defined by its expected
  arguments types.
- Builtin functions can take user-defined functions as parameters.
  - user-defined functions cannot take functions as parameters
  - builtin functions cannot take builtin functions as parameters
    (actually, this is not a limitation since it is trivial to embed a
    builtin function into a user-defined one).

Here is a simple example of builtin function embedding in a user
function.
#+NAME: builtin-function-embedding
#+BEGIN_SRC pugs :exports both :results output
  import math;

  let cosinus: R -> R, x -> cos(x);
  cout << "cosinus(2) = " << cosinus(2) << "\n";
#+END_SRC
Running this example produces
#+results: builtin-function-embedding

Builtin functions are defined when modules are imported, see [[modules]].

** TODO modules<<modules>>


[fn:pugs-def] ~pugs~: Parallel Unstructured Grid Solvers
[fn:MPI-def] ~MPI~: Message Passing Interface
[fn:DSL-def] ~DSL~: Domain Specific Language