Olivier Beaumont

Lionel Eyraud-Dubois

Alena Shilova

Optimal GPU-CPU Offloading Strategies for Deep Neural Network Training

Training Deep Neural Networks is known to be an expensive operation, both in terms of computational cost and memory load. Indeed, during training, all intermediate layer outputs (called activations) computed during the forward phase must be stored until the corresponding gradient has been computed in the backward phase. These memory requirements sometimes prevent to consider larger batch sizes and deeper networks, so that they can limit both convergence speed and accuracy. Recent works have proposed to offload some of the computed forward activations from the memory of the GPU to the memory of the CPU. This requires to determine which activations should be offloaded and when these transfers from and to the memory of the GPU should take place. We prove that this problem is NP-hard in the strong sense, and we propose two heuristics based on relaxations of the problem. We perform extensive experimental evaluation on standard Deep Neural Networks. We compare the performance of our heuristics against previous approaches from the literature, showing that they achieve much better performance in a wide variety of situations.

Introduction

Training for Deep Learning Networks (DNNs) has become a major compute intensive application [START_REF] Szegedy | Going deeper with convolutions[END_REF][START_REF] He | Deep residual learning for image recognition[END_REF][START_REF] Huang | Densely connected convolutional networks[END_REF], typically performed on GPU clusters [START_REF] Goyal | Accurate, large minibatch sgd: Training imagenet in 1 hour[END_REF][START_REF] You | Imagenet training in minutes[END_REF]. The training phase typically involves two traversals of the graph representing the DNN, one in direct order which is called forward propagation and one in reverse order called backward propagation. This incurs high memory usage: the tensors computed during the forward phase, called forward activations, must be kept in memory until the associated backward operation is performed, since they are needed to compute the gradients used to update the weights. Therefore, memory issues become crucial when performing training in DNNs, and the memory limitation of current hardware often prevents data scientists from considering larger models, larger image sizes or larger batch sizes [START_REF] Bulò | In-place activated batchnorm for memory-optimized training of dnns[END_REF][START_REF] Pleiss | Memory-efficient implementation of densenets[END_REF].

For instance, when using ResNet101 with relatively small images of size 224 × 224 with 3 channels and a batch size of 32, the resulting size during training is around 5GB. For applications which require to detect small objects in the images [START_REF] Carranza-Rojas | Going deeper in the automated identification of herbarium specimens[END_REF], the image resolution must be increased, and the memory required for storing activations increases quadratically with the image resolution. The situation is even worse when moving to 3D object recognition [START_REF] Ghadai | Multi-level 3d cnn for learning multi-scale spatial features[END_REF][START_REF] Feng | Gvcnn: Group-view convolutional neural networks for 3d shape recognition[END_REF][START_REF] Su | Multi-view convolutional neural networks for 3d shape recognition[END_REF] or video based DNNs such as 3D-Resnet [START_REF] Hara | Can spatiotemporal 3d cnns retrace the history of 2d cnns and imagenet?[END_REF] or CDC [START_REF] Shou | Cdc: Convolutional-de-convolutional networks for precise temporal action localization in untrimmed videos[END_REF]. In this context, the input consists in a large number of frames (64 for instance), thus inducing a huge memory cost even for small batch sizes. However, even for applications where a small batch size can fit into memory, using larger batch sizes can actually improve the behavior of the training phase. Indeed, batch-normalization [START_REF] Ioffe | Batch normalization: Accelerating deep network training by reducing internal covariate shift[END_REF] is increasingly used in many DNNs, and has been proved to allow faster and more stable training as demonstrated in [START_REF] Santurkar | How does batch normalization help optimization?[END_REF]. Using too small batch sizes degrades these batch statistics [START_REF] Kusumoto | A graph theoretic framework of recomputation algorithms for memory-efficient backpropagation[END_REF] and prevents to benefit from batch normalization techniques.

Many approaches have been proposed in the literature in order to circumvent this memory issue. In this paper, we focus on an offloading approach (also called memory swapping), which consists in reducing memory usage on the GPU by transferring some activations to the memory of the CPU, which is expected to be at least one order of magnitude larger. The corresponding algorithmic question is to determine which activations should be offloaded and when, and also when offloaded activations should be brought back (prefetched) from the memory of the CPU to the memory of the GPU. This approach has been recently considered in [START_REF] Rhu | vdnn: Virtualized deep neural networks for scalable, memory-efficient neural network design[END_REF][START_REF] Garg | Dynamic memory management for gpu-based training of deep neural networks[END_REF], where the authors advocate the general idea and propose several static and dynamic heuristics to decide which activations should be offloaded. In this paper, we provide a deeper analysis of this problem. More specifically, we prove that the general problem, even for sequential models, is strongly NP-complete where only fully integral data transfers are possible and we analyze two relaxations of the problem for which we can derive optimal algorithms. These algorithms can then be used as heuristics for the general problem.

The rest of the paper is organized as follows. In Section 2, we discuss previous works regarding GPU memory offloading, as well as other techniques to reduce memory usage during the training phase. In Section 3, we present the model and notations used throughout the paper, and assess the complexity of the problem. In Section 4, we propose a first relaxation where activations can be partially or completely offloaded into the memory of the CPU, and derive an optimal strategy. In Section 5, we consider the case where partial offloading is not possible, but where communications can be interrupted, and we present a dynamic programming algorithm to find the optimal schedule. In Section 6, we provide experimental results and we assess the efficiency of our heuristics against the previous approach [START_REF] Garg | Dynamic memory management for gpu-based training of deep neural networks[END_REF], before presenting conclusions and perspectives in Section 7.

Related Work

In order to reduce the memory usage of storing the forward activations on a single GPU node, we can identify two kinds of approaches: checkpointing or offloading.

Checkpointing techniques

Checkpointing techniques consist selecting only a few activations that are kept in memory, and then to dynamically recompute the others at runtime. This allows to explore a tradeoff between memory usage and computational cost. This problem has been first widely studied in the context of Automatic Differentiation (AD) [START_REF] Griewank | On automatic differentiation[END_REF]. In the context of AD, computational networks can be seen as long chains, in which the forward activation corresponding to the i-th stage of the chain has to be kept into memory until the i-th backward stage. Checkpointing algorithms consist in determining in advance which forward checkpoints should be kept into memory and which ones should be recomputed from stored checkpoints when performing the backward phase.

Adjoints computation is at the core of many scientific applications, from climate and ocean modeling [START_REF] Adcroft | Mitgcm user manual[END_REF] to oil refinery [START_REF] Brubaker | Engineering Design Optimization using Calculus Level Methods[END_REF]. Many studies have been performed to determine optimal checkpointing strategies for Automatic Differentiation in different contexts. In the case of homogeneous chains, closed form formulas providing the exact position of checkpoints have been proposed [START_REF] Griewank | Algorithm 799: Revolve: an implementation of checkpointing for the reverse or adjoint mode of computational differentiation[END_REF], and in most cases the general algorithmic ingredient is to derive optimal checkpointing strategies using Dynamic Programming [START_REF] Griewank | Algorithm 799: Revolve: an implementation of checkpointing for the reverse or adjoint mode of computational differentiation[END_REF].

The use of checkpointing strategies has recently been advocated for DNN in several papers [START_REF] Gruslys | Memory-efficient backpropagation through time[END_REF][START_REF] Chen | Training deep nets with sublinear memory cost[END_REF][START_REF] Kukreja | Backpropagation for long sequences: beyond memory constraints with constant overheads[END_REF][START_REF] Beaumont | Optimal checkpointing for heterogeneous chains: how to train deep neural networks with limited memory[END_REF][START_REF] Kumar | Efficient rematerialization for deep networks[END_REF][START_REF] Feng | Optimal gradient checkpoint search for arbitrary computation graphs[END_REF][START_REF] Jain | Checkmate: Breaking the memory wall with optimal tensor rematerialization[END_REF]. An implementation of a simple periodic strategy is provided in PyTorch [START_REF] Paszke | Automatic differentiation in pytorch[END_REF][START_REF]Periodic checkpointing in pytorch[END_REF], but it is only optimal for the restricted case of homogeneous chains, whereas DNN models are in general more complicated.

Memory Offloading

Offloading is a potentially complementary approach which consists in offloading some of the forward activations from the memory of the GPU to the memory of the CPU, which is expected to be much larger [START_REF] Rhu | vdnn: Virtualized deep neural networks for scalable, memory-efficient neural network design[END_REF][START_REF] Garg | Dynamic memory management for gpu-based training of deep neural networks[END_REF]. In [START_REF] Rhu | vdnn: Virtualized deep neural networks for scalable, memory-efficient neural network design[END_REF], the authors propose a simple and effective mechanism of memory virtualization, that nevertheless introduces unnecessary idle time by enforcing some synchronization between data transfers and computations of later forward activations. This approach has been later improved in [START_REF] Garg | Dynamic memory management for gpu-based training of deep neural networks[END_REF] by the design of techniques to deal with memory fragmentation. Nevertheless, in both papers, the algorithmic strategies to decide which activations to offload into the memory of the CPU are relatively straightforward. Proposed strategies consist in trying to offload either all activations or only those that correspond to convolutional layers. Indeed, convolutional layers are known to induce a large computational time with respect to their input size, which make them good candidates to overlap offloading and processing. Several follow-up works offer improvements over this first attempt. In order to reduce the overhead incurred by the communications, some authors [START_REF] Rhu | Compressing dma engine: Leveraging activation sparsity for training deep neural networks[END_REF] recommend to add compression to decrease the communication time, while others [START_REF] Kwon | Beyond the memory wall: A case for memory-centric hpc system for deep learning[END_REF] design a memory-centric architecture to help with data transfers. In [START_REF] Meng | Training deeper models by gpu memory optimization on tensorflow[END_REF][START_REF] Le | Tflms: Large model support in tensorflow by graph rewriting[END_REF], the authors implement memory virtualization by manipulating the computational graphs and inserting special operations called swap in and swap out that send the activations in and out of GPU memory. Such an approach can be applied to any arbitrary Computation Graph that represent neural network training graphs. The authors of [START_REF] Le | Tflms: Large model support in tensorflow by graph rewriting[END_REF] improve the candidate selection and prefetching mechanisms by introducing thresholds to filter out different possibilities. Moreover, some works try to combine offloading with other memory optimizing techniques. Memory swapping and memory pooling are implemented together in [START_REF] Zhang | Efficient memory management for gpu-based deep learning systems[END_REF], where candidates for swapping are found by assigning priority scores to all activations. Finally, gradient checkpointing is combined with the simple offloading approach from [START_REF] Rhu | vdnn: Virtualized deep neural networks for scalable, memory-efficient neural network design[END_REF] in [START_REF] Wang | Superneurons: Dynamic gpu memory management for training deep neural networks[END_REF].

F 0 F 1 • • • F L-2 F L-1 F L B 0 B 1 B 2 • • • B L-1 B L x 0 x 1 x 2 x L-2 x L-1 x L x L+1 y L+1 = 1 y L y L-1 y 3 y 2 y 1 y 0 x 0 x 1 x 2 x L-1 x L x 1 x 2 x 3 x L-1 x L x L+1
As a complement to these practical approaches, in this paper we perform the first theoretical analysis of the underlying optimization problem and present both a complexity proof and optimal solutions to two of its relaxations.

Model and Complexity

Computation Model

We consider the training phase of sequential DNNs, as depicted on Figure 1. This training phase consists of two types of computations: forward propagations (F i ) 1≤i≤L and backward propagations (B i ) 1≤i≤L . The forward step F i requires x i as input, and computes x i+1 . The backward step B i requires x i+1 , x i and y i+1 as inputs, and computes y i . The objective of the elementary training phase is to perform the whole computation and to obtain y 0 in the smallest possible time. This computation is performed on a processing device (typically a GPU or TPU) with limited memory M GPU . We denote u Fi the time to process F i , and u Bi the time to process B i . As mentioned before, the training phase is very memory intensive: since they will be needed for the backward phase, all x i values must be stored during the forward phase, and they can only be freed once their corresponding B i operation has been performed.

We will use the following memory model. Each data (x i and y i ) has a given memory usage, denoted respectively by |x i | and |y i |. To perform an operation (either F i or B i ), it is necessary to have all inputs stored into memory, to reserve the memory space to store the output, and to reserve space for the temporary memory usage of the operation, denoted with ex F i for F i and ex B i for B i . For example, running the F i operation requires to have at least

|x i | + |x i+1 | + ex F i memory available.
In order to decrease memory usage, we assume that it is possible to offload some of the data to another memory storage (typically the main memory of the machine). The size of this memory is assumed to be large enough to store all the results and thus is not a constraint; but the speed of data transfers is limited by bandwidth β. The offloaded data can then be prefetched during the backward phase, so that it is available when needed to perform the corresponding backward operation. Such memory hierarchy has been considered by [START_REF] Rhu | vdnn: Virtualized deep neural networks for scalable, memory-efficient neural network design[END_REF][START_REF] Garg | Dynamic memory management for gpu-based training of deep neural networks[END_REF] as well, i.e. there are one GPU with limited memory and one CPU with large enough memory to store all activations of some arbitrary neural network and both are connected with the network with the bandwidth β, which we assume is fully used for any communications. More complicated cases such as multiple GPU and one CPU are out of scope of this paper and they will be left for the future work.

Finally, we make the following assumptions:

(A) y i are not transfered;

(B) Only one transfer can occur at a time;

(C) Transferring x i takes time |x i |/β, i.e. bandwidth is fully used;

(D) Transfers and computations can be completely overlapped, and do not interfere;

(E) x i needs to be stored in memory in its entirety throughout the transfer: during the offloading, the memory is only released after the complete transfer, and during the prefetching, the memory is reserved as soon as the transfer begins.

We can state the decision problem associated to offloading.

Problem 1 (Offloading). Consider a training phase with L operations, with processing times u Fi and u Bi , data sizes |x i | and |y i |, temporary memory usage ex F i and ex B i , where 0 ≤ i ≤ L. Is it possible to perform this computation on a processing device with memory M GPU and bandwidth β between processing device and main memory, with an execution time at most T ?

Preliminary results and lower bound

Proposition 1. For fixed decisions of which data to offload, and in which order transfers should be performed, the best schedule is obtained with a no-wait policy, where each action (computation and data transfers) is performed as early as possible, as soon as data is available and there is enough memory.

Given the activation sizes and the temporary memory usage, it is easy to compute the total amount of used memory (both on the computing device and in the additional memory) during the execution of each operation. We denote by m Fi or m Bi the amount of data required to be stored on both devices to perform F i or B i respectively. Let us additionally denote as M peak the maximum of these values.

m Fi = ex F i + j≤i+1 |x j | m Bi = ex B i + |y i | + |y i+1 | + j≤i+1 |x j | M peak = max max 0≤i≤L m Bi , max 0≤i≤L m Fi
Since any valid schedule must process the operation which achieves the memory peak while using at most M GPU memory on the computing device, the following result holds.

Proposition 2. The amount of data offloaded by any valid schedule is at least M peak -M GPU .

Since any valid schedule must perform all computations, and must transfer at least this amount of data twice (for offloading and prefetching), the following lower bound on the optimal makespan holds true.

Proposition 3. The value LB = max( i u Fi +u Bi , 2 M peak -MGPU β
) is a lower bound on the optimal makespan.

F 0 F 1 • • • F 3m-1 1 F 3m F 3m+1 • • • 1 F 5m-2 F 5m-1 B 0 B 1 • • • B 3m-1 m B 3m B 3m+1 • • • B 5m-2 B 5m-1 u 1 u 2 u 3 u 3m 0 0 V 0 V 0 0 0 0 0 0 0 0 0 0 0 u 1 u 2 u 3 u 3m 0 0 V 0 V u 2 u 3 0 0 V V 0 Figure 2:
The instance of Problem 1 used in the reduction from 3-partition. Activation sizes are indicated on the edges and non-zero execution times of operations are written inside the corresponding nodes.

Complexity Results

Theorem 1. Problem 1 is strongly NP-complete.

Proof. Problem 1 clearly belongs to NP: given the start time of all forward and backward operations, and the set of offloaded data with the corresponding start time of transfers, checking that the schedule satisfies all constraints can be done in linear time. We prove that this problem is strongly NP-hard by a reduction from the 3-partition problem: given a set of integers {u 0 , u 2 , . . . , u 3m-1 } such that i u i = mV , is it possible to partition it into m parts {S 1 , . . . , S m } so that for any j ≤ m, |S j | = 3 and i∈Sj u i = V . This problem is known to be NP-hard in the strong sense.

Given an instance of 3-partition, we consider the following instance of Problem 1, depicted on Figure 2:

• L = 5m, β = V , M GPU = mV , T = 2m; • u Fi = 0 and |x i | = u i for 1 ≤ i < 3m; • u Fi = 1 and |x i | = 0 for i = 3m + 2k, k ∈ [|0, m|]; • u Fi = 0 and |x i | = V for i = 3m + 2k + 1, k ∈ [|0, m|];
• u Bi = 0 and |y i | = 0 for all i, except u B3m = m.

We claim that this instance can be scheduled in time T = 2m if and only if the 3-partition instance is positive.

Let us first assume that there exists a solution to the 3-partition instance, i.e. sets (S j ) 1≤j≤m such that i∈Sj u i = V . We can build a schedule which starts F 3m+2k at time k for 0 ≤ k < m, and executes B 3m from time m to time 2m. At time 0, before the execution of F 3m , the memory usage is exactly mV = i u i . During the execution of F 3m+2k , activations x i for i ∈ S k are transferred. Since β = V , this takes time exactly 1. The memory used at the end of F 3m+2k is thus (m -1)V , which allows to immediately start F 3m+2k+1 . At the end of the forward phase, the memory is filled with m activations of size V . At the beginning of B 3m , the memory is empty: all activations of size u i can be prefetched during the execution of B 3m , allowing to finish the backward phase. This schedule induces no idle time, and finishes in time exactly T = 2m.

Let us now assume that there exists a valid schedule of duration T = 2m, i.e. without any idle time on the processing device. For j < m, let us define the set S j as the indices of the activations whose transfers are included in the execution of F 3m+2j . Since F 3m+1 starts immediately after the end of F 3m , and since memory is only released once the transfer has been completed, the amount of data sent during F 3m is at least V . Since β = V and u F3m = 1, the amount of data is exactly V , thus i∈S0 u i = V . The same argument applies for all j < m, which shows that the sets S j are a valid solution for the 3-partition instance, and completes the proof.

Theorem 1 shows that even when we know which activations should be offloaded, it is difficult to decide the order in which the transfers should be done. Because of this negative complexity result, we study two different relaxations of Problem 1 in the next sections, by relaxing the constraints stating that activations must be sent in entirety before the corresponding memory can be released. This allows to compute optimal solutions in reasonable time, and the resulting algorithms can then be used as heuristic solutions for Problem 1.

Fractional Relaxation

In a first relaxation, let us consider that it is possible to perform partial offloading: any communication can be stopped at any time, and the data that has been transferred up to that time can be released from memory, even if the rest of the activation is still present on the computing device. With this model, it is possible to compute an optimal solution with a greedy algorithm. Let us first prove results about the structure of optimal solutions, and then use that structure to design an optimal greedy algorithm.

Structure of Optimal Solutions

In this section, let us analyze special eager schedules.

• A schedule is said eager if it offloads the first k activations x 0 , x 1 , . . . , x k (where the last one can be partially offloaded).

• A schedule is said ordered if the data is offloaded in order of increasing indices, and prefetched in order of decreasing indices.

Lemma 1. Any valid solution S can be transformed into a eager and ordered solution S with the same makespan.

Proof. Let us denote by M off the amount of activation data offloaded by the schedule S, and let us consider in S the time intervals I off spent offloading data, and the time intervals I fetch spent prefetching data. Let us consider the schedule S in which all operations and data transfers are performed at the same instants as in S, only changing which data is transferred. The first intervals of I off are used to transfer x 0 (since it is possible to stop any communication at any time, using several intervals to transfer x 0 is not a problem), the next ones are used to offload x 1 , and so on, until the amount M off is reached, and similarly for the prefetched data, in reverse order. Clearly S is eager and ordered. Since the x i values become available in the forward phase by order of increasing indices, and are consumed in the backward phase by order of decreasing indices, it is clear that transfers in S are valid: an activation is offloaded only after having been produced, and in the backward phase an activation is prefetched before being used. Furthermore, since transfers occur at the same instants and at the same speed as in S, the memory usage of S is exactly the same as the memory usage of S at any instant. The modified S schedule is thus valid.

Greedy Algorithm

According to this result, we consider only eager and ordered schedules. It is thus sufficient to find the amount of offloaded data which results in the smallest makespan. The next result shows that it is best to offload the least possible amount of data. Lemma 2. Let S and S be no-wait, ordered and eager schedules which offload a quantity of data Q and Q respectively, with Q < Q . Then the makespan of S is not larger than the makespan of S .

Proof. Let us consider the schedule S obtained from S by removing all transfers (offload and prefetch) corresponding to data between Q and Q in the eager order. Since S is valid, all data dependencies are satisfied in S . Let us now prove that S also fulfills the memory constraint.

Consider any time instant t in schedule S , and let us denote by m CPU (t) and m GPU (t) the amount of data stored on the CPU and GPU by S . If m CPU (t) ≤ Q, then the data stored on GPU in S and S are the same, so S is valid at instant t. If m CPU (t) > Q, then the amount of data stored on the GPU in schedule The schedule S is thus a valid, eager and ordered schedule which offloads a quantity of data Q. The schedule S offloads the same data in the same order; since S is no-wait, by Proposition 1 the makespan of S is not larger than the makespan of S , which is equal to the makespan of S .

S is m GPU (t) = m GPU (t) + m CPU (t) -Q. Since S is valid, Q ≥ M peak -M GPU . Furthermore, by definition of M peak , m GPU (t) + m CPU (t) ≤ M peak . Thus, m GPU (t) ≤ m GPU (t) + m CPU (t) + M GPU -M peak ≤ M GPU F 0 F 1 • • • F n-1 1 F n F n+1 F n+2 B 0 B 1 • • • B n-1 1 B n B n+1 B n+2 u 1 u 2 u 3 u n 0 0 V 0 0 0 0 0 0 0 0 0 u 1 u 2 u 3 u n 0 0 V u 2 u 3 0 0 V 0
With these two lemmas, since M peak -M GPU is a lower bound on the amount of data that any schedule has to offload, we can characterize an optimal schedule for this relaxed problem.

Theorem 2. For a given instance, the no-wait, eager, ordered schedule which offloads a quantity M peak -M GPU of data is optimal.

By rounding up the number of offloaded activations, this result provides a heuristic for the original integral problem, that we call Greedy. The Greedy heuristic returns the no-wait, eager, ordered schedule which offloads (entirely) the first k activations, where k is the smallest index such that

i≤k |x i | ≥ M peak -M GPU .
However, it may happen that this Greedy schedule offloads too much data because of the rounding procedure. In the next section, we thus analyze a more sophisticated relaxation in order to obtain a more precise algorithm.

Fractional Communications

Let now consider another formulation of Problem 1, in which an activation must be either entirely offloaded or not offloaded at all. However, it is still allowed to stop a communication at any time and resume it later. In this section, we first prove that this problem is NP-complete in the weak sense, and then propose a pseudo-polynomial optimal algorithm based on Dynamic Programming.

Complexity

Problem 2 (Offloading with interruptions). Consider a training phase with L operations, with processing times u Fi and u Bi , data sizes |x i | and |y i |, temporary memory usage ex F i and ex B i , where 0 ≤ i ≤ L. Is it possible to perform this computation on a processing device with memory M GPU and bandwidth β between the processing device and the main memory, with an execution time at most T , if communications can be interrupted?

Let us first note that Proposition 1 also holds for this problem (it is always better to schedule with a no-wait policy). We can also state a result similar to the one of the fully fractional case. Lemma 3. Any valid solution S can be transformed into an ordered solution S with the same makespan.

The proof is the same as the one of Lemma 1: transforming S using the correct order provides a valid schedule. The result is weaker, because an eager schedule which offloads the same data might not be valid for Problem 2 (the last activation might not be fully offloaded). Let us prove that it is NP-hard by reduction for the 2-Partition problem, which can be stated as: given n positive integers u 1 , u 2 , . . . , u n such that i u i = 2V , is is possible to partition them in two subsets S 1 and S 2 such that i∈S1 u i = i∈S2 u i = V ?

Given an instance of 2-Partition, let us consider an instance I of Problem 2, depicted on Figure 3, with L = n + 3, β = V and M GPU = 2V . All values y i have size 0, and all operations F i and B i take time 0, except for F n and B n whose duration is 1. The input size of the first n operations are given by values u i , and the input size of F n+2 is |x n+2 | = V . Note that I can be computed in polynomial time. We claim that I admits a valid schedule of length T = 2 if and only if the instance of 2-Partition has a solution.

If the 2-partition problem has a solution, then there exist subsets S 1 and S 2 = S \ S 1 such that

i∈S1 u i = V .
It is thus possible to offload all the corresponding activations x i during operation F n (since β = V ), and then to prefetch them during operation B n . This allows to perform F n+1 , F n+2 , B n+2 , B n+1 immediately after F n , since only a quantity V of data from the first activations is stored. Once B n has been performed, all other B i can be performed as well, which results in a schedule of length 2.

On the other hand, let us assume that there exists a schedule of length 2. Let us denote as S 1 the set of indices corresponding to activations that are offloaded in that schedule (remember that each activation is either offloaded completely or not at all), and let Q = i∈S1 u i . Since M peak = 3V , it is clear that Q ≥ V . Since β = V , the makespan of the schedule is at least 2Q V (this is the time it takes to offload and prefetch S 1 ), thus Q ≤ V . This implies that S 1 is a solution of the 2-Partition instance, which completes the proof.

Structure of Optimal Solutions

According to Lemma 3, our objective is now to find the best ordered schedule. In this section, we derive properties of all ordered and no-wait schedules, which will allow to obtain a dynamic programming algorithm in the next section.

Forward and Backward phases

Let us consider any ordered, no-wait schedule S. Let M Fi denote the GPU memory occupied at the end of F i-1 (it should contain all data not offloaded at this instant, plus x i which is the output of F i-1 ). Let ∆ Fi denote the amount of data from x 0 , . . . , x i-1 that S offloads after the end of F i-1 . If this amount is zero, let us denote by Av F the time between the end of the last offload and the end of F i-1 , and let ∆ Fi = -Av F • β. Moreover, let us set ∆ + Fi = max{0, ∆ Fi }. We aim to characterize the delay F i between the end of F i-1 and the start of F i . These notations are depicted in Figure 4.

Let us first remark that since S is a valid schedule, there is enough memory to process F i at some point, which means that ∆ + Fi needs to be large enough,

M Fi -∆ + Fi + ex F i + |x i+1 | ≤ M GPU If M Fi + ex F i + |x i+1 | ≤ M GPU ,
then F i can start immediately after the end of F i-1 , and since S is no-wait, then F i = 0. Otherwise, processing F i can start as soon as enough memory has been released by offloading data at rate β. This yields F i =

M F i +ex F i +|xi+1|-M GPU β
. In summary,

F i = max 0, M Fi + ex F i + |x i+1 | -M GPU β (1)
Let us now derive recursive equations to obtain M Fi+1 and ∆ Fi+1 from M Fi and ∆ Fi . These equations depend on whether x i is offloaded in S.

If x i is offloaded, then the amount of data ready to be offloaded at the end of F i-1 is ∆ + Fi + |x i |. Until the end of F i , the amount of data that can be offloaded is at most ( F i + u Fi )β. Hence we obtain

∆ Fi+1 = ∆ + Fi + |x i | -( F i + u Fi )β (2)
M Fi+1 = M Fi + |x i+1 | -min ∆ + Fi + |x i |, ( F i + u Fi )β . (3)
F i F i-1 F i+1 Execution Communication t memory used M Fi M GPU M Fi+1 slope: -1 β ex F i ex F i + |x i+1 | F i F i+1 ∆ F i β |xi| β ∆ F i+1 β
Figure 4: Notations used in the Forward phase, assuming x i is offloaded and

∆ Fi + |xi| B ≥ F i + u Fi .
If x i is not offloaded, we can write similar equations, except that |x i | is not added to the amount of data to be offloaded. This yields

∆ Fi+1 = ∆ Fi -( F i + u Fi )β (4)
M Fi+1 = M Fi + |x i+1 | -min ∆ + Fi , ( F i + u Fi )β (5)
Let us now derive similar results about the backward phase. We first modify S to process all backward operations and perform all prefetching operations as late as possible without changing the makespan of the schedule. We then define M Bi as the GPU memory occupied right before processing B i-1 (thus it does not take into account the output of B i-1 , which is y i-1 ). Let us also define ∆ Bi as the amount of data from x L , x L-1 , . . . , x i ) that S prefetches before starting B i-1 , and if this amount is zero, then ∆ Bi = -Av B • β, where Av B is the time between the start of B i-1 and the start of the first prefetch operation. Finally, let B i denote the delay between the end of B i and the start of B i-1 .

With the same reasoning as above, we obtain

B i = max 0, M Bi + ex F i + |x i+1 | + |y i+1 | -M GPU β (6)
∆ Bi+1 = |x i | + max 0, ∆ Bi -( B i + u Bi )β if x i is offloaded (7) ∆ Bi+1 = ∆ Bi -( B i + u Bi )β otherwise (8)
Computing M Bi+1 is not necessary, as one can notice that for all i, M Bi -

∆ Bi = |y i | + |x i | + j<i,j not offloaded |x j |, and M Fi -∆ Fi = |x i | + j<i,j not offloaded |x j |. Thus, M Bi -∆ Bi = |y i | + M Fi -∆ Fi
, which allows to compute M Bi once all three other values are known.

Idle time between phases

With these derivations, we have accounted for all idle times in schedule S, except for a possible idle time G between the end of F L and the start of B L . Several situations are possible. If S performs no offloading after F L , and no prefetching before B L (i.e. ∆ F L+1 and ∆ B L+1 ≤ 0), then no idle time occurs and G = 0. If S performs both kind of transfers (i.e. both ∆ F L+1 and ∆ B L+1 are positive), then the idle time corresponds to transferring the total data, G = ∆ F L+1 +∆ B L+1 β . Otherwise, if for example ∆ F L+1 > 0 and ∆ B L+1 ≤ 0, the schedule S can perform offloading during the first backward operations. This is possible only if enough memory is available to perform the operations, so it may result in some idle time if it is necessary to wait until enough data has been offloaded. Let us denote by Av B = -

∆ B L+1 β
the time between the start of B L and the start of the first prefetch operation, and let us set U B j = L i=j+1 u Bi . According to the above derivations, since no prefetching occurs, all the operations performed during this time have no idle time between them. These operations are thus all the B j such that U B j < Av B . We can compute the available memory at any time t after the end of F L , where we pretend that no x i is removed from memory until Av B : M (t) = (M GPU -M F L+1 ) + tβ. At the start of B j , the memory available is M ( G + L i=j+1 u Bi ). Starting B j requires to have at least R B j = ex B j + |y j | + |y j-1 | -i>j+1 |x i | available memory, where we account for the removal of x i s. This provides an upper bound on

G : G ≥ R B j +M F L+1 -M GPU β - L i=j+1
u Bi This is true for all the B j operations that occur during Av B . Another upper bound on G is given by the fact that the offloading must finish before the prefetching starts, G ≥

∆ F L+1 β -Av B .
These are the only constraints on G , hence we have

G = max   0, ∆ F L+1 β -Av B , max j≤L U B j <Av B R B j + M F L+1 -M GPU β -U B j    . ( 9
)
On the other hand, if ∆ B L+1 > 0 and ∆ F L+1 ≤ 0, schedule S can perform prefetching during the last forward operations. Let us denote by Av F = -

∆ F L+1 β
the time between the end of the last offload and the end of F L ,

U F j = L i=j+1 u Fi , and R F j = ex F j -i>j+1 |x i |.
Algorithm 1 Dynamic Programming Algorithm for Fractional Communications 1) and ( 6) 2), ( 3) and ( 7)) 4), ( 5) and ( 8))

Idle i ← HashTable() for 0 ≤ i ≤ L Idle 0 (|x 0 |, 0, 0) = 0 for i = 0, . . . , L do for M Fi , ∆ Fi , ∆ Bi ∈ Idle i do M Bi ← |y i | + ∆ + Bi + M Fi -∆ + Fi if M Fi + ex F i -∆ Fi + |x i+1 | ≤ M GPU and M Bi + ex B i + |x i+1 | + |y i+1 | -∆ Bi ≤ M GPU then Compute F i , B i from equations (
Compute M F , ∆ F , ∆ B if x i is offloaded (equations (
Idle i+1 (M F , ∆ F , ∆ B ) ← min Idle i+1 (M F , ∆ F , ∆ B ), Idle i (M Fi , ∆ Fi , ∆ Bi ) + F i + B i Compute M F , ∆ F , ∆ B if x i is not offloaded (equations (
Idle i+1 (M F , ∆ F , ∆ B ) ← min Idle i+1 (M F , ∆ F , ∆ B ), Idle i (M Fi , ∆ Fi , ∆ Bi ) + F i + B i for M F , ∆ F , ∆ B ∈ Idle L+1 do
Compute G according to equations ( 10) and ( 9)

TotalIdle(M F , ∆ F , ∆ B ) ← Idle L+1 (M F , ∆ F , ∆ B ) + G Get M * F , ∆ * F , ∆ * B which minimizes TotalIdle(M F , ∆ F , ∆ B ) Backtrack in Idle L+1 , .
. . , Idle 0 to obtain optimal offload decisions With considerations similar as above, we obtain

G = max   0, ∆ B L+1 β -Av F , max j≤L U F j <Av F R F j + M B L+1 -M GPU β -U F j    (10)
The results obtained in this section show that in order to compute how future offloading decisions affect the idle time of a schedule, one only needs to know the values of M Fi , ∆ Fi , M Bi , ∆ Bi . The exact decisions of which data from x 0 , . . . , x i-1 has actually been offloaded is not required. This allows to design a Dynamic Programming algorithm to identify the offloading decisions that induce the smallest idle time.

Resulting Algorithm

Dynamic Programming Algorithm

To formalize the dynamic programming algorithm, let us define Idle(i, m, d F , d B ) as the smallest possible sum of idle times between (i) the start of the schedule and the end of F i-1 and (ii) the start of B i-1 and the end of the schedule, for all schedules S such that

M Fi = m, ∆ Fi = d F , ∆ Bi = d B .
Any schedule starts with a memory occupation of |x 0 |, and no idle time, so we can define Idle(0, |x 0 |, 0, 0) = 0, and Idle(0, m, d Once Idle L+1 is computed, TotalIdle can be found by adding the corresponding idle time G between the forward and backward phases. Then, the smallest value in TotalIdle is the smallest possible idle time for any ordered, no-wait schedule. Finally, we can identify which offload decisions have led to this idle time, and then obtain the description of the corresponding schedule.

F , d B ) = ∞
The number of values kept in the hash table can be bounded in the following way: M F and M B are between 0 and M GPU , ∆ F is betweeni u Fi β and i x[i], and ∆ B is fixed once the three other values are specified. The number of possible values is thus

O(M 2 GPU ( i x[i] + u Fi β)), and the complexity of Algorithm 1 is O(LM 2 GPU ( i x[i] + u Fi β))
, which is indeed pseudo-polynomial. This optimal algorithm for the fractional communications model can be turned into a heuristic for the original problem, which we call DynProg. DynProg computes the optimal set of activations for the fractional communications model with Algorithm 1, and outputs the no-wait, ordered schedule which offloads exactly these activations.

Discretization Scheme

To keep its computing time reasonable, we also include in DynProg a discretization scheme by fixing a number S of memory slots (S = 500 is a reasonable value in practice), where each slot has size M S . All memory sizes are then expressed as an integer number of slots, by rounding up if necessary, setting for example ex

F i = ex F i • S M .
The values ex F i , ex B i , and |y i | are rounded up in this way. It is also necessary to discretize the values of u Fi β and u Bi β, but these values need to be rounded down to ensure the feasibility of the produced solution in the original problem. However, since these values are always used as partial sums i j=0 u Fj β, we can perform a more precise rounding procedure: the partial sums are discretized by setting σ F j = S M • i j=0 u Fj β , and then the individual values can be recovered by

u Fi β = σ F i+1 -σ F i .
The same procedure is applied to obtain u Bi β. Finally, we apply special care to the discretization of the |x i | values. It is indeed crucial to obtain values as close as possible to the original values: rounding up directly the values in the same way as we did for the |y i | values ensures that the resulting solution is always feasible, but may require to offload significantly more data to ensure that all the non-offloaded activations fit into memory. On the other hand, rounding down some of the values may result in an unfeasible selection of non-offloaded data. We thus adopt an iterative scheme: we start with an optimistic rounding of the values, and if the resulting solution is not feasible, we increase by one the discretized value |x i | which is below, but closest to |x i |. Since optimal solutions of the dynamic program often choose to offload the first activations, the optimistic rounding is obtained by rounding up the partial sums

Experimental Analysis

This section presents experimental results obtained on three different kinds of networks, whose implementation is available in the torchvision package of PyTorch: ResNet, DenseNet, and Inception v3.

Experimental Setting

We have slightly modified these networks to represent them as linear chains, by grouping each non-linear part of the graph in a virtual layer. We have obtained the values of u F , u B , ex F , ex B , and the sizes of x i and y i by performing measures on sample data. These measurements were performed on a node equipped with a Nvidia Tesla V100-PCIE GPU card with 15.75GB of memory. We also measured the bandwidth β to transfer data using PyTorch from the GPU to the RAM, and obtained around 12.2GB/s.

We use all available depths for ResNet [START_REF] Garg | Dynamic memory management for gpu-based training of deep neural networks[END_REF][START_REF]Periodic checkpointing in pytorch[END_REF]50,101,152) and DenseNet (121, 161, 169 and 201). We use three different image sizes: small images of shape 224 × 224 (which is the default and minimal image size for all models of torchvision), medium images of shape 500 × 500, and large images of shape 1000 × 1000. During the training phase, for higher efficiency, it is classical to process images in batches, where several images are processed independently. For each model and image size, we consider different batch sizes that are powers of 2, starting from the smallest batch size that ensures a reasonable throughput. For each case, we compute schedules with three different algorithms: Greedy (Section 4), DynProg (based on Algorithm 1, see Section 5.3), and Vdnn. AutoSwap [START_REF] Zhang | Efficient memory management for gpu-based deep learning systems[END_REF] is a score-based heuristic uses a weighted average of 4 priority scores to decide which activations should be offloaded in priority. The best weight combination is obtained with Bayesian Optimization. TFLMS [START_REF] Le | Tflms: Large model support in tensorflow by graph rewriting[END_REF] is a heuristic designed for general graphs (not necessarily sequential) in high bandwidth settings, but it does not use any profiling information and thus cannot adapt to the available memory. TFLMS is parameterized with the number of tensors to be offloaded and how many layers in advance the data should be prefetched. We present results with the parameters which, for a given memory size, achieve the lowest makespan among all possible values.

Our Vdnn is a heuristic based on ideas from VDNN++ [START_REF] Garg | Dynamic memory management for gpu-based training of deep neural networks[END_REF]. The original VDNN++ offloads only the input of convolutions, because "CONV layers have a much longer computation latency, being more likely to effectively hide the latency of offload/prefetch", and explores two possibilities: either offload the input of every convolution, or of every other convolution. However, in our setting the layers are not annotated to know which ones corresponds to convolutions. Hence, in Vdnn, we first compute the ratio u F i |xi| for all operations, and then for all possible thresholds, we compute the no-wait, ordered schedule which offloads all the activations whose ratio is above the threshold, and the one which offloads half of them. Vdnn outputs the best schedule out of all these choices.

Note that these experiments cover real scenarios obtained from measurements performed on the plaform described above; however the results presented here are obtained by simulation, in which we compute the schedules for the different heuristics, and estimate their expected duration (and thus the corresponding throughput) based on the measurements. In particular, this allows to consider cases where the memory limit is higher than the available memory on our GPU card.

Representative Results

A representative selection of achieved results is depicted in Figure 5. For each network, we run all algorithms with a memory limit varying from the minimum amount of memory required to run the network, to M peak which allows to process the network with no offloading. In each case, we also compute the lower bound LB (Proposition 3), and the plots show the ratio of the makespan obtained by each algorithm to the lower bound. We observe that both Greedy and DynProg outperform the Vdnn heuristic in all cases, especially in low memory scenarios. Once correctly parameterized, TFLMS is able to obtain optimal makespan for the highest memory limit values. But it is unable to delay forward computations until enough memory is made available through offloading, and thus can not adapt to low memory settings when bandwidth is scarce. AutoSwap often produces the same solution as the Greedy algorithm (for a much higher computational cost), but its performance depends on the random procedure of Bayesian Optimization and is thus very inconsistent. The DynProg approach obtains significantly better performance than Greedy. The difference is small in many cases, except for the DenseNet networks where DynProg is able to consistently obtain almost optimal solutions. The spike that can be observed on these graphs for Greedy and Vdnn correspond to the memory limit M GPU for which both terms of the lower bound LB are almost equal (i.e., the total execution time is very close to the time to transfer M peak -M GPU ). Such cases are significantly more difficult to solve because both criteria need to be optimized carefully.

Overall, DynProg obtains much more stable performance than Vdnn and AutoSwap, and produces solutions over a much wider range than TFLMS. Furthermore, DynProg is able to consistently achieve a ratio below 1.2, which means that its throughput is at least 83% of the highest possible throughput.

We next present the complete set of results, with two different presentations: we first show relative makespan like in the previous plot, which helps in comparing the performance of heuristics. Then we show the throughput (number of images processed by second) obtained for each case, which is the metric of interest for training DNNs.

Complete Results -Ratios

To allow for a better view of the results, the memory limits on the next plots are scaled: 0% corresponds to the amount of memory required to process the chain, and 100% corresponds to the memory peak M peak . The results are shown on Figures 6-12.

Complete results -Througphut

On the next plots, we show raw results, without normalizing with the lower bound. The y axis shows the throughput obtained, which is the number of images processed by seconds: T = batch size makespan . The results are shown on Figures 13141516171819.

Conclusions

In this paper, we address the problem of memory usage during the training phase of Deep Neural Networks. Previous works [START_REF] Rhu | vdnn: Virtualized deep neural networks for scalable, memory-efficient neural network design[END_REF][START_REF] Garg | Dynamic memory management for gpu-based training of deep neural networks[END_REF] have advocated to offload some of the data onto the memory of the CPU, and to prefetch them back when needed. We propose a formal algorithmic model of the corresponding scheduling problem, where the goal is to identify which activations should be offloaded so as to minimize the total execution time. We prove that this problem is NP-hard in the strong sense, and we propose two heuristics based on simpler relaxations of the problem. The Greedy heuristic always offloads the first activations in the network. This very simple technique nevertheless achieves good results in our experimental evaluation. The DynProg algorithm is a more sophisticated approach which takes into account the fact that activations cannot be partially transferred. It allows to obtain better solutions in favorable cases, but achieves worse performance with deeper neural networks, because of the discretization scheme. In any case, both algorithms provide significant improvements over the previous approach. This work opens several promising research directions. A validation with real experiments would allow to confirm the relevance of the assumptions made in the model. Since our theoretical analysis shows that being able to offload activations partially makes the problem much easier, it could be very interesting to assess in which cases this could be technically feasible. Finally, this offloading technique seems to be complementary of the checkpointing approach: some activations can be transferred to the CPU while others can be recomputed. Solving the corresponding algorithmic problem might be challenging, but would certainly yield a significant improvement for training large and deep models.

Figure 1 :

 1 Figure 1: Data dependencies induced the training phase of Sequential Deep Neural Networks.

Figure 3 :

 3 Figure 3: The instance of Problem 2 which corresponds to the 2-partition problem with values u i .

Theorem 3 .

 3 Problem 2 is NP-complete in the weak sense.Proof. Problem 2 clearly belongs to NP, in the same way as Problem 1.

i

  j=0 |x j |, in a similar way as for u Fj β values.

Figure 5 :

 5 Figure 5: Experimental results for image size 224 and batch size 32.

Figure 6 :Figure 7 :Figure 8 :Figure 9 :Figure 10 :Figure 11 :

 67891011 Figure 6: Relative Makespan results for ResNet with small images (224x224).

Figure 12 :Figure 13 :

 1213 Figure 12: Relative Makespan results for Inception.

Figure 14 :

 14 Figure 14: Throughput results for ResNet with medium images (500x500).

Figure 15 :

 15 Figure 15: Throughput results for ResNet with large images (1000x1000).

Figure 16 :

 16 Figure 16: Throughput results for DenseNet with small images (224x224).

Figure 17 :

 17 Figure 17: Throughput results for DenseNet with medium images (500x500).

Figure 18 :

 18 Figure 18: Throughput results for DenseNet with large images (1000x1000).

Figure 19 :

 19 Figure 19: Throughput results for Inception.

  for all other values of m, d F , d B . In order to compute Idle(i, m, d F , d B ) for all i and all relevant values of m, d F , d B , we use hash tables Idle i indexed with (m, d F , d B ), with the understanding that if (m, d F , d B ) is not stored in Idle i , then Idle(i, m, d F , d B ) = ∞. This leads to Algorithm 1, where Idle i values are used to update Idle i+1 values, with two possible cases, either with a schedule that offloads x i , or with a schedule that does not.


M A Sokolovskiy

J Verron

DYNAMICS OF THREE VORTICES IN A TWO-LAYER ROTATING FLUID

Keywords: Mathematics Subject Classification 37J35, 76B47, 76B70 1

, [20] a generalization is given for two-layer rotating fluid in the particular case determined by the conditions of (i) zero total circulation of vortices, and (ii) the equality of the intensities of two vortices. Here, the first of these restrictions is lifted.

The governing equations

We suppose that the vortex motion takes place in an unbounded two-dimensional two-layer medium consisting of liquids with constant densities ρ 1 , ρ 2 (ρ 1 ρ 2 ) in the upper and lower layer, respectively. Suppose that the medium as a whole rotates around an axis perpendicular to the plane coinciding with the 'rigid' (by assumption) top boundary of the upper layer. Again, let one of the vortices be situated in the upper layer and two vortices in the lower layer.

, [20], [40]) u

(1.2)

Equations (1.1) and (1.2) are written in a dimensionless form with the assumption that the layers are of equal thickness (h 1 = h 2 = 1/2). Hereinafter, the superscript is the vortex number, and the subscript is the layer number (the layers are numbered from top to bottom); the numbers of vortices are denoted by Greek letters, and those of layers by Roman letters; (x α j , y α j ), (u α j , v α j ) and κ α j are (x, y)-coordinates, the components of velocity, and the intensity of the vortex α j , 1 respectively; r αβ ij = =

x α i -x β j 2 + y α i -y β j 2 ; γ = 4Ω/ g(ρ 2 -ρ 1 )/ρ * is parameter of stratification; g is acceleration of gravity, ρ * is mean density; K j (R) is the modified Bessel function of the j-th order of argument R. 2 System of equations (1.1)-(1.2) has integral invariants:

Q = 1 2 2 j=1
A j α=1 κ α j -total intensity;

P x = 1 2 2 j=1 A j α=1
x α j κ α j , P y = 1 2 ≡ J q α j , H , ṗα j = -∂H ∂q α j ≡ J p α j , H , j = 1, 2, α = 1, A j (A 1 = 1, A 2 = 2), (1.4) where J is the Jacobian operator, q α j = x α j are generalized coordinates, p α j = y α j κ α j /2 are generalized momenta, and

H = -1 4π 2 j=1 A j α,β=1 α =β κ α j κ β j ln r αβ jj -K 0 γr αβ jj + + A j , A 3-j α, β=1
κ α j κ β 3-j ln r αβ j(3-j) + K 0 γr αβ j(3-j) (1.5) is the Hamiltonian coinciding with the vortex interaction energy. It is easy to verify that invariants M , H and the combination (P x ) 2 +(P y ) 2 are pairwise involutory and hence, as was the case with a homogenous medium [START_REF] Borisov | Poisson Structures and Lie Algebras in Hamiltonian Mechanics[END_REF], [START_REF] Kozlov | Dynamical Systems X. General Theory of Vortices[END_REF], [1], [START_REF] Aref | Integrable, Chaotic and Turbulent Vortex Motion in Two-Dimensional Flows[END_REF], [START_REF] Aref | Three-Vortex Motion with Zero Total Circulation: Addendum[END_REF], [START_REF] Aref | Vortices, Kinematics and Chaos[END_REF], [START_REF] Borisov | Dynamics and Statics of Vortices on a Plane and a Sphere[END_REF], [START_REF] Borisov | Dynamics of Three Vortices on a Plane and a Sphere. -II. General Compact Case[END_REF], [START_REF] Borisov | Dynamics of Three Vortices on a Plane and a Sphere. -III. Noncompact Case. Problem of Collapse and Scattering[END_REF], [START_REF] Borisov | Lie Algebras in Vortex Dynamics and Celestial Mechanics. -IV[END_REF], the problem of three vortices in a two-layer fluid always has a regular solution.

1 Such notation will be used for a vortex with a number of α, concentrated in layer j (in this case, α = 1 for j = 1 and α = 1, 2 for j = 2). 2 Since the argument of the Bessel function depends not only on the distance between the vortices but also on the multiplicative constant γ, hereafter we will use denotations of the form (X, Y ) ≡ γ(x, y), (X α j , Y α j ) ≡ γ(x α j , y α j ).

The specific features of the absolute motion of vortices will be studied below with the help of numerical solution of equations(1.1)-(1.2), supplemented by the appropriate initial conditions. The general properties of the relative motion are conveniently studied based on the qualitative analysis with the use of the so-called trilinear coordinates introduced in [1] and tested in [START_REF] Aref | Three-Vortex Motion with Zero Total Circulation: Addendum[END_REF], [6], [START_REF] Rott | Three-Vortex Motion With Zero Total Circulation[END_REF], [START_REF] Rott | Constrained Three-and Four-Vortex Problems[END_REF] and [START_REF] Sokolovskiy | New Stationary Solutions to the Problem of Three Vortices in a Two-Layer Fluid[END_REF], [START_REF] Gryanik | Dynamics of Baroclinic Vortices with Zero Total Intensity (Hetons)[END_REF], [START_REF] Sokolovskiy | Four-Vortex Motion in the Two Layer Approximation: Integrable Case[END_REF], [START_REF] Sokolovskiy | Dynamics of the Triangular Two-Layer Vortex Structures with Zero Total Intensity[END_REF], where they were applied to particular problems of three and four vortices in a homogeneous and a two-layer fluid, respectively.

Studying the relative and absolute motion of vortices

Phase portraits in trilinear coordinates

Following [1], let us consider a combination of integral invariants (1.3) in the form

L = QM -P x 2 -P y 2 . (2.1)
It is easy to show that L can be expressed in terms of squares of distances between the vortices, which allows this variable to be used for the analysis of the relative motion. The trilinear coordinates (t 1 , t 2 , t 3 ), which at L = 0 satisfy the identity

t 1 + t 2 + t 3 = 3, (2.2)
following from (2.1), in our case are introduced by the relationships

t 1 = 3κ 1 2 κ 2 2 r 12 22 2 L , t 2 = 3κ 1 1 κ 2 2 r 12 12 2 L , t 3 = 3κ 1 1 κ 1 2 r 11 12 2 L (2.3)
and characterize the distances from the coordinate lines t j = 0, (j = 1, 2, 3), constructed on the sides of an equilateral triangle with an altitude of 3, to any point on the plane (see figure 1a). When L = 0, we have

t 1 + t 2 + t 3 = 0 (2.4)
and

t 1 = 3κ 1 2 κ 2 2 r 12 22 2 , t 2 = 3κ 1 1 κ 2 2 r 12 12 2 , t 3 = 3κ 1 1 κ 1 2 r 11 12 2 .
(2.5)

A geometric interpretation of the trilinear coordinates in this case is given in figure 1b.

In the case of a definite problem (when the initial coordinates of vortices and their intensities κ α j are specified) on the plane of variables (t 1 , t 2 , t 3 ), it is necessary to recognize the 'physical region' (PR), within which the triangle inequality holds for the distances between arbitrary three points. In terms of trilinear coordinates, this inequality takes the form

(κ 1 1 t 1 ) 2 + (κ 1 2 t 2 ) 2 + (κ 2 2 t 3 ) 2 2(κ 1 1 κ 1 2 t 1 t 2 + κ 1 1 κ 2 2 t 1 t 3 + κ 1 2 κ 2 2 t 2 t 3 ). (2.6)
Clearly, the points of the boundary of PR where (2.6) becomes an equality are always associated with a continuum of collinear states of the three vortices. Eliminating t 1 from (2.2) and (2.4), we obtain, instead of (2.6), the following, more compact expressions

κ 1 1 + κ 1 2 t 2 + κ 1 1 + κ 2 2 t 3 -3κ 1 1 2 4κ 1 2 κ 2 2 t 2 t 3 (2.7) and κ 1 1 + κ 1 2 t 2 + κ 1 1 + κ 2 2 t 3 2 4κ 1 2 κ 2 2 t 2 t 3 (2.8)
for the cases L = 0 and L = 0, respectively. Isolines of the Hamiltonian (1.5), expressed in terms of trilinear coordinates are phase portraits of the relative motion of the vortex structure within the PR.

In the case of κ 1 2 = κ 2 2 = κ; κ in figure 4, by relationships

γ 2 L = -5, µ = -2.5 -(a); γ 2 L = 2, µ = 1 -(b); γ 2 L = 2, µ = 2 -(c
t 2 t 3 9 4 , t 1 < 0; t 2 , t 3 0; γ 2 L = -2 t 1 6; t 2 , t 2 < 0; γ 2 L = 2 at µ = -1; (2.14)
and in figures 5a and 5b (where L = 0):

t 2 + t 3 -9 2 16t 2 t 3 , t 1 0; t 2 , t 3 < 0 at µ = -1.5; t 2 + t 3 + 3 2 16t 2 t 3 , t 1 0; t 2 , t 3 < 0 at µ = -0.5. (2.15)
Note on choreographies. Of great importance in the celestial mechanics and the theory of vortices are closed trajectories, which correspond to periodic solutions of the problem. Running a few steps forward, we note that, in the interpretation of the relative motion of vortices with the use of the trilinear coordinates, periodic motions are clearly distinguishable, because the respective phase curves begin and end on the PR boundary. The respective trajectories of vortices (or material bodies) are referred to as choreographies [START_REF] Simó | New families of solutions in N -body problem[END_REF], [START_REF] Simó | Dynamical Properties of the Figure Eight Solution of the Three-Body Problem[END_REF], [START_REF] Chenciner | Simple Choreographic Motions of N Bodies: a Preliminary Study[END_REF], [START_REF] Borisov | Mathematical Methods of Vortex Structure Dynamics. In: Fundamental and Applied Problems in the Theory of Vortices[END_REF].

Although the solution of the three-body problem in the form of a closed circular trajectory was first obtained by Lagrange as far back as 1772 [START_REF] Chenciner | Simple Choreographic Motions of N Bodies: a Preliminary Study[END_REF], and later, starting from works of Hénon [START_REF] Hénon | Family of Periodic Solutions of the Planar Three-Body Problem and Their Stability[END_REF] and Moore [START_REF] Moore | Braids in Classical Dynamics[END_REF], the class of available solutions (for problems of celestial mechanics with a uniform Newtonian potential) with closed trajectories had extended significantly, the notion of choreography was coined by Simó [START_REF] Simó | New families of solutions in N -body problem[END_REF] quite recently (2002).

Choreographies can be • absolute when the trajectories are closed in a fixed coordinate system;

• relative when the trajectories become closed in a coordinate system rotating uniformly around a centre of rotation (or moving with a constant translational velocity, if the centre of rotation lies at infinity); • simple when all the bodies (or vortices) move along the same trajectory; and • complex when at least one body moves along a trajectory of its own.

First attempts to apply the notion of choreography to vortex dynamics in a homogenous fluid on a plane and a sphere were made by Borisov, Mamaev, and Kilin in [START_REF] Borisov | Mathematical Methods of Vortex Structure Dynamics. In: Fundamental and Applied Problems in the Theory of Vortices[END_REF], [START_REF] Borisov | Absolute and Relative Choreographies in the Problem of Point Vortices Moving on a Plane[END_REF].

Taking into account the two-layer character of the vortex structures considered in this study, the problem of distinguishing choreographies becomes even more important. As will be shown below, there can be situations in twolayer media when the vortices belonging to different layers have common trajectories. We will not study choreographies systematically but will repeatedly mention the cases of existence of this special class of vortex motions.

The patterns of phase portraits (figures 2-5) demonstrate a wide diversity of their topological properties. A detailed discussion of some of them given below in connection with specific problems will be accompanied by demonstrations of absolute trajectories of vortices and their relative choreographies. However, first we will obtain the conditions of existence of singular points on the phase plane, which correspond to stationary solutions of the equations of motion (1.1)-(1.2).

Analysis of stationary states

Judging from figures 2-5, we can conclude that two types of singular points -stable elliptic and unstable hyperbolic -can appear in the problems of the class considered here. The former always belong to the boundary of PR, and the latter coincide with the self-intersection points of separatrices. Both elliptic and hyperbolic singularities are obviously associated with rotational motions of rigid configurations of three vortices -collinear in the former case, and isosceles-triangular in the latter case. These rotations take place around the vorticity center of the vortex structure, which has the coordinates The boundary elliptic singular points are also divided into polar -in the upper and (or) lower parts of the PR for finite (infinite) PR -and side points. If for a polar point we have t 1 = 0, it corresponds to a structure comprising an upper-layer vortex and theoretically merged lower-layer vortices; otherwise, if t 1 = 0, this point corresponds to an ordinary or inverse roundabout. Of particular interest is the existence of side elliptic points, which are phase prototypes of eccentric roundabouts3 . We have

(X c , Y c ) = γP x Q , γP y Q . (2.16)
(X c , Y c ) =        2(R -B) µ + 2
, 0 for side elliptical singular points; 0, 2B sin ϕ µ + 2 for hyperbolic singular points.

(2.17)

The meaning of variables B, R and ϕ, where ϕ ∈ [0; π/2), is clarified by figures 6a, 7a, 9a and 10a, where, without loss of generality, it is assumed that, in the initial moment, the upper layer vortex 1 1 localizes in the origin of coordinates (X, Y ) = (0, 0), and (i) the collinear configuration belongs to the X-axis, and

(ii) the triangular configuration is symmetrical about the Y axis.

For collinear configurations, we obviously have

X c 0 at µ < -2, B R or µ > -2, B R, X c 0 at µ > -2, B R or µ < -2, B R, (2.18)
and for triangular configurations,

Y c 0 at µ > -2, Y c 0 at µ < -2. (2.19)
At µ → -2 (when the total intensity equals zero) the vorticity center tends to an infinite point 4 .

Let us consider some interesting particular combinations of coordinates (2.17).

• For the collinear structure, the vorticity center belongs to the center of the segment between vortices 1 2 and 2 2 at µ = 0; -coincides with vortex 2 2 at µ = 0 and

B = 2R(1 + µ) µ ; (2.

20)

an axisymmetric structure for which the center of rotation does not coincide with the central vortex. In this work, the latter definition is generalized to the asymmetric case. 4 This case is studied in [START_REF] Sokolovskiy | Four-Vortex Motion in the Two Layer Approximation: Integrable Case[END_REF].

coincides with the center of the segment between vortices 1 1 and 2 2 at µ = 2 and

B = 2Rµ µ -2 .
(2.21)

• For the triangular structure, the vorticity center lies within the triangle formed by the vortices at µ > 0;

lies outside this triangle at µ < 0, and, as mentioned above, at µ = -2 -it moves to infinity;

coincides with the center of the circle circumscribed around this triangle at µ = -2 cos 2ϕ; the center of this circle, in its turn, lies within the triangle at ϕ > π/4 (when 0 < µ < 2) and outside the triangle at ϕ < π/4 (when -2 < µ < 0).

We obtain dispersion equations describing the conditions of existence of steady-state solutions either of collinear type in the form of an eccentric roundabout or in the form of an isosceles triangle.

Eccentric roundabout

Using equations (1.1)-(1.2), we readily obtain the dispersion equation for a steady state configuration in the form of a uniformly rotating collinear structure:

F (B, R; µ) = 1 2 + 2R(1 + µ) B(2R -B) + K 1 (2R)- (2R + Bµ)K 1 (2R -B) + [2R(1 + µ) -Bµ]K 1 (B) 2(R -B) = 0. (2.22)
Clearly, F is a monotonically decreasing function of both B and R, and

F ∼    1/ε 1 > 0, ε 1 = min(B, R) at B, R ≪ 1; (1 + µ) 1/ε 2 -K 1 (ε 2 ) , ε 2 = B ≪ 1 at R ≫ 1; (1 + µ) 1/ε 3 -K 1 (ε 3 ) , ε 3 = 2R -B ≪ 1 at B → 2R, R ≫ 1.
(2.23)

Expressions in square brackets in (2.23) are always positive; therefore, the necessary condition for the existence of solution of equation (2.22) can be written as the inequality µ < -1.

(2.24)

Thus, the coordinates of the center of rotation for stationary states in the form of an eccentric roundabout correspond to a narrower interval of parameter µ than that specified by (2.18) for arbitrary collinear configurations. Namely:

X c 0 at µ < -2, B R or -1 > µ > -2, B R, X c 0 at -1 > µ > -2, B R or µ < -2, B R. (2.25)
The angular velocity of the rectilinear vortical structure with respect to the vorticity center (2.25) is

ω = γκ(µ + 2) 4π(2R + Bµ) B + 2Rµ 2BR -µK 1 (B) + K 1 (2R) . (2.26)
The case µ = -2 (i. e. Q = 0), when the angular velocity becomes zero, and the rotation center (2.17) moves to an infinite point, was studied in [START_REF] Sokolovskiy | Four-Vortex Motion in the Two Layer Approximation: Integrable Case[END_REF].

Figures 6 and7 demonstrate the character of solution to equation (2.22) at µ = -2.5 and µ = = -1.5, respectively, within a wide range of values of integral L. In particular, they show the abovementioned property that the sign of L is constant at µ -2 and can change at 0 > µ > -2. Note that curve X c (R) is symmetrical with respect to the R axis, and B(R) is symmetrical with respect to straight line B = R, and with growing R, the value of B asymptotically approaches 0 or 2R at R > B and R < B, respectively -see (2.23). To illustrate these effects, auxiliary half-rays B = R and B = 2R are shown in the figures by dotted lines. Obviously, at R → ∞, the angular velocity tends to zero.

The motion of vortices in the cases when both the dispersion equation (2.22) and the relationships (2.20) or (2.21) are valid is illustrated in figure 8. Here and in all subsequent figures, where vortex trajectories (or choreographies) are shown, triangular markers denote the instantaneous positions of the upper-layer vortex, and circle and box markers denote the positions of the lower-layer vortices. The rotation in figure 8a is cyclonic, and the center of vorticity coincides with the vortex 2 2 , while the rotation in figure 8b is anticyclonic and the trajectories of vortices 1 1 and 1 2 coincide (the calculation, however, was terminated before they begin to overlap). Clearly, the collinearity condition for configuration is valid everywhere in this case, but the segments passing through three vortices are shown in the figure only for two time moments for each experiment, i.e., the initial (when the vortices lie on the X-axis) and the final (target).

Triangular stationary structure

Similarly, we can obtain relationships between the geometrical parameters of a stationary unstable vortex structure in the form of a uniformly rotating isosceles triangle (figures 9 and 10), which corresponds to a fixed hyperbolic singular point of the phase portrait.

The dispersion equation relating the side length of triangle B and the angle at its base ϕ has the form 1 + 2B cos ϕK 1 (2B cos ϕ)

cos 2 ϕ = -4(1 + µ) 1 -BK 1 (B) . ( 2

.27)

Since the expression in the left-hand part of (2.27) is always greater than unity and 0 < 1-BK 1 (B) < 1, the necessary condition for the existence of solution to the dispersion equation takes the form µ < -5/4.

(2.28)

It is worth mentioning here that the comparison of conditions (2.24) and (2.28) yields the explanation of the fact that at µ = -1 (figure 4), there are no side elliptic singular points, while there is a hyperbolic singular point. The asymptotic properties of modified Bessel functions [START_REF] Prudnikov | Integrals and Series. V. 2. Special Functions[END_REF] determine the restrictions on the angle at the base of the triangle at B ≫ 1:

arccos ϕ < - 1 4(1 + µ) . (2.29)
clearly, the values of parameter µ should lie on the half-line (2.28) on which the expression under the radical in (2.29) is always positive. Formula for the angular velocity of the triangular structure relative to the vorticity center (0, Y c ) takes the form

ω = γκ(2 + µ) 4πB 2 1 -BK 1 (B) . (2.30)
If κ > 0, it becomes ω 0 at µ -2, and, in view of (2.28), ω 0 at -5/4 > µ -2. Analysis of the phase portraits in figures 2-5 shows that, at least in this (relatively wide) interval of external parameters, no more than three qualitatively different types of motion are possible:

• Type {1} -double trapping, when all three vortices are involved in joint rotational motion in the direction induced either by their resultant intensity or by the intensity of the 'strongest' vortex. This motion can be conventionally determined by the scheme:

1 1 1 2 2 2 .
(2.33)

• Type {2} -simple trapping, when one of the lower-layer vortices combines with the upper-layer vortex and one of the following formulas is valid:

1 1 1 2 + 2 2 or 1 1 2 2 + 1 2 .
(2.34)

• Type {3} -the interaction between the two lower-layer vortices predominates, and the motion is determined by the scheme 1 2

2 2 + 1 1 . (2.35)

Absolute motion of vortices. Choreographies

Analysis of the relative motion of vortex structures is known to be unable to yield a complete picture of their actual dynamics. In this paragraph, in addition to the above-considered specific cases of vortex motion trajectories under stationary conditions, we will discuss in greater detail the features of the absolute motion of triangular vortex structures in combination with the qualitative analysis based on some of figures 2-5. Characteristic examples of trajectories of all three types of absolute motion and the respective relative choreographies at µ = -2.5 are given in figure 11, 5 where, obviously, the choreographies are plotted in the system of coordinates uniformly rotating around the vorticity center in the anticyclonic direction. It is worth mentioning that, on the relative choreographies of motions of types {1} and {3}, the lower-layer vortices 1 2 and 2 2 move along the same trajectory, always remaining on the opposite ends of its 'diameter'. However, in the first case, they (in a moving coordinate system) move along closed trajectories in the clockwise direction induced by the upper-layer vortex; in the second case, because of the predominance of intralayer interaction, they move in the opposite direction (formulas (2.33), (2.35)). For motions of type {1}, vortex 1 1 stays within the closed figure described by the lower-layer vortices, and for motions of type {3}, outside this figure.

Figures 12a-13a illustrate in more detail the interaction of type {2} for the cases when the image point of the phase plane coincides with the side elliptic point and we have a stable rotation of a rigid collinear structure, i.e., an eccentric roundabout. Figures 12b-13b demonstrate the unstable character of the triangular construction (the respective phase trajectory lies in the vicinity of the separatrix): the vortices alternately either occupy the vertices of the isosceles triangle or lie along a straight line, as the image point moving along the phase curve lies either in a vicinity of the hyperbolic singularity or on the boundary of PR, respectively. It is worth mentioning that the relative choreographies in figure 12 have been obtained in the coordinate system rotating around the vorticity center with a negative angular velocity as µ < -2 (see the formula for angular velocity (2.26)), and those in figure 13, where µ > -2, with a positive velocity 6 . Let us consider the motion at µ = 1, when the intensities of all three vortices are equal. The behavior of analogous vortex structures in a homogeneous fluid was exhaustively studied in [START_REF] Meleshko | Dynamics of Vortex Structures[END_REF], [START_REF] Borisov | Mathematical Methods of Vortex Structure Dynamics. In: Fundamental and Applied Problems in the Theory of Vortices[END_REF], [START_REF] Tavantzis | The Dynamics of Three Vortices Revisited[END_REF].

Isolines of the Hamiltonian in trilinear coordinates for this case are shown in figure 2b. Note that in this simplest case, the motions of vortices are also relatively simple and all belong to a single class {3}, where the lower-layer vortices are involved in rotation around one another, while the vortex 1 1 can lie either within or outside the construction formed by them. The phase portrait shows all motions to be periodic. Moreover, twice within each period, every evolving triangular structure degenerates into a collinear configuration and once takes the form of an isosceles triangle. A particular case of the collinear state can be the one when the upper-layer vortex lies precisely above one of the lower-layer vortices in the tangency points of the PR boundary with the lateral sides of the coordinate triangle: t = (1.5, 0, 1.5) and t = (1.5, 1.5, 0). The portrait contains two singular elliptic points. The upper point corresponds to a stable stationary state -a roundabout, when the vortices with subscripts 1 2 and 2 2 , located symmetrically about the immovable central vortex of the upper layer, uniformly rotate along a circumference. The lower point corresponds to a structure in which merged lower-layer vortices rotate around a common center with a theoretically infinite angular velocity.

Figures 14.1.a-14.1.d illustrate the set of possible motions of such system in a sequence when passing from one experiment to another corresponds to the motion of the image point in the phase space along the boundary of PR from the upper polar singularity toward the lower polar singularity.

Figure 14.1.a shows a realization of a stationary roundabout. This is an example of absolute choreography. Obviously, the relative choreography (14.2.a) in the system of coordinates rotating with an angular velocity of peripheral vortices, is determined by three fixed points. shows the initial stage of the evolution of a vortex structure in the case when the vortices 1 1 and 2 2 coincide in the initial moment. The behavior of the system is such that the vortex periodically and alternatively coincides with the upper-layer vortices. This can be clearly seen in both the figures of relative trajectories (14.1.b) and relative choreographies (14.2.b). Figures 14.1.c and 14.1.d give the trajectories of absolute motion for the situations where vortex 1 1 all the time stays out of the structure 1 2 -2 2 . In the former case, the upper-layer vortices periodically coincide with the vorticity center of the system, and in the latter case, they alternatively change places always remaining on the same side with respect to this center.

Clearly, as the initial point of the phase trajectories moves farther toward the lower pole, the trajectories of the lower-layer vortices will describes a figure more and more similar to a circle with a radius tending to zero. This is also confirmed by the fact that the phase trajectories in the lower part of the figure are practically horizontal lines.

It should be mentioned that most examples considered above (figures 8, 11, 12, 14) represent the cases of compact PR. Now let us consider characteristic features of the absolute motions allowing for infinite trajectories. Suppose that µ = -1. The respective phase portraits are given in figure 4. Figure 15 gives mostly the initial parts of trajectories of absolute motion of vortices at L < 0. Experiments 1 and 2, demonstrating a realization of regime {2}, start from collinear initial configu-rations. They show that mutual capture of vortices 1 1 and 1 2 takes place via the formation of a two-layer pair moving to infinity, whereas vortex 2 2 at t → ∞ takes some fixed position. In these figures, segments connect the positions of vortices in the initial and final moments of calculations.

Two more examples of infinite motions are given by experiments 6 (type {2}) and 7 (type {1}). Here, the initial structures are triangles (not isosceles), such that the distances r 21 22 (or trilinear coordinates t 1 ) for these two cases are identical. In this case, obviously, coordinates Y 2 2 are chosen based on the conservation condition of invariant L. The configurations marked out in the first fragment are: a) initial, b) in the form of an isosceles triangle (corresponds to the image point lying on the line t 2 = t 3 ), c) the configuration antisymmetric to the initial one (in which the sides of the triangle with the lengths of r 12 12 and r 11 12 change places), and d) final. The configurations in the second fragment are the same with the exception that in item b), we have a collinear configuration (corresponding to the position of the image point on the PR boundary), and in item c), the configuration is symmetrical. Clearly, the possible class of vortex structures of motions of type {2} includes the collinear state and does not include the configuration in the form of an isosceles triangle; the situation with motions of type {1} is inverse.

Experiments 3, 4 and 5 demonstrate trajectories of finite type {3} with a dominating character of intralayer interaction between lower-layer vortices 1 2 and 2 27 . Fragments 3-5 present trajectories in time intervals of T /2, 2T and 188T with marked out instantaneous configurations with intervals of T /8, T /2 and 12T , respectively (T is the period of relative motion). Thus, motions of type {3}, as before, are represented by rotations along loop-like trajectories of cyclonic lower-layer vortices around some variable center. After each quarter of period, the vortex structure forms, in turn, either a collinear configuration or a configuration in the form of an isosceles triangle with successive rearrangement of vortices 1 2 and 2 2 . Figure 16 shows trajectories of absolute finite character (left) and the respective choreographies (right). Lines (a) and (c) correspond to figures 15.3 and 15.4, respectively. In both cases, the calculation time is extended (as compared with the previous figure) and amounts to about 6T . The origin of coordinates coincides with the vorticity center of the vortex structure. During the construction of choreographies, the coordinate system was subjected to a rotational motion -counterclockwise in figure 16a and clockwise in figure 16c. It is clear that in this case, an intermediate state of the systems must exist when no rotation of the coordinate system is required and, thus, an absolute choreography is observed. This is the situation shown in figure 16b (the parameter values are given in the capture). The markers and the segments connecting them on the choreographies show the collinear and triangular positions of vortices with a time step of a quarter of the period.

Thus, in this case, as in the four-vortex problem with P = 0 and M = 0 [START_REF] Gryanik | Dynamics of Baroclinic Vortices with Zero Total Intensity (Hetons)[END_REF], one of the phase trajectories in the class of finite motions {3} is associated with an absolute choreography. It should be noted that this periodic solution is stable: the system of vortices in figure 16b has accomplished more than 20 cycles.

Conclusions

Possible motions of a system of three vortices were studied under the assumption that one vortex lies in the upper layer of a two-layer fluid and the other two vortices, which have equal intensities, lie in the lower layer. The associated dynamic system was studied in a wide range of external parameters. The relative motions were classified, in particular, dispersion relationships, determining the characteristics of stationary states, were obtained. A new stable stationary solutioneccentric roundabout was obtained; it generalizes the solution of triton type [START_REF] Sokolovskiy | New Stationary Solutions to the Problem of Three Vortices in a Two-Layer Fluid[END_REF], [START_REF] Sokolovskiy | Dynamics of the Triangular Two-Layer Vortex Structures with Zero Total Intensity[END_REF] to the case of nonzero total intensity. The existence of relative and absolute choreographies for a system of three vortices was established. Numerically calculated absolute trajectories of vortices are exemplified. Clearly, these examples fail to reflect the entire spectrum of possible motions for the three-vortex system considered. However, the examples give an idea of their character. Additionally, the analysis of the presented phase portraits extends our knowledge about the properties of this dynamic system. In particular, the upper part of figure 4 (where γ 2 L = 2 and µ = -1) shows that the respective symmetrical roundabout is unstable (it corresponds to the bottom point of PR). The problem of stability of axisymmetric structures in the presence of a central vortex (in the manner of works [START_REF] Aref | Point Vortex Motions with a Center of Symmetry[END_REF], [START_REF] Borisov | Mathematical Methods of Vortex Structure Dynamics. In: Fundamental and Applied Problems in the Theory of Vortices[END_REF], [START_REF] Gryanik | About Theoretical Models of Localized Quasi-Geostrophic Vortices in the Atmosphere and Ocean[END_REF], [START_REF] Kurakin | Stability of Stationary Rotation of a Regular Polygon[END_REF], [START_REF] Bauer | Stability of Rectilinear Geostrophic Vortices in Stationary Equilibrium[END_REF], [START_REF] Cabral | Stability of Relative Equilibria in the Problem on N + 1 Vortices[END_REF], [START_REF] Kurakin | The Stability of Stationary Rotation of a Regular Vortex Polygon[END_REF], [START_REF] Morikawa | Interacting Motion of Rectilinear Geostrophic Vortices[END_REF], but applied to a two-layer fluid medium) will be a subject of a separate paper.

The obtained results yield the explanation of some mechanisms of vortex interaction in a stratified rotating fluid and can be of use in the interpretation of specific features of real vortex structure trajectories in the ocean and atmosphere.

γr α1 21 K 1 γr α1 21 .

 2121

3 )

 3 System of equations (1.1)-(1.2) can be written in the Hamiltonian form qα j = ∂H ∂p α j

Fig. 1 .

 1 Fig. 1. Scheme of trilinear coordinates and axes orientation at L = 0 -(a) and L = 0 -(b).

Fig. 3 .

 3 Fig. 3. Phase portraits of three-vortex problems in PR (2.13) at µ = -1.5 and L < 0 -in the lower part, L > 0 -in the upper part: γ 2 |L| = 2 -(a); in the lower part γ 2 |L| = 0.65 -(b). Figures mark the points in PR that correspond to the initial conditions for the experiments the results of which are given in figure 13.

Fig. 4 .

 4 Fig. 4. Phase portraits of three-vortex problems in PR (2.14) at µ = -1 and γ 2 L = -2 -in the lower part, γ 2 L = 2 -in the upper part. Figures mark the points in PR that correspond to the initial conditions for the experiments the results of which are given in figure 15.

Fig. 5 .

 5 Fig. 5. Phase portraits of three-vortex problems in PR (2.15) at L = 0: µ = -1.5 -(a) and µ = -0.5 -(b).

Fig. 6 .

 6 Fig. 6. Scheme of the initial position of vortices for an eccentric roundabout at B > R (top) and at B < R (bottom) -(a); dispersion curve B(R), satisfying (2.22), the curve of coordinates of the vorticity center Xc(R), the angular velocity (2.26) of rotation of the vortex construction ω(R) and the value of the integral invariant L(B) (2.9) at µ = -2.5 -(b). The tilted dashed lines have equations B = R and B = 2R. The branches of curves B(R) and Xc(R) for the case B > R are shown by the full line, and those for B < R by the dashed line. The points of intersection of the curves with the vertical dashed line correspond to the parameters of numerical experiments at γ 2 L = -5, shown in figure 12a.

Fig. 7 .

 7 Fig. 7. The same as in figure 6 but for µ = -1.5. The points of intersection of the curves with the vertical dashed line correspond to the parameters of numerical experiments at γ 2 L = -2, shown in figure 13a.

Fig. 8 . 1 1✂Fig. 9 .

 819 Fig. 8. Trajectories of absolute motion of vortex structureseccentric roundabouts -in the case when the conditions (2.20) -(a) and (2.21) -(b) are satisfied. Parameter values: µ = -1.5757, B = 0.7162, R = 1.0526, L = 0.5746 -(a); µ = -2.1598, B = 0.8551, R = 0.8235, L = -0.2211 -(b). The full curves represent the upper-layer vortex trajectories ✁ 1 1✂ , and the dashed line, those of the lower-layer vortices

Fig. 10 .

 10 Fig. 10. The same as in figure 9, but at µ = -1.5. The vertical dashed line reflects the parameters of experiments in figure 13b.

Fig. 13 .

 13 Fig. 13. Trajectories of the absolute motion -1 and the respective choreographies -2 for motions of type {2} at γ 2 L = -2, µ = -1.5 and the initial conditions with X 2 2 = 0.1949, t = (-10.2416; 0.0855; 13.1561) -(a); X 2 2 = -0.2338, t = (-2.6901; 0.1230; 5.5671) -(b). Rectilinear segments in figure (a) show the collinear positions of vortices in the initial and final time moments in the calculation, and in figure (b), at intervals of a quarter of period (the configurations are, alternatively, collinear or in the form of an isosceles triangle). The trilinear coordinates for cases (a) and (b) correspond to points with numbers 1, 2 in figure 3a.

  Figure 14.1.b

Fig. 14 . 1 -

 141 Fig. 14. 1 -Absolute trajectories of three vortices in the case µ = 1 and L = 2. 2 -The respective relative choreographies in the coordinate system rotating in the cyclonic direction about the vorticity center of the system (its coordinates are marked by crosses). The initial conditions correspond to the following parameter values: (a) -X 2 2 = -1/ √ 3, t = (2.0, 0.5, 0.5) -point 1 in figure 2b; (b) -X 2 2 = 0, t = (1.5, 0.0, 1.5) -point 2 in figure 2b; (c) -X 2 2 = 1/ √ 3, t = (0.5, 0.5, 2.0) -point 3 in figure 2b; (d) -X 2 2 = 0.75, t = (0.2227, 0.8438, 1.9335) -point 4 in figure 2b.

Fig. 15 .

 15 Fig. 15. Trajectories of absolute motion at µ = -1, γ 2 L = -2 and the initial conditions corresponding to the enumerated points on the phase plane (figure 3b) with the following parameters: 1 -X 2 2 = 0.45, t = (-4.7112; 7.4074; 0.3038); 2 -X 2 2 = 0.4810, t = (-3.8304; 6.4824; 0.3770); 3 -X 2 2 = 0.4811, t = (-3.3750; 6.0000; 0.3750); 4 -X 2 2 = = 0.5, t = (-4.7112; 7.4074; 0.3038); 5 -X 2 2 = 0.9, t = (-0.0669; 1.8519; 1.2150); 6 -X 2 2 = 2.5, t = = (-15.0000; 8.6250; 9.3950); 7 -X 2 2 = 3, t = (-15.0000; 4.5000; 13.5000).

Fig. 16 .

 16 Fig. 16. Absolute trajectories of three vortices (left) and the respective choreographies (right) for motions of the type {3} at collinear initial state and the following parameter values: µ = -1, γ 2 L = -2 and (a) -X 2 2 = = 0.4811, t = (-3.3750; 6.0000; 0.3750); (b) -X 2 2 = 0.4833, t = (-3.7722; 6.4218; 0.3504); (c) -X 2 2 = 0.5, t = = (-4.7112; 7.4074; 0.3038).

The terms ordinary and inverse roundabouts are introduced in[START_REF] Sokolovskiy | Four-Vortex Motion in the Two Layer Approximation: Integrable Case[END_REF] to define rigid axisymmetric triangular structures for which the peripheral vortices of the lower layer rotate with respect to the central vortex in the upper layer in the direction induced by this vortex or in the direction determined by the intralayer interaction, respectively. The ordinary roundabout is a two-layer analogue of the so-called tripolar vortices or tripoles in a homogeneous fluid[START_REF] Beckers | The Observation of a Triangular Vortex in a Rotating Fluid[END_REF],[START_REF] Carnevale | Emergence and Evolution of Triangular Vortices[END_REF],[START_REF] Van Heijst | Tripolar Vortices in a Rotating Fluid[END_REF],[START_REF] Van Heijst | Laboratory Experiments on the Tripolar Vortex in a Rotating Fluid[END_REF],[START_REF] Kloosterziel | On Tripolar Vortices[END_REF],[START_REF] Polvani | The Tripole: A New Coherent Vortex Structure of Incompressible Two-Dimensional Flows[END_REF],[START_REF] Velasco Fuentes | Unsteady Behaviour of a Topography-Modulated Tripole[END_REF]; inverse roundabout has a purely baroclinic nature. The term eccentric roundabout is assigned in[START_REF] Sokolovskiy | Four-Vortex Motion in the Two Layer Approximation: Integrable Case[END_REF] to

Hereafter, in the captions to the figures where absolute trajectories are shown, we will use denotations of the form t = (t1, t2, t3) to specify the initial coordinates of the respective image points on the phase plane.

See also the behavior of curves ω(B) in figures 9 and 10.

Note that notwithstanding the fact that the initial conditions in experiments 2 and 3 are very similar (the markers with numbers 2 and 3 in figure4are practically indistinguishable), the respective motions belong to different types since the image points on the phase plane belong to trajectories lying in a vicinity of the separatrix but on different sides of it.

Acknowledgement

We are grateful to our colleagues V. M. Gryanik, A. V. Borisov and I. S. Mamaev for fruitful discussions of the results of this study. This work was financially supported by the Russian Foundation for Basic Research (project ✄ 04-05-64367).

Classification of motions of triangular vortical structures

In this paragraph, we will discuss the features of different types of motion of vortices based, on the one hand, on the classification provided by phase portraits in figures 2-5, and, on the other hand, on studying the behavior of trajectories of the absolute motion of vortices and their relative choreographies by numerically solving the equations of motion (1.1)-(1.2).

The choice of the initial conditions in the numerical experiments coordinated with the analysis of phase portraits will be based on the following considerations.

(i) If the phase curves of the class of motions considered reach the boundary of PR and, thus, the set of possible vortex configurations contains a collinear structure, we will suppose that all the three vortices belong to the X-axis, the vortex 1 1 localizes in the origin, and the coordinates of vortices 1 2 and 2 2 satisfy the conservation condition of the invariant L, i. e., at t = 0,

(2.31) (ii) If the phase trajectories do not reach the boundary of PR, we will assume that at t = 0, only vortices 1 1 and 1 2 belong to the X-axis, and the vortex 2 2 does not belong to it:

(2.32)

Bearing in mind (2.31)-(2.32), when specifying the parameters of the initial state in the captions to figures showing vortex trajectories, we will show only the values of the X 2 2 coordinate in the first case, and the pair of coordinates (X 2 2 , Y 2 2 ) in the second case.