Performance Tips

Any code for which performance is important should be placed in a function. Code inside functions is usually much faster than top-level code, due to the way the Julia compiler works.

Functions are important not only in terms of performance: they are more reusable and testable, and they make it easier to understand what actions are being performed and what the input and output data are. The Julia style guide also recommends write functions, not just scripts.

Functions should accept arguments, and not work directly with global variables (see the next paragraph).

The value of an untyped global variable can change at any time, which can lead to a change in its type. Therefore, it is difficult for the compiler to optimize code that uses global variables. This also applies to variables whose values are types, that is, to type aliases at the global level. Variables should be local or passed to functions as arguments where possible.

Global names are often constants, and if you declare them in this way, you can significantly improve performance.:

If it is known that the type of the global variable will not change., it should be annotated.

The use of untyped global variables can be optimized by annotating their types at the time of use.:

From the point of view of style, it is better to pass arguments to functions. This improves the reusability of the code and allows for a better understanding of the input and output data.

In the REPL session, this code is:

equivalent to the following:

Therefore, all the performance-related features described above are applicable.

A macro is useful for measuring performance. @time. Here we will repeat the above example with a global variable, but this time we will remove the type annotation.:

On the first call ('@time sum_global()'), the function is compiled. (If there is a macro during the session @time has not been used yet, and the functions necessary for timing will also be compiled.) You should not take the results of this execution seriously. Note that for the second execution, in addition to the time, it is also reported that a significant amount of memory has been allocated. In this case, the sum of all the elements in a vector of 64-bit floating-point numbers is simply calculated, so there is no need to allocate memory (in the heap).

It should be clarified that the macro @time reports the allocation of memory exactly in the _ period_. Usually, such allocation is required either for mutable objects, or for creating or enlarging variable-size containers (for example, Array, Dict, strings, or type-unstable objects whose type is known only at runtime). Allocating (or deallocating) such memory blocks may require a resource-intensive function call for libc (for example, malloc in C) and tracking by the garbage collector. On the contrary, storing immutable values such as numbers (except bignum), tuples, and immutable structures (struct) requires much less resources (for example, they can be stored on the stack or in CPU register memory), so you don’t have to worry about performance when allocating memory for them.

Unexpected memory allocation is almost always a sign of some kind of code problem, usually related to type stability or the creation of many small temporary arrays. Therefore, even if you do not take into account the memory consumption, there is a high probability that the code generated for the function is far from optimal. Pay attention to these symptoms and follow the advice below.

In this particular case, memory is allocated due to the use of an unstable global variable x'. If you instead pass `x to the function as an argument, no more memory will be allocated (some of the allocated memory is lower due to the fact that the macro @time is executed in the global area). After the first call, this option will be much faster.:

The allocation of 1 in the output is due to the fact that the macro @time itself is executed in the global area. If the time is counted in the function itself, then no memory will be allocated at all.:

In some situations, a function may need to allocate memory during operation, and this complicates the picture. In such cases, you can use one of the following tools to diagnose problems, or write a version of the function in which memory allocation and algorithmic aspects are isolated from each other (see section Pre-allocation of memory for output data).

The Julia language and its ecosystem of packages offer tools that can help diagnose problems and improve code performance.:

When working with parameterized types, including arrays, it is best to avoid parameterization with abstract types whenever possible.

Consider the following code:

Since a is an array of an abstract type Real, it should be possible to store any value of Real' in it. Since `Real objects can be of arbitrary size and structure, the variable a must be represented by an array of pointers to individual Real objects in memory. However, if only one type of numbers could be stored in a, for example Float64, the storage efficiency would increase:

When assigning numbers to the variable a, they will now be converted to the Float64 type, and a will be stored as a continuous block of 64-bit floating-point values that can be efficiently operated on.

If you can’t avoid using containers with abstract value types, sometimes it’s better to parameterize them with the Any type so that the types aren’t checked at runtime. For example, IdDict{Any, Any} will run faster than IdDict{Type, Vector}.

See also the discussion in the section Parametric types.

In many languages with optional type declarations, adding declarations is the main way to improve code performance. In Julia, this is not the case. In Julia, the compiler usually knows the types of all arguments to functions, local variables, and expressions. However, there are a number of special situations in which ads can be useful.

If the function is implemented as several small definitions, the compiler can directly call the most appropriate code or even embed it.

Here is an example of a "composite function" that should actually be implemented as multiple definitions.:

This option would be more concise and effective.:

However, it should be noted that the compiler is very effective at optimizing dead-end branches in the code, as in the example of `mynorm'.

It is desirable that the function always returns a value of the same type whenever possible. Consider the following definition:

Although it seems quite harmless, the problem is that 0 is an integer (like Int), and x can be of any type. Thus, depending on the value of x, this function can return a value of one of two types. This behavior is acceptable, and in some cases may be desirable. However, the situation can be easily fixed as follows:

There is also a function oneunit and a more general function 'oftype(x, y), which returns the argument `y, converted to type `x'.

A similar problem of "type stability" is relevant for variables that are reused inside a function.:

The local variable x is an integer at first, but after one iteration of the loop it becomes a floating-point number (as a result of applying the operator /). Therefore, it is more difficult for the compiler to optimize the loop body. There are several ways to solve this problem.:

Many functions are built according to a certain pattern: first, some preparatory work is done, and then many iterations are repeated to perform basic calculations. If possible, it is advisable to separate these basic calculations into separate functions. For example, the following fictional function returns an array of random type:

It would be correct to write it like this:

The Julia compiler specializes the code taking into account the types of arguments at function boundaries, therefore, in the initial implementation, the type a is unknown during the execution of the loop (since it is randomly selected). The second version is likely to run faster, as the inner loop can be recompiled as part of fill_twos! for different types of a.

The second form is also more correct in terms of style and makes code easier to reuse.

This template is repeatedly used in the Base Julia module. For example, you can familiarize yourself with the functions vcat and hcat in abstractarray.jl , as well as with the function fill!, which could be used instead of writing your own implementation of fill_twos!.

Functions such as strange_twos are used when working with data of an undefined type, for example, loaded from an input file that may contain integers, floating-point numbers, strings, or something else.

Let’s say you want to create an N-dimensional array with a size of 3 on each axis. Such an array can be created as follows:

This approach is very effective: the compiler can determine that A is an Array'.{Float64,2}, since it knows the type of value for initialization (5.0::Float64) and the dimension ((3, 3)::NTuple{2,Int}It follows that the compiler can generate very efficient code for further use of the array `A in the same function.

But now let’s assume that you need to write a function that creates a 3×3×…​ array of arbitrary dimension. The following implementation suggests itself:

It works, but (as you can see for yourself with the help of @code_warntype array3(5.0, 2)the problem is that the output type cannot be deduced: the argument `N is an Int type value, and the type inference engine cannot predict its value. This means that the code that uses the result of this function must be conservative: the type must be checked every time A is accessed. This code will run very slowly.

One of the very effective ways to solve such problems is to use removing barriers between functions. However, in some cases, it may be necessary to completely eliminate type instability. One approach may be to pass the dimension as a parameter, for example by using Val{T}() (see the section Types-values):

Julia has a special version of ntuple that accepts an instance of Val{::Int} as the second parameter; by passing N as the type parameter, you tell its "value" to the compiler. Therefore, this version of array3 allows the compiler to predict the return type.

However, using this technique may conceal unexpected subtleties. For example, calling array3 from a function like the following would be useless:

The same problem arises here: the compiler cannot determine n, so it does not know the type of Val(n). Incorrect use of Val can often lead to a decrease in performance. (The use of the above code template is advisable only if the Val is correctly combined with the barrier method between functions, which makes the core function more efficient.)

Here is an example of the correct use of Val.

In this example, N is passed as a parameter, so its "value" is known to the compiler. In fact, Val(T) only works if T is hard-coded, is a literal (Val(3)), or is already specified in the type scope.

Having appreciated the advantages of multiple dispatching, a programmer can get carried away and start applying it in all situations in a row. For example, it can use the following structure to store information:

and then perform dispatching of objects: Car{:Honda,:Accord}(year, args...).

This may be advisable if one of the following conditions is met.

If the latter is true, a function that processes such a homogeneous array can be effectively specialized: The type of each element is known to Julia in advance (all objects in the container have the same specific type), so when compiling the function, Julia can determine suitable method calls (so this will not need to be done at runtime), thereby generating efficient code to process the entire list.

If these conditions are not met, most likely, there will be no benefit. Moreover, the "combinatorial explosion of types" caused by this approach will have the opposite effect. If the type items[i+1] differs from the type item[i], the Julia environment will have to determine the type at runtime, search for the appropriate method in the method tables, determine the appropriate method (by crossing types), establish whether its JIT compilation has already been performed (and perform it, if this is not the case), and then execute the call. In fact, you are instructing the type system and the JIT compilation mechanism to perform an action in the code similar to the branching operator or dictionary search.

A number of test results comparing type dispatch, dictionary search, and the branch operator can be found in https://groups.google.com/forum /#!msg/julia-users/jUMu9A3QKQQ/qjgVWr7vAwAJ[forum discussion].

Perhaps the compile-time effect is even more serious than the runtime effect: Julia compiles specialized functions for each individual Car' object.{Make, Model}. If there are hundreds or thousands of types, then for each function that accepts such an object as a parameter (its own custom function get_year', the universal function `push! from the Base Julia module or any other), hundreds or thousands of variants will be compiled. Each of them increases the cache size of the compiled code, the length of the internal lists of methods, and so on. Excessive enthusiasm for values as parameters can easily lead to resource exhaustion.

Julia multidimensional arrays are stored in memory in columns, that is, they are placed column by column. This can be checked using the vec function or the syntactic construct [:], as shown below (note that the array elements follow in the order [1 3 2 4], not [1 2 3 4]).

This array ordering convention is accepted in many languages, including Fortran, Matlab, R, and others. An alternative to column-based memory allocation is row-based; this convention is accepted in C and Python (numpy), among others. The order in which arrays are placed in memory can seriously affect performance when iterating through them. Remember a simple rule: if an array is stored in memory in columns, the first index changes the fastest. This means that iteration will be faster if the index used in the deepest nested loop is the first in the slice expression. Keep in mind that when indexing an array using :, an implicit loop is performed with alternate access to each element in a certain dimension. For example, columns are retrieved faster than rows.

Consider the following artificial example. Suppose we need to write a function that accepts an object Vector and returns a square matrix Matrix, the rows or columns of which are filled with copies of the input vector. Let’s say it doesn’t matter to us whether rows or columns will be filled with these copies (perhaps the corresponding behavior can be configured in another part of the code). There are at least four ways to implement this (in addition to the recommended built-in function call repeat):

Now we will measure the execution time of each of these functions by entering a vector of size 10000 by 1:

Note that the 'copy_cols` function is much faster than copy_rows'. This is quite expected, because the `copy_cols function takes into account the order in which the Matrix object is placed in memory by columns and fills it in one column at a time. In addition, the 'copy_col_row` function is much faster than copy_row_col because it follows the rule that the first element in the slice expression should be used in the deepest nested loop.

If a function returns an Array or another complex type of object, it may need to allocate memory. Unfortunately, memory allocation and the reverse operation, garbage collection, are often serious bottlenecks.

Sometimes you can avoid having to allocate memory each time a function is called by pre-allocating memory for the output data. Let’s look at the simplest example.

with

Time measurement results:

There are other advantages to preallocating memory. For example, the caller can control the output type of the algorithm. In the example above, if desired, you could pass SubArray instead Array.

Abusing preallocation of memory can make the code less elegant, so it may be necessary to evaluate the appropriateness of its use and measure performance. However, for "vectorized" functions (executed element-by-element), you can use the convenient syntax x.= f.(y) for in-place operations with combined loops without temporary arrays (see the section on dot syntax for vectorizing functions).

Some subtypes Number, such as BigInt or BigFloat, can be implemented as types mutable struct or have mutable components. The arithmetic interfaces in Julia Base usually choose convenience over efficiency in such cases, so their simplest use can lead to suboptimal performance. Package Abstractions https://juliahub.com/ui/Packages/General/MutableArithmetics ['MutableArithmetics'], on the other hand, allows you to use the mutability of such types to write fast code that allocates the necessary volume. 'MutableArithmetics` can also be used to explicitly copy values of mutable arithmetic types when necessary. 'MutableArithmetics` is a custom package and is not associated with the Julia project.

Julia has a special dot syntax, which converts any scalar function into a "vectorized" function call, and any operator into a "vectorized" operator. At the same time, it has a special property: nested "point calls" are combined at the syntactic level into a single loop without allocating memory for temporary arrays. When using .= and similar assignment operators, the result can also be stored in place in a pre-stored array (see above).

In the context of linear algebra, this means that although operations such as vector + vector and vector * scalar' are defined, it may be better to use `vector' instead.+ vector and vector .* scalar, since the resulting loops can be combined with neighboring calculations. For example, consider two functions:

The functions f and fdot' perform the same calculations. However, the function `fdot' (defined using a macro @.`) is performed much faster when applied to an array:

That is, fdot(x) is ten times faster and allocates six times less memory than f(x), because each operation of * and + in f(x)' allocates memory for a new temporary array and it is executed in a separate loop. In this example, the speed of the function `f.(x) is the same as the function fdot(x), but in many situations it is more convenient to simply add a few dots to expressions than to define a separate function for each vectorized operation.

The above-mentioned combination of point loops allows you to express high-performance operations in concise and idiomatic code. However, it is important to remember that the combined operation will be calculated at each iteration of the translation. This means that in some situations, especially in the presence of composite or multidimensional translations, an expression with dot calls may evaluate the function more times than expected. For example, we want to build a random matrix whose rows have a Euclidean norm equal to one. We can write something like the following:

It will work. However, such an expression will actually recalculate sqrt(d[i]) for each element in the string x[i, :], which means that many more square roots are calculated than necessary. To find out exactly which indexes will be iterated on, you can call Broadcast.combine_axes for the arguments of the combined expression. A tuple of ranges will be returned, the entries of which correspond to the axes of the iteration. The product of the lengths of these ranges will be the total number of calls to the combined operation.

It follows that if some components of the translation expression are constant along the axis, for example, sqrt in the second dimension in the previous example, then performance can be improved by forcibly "separating" these components, that is, pre-allocating the result of the translation operation, reusing the cached value along its constant axis. Some of these possible approaches use temporary variables, enclose the components of a dot expression in an identity, or use an equivalent internally vectorized (but not combined) function.

In each of these options, an approximately three-fold acceleration is achieved by highlighting. For large broadcast objects, such acceleration can be asymptotically very large.

In Julia, an array slice expression, such as array[1:5, :], creates a copy of the data (but not if it is used on the left side of the assignment: array[1:5,:] =... performs assignment in place of the corresponding part of the array array). If a lot of operations are performed with a slice, this can be advantageous from a performance point of view, since it is more efficient to work with a small continuous copy than to access the source array by indexes. On the other hand, if a small number of simple operations are performed with a slice, the cost of memory allocation and copying can be significant.

An alternative is to create a "representation" of the array, that is, an array object (SubArray) that references the data of the original array in place without creating a copy. (When writing to the view, the data in the original array also changes.) You can do this for individual slices by calling the function view, or, more simply, for the entire expression or block of code, specifying @views before the expression. For example:

Note that the fview version of the function is three times faster and allocates less memory.

Arrays are stored in memory as contiguous blocks, which makes it possible to vectorize the CPU and reduce the number of memory accesses due to caching. For the same reasons, arrays are recommended to be accessed by columns (see above). Random access patterns and disjointed representations can seriously slow down array operations due to inconsistent memory access.

Copying data accessed randomly into a continuous array before re-accessing it can lead to significant speed gains, as in the example below. In this case, the matrix elements are accessed at random indexes, and then multiplied. Copying to regular arrays speeds up multiplication, despite the cost of copying and allocating memory.

Provided that there is enough memory, the cost of copying a representation to an array is more than offset by the speed gain due to repeated multiplication in a continuous array.

If the calculations involve many small (< 100 elements) arrays of a fixed size (which is known before execution), then you can try using https://github.com/JuliaArrays/StaticArrays.jl [the StaticArrays.jl package]. It allows you to represent such arrays in such a way that unnecessary memory allocations in the heap are not required, and the compiler will be able to specialize the code taking into account the size of the array, for example, by fully deploying vector operations (eliminating loops) and storing elements in CPU registers.

For example, in two-dimensional geometric calculations, there may be many operations with two-component vectors. By using the SVector' type from the StaticArrays package.convenient vector notation and operations such as `norm(3v - w) become available to jl in relation to the vectors v and w, and the compiler can expand the code into a minimal calculation equivalent to @inbounds hypot(3v[1]-w[1], 3v[2]-w[2]).

When writing data to a file (or to another I/O device), building intermediate lines entails additional costs. Instead of:

use this one:

In the first version of the code, the string is first constructed and then written to a file, while in the second version the values are written to the file directly. Also, note that in some cases, string interpolation code may be less readable. Compare:

with this option:

When executing a remote function in parallel mode, the following code:

it will be faster than this:

The first option involves a network access cycle for each workflow, and in the second, two network calls occur: the first is performed by a macro @spawnat, and the second one is a method fetch (or even wait). In addition, the methods fetch and wait are executed sequentially, which generally reduces performance.

The outdated function performs an internal search for a single output of the corresponding warning. This additional search operation can significantly slow down the work, so all cases of using outdated functions should be adjusted according to the instructions in the warnings.

There are a number of less important points that should be considered when working with continuous internal loops.

Sometimes you can improve optimization by declaring certain program properties.

The standard 1:n idiom for indexing an AbstractArray array is unsafe if non-standard indexing is applied in the array, and may lead to a segmentation fault if boundary checking is disabled. Use LinearIndices(x) or eachindex(x) instead (see also the chapter Arrays with custom indexes).

Here is an example of using the @inbounds and @simd annotations at the same time (in this case, @noinline is used so that the optimizer does not disrupt the testing process):

On a computer with a 2.4GHz Intel Core i5 processor, the result will be as follows:

(Performance is measured in GFlop/sec, and the higher the score, the better.)

Here is an example of using all three annotations. This program first calculates the final difference of a one-dimensional array, and then the L2 norm for the result.:

On a computer with an Intel Core i7 processor running at 2.7 GHz, the result will be as follows:

Here, the --math-mode=ieee parameter disables the @fastmath macro so that the results can be compared.

In this case, using @fastmath increases the speed by about 3.7 times. Such an increase is unusual - as a rule, it happens less. (In this particular example, the working set of test data is quite small and fits completely in the processor’s L1 cache, so delays in accessing memory do not matter and the time spent on the calculation is primarily determined by the CPU load. In most real-world programs, the situation is different.) In addition, in this case, this optimization does not affect the result, but usually the result is slightly different. In some cases, especially when using numerically unstable algorithms, the result may vary greatly.

The @fastmath annotation rearranges floating-point expressions, for example, by changing the order of calculations or assuming that some special cases (inf, nan) are impossible. In this case (on this particular computer), the main difference is that the expression 1 /(2*dx) in the 'derive` function is output from the loop (calculated outside it), as if we had written idx = 1 /(2*dx)'. Inside the loop, the expression `+…​ / (2*dx)+ is converted to ... * idx, which is calculated much faster. Of course, both the actual optimization applied by the compiler and the resulting speed gain largely depend on the hardware. You can use the Julia function to find out how the generated code has changed. code_native.

Note that the macro @fastmath also assumes that there are no NaN values during the calculation, which can lead to unexpected behavior.:

Subnormal numbers, previously called https://en.wikipedia.org/wiki/Denormal_number [denormalized], are useful in many situations, but on some hardware they lead to lower performance. Challenge 'set_zero_subnormals(true) allows floating-point operations to interpret subnormal input or output values as zero, which can improve performance on some hardware. Challenge 'set_zero_subnormals(false) strictly prescribes the processing of subnormal values according to the IEEE standard.

The following is an example of the significant impact of subnormal numbers on performance on some hardware.

The result will be something like this:

Note that each even iteration runs significantly faster.

In this example, many subnormal numbers are generated because the values in a form an exponentially decreasing curve that gradually flattens out.

Treating subnormal numbers as zeros should be used with caution, because it can disrupt the operation of some expressions, for example, x-y == 0 implies x == y:

In some applications, an alternative to zeroing out subnormal numbers is to introduce a little noise. For example, instead of initializing object a with zeros, initialize it as follows:

The macro @code_warntype' (or its variant as a function `code_warntype) can sometimes be useful for diagnosing type problems. Here is an example:

Interpreting the macro output @code_warntype', as well as related macros @code_lowered`, @code_typed, @code_llvm and @code_native requires some experience. Your code is presented in a form that has been significantly redesigned during the generation of compiled machine code. Most expressions have type annotations in the form of ::T (where T can be, for example, Float64). The most important feature of the macro '@code_warntype' is that non-specific types are displayed in red. Since this document is written in Markdown format, which does not imply color highlighting, the red text in it is indicated in capital letters.

The first line shows the output return type of the function — Body::Float64. The following lines represent the body of the function f in the SSA Julia intermediate representation format. Numbered fields are labels representing the goals of transitions (via goto) in the code. If you look at the function body, you can see that pos is called first and the return type Union is output.{Float64, Int64}`. It is displayed in capital letters, as it is not specific. This means that it is impossible to accurately determine the returned type pos based on the input types. However, the result of the operation y*x is of type Float64 regardless of whether y is of type Float64 or Int64'. As a result, the result is `f(x::Float64) will not be unstable in type, even if some intermediate calculations are such.

How you use this information is up to you. Obviously, it would be best to fix the pos so as to ensure type stability. In this case, all the variables in f will be specific and performance will be optimal. However, in certain situations, such temporal instability of types does not really matter. For example, if pos is never used in isolation, the fact that the output type is f is stable (for input values of type Float64), will protect further code from the effects of type instability. This is especially important in cases where it is difficult or impossible to eliminate type instability. In such situations, the above tips (for example, about adding type annotations or separating functions) are the best way to contain the effects of type instability. Also keep in mind that even the Base Julia module has functions with type instability. For example, the function findfirst returns the index of the position in the array where the key is found, or nothing if the key is not found — a clear case of type instability. To make it easier to identify cases of type instability that may be relevant, Union' unions containing `missing or `nothing' are highlighted in yellow rather than red.

Below are examples that can help you figure out how to interpret expressions marked as containing non-finite types.:

Consider the following example, which defines an internal function.

The 'abmult` function returns the f function, which multiplies its argument by the absolute value of r'. The internal function assigned to `f' is called a closure. Internal functions are also used in the language for `do blocks and generator expressions.

This style of writing code presents performance challenges. When converting the above code into lower-level instructions, the analyzer rearranges it by extracting an internal function and placing it in a separate code block. "Captured" variables such as r, which are shared between internal functions and the external scope, are also extracted to a dedicated "container" on the heap, accessible by both internal and external functions. The reason is that according to the language specification, the variable r in the inner scope must be identical to r in the outer scope even after changing the r in the outer scope (or in another inner function).

The analyzer was mentioned in the previous paragraph. This refers to the compilation stage where the module containing abmult is loaded for the first time, as opposed to the later stage when it is first called. The analyzer does not know that Int is a fixed type or that the operator r = -r converts the value of Int to another value `Int'. The magic of type inference happens at a later stage of compilation.

Thus, the analyzer does not know that r has a fixed type (Int) or that the value of r does not change after creating the internal function (and that, therefore, the container is not needed). Therefore, the analyzer creates code for a container containing an object of an abstract type, such as Any, which requires type dispatching for each instance of r at runtime. You can check this by applying the macro @code_warntype to the above function. Both packaging and dispatching types at runtime can reduce performance.

If the captured variables are used in a part of the code for which performance is extremely important, then you can use the tips below. First, if it is known that the type of the captured variable does not change, it can be explicitly declared using the type annotation (applied to the variable itself, and not in the right part):

Type annotation partially compensates for the performance loss due to capture, since the analyzer can associate a specific type with an object in the container. Further, if the captured variable does not need to be packaged at all (since it will not be repackaged after the closure is created), this can be indicated using the let block as follows.

The 'let` block creates a new variable r, whose scope is limited by an internal function. The second technique allows you to fully restore performance in the presence of captured variables. Keep in mind: this aspect of the compiler is evolving rapidly, and it is possible that in future versions it will not be necessary to annotate types in such detail to ensure performance. In the meantime, it is possible to automate the addition of let statements, as in abmult3, using packages from community members such as https://github.com/c42f/FastClosures.jl [FastClosures].

The information in this section relates to Julia’s multithreaded code, which performs linear algebra operations in each thread. Indeed, these linear algebra operations involve the use of BLAS/LAPACK calls, which themselves are multithreaded. In this case, it is necessary to ensure that the cores are not overloaded due to two different types of multithreading.

Julia compiles and uses its own copy of OpenBLAS for linear algebra with the number of threads controlled by the environment variable OPENBLAS_NUM_THREADS'. It can be set as a command line parameter when starting Julia, or changed during a Julia session using `BLAS.set_num_threads(N) (the BLAS submodule is exported using `using LinearAlgebra'). Its current value can be accessed via `BLAS.get_num_threads()'.

If the user does not specify anything, Julia tries to choose a reasonable value for the number of OpenBLAS streams (for example, based on the platform, Julia version, etc.). However, it is usually recommended to check and set the value manually. OpenBLAS operates as follows.

When running Julia in multithreaded mode with JULIA_NUM_THREADS=X, it is usually recommended to set OPENBLAS_NUM_THREADS=1. Given the behavior described above, increasing the number of BLAS threads to N>1 can easily lead to performance degradation, especially when N<<X. However, this is just a rule of thumb, and the best way to set each number of threads is to experiment in a specific existing application.

As an alternative to OpenBLAS, there are several other backends that can help with linear algebra performance. Notable examples are https://github.com/JuliaLinearAlgebra/MKL.jl [MKL.jl] and https://github.com/JuliaMath/AppleAccelerate.jl [AppleAccelerate.jl].

These are external packages, so we won’t discuss them in detail here. See the relevant documentation (especially because their behavior regarding multithreading differs from that of OpenBLAS).