Add an equation to the DAE system to be solved in the phase space spanned by `PHASE_SPACE`

.

Assign an expression to a variable, a vector or a matrix.

Catastrophically abort the execution and quit wasora.

Whenever the instruction `ABORT`

is executed, wasora quits without closing files or unlocking shared memory objects. The objective of this instruction is, as illustrated in the examples, either to debug complex input files and check the values of certain variables or to conditionally abort the execution using `IF`

clauses.

Define a scalar alias of an already-defined indentifier.

The existing object can be a variable, a vector element or a matrix element. In the first case, the name of the variable should be given as the existing object. In the second case, to alias the second element of vector `v`

to the new name `new`

, `v(2)`

should be given as the existing object. In the third case, to alias second element (2,3) of matrix `M`

to the new name `new`

, `M(2,3)`

should be given as the existing object.

Call a previously dynamically-loaded user-provided routine.

Explicitly close an already-`OPEN`

ed file.

Mark a scalar variable, vector or matrix as a constant.

Give a default value for an optional commandline argument.

If a `$n`

construction is found in the input file but the commandline argument was not given, the default behavior is to fail complaining that an extra argument has to be given in the commandline. With this keyword, a default value can be assigned if no argument is given, thus avoiding the failure and making the argument optional.

Explicitly mark variables, vectors or matrices as “differential” to compute intial conditions of DAE systems.

Ask wasora not to evaluate assignments at parse time.

Define a file, either as input or as output, for further usage.

Fit a function of one or more arguments to a set of data.

The function with the data has to be point-wise defined. The function to be fitted hast to be parametrized with at least one of the variables provided after the `VIA`

keyword. Only the names of the functions have to be given. Both functions have to have the same number of arguments. The initial guess of the solution is given by the initial value of the variables listed in the `VIA`

keyword. Analytical expressions for the gradient of the function to be fitted with respect to the parameters to be fitted can be optionally given with the `GRADIENT`

keyword. If none is provided, the gradient will be computed numerically using finite differences. A range over which the residuals are to be minimized can be given with `RANGE_MIN`

and `RANGE_MAX`

. For multidimensional fits, the range is an hypercube. If no range is given, all the definition points of the function witht the data are used for the fit. Convergence can be controlled by given the relative and absolute tolreances with `DELTAEPSREL`

(default 1e-4) and `DELTAEPSABS`

(default 1e-6), and with the maximum number of iterations `MAX_ITER`

(default 100). If the optional keyword `VERBOSE`

is given, some data of the intermediate steps is written in the standard output.

Define a function of one or more variables.

The number of variables n is given by the number of arguments given between parenthesis after the function name. The arguments are defined as new variables if they had not been already defined as variables. If the function is given as an algebraic expression, the short-hand operator `:=`

can be used. That is to say, `FUNCTION f(x) = x^2`

is equivalent to `f(x) := x^2`

. If a `FILE_PATH`

is given, an ASCII file containing at least n+1 columns is expected. By default, the first n columns are the values of the arguments and the last column is the value of the function at those points. The order of the columns can be changed with the keyword `COLUMNS`

, which expects n+1 expressions corresponding to the column numbers. A function of type `ROUTINE`

calls an already-defined user-provided routine using the `CALL`

keyword and passes the values of the variables in each required evaluation as a `double *`

argument. If `MESH`

is given, the definition points are the nodes or the cells of the mesh. The function arguments should be (x), (x,y) or (x,y,z) matching the dimension the mesh. If the keyword `DATA`

is used, a new empty vector of the appropriate size is defined. The elements of this new vector can be assigned to the values of the function at the i-th node or cell. If the keyword `VECTOR`

is used, the values of the dependent variable are taken to be the values of the already-existing vector. Note that this vector should have the size of the number of nodes or cells the mesh has, depending on whether `NODES`

or `CELLS`

is given. If `VECTOR_DATA`

is given, a set of n+1 vectors of the same size is expected. The first n+1 correspond to the arguments and the last one is the function value. Interpolation schemes can be given for either one or multi-dimensional functions with `INTERPOLATION`

. Available schemes for n=1 are:

- linear
- polynomial, the grade is equal to the number of data minus one
- spline, cubic (needs at least 3 points)
- spline_periodic
- akima (needs at least 5 points)
- akima_periodic (needs at least 5 points)
- steffen, always-monotonic splines-like (available only with GSL >= 2.0)

Default interpolation scheme for one-dimensional functions is `(*gsl_interp_linear)`

.

Available schemes for n>1 are:

- nearest, f(\vec{x}) is equal to the value of the closest definition point
- shepard, inverse distance weighted average definition points (might lead to inefficient evaluation)
- shepard_kd, average of definition points within a kd-tree (more efficient evaluation provided
`SHEPARD_RADIUS`

is set to a proper value) - bilinear, only available if the definition points configure an structured hypercube-like grid. If n>3,
`SIZES`

should be given.

For n>1, if the euclidean distance between the arguments and the definition points is smaller than `INTERPOLATION_THRESHOLD`

, the definition point is returned and no interpolation is performed. Default value is square root of `9.5367431640625e-07`

. The initial radius of points to take into account in `shepard_kd`

is given by `SHEPARD_RADIUS`

. If no points are found, the radius is double until at least one definition point is found. The radius is doubled until at least one point is found. Default is `1.0`

. The exponent of the `shepard`

method is given by `SHEPARD_EXPONENT`

. Default is `2`

. When requesting `bilinear`

interpolation for n>3, the number of definition points for each argument variable has to be given with `SIZES`

, and wether the definition data is sorted with the first argument changing first (`X_INCREASES_FIRST`

evaluating to non-zero) or with the last argument changing first (zero). The function can be pointwise-defined inline in the input using `DATA`

. This should be the last keyword of the line, followed by N=k\cdot (n+1) expresions giving k definition points: n arguments and the value of the function. Multiline continuation using brackets `{`

and `}`

can be used for a clean data organization. See the examples.

Record the time history of a variable as a function of time.

Begin a conditional block.

```
IF expr
<block_of_instructions_if_expr_is_true>
[ ELSE ]
[block_of_instructions_if_expr_is_false]
ENDIF
```

Define whether implicit declaration of variables is allowed or not.

By default, wasora allows variables (but not vectors nor matrices) to be implicitly declared. To avoid introducing errors due to typos, explicit declaration of variables can be forced by giving `IMPLICIT NONE`

. Whether implicit declaration is allowed or explicit declaration is required depends on the last `IMPLICIT`

keyword given, which by default is `ALLOWED`

.

Include another wasora input file.

Includes the input file located in the string `file_path`

at the current location. The effect is the same as copying and pasting the contents of the included file at the location of the `INCLUDE`

keyword. The path can be relative or absolute. Note, however, that when including files inside `IF`

blocks that instructions are conditionally-executed but all definitions (such as function definitions) are processed at parse-time independently from the evaluation of the conditional. The optional `FROM`

and `TO`

keywords can be used to include only portions of a file.

Define how initial conditions of DAE problems are computed.

In DAE problems, initial conditions may be either:

- equal to the provided expressions (
`AS_PROVIDED`

) - the derivatives computed from the provided phase-space variables (
`FROM_VARIABLES`

) - the phase-space variables computed from the provided derivatives (
`FROM_DERIVATIVES`

)

In the first case, it is up to the user to fulfill the DAE system at t = 0. If the residuals are not small enough, a convergence error will occur. The `FROM_VARIABLES`

option means calling IDA’s `IDACalcIC`

routine with the parameter `IDA_YA_YDP_INIT`

. The `FROM_DERIVATIVES`

option means calling IDA’s `IDACalcIC`

routine with the parameter IDA_Y_INIT. Wasora should be able to automatically detect which variables in phase-space are differential and which are purely algebraic. However, the `DIFFERENTIAL`

keyword may be used to explicitly define them. See the (SUNDIALS documentation)[https://computation.llnl.gov/casc/sundials/documentation/ida_guide.pdf] for further information.

Load a wasora plug-in from a dynamic shared object.

A wasora plugin in the form of a dynamic shared object (i.e. `.so`

) can be loaded either with the `LOAD_PLUGIN`

keyword or from the command line with the `-p`

option. Either a file path or a plugin name can be given. The following locations are tried:

- the current directory
`./`

- the parent directory
`../`

- the user’s
`LD_LIBRARY_PATH`

- the cache file
`/etc/ld.so.cache`

- the directories
`/lib`

,`/usr/lib`

,`/usr/local/lib`

If a wasora plugin was compiled and installed following the `make install`

procedure, the plugin should be loaded by just passing the name to `LOAD_PLUGIN`

.

Load one or more routines from a dynamic shared object.

Call the `m4`

macro processor with definitions from wasora variables or expressions.

Define a matrix.

Find the combination of arguments that give a (relative) minimum of a function, i.e. run an optimization problem.

Systematically sweep a zone of the parameter space, i.e. perform a parametric run.

Define which variables, vectors and/or matrices belong to the phase space of the DAE system to be solved.

Print plain-text and/or formatted data to the standard output or into an output file.

Each argument `object`

that is not a keyword is expected to be part of the output, can be either a matrix, a vector, an scalar algebraic expression. If the given object cannot be solved into a valid matrix, vector or expression, it is treated as a string literal if `IMPLICIT`

is `ALLOWED`

, otherwise a parser error is raised. To explicitly interpret an object as a literal string even if it resolves to a valid numerical expression, it should be prefixed with the `TEXT`

keyword. Hashes `#`

appearing literal in text strings have to be quoted to prevent the parser to treat them as comments within the wasora input file and thus ignoring the rest of the line. Whenever an argument starts with a porcentage sign `%`

, it is treated as a C `printf`

-compatible format definition and all the objects that follow it are printed using the given format until a new format definition is found. The objects are treated as double-precision floating point numbers, so only floating point formats should be given. The default format is `"%g"`

. Matrices, vectors, scalar expressions, format modifiers and string literals can be given in any desired order, and are processed from left to right. Vectors are printed element-by-element in a single row. See `PRINT_VECTOR`

to print vectors column-wise. Matrices are printed element-by-element in a single line using row-major ordering if mixed with other objects but in the natural row and column fashion if it is the only given object. If the `FILE`

keyword is not provided, default is to write to stdout. If the `NONEWLINE`

keyword is not provided, default is to write a newline ‘’ character after all the objects are processed. The `SEP`

keywords expects a string used to separate printed objects, the default is a tab ‘DEFAULT_PRINT_SEPARATOR’ character. Use the `NOSEP`

keyword to define an empty string as object separator. If the `HEADER`

keyword is given, a single line containing the literal text given for each object is printed at the very first time the `PRINT`

instruction is processed, starting with a hash `#`

character. If the `SKIP_STEP`

(`SKIP_STATIC_STEP`

)keyword is given, the instruction is processed only every the number of transient (static) steps that results in evaluating the expression, which may not be constant. By default the `PRINT`

instruction is evaluated every step. The `SKIP_HEADER_STEP`

keyword works similarly for the optional `HEADER`

but by default it is only printed once. The `SKIP_TIME`

keyword use time advancements to choose how to skip printing and may be useful for non-constant time-step problems.

Print one or more functions as a table of values of dependent and independent variables.

Print the elements of one or more vectors.

Read data (variables, vectors o matrices) from files or shared-memory segments.

Perform either a wait or a post operation on a named shared semaphore.

Execute a shell command.

Solve a non-linear system of n equations with n unknowns.

Force transient problems to pass through specific instants of time.

The time step `dt`

will be reduced whenever the distance between the current time `t`

and the next expression in the list is greater than `dt`

so as to force `t`

to coincide with the expressions given. The list of expresssions should evaluate to a sorted list of values.

Define one or more scalar variables.

Define a vector.

Sort the elements of a vector into ascending numerical order.

Write data (variables, vectors o matrices) to files or shared-memory segments. See the `READ`

keyword for usage details.

Flag that indicates whether the overall calculation is over.

Flag that indicates whether the parametric, optimization of fit calculation is over or not. It is set to true (i.e. \neq 0) by wasora whenever the outer calculation is considered to be finished, which can be that the parametric calculation swept the desired parameter space or that the optimization algorithm reached the desired convergence criteria. If the user sets it to true, the current step is marked as the last outer step and the transient calculation ends after finishing the step.

Flag that indicates whether the static calculation is over or not. It is set to true (i.e. \neq 0) by wasora if `step_static`

\ge `static_steps`

. If the user sets it to true, the current step is marked as the last static step and the static calculation ends after finishing the step.

Flag that indicates whether the transient calculation is over or not. It is set to true (i.e. \neq 0) by wasora if `t`

\ge `end_time`

. If the user sets it to true, the current step is marked as the last transient step and the transient calculation ends after finishing the step.

Actual value of the time step for transient calculations. When solving DAE systems, this variable is set by wasora. It can be written by the user for example by importing it from another transient code by means of shared-memory objects. Care should be taken when solving DAE systems and overwriting `t`

. Default value is 1/16, which is a power of two and roundoff errors are thus reduced.

Final time of the transient calculation, to be set by the user. The default value is zero, meaning no transient calculation.

Dummy index, used mainly in vector and matrix row subindex expressions.

A very big positive number, which can be used as `end_time = infinite`

or to define improper integrals with infinite limits. Default is 2^{50} \approx 1 \times 10^{15}.

Flag that indicates if the current step is the initial step of an optimization of fit run.

Flag that indicates if wasora is solving the iterative static calculation. Flag that indicates if wasora is in the first step of the iterative static calculation. Flag that indicates if wasora is in the last step of the iterative static calculation.

Flag that indicates if wasora is solving transient calculation.

Flag that indicates if wasora is in the first step of the transient calculation.

Flag that indicates if wasora is in the last step of the transient calculation.

Dummy index, used mainly in matrix column subindex expressions.

Maximum bound for the time step that wasora should take when solving DAE systems.

Minimum bound for the time step that wasora should take when solving DAE systems.

The number of online available cores, as returned by `sysconf(_SC_NPROCESSORS_ONLN)`

. This value can be used in the `MAX_DAUGHTERS`

expression of the `PARAMETRIC`

keyword (i.e `ncores/2`

).

This should be set to a mask that indicates how to proceed if an error ir raised in any routine of the GNU Scientific Library.

This should be set to a mask that indicates how to proceed if an error ir raised in any routine of the SUNDIALS IDA Library.

This should be set to a mask that indicates how to proceed if Not-A-Number signal (such as a division by zero) is generated when evaluating any expression within wasora.

A double-precision floating point representaion of the number \pi, equal to `math.h`

’s `M_PI`

constant.

The UNIX process id of wasora (or the plugin).

If this variable is not zero, then the transient problem is run trying to syncrhonize the problem time with realtime, up to a scale given. For example, if the scale is set to one, then wasora will advance the problem time at the same pace that the real wall time advances. If set to two, wasora’s time wil advance twice as fast as real time, and so on. If the calculation time is slower than real time modified by the scale, this variable has no effect on the overall behavior and execution will proceed as quick as possible with no delays.

Maximum allowed relative error for the solution of DAE systems. Default value is is 1 \times 10^{-6}. If a fine per-variable error control is needed, special vector `abs_error`

should be used.

Number of steps that ought to be taken during the static calculation, to be set by the user. The default value is one, meaning only one static step.

Indicates the current step number of the iterative outer calculation (parametric, optimization or fit). Indicates the current step number of the iterative inner calculation (optimization or fit).

Indicates the current step number of the iterative static calculation.

Indicates the current step number of the transient static calculation.

Actual value of the time for transient calculations. This variable is set by wasora, but can be written by the user for example by importing it from another transient code by means of shared-memory objects. Care should be taken when solving DAE systems and overwriting `t`

.

A very small positive number, which is taken to avoid roundoff errors when comparing floating point numbers such as replacing a \leq a_\text{max} with a < a_\text{max} + `zero`

. Default is (1/2)^{-50} \approx 9\times 10^{-16} .

Number of cells of the unstructured grid. This number is the actual quantity of volumetric elements in which the domain was discretized.

Number of total elements of the unstructured grid. This number include those surface elements that belong to boundary physical entities.

Number of nodes of the unstructured grid.

Returns the absolute value of the argument x.

Computes arc in radians whose cosine is equal to the argument x. A NaN error is raised if |x|>1.

Computes arc in radians whose sine is equal to the argument x. A NaN error is raised if |x|>1.

Computes, in radians, the arc tangent of the argument x.

Computes, in radians, the arc tangent of quotient y/x, using the signs of the two arguments to determine the quadrant of the result, which is in the range [-\pi,\pi].

Returns the smallest integral value not less than the argument x.

Returns the value of a certain clock in seconds measured from a certain (but specific) milestone. The kind of clock and the initial milestone depends on the optional flag f. It defaults to zero, meaning wall time since the UNIX Epoch. The list and the meanings of the other available values for f can be checked in the `clock_gettime (2)`

system call manual page.

Computes the cosine of the argument x, where x is in radians. A cosine wave can be generated by passing as the argument x a linear function of time such as \omega t+\phi, where \omega controls the frequency of the wave and \phi controls its phase.

Computes the hyperbolic cosine of the argument x, where x is in radians.

Computes the time derivative of the signal x using the difference between the value of the signal in the previous time step and the actual value divided by the time step. For t=0, the return value is zero. Unlike the functional `derivative`

, this function works with expressions and not with functions. Therefore the argument x may be for example an expression involving a variable that may be read from a shared-memory object, whose time derivative cannot be computed with `derivative`

.

Filters the first argument x with a deadband centered at zero with an amplitude given by the second argument a.

Checks if the two first expressions a and b are equal, up to the tolerance given by the third optional argument \epsilon. If either |a|>1 or |b|>1, the arguments are compared using GSL’s `gsl_fcmp`

, otherwise the absolute value of their difference is compared against \epsilon. This function returns zero if the arguments are not equal and one otherwise. Default value for \epsilon = 10^{-16}.

Computes the exponential function the argument x, i.e. the base of the natural logarithms raised to the x-th power.

Computes the first exponential integral function of the argument x. If x equals zero, a NaN error is issued.

Computes the second exponential integral function of the argument x.

Computes the third exponential integral function of the argument x.

Computes the n-th exponential integral function of the argument x. If n equals zero or one and x zero, a NaN error is issued.

Returns the largest integral value not greater than the argument x.

Computes the zero-centered Heaviside step function of the argument x. If the optional second argument \epsilon is provided, the discontinuous step at x=0 is replaced by a ramp starting at x=0 and finishing at x=\epsilon.

Performs a conditional testing of the first argument a, and returns either the second optional argument b if a is different from zero or the third optional argument c if a evaluates to zero. The comparison of the condition a with zero is performed within the precision given by the optional fourth argument \epsilon. If the second argument c is not given and a is not zero, the function returns one. If the third argument c is not given and a is zero, the function returns zero. The default precision is \epsilon = 10^{-16}. Even though `if`

is a logical operation, all the arguments and the returned value are double-precision floating point numbers.

Computes the time integral of the signal x using the trapezoidal rule using the value of the signal in the previous time step and the current value. At t = 0 the integral is initialized to zero. Unlike the functional `integral`

, this function works with expressions and not with functions. Therefore the argument x may be for example an expression involving a variable that may be read from a shared-memory object, whose time integral cannot be computed with `integral`

.

Idem as `integral_dt`

but uses the backward Euler rule to update the integral value. This function is provided in case this particular way of approximating time integrals is needed.

Returns one if the argument x rounded to the nearest integer is even.

Returns true if the argument~x is in the interval~[a,b), i.e. including~a but excluding~b.

Returns one if the argument x rounded to the nearest integer is odd.

Computes the regular cylindrical Bessel function of zeroth order evaluated at the argument x.

Filters the first argument x(t) with a first-order lag of characteristic time \tau, i.e. this function applies the transfer function !bt [ G(s) = ] !et to the time-dependent signal x(t), by assuming that it is constant during the time interval [t-\Delta t,t] and using the analytical solution of the differential equation for that case at t = \Delta t with the initial condition y(0) = y(t-\Delta t).

Filters the first argument x(t) with a first-order lag of characteristic time \tau, i.e. this function applies the transfer function !bt [ G(s) = ] !et to the time-dependent signal x(t) by using the bilinear transformation formula.

Filters the first argument x(t) with a first-order lag of characteristic time \tau, i.e. this function applies the transfer function !bt [ G(s) = ] !et to the time-dependent signal x(t) by using the Euler forward rule.

Returns the value the signal x had in the previous time step. This function is equivalent to the Z-transform operator “delay” denoted by z^{-1}\left[x\right]. For t=0 the function returns the actual value of x. The optional flag p should be set to one if the reference to `last`

is done in an assignment over a variable that already appears insi expression x. See example number 2.

Limits the first argument x to the interval [a,b]. The second argument a should be less than the third argument b.

Limits the value of the first argument x(t) so to that its time derivative is bounded to the interval [a,b]. The second argument a should be less than the third argument b.

Computes the natural logarithm of the argument x. If x is zero or negative, a NaN error is issued.

Returns the integer index i of the maximum of the arguments x_i provided. Currently only maximum of ten arguments can be provided.

Returns the integer index i of the minimum of the arguments x_i provided. Currently only maximum of ten arguments can be provided.

Returns the maximum of the arguments x_i provided. Currently only maximum of ten arguments can be provided.

Returns the minimum of the arguments x_i provided. Currently only maximum of ten arguments can be provided.

Returns the remainder of the division between the first argument a and the second b. Both arguments may be non-integral.

Returns one if the first argument x is zero and zero otherwise. The second optional argument \epsilon gives the precision of the “zero” evaluation. If not given, default is \epsilon = 10^{-16}.

Returns a random real number uniformly distributed between the first real argument x_1 and the second one x_2. If the third integer argument s is given, it is used as the seed and thus repetitive sequences can be obtained. If no seed is provided, the current time (in seconds) plus the internal address of the expression is used. Therefore, two successive calls to the function without seed (hopefully) do not give the same result. This function uses a second-order multiple recursive generator described by Knuth in Seminumerical Algorithms, 3rd Ed., Section 3.6.

Returns a random real number with a Gaussian distribution with a mean equal to the first argument x_1 and a standard deviation equatl to the second one x_2. If the third integer argument s is given, it is used as the seed and thus repetitive sequences can be obtained. If no seed is provided, the current time (in seconds) plus the internal address of the expression is used. Therefore, two successive calls to the function without seed (hopefully) do not give the same result. This function uses a second-order multiple recursive generator described by Knuth in Seminumerical Algorithms, 3rd Ed., Section 3.6.

Rounds the argument x to the nearest integer. Halfway cases are rounded away from zero.

Computes a sawtooth wave betwen zero and one with a period equal to one. As with the sine wave, a sawtooh wave can be generated by passing as the argument x a linear function of time such as \omega t+\phi, where \omega controls the frequency of the wave and \phi controls its phase.

Returns minus one, zero or plus one depending on the sign of the first argument x. The second optional argument \epsilon gives the precision of the “zero” evaluation. If not given, default is \epsilon = 10^{-16}.

Computes the sine of the argument x, where x is in radians. A sine wave can be generated by passing as the argument x a linear function of time such as \omega t+\phi, where \omega controls the frequency of the wave and \phi controls its phase.

Computes the hyperbolic sine of the argument x, where x is in radians.

Computes the positive square root of the argument x. If x is negative, a NaN error is issued.

Computes a square function betwen zero and one with a period equal to one. The output is one for 0 < x < 1/2 and goes to zero for 1/2 < x < 1. As with the sine wave, a square wave can be generated by passing as the argument x a linear function of time such as \omega t+\phi, where \omega controls the frequency of the wave and \phi controls its phase.

Computes the tangent of the argument x, where x is in radians.

Computes the hyperbolic tangent of the argument x, where x is in radians.

Returns one if the first argument x is greater than the threshold given by the second argument a, and zero otherwise. If the optional third argument b is provided, an hysteresis of width b is needed in order to reset the function value. Default is no hysteresis, i.e. b=0.

Returns one if the first argument x is less than the threshold given by the second argument a, and zero otherwise. If the optional third argument b is provided, an hysteresis of width b is needed in order to reset the function value. Default is no hysteresis, i.e. b=0.

Computes a triangular wave betwen zero and one with a period equal to one. As with the sine wave, a triangular wave can be generated by passing as the argument x a linear function of time such as \omega t+\phi, where $$ controls the frequency of the wave and \phi controls its phase.

Computes the derivative of the expression f(x) given in the first argument with respect to the variable x given in the second argument at the point x=a given in the third argument using an adaptive scheme. The fourth optional argument h is the initial width of the range the adaptive derivation method starts with. The fifth optional argument p is a flag that indicates whether a backward (p < 0), centered (p = 0) or forward (p > 0) stencil is to be used. This functional calls the GSL functions `gsl_deriv_central`

or `gsl_deriv_forward`

according to the indicated flag p. Defaults are h = (1/2)^{-10} \approx 9.8 \times 10^{-4} and p = 0.

Finds the value of the variable x given in the second argument which makes the expression f(x) given in the first argument to take local a minimum in the in the range [a,b] given by the third and fourth arguments. If there are many local minima, the one that is closest to (a+b)/2 is returned. The optional fifth argument \epsilon gives a relative tolerance for testing convergence, corresponding to GSL `epsrel`

(note that `epsabs`

is set also to \epsilon). The sixth optional argument is an integer which indicates the algorithm to use: 0 (default) is `quad_golden`

, 1 is `brent`

and 2 is `goldensection`

. See the GSL documentation for further information on the algorithms. The seventh optional argument p is a flag that indicates how to proceed if there is no local minimum in the range [a,b]. If p = 0 (default), a is returned if f(a) < f(b) and b otherwise. If p = 1 then the local minimum algorimth is tried nevertheless. Default is \epsilon = (1/2)^{-20} \approx 9.6\times 10^{-7}.

Computes the integral of the expression f(x) given in the first argument with respect to variable x given in the second argument over the interval [a,b] given in the third and fourth arguments respectively using a non-adaptive procedure which uses fixed Gauss-Kronrod-Patterson abscissae to sample the integrand at a maximum of 87 points. It is provided for fast integration of smooth functions. The algorithm applies the Gauss-Kronrod 10-point, 21-point, 43-point and 87-point integration rules in succession until an estimate of the integral is achieved within the relative tolerance given in the fifth optional argument \epsilon It correspondes to GSL’s `epsrel`

parameter (`epsabs`

is set to zero).

The rules are designed in such a way that each rule uses all the results of its predecessors, in order to minimize the total number of function evaluations. Defaults are \epsilon = (1/2)^{-10} \approx 10^{-3}. See GSL reference for further information.

Computes the integral of the expression f(x) given in the first argument with respect to variable x given in the second argument over the interval [a,b] given in the third and fourth arguments respectively using the n-point Gauss-Legendre rule, where n is given in the optional fourth argument. It is provided for fast integration of smooth functions with known polynomic order (it is exact for polynomials of order 2n-1). This functional calls GSL function `gsl_integration_glfixedp`

. Default is n = 12. See GSL reference for further information.

Computes the integral of the expression f(x) given in the first argument with respect to variable x given in the second argument over the interval [a,b] given in the third and fourth arguments respectively using an adaptive scheme, in which the domain is divided into a number of maximum number of subintervals and a fixed-point Gauss-Kronrod-Patterson scheme is applied to each quadrature subinterval. Based on an estimation of the error commited, one or more of these subintervals may be split to repeat the numerical integration alogorithm with a refined division. The fifth optional argument \epsilon is is a relative tolerance used to check for convergence. It correspondes to GSL’s `epsrel`

parameter (`epsabs`

is set to zero). The sixth optional argument 1\leq k \le 6 is an integer key that indicates the integration rule to apply in each interval. It corresponds to GSL’s parameter `key`

. The seventh optional argument gives the maximum number of subdivisions, which defaults to 1024. If the integration interval [a,b] if finite, this functional calls the GSL function `gsl_integration_qag`

. If a is less that minus the internal variable `infinite`

, b is greater that `infinite`

or both conditions hold, GSL functions `gsl_integration_qagil`

, `gsl_integration_qagiu`

or `gsl_integration_qagi`

are called. The condition of finiteness of a fixed range [a,b] can thus be changed by modifying the internal variable `infinite`

. Defaults are \epsilon = (1/2)^{-10} \approx 10^{-3} and k=3. The maximum numbers of subintervals is limited to 1024. Due to the adaptivity nature of the integration method, this function gives good results with arbitrary integrands, even for infinite and semi-infinite integration ranges. However, for certain integrands, the adaptive algorithm may be too expensive or even fail to converge. In these cases, non-adaptive quadrature functionals ought to be used instead. See GSL reference for further information.

Computes product of the N=b-a expressions f(i) given in the first argument by varying the variable~i given in the second argument between~a given in the third argument and~b given in the fourth argument,~i = a, a+1, \dots ,b-1,b.

Computes the value of the variable x given in the second argument which makes the expression f(x) given in the first argument to be equal to zero by using a root bracketing algorithm. The root should be in the range [a,b] given by the third and fourth arguments. The optional fifth argument \epsilon gives a relative tolerance for testing convergence, corresponding to GSL `epsrel`

(note that `epsabs`

is set also to \epsilon). The sixth optional argument is an integer which indicates the algorithm to use: 0 (default) is `brent`

, 1 is `falsepos`

and 2 is `bisection`

. See the GSL documentation for further information on the algorithms. The seventh optional argument p is a flag that indicates how to proceed if the sign of f(a) is equal to the sign of f(b). If p=0 (default) an error is raised, otherwise it is not. If more than one root is contained in the specified range, the first one to be found is returned. The initial guess is x_0 = (a+b)/2. If no roots are contained in the range and p \neq 0, the returned value can be any value. Default is \epsilon = (1/2)^{-10} \approx 10^{3}.

Computes sum of the N=b-a expressions f_i given in the first argument by varying the variable i given in the second argument between a given in the third argument and b given in the fourth argument, i=a,a+1,\dots,b-1,b.

Computes the dot product between vectors \vec{a} and \vec{b}, which should have the same size.

Returns the biggest element of vector \vec{b}, taking into account its sign (i.e. 1 > -2).

Returns the index of the biggest element of vector \vec{b}, taking into account its sign (i.e. 2 > -1).

Returns the smallest element of vector \vec{b}, taking into account its sign (i.e. -2 < 1).

Returns the index of the smallest element of vector \vec{b}, taking into account its sign (i.e. -2 < 1).

Computes euclidean norm of vector \vec{b}. Other norms can be computed explicitly using the `sum`

functional, as illustrated in the example.

Returns the size of vector \vec{b}.

Computes the sum of all the components of vector \vec{b}.