François Pacaud

Michel Schanen

Sungho Shin

Daniel Adrian Maldonado

Mihai Anitescu

Parallel Interior-Point Solver for Block-Structured Nonlinear Programs on SIMD/GPU Architectures

We investigate how to port the standard interior-point method to new exascale architectures for block-structured nonlinear programs with state equations. Computationally, we decompose the interior-point algorithm into two successive operations: the evaluation of the derivatives and the solution of the associated Karush-Kuhn-Tucker (KKT) linear system. Our method accelerates both operations using two levels of parallelism. First, we distribute the computations on multiple processes using coarse parallelism. Second, each process uses a SIMD/GPU accelerator locally to accelerate the operations using fine-grained parallelism. The KKT system is reduced by eliminating the inequalities and the state variables from the corresponding equations, to a dense matrix encoding the sensitivities of the problem's degrees of freedom, drastically minimizing the memory exchange. We demonstrate the method's capability on the supercomputer Polaris, a testbed for the future exascale Aurora system. Each node is equipped with four GPUs, a setup amenable to our two-level approach. Our experiments on the stochastic optimal power flow problem show that the method can achieve a 50x speed-up compared to the state-of-the-art method.

Introduction

Solving complex engineering problems often resorts to the solution of large-scale blockstructured nonlinear programs. As such, there has been a long interest in designing efficient nonlinear optimization algorithms, particularly by using parallel computing. Parallelism can happen at two levels. At first, coarse parallelism splits the program into large computational chunks, usually dispatched to multiple processors using a message-passing interface in distributed memory. In this paradigm, the parallel algorithm is designed to minimize the communication between the different processes. In a complementary direction, fine-grained parallelism breaks down the program into small tasks, fast to compute in shared memory. This method requires a large number of processors to be efficient, and it is usually better on SIMD architectures with low communication overhead, as provided by Graphical Processing Units (GPUs). In the mathematical optimization community, coarse parallelism has traditionally been used to solve large-scale block-structured optimization problems, as encountered in dynamic or stochastic nonlinear programs. On the contrary, fine-grained parallelism has gained attraction only recently, with the renewed interests for machine learning applications and stochastic gradient algorithms. In this work, we combine coarse and fine-grained parallelism to solve block-structured nonlinear problems on new exascale architectures, where the solution algorithm is streamlined on different GPUs using CUDA-aware MPI.

Literature review

In his pioneering work [START_REF] Schnabel | Parallel computing in optimization[END_REF][START_REF] Schnabel | A view of the limitations, opportunities, and challenges in parallel nonlinear optimization[END_REF], Robert Schnabel identified three practical approaches to run optimization algorithms in parallel: (i) parallelize the function evaluations; (ii) parallelize the linear algebra; and (iii) parallelize the optimization algorithm itself.

The first attempt to parallelize the evaluations has been to streamline the computation of the derivatives using finite-differences [START_REF] Nash | A general-purpose parallel algorithm for unconstrained optimization[END_REF]. Soon, it has been noted that parallelizing the forward pass in automatic differentiation (AD) is also straightforward, provided that we can propagate the tangents (encoding the first-order sensitivity) in parallel [START_REF] Hovland | Automatic differentiation of parallel programs[END_REF]. Unfortunately, doing the same in the reverse pass is not trivial, as adjoining a mutable code leads to race conditions (e.g., every read becomes a write operation). This has led to extensive research on adapting automatic differentiation to parallel environments [START_REF] Bücker | Bringing together automatic differentiation and OpenMP[END_REF][START_REF] Hovland | Automatic differentiation for message-passing parallel programs[END_REF][START_REF] Moses | Reverse-mode automatic differentiation and optimization of GPU kernels via Enzyme[END_REF]. Now, most state-of-the-art differentiable tools employ a Domain Specific Language (DSL) constraining the user to specific differentiable operations. In particular, this approach has been adopted mainly in machine learning, leading to the development of fast AD libraries efficiently generating the derivatives efficiently on hardware accelerators such as GPUs or TPUs [START_REF] Bradbury | JAX: composable transformations of Python+NumPy programs[END_REF][START_REF] Paszke | Pytorch: An imperative style, high-performance deep learning library[END_REF].

The parallelization of linear algebra is usually more involved, as most large-scale optimization methods fall back on the solution of sparse indefinite Karush-Kuhn-Tucker (KKT) systems [START_REF] Nocedal | Numerical optimization[END_REF]. In the 1980s, preliminary results were obtained by running iterative methods in parallel, using block-Krylov [START_REF] Saad | On the rates of convergence of the Lanczos and the block-Lanczos methods[END_REF] or block-truncated Newton methods [START_REF] Nash | Block truncated-Newton methods for parallel optimization[END_REF]. However, block iterative algorithms are quickly limited by the lack of generic preconditioners for KKT systems. The 1990s witnessed the emergence of the interiorpoint methods (IPM), together with the development of large-scale sparse direct linear solvers [START_REF] Duff | MA57-a code for the solution of sparse symmetric definite and indefinite systems[END_REF][START_REF] Schenk | Solving unsymmetric sparse systems of linear equations with PARDISO[END_REF]. In IPM, a significant portion of the time is spent solving a sequence of (indefinite) KKT systems, hence the method directly benefits from efficient sparse linear solvers able to run in parallel [START_REF] Amestoy | Multifrontal parallel distributed symmetric and unsymmetric solvers[END_REF][START_REF] Duff | Developments and trends in the parallel solution of linear systems[END_REF]. In the 2000s, it was shown that, for blockstructured optimization problems as we consider here, the layout of the optimization problem can be exploited further in a Schur complement approach to solve the Newton step in parallel [START_REF] Birge | Computing block-angular Karmarkar projections with applications to stochastic programming[END_REF][START_REF] Choi | Exploiting special structure in a primal-dual pathfollowing algorithm[END_REF][START_REF] Gondzio | Parallel interior-point solver for structured linear programs[END_REF][START_REF] Jessup | Parallel factorization of structured matrices arising in stochastic programming[END_REF][START_REF] Petra | An augmented incomplete factorization approach for computing the Schur complement in stochastic optimization[END_REF][START_REF] Word | Efficient parallel solution of large-scale nonlinear dynamic optimization problems[END_REF][START_REF] Zavala | Interior-point decomposition approaches for parallel solution of large-scale nonlinear parameter estimation problems[END_REF]. These developments led to the development of mature decomposition-based parallel nonlinear solvers for scenario-based problems in the 2010s [START_REF] Chiang | Structured nonconvex optimization of large-scale energy systems using PIPS-NLP[END_REF][START_REF] Gondzio | Exploiting structure in parallel implementation of interior point methods for optimization[END_REF][START_REF] Rodriguez | Scalable parallel nonlinear optimization with PyNumero and Parapint[END_REF][START_REF] Shin | Graph-based modeling and decomposition of energy infrastructures[END_REF][START_REF] Zhu | Exploiting modern computing architectures for efficient large-scale nonlinear programming[END_REF].

Eventually, running an optimization algorithm fully in parallel generally requires a subtle combination of (i) and (ii), often devolving to a software engineering problem. The challenge is to evaluate the derivatives and solve the resulting KKT system each in parallel; all this while minimizing the communication between the different processes. This has led to the development of different prototypes for MPI-parallel modelers [START_REF] Colombo | A structureconveying modelling language for mathematical and stochastic programming[END_REF][START_REF] Huchette | Parallel algebraic modeling for stochastic optimization[END_REF][START_REF] Rodriguez | Scalable parallel nonlinear optimization with PyNumero and Parapint[END_REF][START_REF] Watson | PySP: modeling and solving stochastic programs in python[END_REF], most of them extending a specific AD backend [START_REF] Bussieck | General algebraic modeling system (GAMS)[END_REF][START_REF] Dunning | JuMP: A modeling language for mathematical optimization[END_REF][START_REF] Fourer | A modeling language for mathematical programming[END_REF]. Such approaches have been successfully applied to solve large-scale block-structured nonlinear problems, as encountered in stochastic programming and dynamic optimization.

Contributions

In this article, we introduce a new parallel algorithm to solve block-structured nonlinear programs involving state equations on exascale supercomputers. Our algorithm uses the parallel interior-point solver MadNLP [START_REF] Shin | Graph-based modeling and decomposition of energy infrastructures[END_REF], using two layers of parallelism to streamline both the evaluation of the derivatives and the solution of the KKT system. This framework targets new exascale supercomputers, where each node is assigned to multiple GPUs connected with a unified memory (designed to have fast memory exchange between the different GPUs).

We demonstrate the capability of the algorithm on scenario-based power flow problems (block-OPF), here formulated as two-stage stochastic nonlinear programs. The scenarios can be stochastic or represent contingencies (which can be interpreted as stochastic outcomes with uniform distribution), as is the case of the very widely used security-constrained AC optimal power flow (SC-ACOPF) problem [START_REF] Capitanescu | State-of-the-art, challenges, and future trends in security constrained optimal power flow[END_REF]. SC-ACOPF is one of the core analyses undertaken in the planning, operational planning, and realtime operation of transmission systems [START_REF] Capitanescu | State-of-the-art, challenges, and future trends in security constrained optimal power flow[END_REF]. SC-ACOPF is run several times a day by many operators in the US and the world. For brevity, we will refer to such problems as stochastic.

The block structure of such problems is given by the different scenarios associated with the stochastic problem, leading to potential parallelism in both the evaluation of the derivatives and the solution of the resulting block-angular KKT system. The parallel solution of the block-OPF problem with a Schur complement approach has been studied extensively both with PIPS-NLP [START_REF] Chiang | Structured nonconvex optimization of large-scale energy systems using PIPS-NLP[END_REF][START_REF] Schanen | Toward multiperiod ac-based contingency constrained optimal power flow at large scale[END_REF] (multiprocessing) and with Beltistos [START_REF] Kardoš | Two-level parallel augmented Schur complement interior-point algorithms for the solution of security constrained optimal power flow problems[END_REF][START_REF] Kardoš | Beltistos: A robust interior point method for large-scale optimal power flow problems[END_REF] (multiprocessing + factorization of the dense Schur complement on the GPU). Compared to the state-of-the-art solver Beltistos, our approach carries out almost all computation on the GPUs including a global CUDA-aware MPI reduction, from the evaluation of the derivatives to the assembling of the Schur complement. We test our implementation on the pre-exascale supercomputer Polaris, where each node is equipped with 4 A100 GPUs, and we solve block-OPF problems with up to 9,251 nodes.

Problem statement

In systems engineering, it is common to encounter optimization problems with relatively few degrees of freedom -"controls". Then, the goal is to appropriately fix the values for the degrees of freedom, e.g., by minimizing a given operational cost while satisfying the physical equations of the problem. In that context, the internal state of the system is described by a state variable x ∈ R nx , whose values depend on the current controls u ∈ R nu associated with the problem's degrees of freedom. If the problem is well-posed, this translates to the state equation g(x, u) = 0, where the function g exhibits the physical structure of the problem (e.g., a differential equation encoding a dynamics, or a nonlinear network flow associated with static balance equations). When the system faces uncertainties, it is often appropriate to choose a control u feasible under a finite set of conditions (or scenarios). That is, the control u must satisfy N different state equations

g i (x i , u) = 0 for all i = 1, • • • , N, (1)
where the state x i now depend on the current scenario i. The variables x i can be assimilated into a recourse variable. The N functions g 1 , • • • , g N define the block structure of the problem.

Block-structured nonlinear programs

In addition to satisfying the N state equations (1), we aim at minimizing the average operating costs on the N different scenarios. The corresponding problem formulates as a two-stage nonlinear program, which, in our case, is a nonlinear program with partially separable structure [START_REF] Demiguel | On decomposition methods for a class of partially separable nonlinear programs[END_REF]:

min x1,••• ,xN , u N i=1 f i (x i , u) s.t.      x i ≥ 0 , u ≥ 0 g i (x i , u) = 0 , h i (x i , u) ≤ 0 , ∀i = 1, • • • , N , (2)
with

f i : R nx × R nu → R, g i : R nx × R nu → R nx , h i : R nx × R nu → R m smooth
functions encoding the objective, the state equations, and the operational constraints, respectively. We note that the number of variables (N × n x + n u ) and constraints (N × (m + n x )) are linearly proportional to the number of blocks N .

In addition, if we introduce local control variables u 1 , • • • , u N with the additional coupling constraint u 1 = • • • = u N = u, we get a problem with a separable structure, solvable using the primal decomposition method; at the expense of increasing the search space [START_REF] Demiguel | On decomposition methods for a class of partially separable nonlinear programs[END_REF][START_REF] Ruszczynski | Interior point methods in stochastic programming[END_REF].

By introducing slack variables s 1 , • • • , s N , we rewrite (2) in standard form:

min x1,••• ,xN , s1,••• ,sN , u N i=1 f i (x i , u) s.t.      u ≥ 0 , x i ≥ 0 , s i ≥ 0 g i (x i , u) = 0 , h i (x i , u) + s i = 0 , ∀i = 1, • • • , N . (3)
We define y i ∈ R nx the multipliers (or adjoints) associated to the equality constraints g i (x i , u) = 0, z i ∈ R m the multipliers associated to the operational constraints h i (x i , u) + s i = 0, as well as λ, κ i , ν i the three multipliers associated to the respective bound constraints u ≥ 0, x i ≥ 0, s i ≥ 0. The Lagrangian associated to (3) is:

L(x, u, s; y, z, λ, µ, ν) := N i=1 f i (x i , u) + y i g i (x i , u) + z i h i (x i , u) + s i -κ i x i -ν i s i -λu , (4) with x := (x 1 , • • • , x N ), s := (s 1 , • • • , s N ), y := (y 1 , • • • , y N ), z := (z 1 , • • • , z N ).
To simplify the notations, we define the extended objective function and the extended constraints:

f (x, u) := i=1 f i (x i , u) , g(x, u) :=    g 1 (x 1 , u) . . . g N (x N , u)    , h(x, u) :=    h 1 (x 1 , u) . . . h N (x N , u)    .
We assume the functions f, g, h are twice differentiable. We denote

H = ∂ (x,u) h(x, u) ∈ R N m×(N nx+nu)
Jacobian of the inequality cons.

G = ∂ (x,u) g(x, u) ∈ R N nx×(N nx+nu)
Jacobian of the equality cons.

W = ∇ 2 (x,u) L(x, u, s; •) ∈ R (N nx+nu)×(N nx+nu)
Hessian of Lagrangian.

Interior-point method

The interior-point method (IPM) [START_REF] Nocedal | Numerical optimization[END_REF]Chapter 19] is a classical approach to solve (3).

KKT system

The Karush-Kuhn-Tucker (KKT) equations associated to (3) can be expressed as

∇ x f i + (G i x ) y i + (H i x ) z i -κ i = 0, ∀i = 1, • • • , N (5a)
N i=1 ∇ u f i + (G i u ) y i + (H i u ) z i -λ = 0, (coupling) (5b)
z i -ν i = 0, ∀i = 1, • • • , N (5c) g i (x i , u) = 0, ∀i = 1, • • • , N (5d) h i (x i , u) + s i = 0, ∀i = 1, • • • , N (5e)
X i κ i = 0, (x i , κ i ) ≥ 0, ∀i = 1, • • • , N (5f)
S i ν i = 0, (s i , ν i ) ≥ 0, ∀i = 1, • • • , N (5g)
U λ = 0, (u, λ) ≥ 0, (5h)
where U = diag(u), X i = diag(x i ), S i = diag(s i ).

The interior-point method uses a homotopy parameter µ > 0 to replace the complementarity constraints (5f)-(5g)-(5h) by the smooth approximations: X i κ i = µe nx , S i ν i = µe m , U λ = µe nu (e n being the vector of all ones of dimension n). The resulting (smooth) system of nonlinear equations can be solved iteratively using Newton method, where at each iteration, the descent direction is updated by solving the following augmented linear system:

    W + Σ p 0 G H 0 Σ s 0 I G 0 0 0 H I 0 0         p d p s p y p z     = -     r 1 r 2 r 3 r 4     (6)
with

r 1 = ∇ x f + G x y + H x z -µX -1 e nx ∇ u f + G u y + H u z -µU -1 e nu , r 2 = z -µS -1 e m , r 3 = g(x, u), r 4 = h(x, u) + s. The primal descent direction p d decomposes as p d = (p x1 , • • • , p xN , p u ).

Block angular structure

The linear system ( 6) is sparse and symmetric indefinite, and can be factorized using the Bunch-Kaufman algorithm. However, it is often beneficial to exploit its blockangular structure. Indeed, both the Hessian of the Lagrangian and the Jacobians have a block-angular structure, given as

W =      W x1x1 W x1u . . . . . . W xN xN W xN u W ux1 . . . W uxN W uu      , G =    G 1 x1 G 1 u . . . . . . G N xN G N u    .
By reordering the linear system (6), we can expose the block-angular structure of the KKT system as:

     A 1 B 1 . . . . . . A N B N B 1 . . . B N A 0      (7) with A 0 = W uu , A i =     W xixi + Σ xi 0 G xi H xi 0 Σ si 0 I G xi 0 0 0 H xi I 0 0     , B i =   W xiu (G i u ) (H i u )   .
The block-angular structure ( 7) can be exploited to solve the KKT linear system in parallel using a Schur complement approach. In that case, the submatrices A i can be factorized independently to assemble the Schur complement in parallel [START_REF] Chiang | Structured nonconvex optimization of large-scale energy systems using PIPS-NLP[END_REF].

Condensation and reduction

Instead of reordering the augmented KKT system (6) as a block angular matrix [START_REF] Capitanescu | State-of-the-art, challenges, and future trends in security constrained optimal power flow[END_REF], we propose an alternative approach based on successive condensation and reduction of the KKT system, following the method introduced in [START_REF] Pacaud | Accelerating condensed interior-point methods on SIMD/GPU architectures[END_REF]. If the structure is welldefined, we show that we can condense the KKT system (6) to a dense matrix with size n u × n u in two steps: first, by removing the inequality constraints in [START_REF] Byrne | MPI. jl: Julia bindings for the Message Passing Interface[END_REF], then by exploiting the structure of the equality constraints to reduce the condensed system to a dense matrix. The condensation and reduction steps are illustrated in Figure 1.

Condensation step

The condensation step allows reducing the size of the KKT system drastically if the number of inequality constraints is large1 .

Proposition 2.1 (Condensed KKT system). The linear system (6) is equivalent to

K + Σ p G G 0 p d p y = - r 1 + H (Σ s r 4 -r 2 ) r 3 , (8)
where descent directions p s and p z are recovered as

K ∈ R (N nx+nu)×(N nx+nu) is the condensed matrix K := W + H Σ s H. The
p z = Σ s Hp d + r 4 -r 2 , p s = -Σ -1 s r 2 + p z . (9)
Proof. See [31, Theorem 2.2].

The condensed matrix K inherits the block-angular structure of the Hessian of the Lagrangian W .

Proposition 2.2. The condensed matrix K = W +H Σ s H has a block-angular structure, given as

K =      K x1x1 K x1u . . . . . . K xN xN K xN u K ux1 . . . K uxN K uu      (10)
where we have defined the condensed blocks

K xixi := W xixi + (H i xi ) Σ si H i xi , K uxi := W uxi + (H i u ) Σ si H i xi and K uu := W uu + N i=1 (H i u ) Σ si H i u .
Proof. This is proved by induction.

Reduction step

In addition, we can exploit the structure of the equality constraints g 1 , • • • , g N to further reduce the size of the linear system (8) down to a dense matrix with size n u × n u . Equation ( 10) exhibits the structure w.r.t. the state x and the control u, we rewrite as such the condensed KKT system (8) as

            K x1x1 K x1u (G 1 x1 ) . . . . . . . . . K xN xN K xN u (G N xN ) K ux1 . . . K uxN K uu (G 1 u ) . . . (G 1 u ) G 1 x1 G 1 u . . . . . . G N xN G N u                         p x1 . . . p xN p u p 1 y . . . p N y             = -            r1 1 . . . rN 1 r2 r1 3 . . . rN 3           
, where we have renamed the right-hand-side in (8) as r.

Proposition 2.3 (Reduction). Assume that for all i = 1, • • • , N the Jacobian matrices G i x ∈ R nx×nx are invertible. Then the linear system (8) is equivalent to

Kuu p u = -r 2 + N i=1 (G i u ) (G i x ) -ri 1 + K uxi -(G i u ) (G i x ) -K xixi (G i x ) -1 ri 3 ( 11
)
with Kuu := Z KZ and Z ∈ R (nu+N nx)×nu is the reduction operator defined as

Z =      -(G 1 x ) -1 G 1 u . . . -(G N x ) -1 G N u I      . ( 12
)
The descent directions p x and p y are recovered as

p i x = -(G i x ) -1 ri 3 + G i u p u p i y = -(G i x ) -ri 1 + K xixi p i x + K xiu p u . (13)
Proof. See [31, Theorem 2.1].

The reduction [START_REF] Demiguel | On decomposition methods for a class of partially separable nonlinear programs[END_REF] is equivalent to a Schur complement approach applied to the condensed KKT system [START_REF] Chiang | Structured nonconvex optimization of large-scale energy systems using PIPS-NLP[END_REF]. In Proposition (2.1), we have shown that the condensed matrix K has a block-angular structure. The associated condensed KKT system [START_REF] Chiang | Structured nonconvex optimization of large-scale energy systems using PIPS-NLP[END_REF] is also inheriting a block-angular structure in the form of [START_REF] Capitanescu | State-of-the-art, challenges, and future trends in security constrained optimal power flow[END_REF], where the blocks are given by

A 0 = K uu , A i = K xixi (G i x ) G i x 0 , B i = K xiu G i u . ( 14
) Proposition 2.4. Assume that for each i = 1, • • • , N the Jacobian G i x is invertible. Let S uu = A 0 -N i=1 B i A -1
i B i be the Schur complement associated to the blockangular system (7) with the matrices (A i , B i ) defined in [START_REF] Dunning | JuMP: A modeling language for mathematical optimization[END_REF]. Then, the Schur complement S uu is equal to the reduced matrix K uu defined in [START_REF] Demiguel | On decomposition methods for a class of partially separable nonlinear programs[END_REF]: S uu = Z KZ.

Proof. First, note that if the Jacobian G i

x is invertible, then the block matrix A i defined in ( 14) is also invertible, with

A -1 i = 0 (G i x ) -1 (G i x ) - -(G i x ) -K xixi (G i x ) -1 . (15)
Using ( 14)-( 15), we expand the expression of the terms in the sum constituting the Schur complement S uu :

B i A -1 i B i = K uxi (G i u ) 0 (G i x ) -1 (G i x ) - -(G i x ) -K xixi (G i x ) -1 K xiu G i u , = (G i u ) (G i x ) -K xiu + K uxi (G i x ) -1 (G i u ) -(G i u ) (G i x ) -K xixi (G i x ) -1 G i u .
Hence, the Schur complement

S uu = A 0 -N i=1 B i A -1 i B i expands as S uu = K uu - N i=1 (G i u ) (G i x ) -K xiu + K uxi (G i x ) -1 (G i u ) -(G i u ) (G i x ) -K xixi (G i x ) -1 G i u = Z KZ .
We recover the expression of the reduced matrix K uu in Proposition 2.3.

Discussion

Hence, we can interpret the reduction step as a Schur complement approach. Forming the Schur complement has always been the bottleneck when solving distributed block angular problems in parallel [START_REF] Chiang | Structured nonconvex optimization of large-scale energy systems using PIPS-NLP[END_REF][START_REF] Lubin | Scalable stochastic optimization of complex energy systems[END_REF]. Its reduction operation involves large memory transfers between the processes, with the number of transfers being on the order of O(log(p)), where p is the number of processes. Due to the quasi-shared memory architecture on GPUs, the reduction can be implemented efficiently [START_REF] Pacaud | Accelerating condensed interior-point methods on SIMD/GPU architectures[END_REF]. In the next section, we propose to extend [START_REF] Pacaud | Accelerating condensed interior-point methods on SIMD/GPU architectures[END_REF] to assemble the reduced matrix K uu using two levels of parallelism, using both MPI and CUDA, thus reducing the reliance on distributed memory.

Parallel implementation

In the previous section, we have detailed the structure of block-angular nonlinear programs and presented the condensation and reduction steps for the KKT system. The loose coupling between the blocks is favorable for parallelizing the evaluation of the derivatives and the solution of the block-angular KKT system. Globally, we can distribute the computation on different processes using MPI (coarse parallelism). Locally, we can further streamline the computation using GPU accelerators (fine-grained parallelism). This paradigm, with its two levels of parallelism, is directly in line with what is currently offered by the new exascale architectures, where each node has 4 to 8 GPUs, all sharing a unified memory for fast communication. We present in §3.1 how we streamline the evaluation of the model using automatic differentiation, and in §3.2 how we parallelize the solution of the KKT system.

Parallel automatic-differentiation

First, we present how to evaluate the model in parallel using automatic differentiation [START_REF] Griewank | Evaluating derivatives: principles and techniques of algorithmic differentiation[END_REF]. We illustrate the procedure in Figure 2. The goal of the algorithm is to streamline the evaluation of the N scenarios on N/M GPUs, M being the number of scenarios evaluated locally on each GPU (we suppose here that N is a multiple of M ).

root g1 , • • • , g4 g 1 3 , • • • , g 1 6 g 1 g 4 • • • g 13 g 16 • • • • • • • • •

Local parallelism

The first level of parallelism streamlines the evaluation of the model on SIMD/GPU devices. We have designed our implementation to run entirely on the GPU device, to avoid any data transfer between the host and the device.

3.1.1.1. Block evaluation. We suppose that the nonlinear functions (f i , g i , h i ) share the same structure, its expressions yielding the same Abstract Syntax Tree (AST) for all i = 1, • • • , M . We illustrate the block evaluation on a simple abstract tree, but the reasoning extends to more complicated structures. We suppose that for all i, the functions f i , g i , h i depend linearly on a nonlinear basis matrix ψ : R nx × R nu → R nb : that is, there exists three sparse matrices L f , L g , L h such that

f i (x i , u) = L f ψ(x i , u) , g i (x i , u) = L g ψ(x i , u) , h i (x i , u) = L h ψ(x i , u) . (16)
Suppose we aim to evaluate the M functions g 1 , • • • , g M in batch for the states x 1 , • • • , x M . The structure ( 16) is directly amenable for SIMD evaluation. We denote by

X M = (x 1 , • • • , x M ) ∈ R nx×M
the dense matrix obtained by concatenating the M states together. By using a proper GPU kernel or a parallel modeler, we can evaluate the basis in a SIMD fashion and build the matrix Ψ(X M , u)

:= ψ(x 1 , u), • • • , ψ(x M , u) ∈ R nb×M . Then, evaluating the functions g 1 , • • • , g M simul-
taneously translates to the evaluation of one SpMM product:

g 1 (x 1 , u), • • • , g M (x M , u) = L g Ψ(X M , u) ∈ R nx×M . ( 17
)
The total memory required in the two successive operations is O((n x + n b ) × M ), and depends linearly on the number of blocks M . We note the SpMM operations are generally implemented efficiently in the vendor library (cusparse for CUDA, rocSPARSE for AMDGPU).

3.1.1.2. First-order derivatives. Suppose that for a given i we have a differentiable implementation gb i : R nd → R nx associated to the function g i . We aim to evaluate the Jacobian-matrix products (∇g i )D for p tangents encoded in a matrix D ∈ R nd×p using forward-mode AD and operator overloading. This operation translates to propagating forward a vector of dual numbers. Denoting by d ∈ D nd p the dual number encoding the p tangents stored in D, evaluating (∇g i )D simply amounts to call gb i (d) and extract the results in the dual numbers returned as a result. As G i is sparse, we can apply the technique of Jacobian coloring [START_REF] Griewank | Evaluating derivatives: principles and techniques of algorithmic differentiation[END_REF] to compress the independent columns of the sparse matrix G i and reduces the number of required seeding tangents p needed to evaluate the full Jacobian.

Suppose now we want to evaluate the sparse Jacobians G 1 , • • • , G M in batch. As the functions g i are based on the same AST, their respective Jacobians G 1 , • • • , G M are sharing the same sparsity pattern. By seeding a matrix of dual numbers

D M = (d 1 , • • • , d M ) ∈ D nd×M p
, we can use the same operation as [START_REF] Gondzio | Parallel interior-point solver for structured linear programs[END_REF] to streamline the evaluation of the M Jacobian-vector products using the SIMD kernel Ψ(•) and SpMM operations:

gb 1 (d 1 ), • • • , gb M (d M ) := L g Ψ(D M ) ∈ D nx×M p . ( 18
)
Once the results are evaluated, it remains to uncompress the dual outputs to build the M sparse Jacobians G 1 , • • • , G M . Hence, we can streamline the evaluation of the Jacobian along with the number of tangents p and the number of blocks M . This comes at the expense of increasing memory usage to O((n x + n b + n d ) × M × p) (to store the dual matrices associated to the input, the intermediate basis Ψ and the output).

3.1.1.3. Second-order derivatives. The evaluation of the second-order derivatives follows the same procedure, using forward-over-reverse AD. For each i, we suppose available an adjoint function adj gb i : R nd × R nx → R nd which for any primal x ∈ R nd and adjoint y ∈ R nx evaluates the Jacobian-transpose vector product (G i (x)) y (reverse-mode). Using forward-mode AD on top of adj gb i , we can compute the second-order derivatives y ∇ 2 g i (x)V for p directions V by calling adj gb i (x, y).

Using Hessian coloring, we can compress the independent columns of the sparse matrix y ∇ 2 g(x) and reduce the number of seeding tangents p required to evaluate the full Hessian. We note that in general obtaining an adjoint adj gb i running in parallel is nontrivial due to potential race conditions incurred by the control flow reversal of the original code.

Computing the Hessian y ∇ 2 g i (x) in parallel for i = 1, • • • M amounts to defining two matrices of dual numbers

X M = (X 1 , • • • , X M ) ∈ D nd×M p , Y M = (y 1 , • • • , y M ) ∈ D nx×M p and evaluate ∇Ψ(X M ) L g Y M .
The dual outputs are uncompressed to build the M sparse Hessians (as the sparsity pattern of the Hessians is different than those of the Jacobians, the matrix X M employed here is different than the one used in ( 18)). The total memory required to store the duals is O((2n x + n d + n b ) × M × p). For more details, we refer to the vector forward mode as described in [START_REF] Griewank | Evaluating derivatives: principles and techniques of algorithmic differentiation[END_REF].

Global parallelism

Now, if we have several GPUs at our disposal, we can push the parallelism further by distributing the evaluations using multiprocessing and a Message Passing Interface (MPI) library. Coming back at our original problem (2), we illustrate in Figure 2 how to dispatch the evaluation of the N nonlinear constraints g 1 , • • • , g N (the same reasoning applies to the objectives f 1 , • • • , f N and the inequality constraints h 1 , • • • , h N ). We use the streamlined implementation described in the previous subsection to evaluate the constraints in a batch of size M : the first GPU evaluates the constraints g 1 , • • • , g M , the second GPU evaluates g M +1 , • • • , g 2M , and so on. In total, the evaluation of the N constraints requires N/M GPUs (if M = 1, each GPU evaluate one constraint; if M = N , we use only one GPU evaluating all the constraints).

The implementation has been designed to minimize the communication between the different processes: each batch g 1 , • • • , g M stores the data it needs locally, the only data exchange with the other processes being the vector of input and the vector of output. In addition, we will see in the next section we do not have to transfer the firstand second-order information if a parallel linear solver is being used.

Parallel KKT solver

By exploiting the block-angular structure of the KKT system, we can solve the Newton step in parallel using a Schur complement approach. The challenge lies in the computation of the Schur complement matrix

S = A 0 -N i=1 B i A -1 i B i . Each prod- uct B i A -1
i B i requires the factorization of the matrix A i and the solution of a linear system with multiple (sparse) right-hand-side A -1 i B i . State-of-the-art methods are evaluating the Schur complement using an incomplete augmented factorization applied on the auxiliary matrix A i B i B i 0 , as currently implemented in the Pardiso linear solver [START_REF] Petra | An augmented incomplete factorization approach for computing the Schur complement in stochastic optimization[END_REF]. Here, we use an alternative approach building on the reduced KKT system §2.3.2 (equivalent to the Schur complement approach). As the reduction can be streamlined on GPU accelerators [START_REF] Pacaud | Accelerating condensed interior-point methods on SIMD/GPU architectures[END_REF], this approach can assemble the Schur complement in parallel using CUDA-aware MPI. We illustrate the parallel computation of the Schur complement in Figure 3.

G 1:M x =    G 1 x1 . . . G M xM    , G 1:M u =    G 1 u . . . G M u    , H 1:M =    H 1 x1 H 1 u . . . . . . H M xM H M u    ,
and sparse Hessian

W 1:M =      W x1x1 W x1u . . . . . . W xM xM W xM u W ux1 . . . W uxM W uu      .
Once the sparse matrices are obtained, we recover the condensed matrix x Q = LU (P , Q being two permutation matrices), computing (G 1:M

K 1:M = W 1:M + (H 1:M ) Σ(H 1:M ) (Proposition
x ) -1 b translates to two backsolves (SpSV) and two matrix-vector multiplications (SpMV), as (G

1:M x ) -1 b = QU -1 L -1 P b.
Local reduction. Once the sparse matrices are built, we evaluate locally the reduced matrix K 1:M uu on the GPU, using div(n u , n batch ) + 1 matrix-matrix product K 1:M uu V (with V ∈ R nu×nbatch a dense matrix encoding n batch vectors of the Cartesian basis of R nu ). The evaluation of one batched matrix-matrix product

K 1:M uu V = (Z K 1:M Z)V proceeds in three steps (1) Solve T x = -(G 1:M x ) -1 (G 1:M u V ). (2) Evaluate L x L u := K 1:M xx K 1:M xu K 1:M ux K 1:M uu T x V .
(

) Set K 1:M uu V = L u -G 1:M u (G 1:M x ) -L x . 3
In total, we need 2 SpSM and 3 SpMM operations in the first step, 1 SpMM in the second step, and 2 SpSM and 3 SpMM operations in the third step, giving a total of 4 SpSM and 7 SpMM operations. More than the computation, the reduction is limited by the memory, as we have to store the three buffers L x , T x , T u with a total size of (2M × n x + n u ) × n batch . If n x is too large, it is in our interest to reduce M (by using more GPUs) or to reduce n batch (at the expense of computing more matrix-matrix product

K 1:M uu V ).
Global reduction. Once we obtain the locally reduced matrices

K nM +1:(n+1)M uu for n = 0, • • • , N/M -1, we can assemble the global reduced matrix K uu = N/M -1 n=0 K nM +1:(n+1)M uu
using one all reduce (MPI Allreduce) operation. The size of the reduced matrix K uu is n u × n u , hence limiting the memory transfer required in the algorithm.

Discussion

We have presented a practical way to assemble the Schur complement on multi-GPU architectures. The parallelism occurs both at the local level (SIMD evaluations on the GPUs) and at the global level (distributed computation with MPI). The algorithm has the advantage of assembling the sparse Jacobians and Hessians only locally, as the reduction occurs before proceeding to the memory transfer with MPI Allreduce. The reduced matrix has a dimension n u × n u , which compresses the memory transfer significantly if the number of degrees of freedom n u is small. However, this comes at the expense of storing a vector of dual numbers (whose memory is linearly proportional to the number of blocks M evaluated locally and the number of tangents p being employed to evaluate the sparse derivatives) and additional buffers in the reduction algorithm. In the next section, we will test an implementation of the algorithm on CUDA GPUs, and show that the algorithm is practical.

Numerical results

We demonstrate the capabilities of the algorithm we introduced in Section §3 on the supercomputer Polaris, using CUDA-aware MPI to dispatch the solution on multiple GPUs. We present in §4.1 the stochastic optimal power flow problem, and give in §4.2 detailed assessments of the algorithms we have introduced earlier in §3. Eventually, we present in §4.3 a benchmark comparing our parallel solution algorithm with a state-of-the-art solution method running on the CPU.

Settings

Case study: the block-structured optimal power flow

The stochastic optimal power flow problem aims at finding an optimal dispatch for the generators u. The solution u should minimize the operational costs while satisfying the physical constraints (power flow equations g(x, u) = 0, here playing the role of the state equations) and operational constraints (line flow constraints h(x, u) ≤ 0) on a given set of scenarios. Each scenario is assigned given load parameters (energy demands) and potential contingencies (line tripping). The values of the state x depend on the local scenario we are in, the state x being the recourse variable in our case. As such, the problem has a partially separable structure as introduced in Problem (2), the control u being shared across all scenarios. We refer to [START_REF] Capitanescu | State-of-the-art, challenges, and future trends in security constrained optimal power flow[END_REF] for the original presentation of the stochastic optimal power flow problem and to [START_REF] Chiang | Structured nonconvex optimization of large-scale energy systems using PIPS-NLP[END_REF][START_REF] Kardoš | Two-level parallel augmented Schur complement interior-point algorithms for the solution of security constrained optimal power flow problems[END_REF][START_REF] Kardoš | Structure-exploiting interior point methods[END_REF][START_REF] Lubin | Scalable stochastic optimization of complex energy systems[END_REF] for practical algorithms solving the stochastic optimal power flow problem (some also focus on the multistage setting, which is not covered in this article). For our benchmark, we look at reference instances provided by MATPOWER [START_REF] Zimmerman | MATPOWER: Steady-state operations, planning, and analysis tools for power systems research and education[END_REF], whose characteristics are detailed in Table 1. We recall that in our case, the size of the Schur complement matrix Kuu is given by the number of controls n u .

Implementation

The algorithm has been implemented entirely in Julia 1.8. The Schur complement approach has been developed as an extension of the nonlinear optimization solver MadNLP [START_REF] Shin | Graph-based modeling and decomposition of energy infrastructures[END_REF], using CUDA-aware MPI as provided in [START_REF] Byrne | MPI. jl: Julia bindings for the Message Passing Interface[END_REF]. We have used the package ExaPF as a nonlinear modeler for the optimal power flow problem. All the results presented here have been generated on the supercomputer Polaris equipped with a total of 560 nodes, each node having with 1 CPU and 4 A100 GPUs.

Assessment of the parallel implementation

Assessing the performance of the parallel automatic differentation

We first assess the performance of the parallel automatic differentiation we introduced in §3.1 in a multi-GPU setting. We compare the performance we obtain with a CPU implementation. We use case1354pegase as a representative instance, and display the time spent in the automatic differentiation as we increase the total number of scenarios N . The results are displayed in Figure 4.

We observe that the computation time depends linearly on the number of scenarios, as expected. For N = 8, it is not worthwhile dispatching the evaluation on multiple GPUs as the problem is small enough to be evaluated on a single GPU. For N = 512, the evaluation time is 12.3s on the CPU, compared to 0.50, 0.41, 0.31, and 0.28s using 1, 2, 4 and 8 GPUs, respectively. Hence, we get a 40x speed-up when evaluating the derivatives in a multi-GPU setting, and it is not worthwhile to use more than 4 GPUs (one node).

Assessing the performance of the parallel KKT solver

We proceed to the same performance analysis to assess the performance of the parallel KKT solver detailed in §3.2. We compare the time required to evaluate the full solution of the KKT system afresh (including reduction time, factorization time and backsolve time) on case1354pegase as we increase the number of scenarios N . As a reference, we give the time taken by the sparse linear solvers HSL MA27 (single-threaded) and HSL MA57 (multi-threaded). The results are displayed in Figure 5.

On the left, we display the evolution of the time spent in the linear solver as we increase the number of scenarios. For N = 512, we observe that we get a linear speedup as we increase the number of GPUs: using 8 GPUs, the parallel KKT solver is 40x faster than using HSL MA27 on the CPU. Interestingly, we observe that HSL MA57 is not faster than HSL MA27, despite being multithreaded. This is consistent with the observation made in [START_REF] Tasseff | Exploring benefits of linear solver parallelism on modern nonlinear optimization applications[END_REF], and illustrates the difficulty of parallelizing effectively the sparse LDL factorization (Bunch-Kaufman). On the right, we display a performance profile detailing the time spent in MA27 and the parallel KKT solver on case1354pegase with N = 512 scenarios. We observe that most of the time in HSL MA27 is spent on factorizing the sparse augmented KKT system [START_REF] Byrne | MPI. jl: Julia bindings for the Message Passing Interface[END_REF]. On the other side, the factorization of the dense reduced matrix Kuu is trivial using LAPACK on the GPU; the bottleneck in the parallel KKT solver is the reduction algorithm itself. Fortunately, the reduction algorithm can run in parallel: we get a linear speed-up as we increase the number of GPUs used in the reduction algorithm.

Assessing the memory consumption

We have observed in §3.1 that the total memory required to store the duals is O((2n x + n d + n b ) × M × p), with M being the number of scenarios stored locally (M = N on 1 GPU, M = N/2 on 2 GPUs) and p the number of tangents. We display in Table 2 the memory taken by the automatic differentiation backend and by the parallel KKT solver for case1354pegase as we increase the number of scenarios N . We note that storing the duals is expensive in terms of memory, with up to 10.9GB for N = 512 on one GPU (as a reference, each NVIDIA A100 GPU on Polaris has 40GB of memory available). By evaluating the model on different processes with MPI, we can split the memory consumption on the different GPUs we are using, leading to better use of the resource at our disposal. First, we are interested in the scaling of the parallel algorithm in relation to the total number of scenarios N . We consider the case118 instance, and increase the number of scenarios N from 8 up to 2,048. For each N , we solve the block-structured OPF problem with MadNLP using our parallel KKT solver, and we compare with the performance we obtained with HSL MA27. The results are displayed in Figure 6. We observe that the solver HSL MA27 is initially faster than our parallel KKT solver, as the problem is too small to benefit from parallelism. However, as soon as N ≥ 16 the parallel KKT solver becomes competitive with HSL MA27. The relative performance is improving as we increase the number of scenarios N : for N = 512, we get a 68x speedup when using 8 GPUs, compared to the reference HSL MA27 (10.4s versus 712s). Interestingly, using 2 nodes (=8 GPUs) does not lead to any speed-up compared to a single node (=4 GPUs) if N ≤ 256; this setting is attractive only when the size of the problem becomes sufficiently large (N ≥ 1024) to compensate for the additional memory exchange.

Assessing the parallel performance w.r.t. the size of the problem

Second, we increase the size of the problems. We set a fixed number of scenarios N = 8, and look at the time to solution for case1354pegase, case2869pegase and case9241pegase. We detail the respective dimension of each problem in Table 3. We display the results in Figure 7, and give the detailed benchmark in Table 4. On the left (a), we display the total time required to find the solution of the three instances as a function of the number of GPUs; on the right (b), we show the performance profile associated to case9241pegase. In (a), we observe that overall the parallel algorithm is faster than the CPU implementation. The parallel algorithm scales well as we increase the number of GPUs we are using, the parallel algorithm being 35x faster than the reference when using 8 GPUs to solve case9241pegase. In (b), we detail the time

Assessing the parallel performance on a very large-scale instance

We finish our numerical experiments by solving a very large-scale instance: case1354pegase with N = 512 scenarios. The dimension of the resulting optimization problem is displayed in Table 3: the problem has more than 1 million variables, and 4 millions constraints. We solve this instance on resp. 1 node, 2, 4 and 8 nodes (resp. 4, 8, 16 and 32 GPUs). The results are displayed in Figure 8. We observe that the scaling is almost perfect when we use 2 nodes (8 GPUs) instead of a single node (4 GPUs) but we do not observe the same behavior when we increase the number of nodes to 4 and 8. On that instance, the gain we get when using 8 nodes (32 GPUs) is marginal compared to when using 4 nodes (16 GPUs): the solving time only decreases from 67s to 58s. This corroborate our observations: it is better to pack all the computation on a single node to use four A100 GPUs connected together via unified memory (NVLINK has a transfer rate of 600GB/s). When we have to use more than 2 nodes, the memory transfers are more involved as they have to pass through the network of the supercomputer.

Conclusion

We show promising results for leveraging massively parallel SIMD architectures like GPUs for block-structured nonlinear programs. The parallelism is applied to both the derivative evaluation and the solution of the KKT linear system. The main operation in the KKT algorithm is the assembling of the Schur complement, the factorization of the dense Schur complement being fast to carry on the GPU. At all levels, the method benefits significantly from the massive parallelism, achieving a speedup of around 40 for the derivatives compared to a sequential CPU implementation. The speedup is very application dependent, not least on the Hessian coloring and the problem's structure. The assembling of the Schur complement is bottlenecked by a distributed reduction operation bound by the interconnect's latency and throughput between GPUs. Current, so-called super nodes with multiple GPUs connected via fast networks like NVLINK greatly accelerate this operation. Lastly, our method is limited by the memory capacity of the GPU accelerators as it grows linearly with the number of problem blocks. In the context of ACOPF we are confident that upcoming GPUs will provide enough memory to solve a large number of scenarios in parallel, even for the largest grid instances (e.g., Eastern Interconnection with 70,000 nodes).

With the upcoming release of the Aurora supercomputer, these SIMD architectures will allow new science in regimes that were impossible with previous CPU architectures.

Figure 1 .

 1 Figure 1.. Successive reductions for a block-structured nonlinear problem with N = 3: Augmented system (6), Condensed system (8), Reduced system (11).

Figure 2 .

 2 Figure 2.. Parallel evaluation of the derivatives for g 1 , • • • g N on 4 GPUs: we have a total of N = 16 scenarios, each GPU evaluating M = 16/4 = 4 scenarios locally.

Figure 3 ..

 3 Figure 3.. Parallel computation of the Schur complement.

2 . 1 )

 21 using one SpGEMM operation and we factorize the matrix G 1:M x using a sparse LU factorization (potentially running in batch as the matrices G 1 x1 , • • • , G M xM are sharing the same sparsity pattern). Once the matrix G 1:M x is factorized as P G 1:M

Figure 4 .

 4 Figure 4.. Time spent to evaluate the model and its derivatives with automatic differentiation.

Figure 5 .

 5 Figure 5.. Time spent to solve the KKT system for case1354pegase.

Figure 7 .

 7 Figure 7.. For a fixed number of scenarios N = 8, (a) total time spent solving the block-OPF case1354pegase, case2869pegase and case9241pegase with MadNLP (b) performance profile for case9241pegase with varying number of GPUs.

  1354pegase

Figure 8 .

 8 Figure 8.. Solving case1354pegase with N = 512

Table 1 .

 1 . MATPOWER instances used in the benchmark.

	Name	#bus #lines #gen	n x	n u
	case118	118	186	54	181	107
	case1354pegase 1,354	1,991	260	2,447	519
	case2869pegase 2,869	4,582	510	5,227 1,019
	case9241pegase 9,241 16,049 1,445 17,036 2,889

Table 2 .

 2 . Memory consumption in MB 4.3. Parallel solution of the block-structured OPF problemWe analyze the parallel performance of our implementation on block-structured OPF problems.

		1 GPU		2 GPUs
	N	AD KKT solver	AD KKT solver
	8	171.1	92.3	85.5	48.1
	16	342.2	181.5	171.1	93.1
	32	684.3	360.0	342.2	183.2
	64	1,368.7	716.8	684.3	363.2
	128	2,737.3	1,430.5 1,368.7	723.4
	256	5,474.7	2,858.0 2,737.3	1,443.6
	512 10,949.3	5,712.8 5,474.7	2,884.1

4.3.1. Assessing the parallel performance w.r.t. the number of scenarios

  Figure 6.. Time to solve the block-structured OPF problem case118 as a function of the number of scenarios N .spent in the different operations for case9241pegase: the time spent to factorize the Schur complement with Lapack (using cusolve) is constant as the size of the Schur complement remains the same as we increase the number of GPUs. We observe that the time spent in the AD decreases linearly with the number of GPUs exploited, but the relative time spent in AD is negligible (less than 5% of the total time). Most of the time is spent in the parallel reduction, as discussed earlier in §4.2.2.

	Time to solution [s]	10 0 10 1 10 2 10 3	ma27 (CPU) polaris (1 GPUs) polaris (2 GPUs) polaris (4 GPUs) polaris (8 GPUs)	case118
			8 16 32 64 128 256 512 1,024 2,048 N scenarios
				N		nvar	ncon Kuu (mb)
			1354pegase	8	20,095	53,520	2.1
			2869pegase	8	42,835	119,216	7.9
			9241pegase	8	139,177	404,640	63.7
			1354pegase 512 1,253,383 4,425,280	2.1

Table 3 .

 3 . Dimension of the instances we have used in our benchmark.

Table 4

 4

						2869pegase			9241pegase
		#it AD KKT Tot. #it AD KKT Tot. #it	AD KKT	Tot.
	CPU	44 2.6	4.2	7.0	77 11.9	27.4 40.3 136 205.6 771.8 984.1
	1 GPU	44 0.3	1.8	2.1	93	1.1	11.7 12.8	98	5.5 112.3 117.8
	2 GPUs	44 0.3	1.1	1.4	93	0.8	7.4	8.2	98	3.4	56.8	60.2
	4 GPUs	44 0.3	1.0	1.3	93	0.8	5.7	6.5	98	2.3	35.8	38.1
	8 GPUs	44 0.2	1.0	1.2	93	0.6	5.1	5.7	98	1.4	26.4	27.7

.. Detailed results

It is equivalent to the normal equations in linear programming[START_REF] Nocedal | Numerical optimization[END_REF] Chapter 16, p.412]

Acknowledgment This material was based upon work supported by the U.S. Department of Energy, Office of Science, Office of Advanced Scientific Computing Research (ASCR) under Contract DE-AC02-06CH11347 and by NSF through award CNS-1545046. The authors gratefully acknowledge the funding support from the Applied Mathematics Program within the U.S. Department of Energy's (DOE) Office of Advanced Scientific Computing Research (ASCR) as part of the project ExaSGD. This research used resources of the Argonne Leadership Computing Facility, which is a DOE Office of Science User Facility supported under Contract DE-AC02-06CH11357.


Variations étymologiques sur l'étymon sino-tibétain 'étau, pinces'

Keywords: étymologie, changement sémantique, linguistique historique, rgyalronguique, sino-tibétain 2 L'étymon en birmo-rgyalronguique, tibétain et chinois Dans Rgya. oriental Japhug ta-mɢom

published or not. The documents may come

Variations étymologiques sur l'étymon sino-tibétain 'étau, pinces'

Yunfan

Introduction

La linguistique historique du sino-tibétain est encore dans son enfance. Bien qu'il y ait des travaux à visée générale, par exemple, sur la phylogénie sino-tibétaine (Zhang et al. 2019a ;[START_REF] Sagart | Dated language phylogenies shed light on the ancestry of Sino-Tibetan[END_REF], ou sur la reconstruction du chinois archaïque [START_REF] Baxter | Old Chinese : A new reconstruction[END_REF], les études détaillées se concentrant sur l'analyse d'étymologie individuelle ne sont pas suffisantes. 1 Cet article s'intéresse à un étymon sino-tibétain qui a subi plusieurs changements sémantiques. Cet étymon, ayant les sens 'étau/pince ', 'molaire', 'mâchoire', 'mâcher' et 'vulve, vagin', se trouve dans au moins 60 variétés dans la famille. Dans cet article, je me concentre d'abord sur ses réflexes dans les langues birmo-rgyalronguiques, le tibétain et le chinois (la section 2), ensuite, je présente de façon panoramique les réflexes dans les autres langues sino-tibétaines (la section 3), avant de passer à la conclusion avec une analyse qualitative sur la distribution géographique de tous les sens et colexifications (la section 4). 2zgám, en 'brongrdzong zgəmʔ et en guanyinqiao zgəḿ, toutefois, reste inexpliqué (Jacques 2014 : 200)

'Mâchoire'

Le chinois archaïque préserve l'étymon sous le sens 'mâchoire' : 頷 *[ɢ]ˤ[ə]mʔ (Zhang et al. 2019b). Bien que Schuessler (2006 : 270-271) Le choyul de Pubarong est la seule variété présentant ce sémantisme. Dans le choyul de Xinlong attesté par Huang and Dai (1992), 'vagin, vulve' se présente avec la forme sue⁵⁵, qui pourrait remonter à *sok,9 ce qui est probablement cognat avec le birman cok 'vagin, vulve' et le chang suk 'vagin, vulve'.10 Le choyul de Xinlong sue⁵⁵ est donc un mot hérité, et le choyul de Pubarong nɢów 'vagin, vulve' innové, selon notre analyse, de l'étymon sino-tibétain signifiant 'étau, pinces'.

L'étymon dans d'autres langues sino-tibétaines

L'étymon en question se trouve dans des variétés autre que le birmo-rgyalronguique, le tibétain et le chinois.

Réflexes à sens unique

'Étau, pinces' + 'mâcher'

Le kulung colexifie kamta 'étau, pinces' et le verbe kʰam-u 'mâcher'. L'innovation lexicale en choyul, où un changement récent d''étau, pinces' à 'vagin, vulve' a eu lieu, semble être un cas particulier. Ce changment sémantique est unique dans les langues sino-tibétaines. 11 Comme le montre sur la carte, le choyul se situe entre les langues rgyalronguiques (avec le sens hérité, 'étau, pince') et le minyag (avec un sens innové, 'molaire'). Par conséquent, nous pouvons imaginer deux chemins de changement sémantiques :

'Dent, molaire' + 'mâchoire'

1. 'étau, pince' → 'vagin, vulve' 2. 'étau, pince' → 'molaire' → 'vagin, vulve' Gong (2018a : 336-340) rapporte l'histoire de Vagina dentata circulant récemment dans la région rgyalronguique. La version documentée par Gong (2018a), qui est en zbu de Ngyaltsu, se relève d'une histoire sur un couple, où le mari fit une farce à leur voisin qui eut une liaison avec la femme, en disant que des dents poussèrent depuis son vagin. Quand le voisin se couchait avec la femme la même nuit, il poussa son genou que saisit ensuite la femme avec la main. Le voisin crut que les doigts de la femme étaient des dents qui poussèrent de la vulve. Selon Gong (2018a), cette histoire est plutôt d'origine d'Asie du nord-est, et s'est répdandue au sud-ouest. Il se peut que cette histoire soit liée au changement sémantique en question en choyul. Ce changement sémantique est typologiquement rare, le sens 'vagin, vulve' est sémantiquement distant d''étau, pinces' ou de 'molaire'. Il ne semble pas qu'un changement sémantique comme celui-ci puisse facilement apparaître. Il est plus probable qu'il y ait une motivation externe liant déjà 'vagin, vulve' à 'étau, pinces' ou à ''molaire', déclenchant le changement sémantique, et l'histoire de Vagina dentata peut très bien être cette motivation externe. Cette hypothèse implique que le second chemin de changement sémantique, via une étape intermédiaire, soit 'molaire', est préférée. Néanmoins, nous ne disposons pas d'évidence sur la date exacte où Vagina dentata entra dans la région rgyalronguique. Si cette histoire s'est introduite dans la région plus tard que le changement sémantique, notre hypothèse sera refutée.

  PST-faire.tomber.II.1SG 'J'ai fait tomber l'or que j'ai mis dans le pli de mes vêtements dans l'eau.' un temple avec une tour dans celle-ci. On a mis de beaux dessins de la résidence royales où habitent des personnes immortelles. (Leilin 04.34A.crimes par lesquels il est légal de répudier une femme, le premier est la jalousie.(Jacques 2007 : 99)

FIG. 1 .FIG. 2 .

 12 FIG. 1. Distribution des différents sens de l'étymon en question (sens uniques)

2.2 'Dent, molaire'

  , mais pourrait être lié à la préinitiale nasale manifestée par le tangoute 𘖩 Correspondances uvulaires entre le japhug et le khroskyabs La forme sqɔ en geshiza est incertaine. La rime -ɔ en geshiza correspond régulièrement à *-Vk/*-Vq dans d'autres langues apparentées. Par exemple, le geshiza ɕaeŋɔ 'hameçon' correspond au japhug ʑɴɢoʁ 'accrocher' et au siyuewu sŋǽʁ 'accrocher', et le geshiza mtsʰɔ 'tamiser' correspond au japhug tsʰaʁ 'crible', etc. Cette rime n'est jamais attestée comme correspondance à *-Vm ailleurs. Une possibilité est que la forme sqɔ soit un ancien état lié de la forme attendue †sqo, 4 puisque la voyelle -ɔ apparait dans certains états liés en geshiza : ɕua 'nuit' vs ɕɔ-tɕin (nuit-milieu) 'minuit', va 'cochon' vs vɔ-sʰu (cochon-rectum) Gong 2003 et N.Coʶ avec une initiale prénasalisée par Gong 2020 ; 2021) correspond régulièrement à *-Vm dans d'autres langues, tout comme quelques autres rimes reconstruites avec la même voyelle -o. Voir le tableau 3. 𗀫En tangoute, on ne trouve pas de trace correspondant aux préinitiales *sattestées dans d'autres langues ouest-rgyalronguiques et en tibétain, sinon on s'attendrait une voyelle dite "tendue" : †ɴqọʶ. la préinitiale rhotique r-, qui est probablement un préfixe pour les parties du corps, comparons le tibétain ང་ r-kaŋ 'jambe' et ་ rna 'nez'(Hill 2014 : 630). Ce préfixe est préservé dans rŋá 'visage', rqê 'gorge' et rŋî 'genou' en guanyinqiao. C'est également la raison pour laquelle on ne trouve pas de rdans le guanyinqiao zgəḿ (<

	4787 uʶ² (voie-LOC) 'dans la rue', la Japhug Khroskyabs (Wobzi) Khroskyabs (Siyuewu) Sens ɴqoʶ² 'étau'. ɴɢlɯt glǽ glód se casser ʑɴɢoʁ sŋá sŋǽʁ accrocher ɴɢu qú ɴqúz lâche La rime du tangoute 𘖩 Japhug Siyuewu Tangoute No. rime Sens jpum pəḿ 𘆋 1805 wọʶ¹ 73 épais tɤ-jpɣom rpʰəḿ 𗎈 4053 wọʶ¹ 73 glace kɯɕnom snəm 𘌇 5990 73 épis Minyag (Jiagenba) Minyag (Pengbuxi) Njorogs Japhug Sens sɔ́sɔ̠ ¹ xsôm χsɯm susceptibles d'être attachées. Par exemple, dans 𗵘 0020 tɕa¹-𗅁 2983 trois dzɔ́dzɔ̠ ¹ dzâm ndzom pont ndɔ́ndɔ̠ ¹ tʰôm tɤ-mtʰɯm viande TAB. 5. Correspondances de -ɔ en minyag forme 𗵘 tɕa¹ 'voie' est apparentée à l'état lié du wobzi tɕʰae-'voie'. On s'attendrait †tɕi² 0020 en tangoute s'il s'agissait d'une forme indépendante, comme celle en wobzi, tɕʰî 'voie'. Le tangoute 𗋕 2019 tʰa¹ semble être la seule forme avec le sens 'ceci' devant 𗅁 2983 uʶ² 'LOC'. nʕọʶ² tɤ-snom snəḿ 𘚋 0549 nʕọʶ¹ 73 soeur χsɯm xsəm 𘕕 5865 73 trois Les réflexes en guanyinqiao rqəm-et en tangoute 𘟗 0039 Les autres formes apparentées, telles que 𗣜 2173 tʰu² 'ceci' et 𘌽 5354 tʰɨ² 'ceci', n'y sont jamais qowʶr² (< prétangoute *rqVʶm²) attestées, voir l'exemple (2). L'exemple 2a montre 𗋕 2019 tʰa¹ 'DEM' suivi du locatif 𗅁 2983 uʶ², sọʶ¹ tɯ-rnom ʁrnəm 𗅾 2915 noʶ¹ 51 côte n-pʰəḿ 𗀫 0014 mpʰoʶ¹ 54 couvrir Le sens 'dent, molaire' est attesté de façon certaine en minyag et en tangoute, comme le montre le tableau 4. et l'exemple 2b montre l'indépendance syntaxique de 𘌽 5354 tʰɨ² 'ceci', en tant que sujet de la copule 𗨻 2226 weʶ² 'être'. dans.le.sein-LOC TAB. 2. 'saucisson'. (1) a. Tangoute TAB. 3. Correspondances des rimes tangoutes avec la voyelle -o 𗫟 3163 pʰowʶ¹-𗅁 2983 uʶ² 𗋚 2590 wɨ²-𗶠
	Branche Rgyal. Sud-Ouest 𗶱 4479 tɕọ¹	0014 ɕʕowʶ¹ Minyag (Jiagenba) Variété 4763 Minyag (Pengbuxi) servir 𗗭 5984 ntɕow²	mpʰoʶ¹ 'couvrir', en particulier, Forme rassembler həgɔ́'molaire' χəɢɔ¹ 'molaire' servir
	partage le même numéro de rime avec 𘖩 4787 Rgyal. occidental Tangoute 𗅌 3498 ʁoʶr² sinueux Khroskyabs (Guanyinqiao) rqəmɕɣî 'molaire' ɴqoʶ² 'étau', No. 54. 𘟗 0039 qowʶr² 'dent' 𘌁 nʁowʶr¹ courbé 5428 𗀫 0014 mpʰoʶ¹ couvrir 𗫟
	TAB. 4. 'Molaire, dent'
	Dans les deux variétés du minyag et le guanyinqiao, les formes sont des composés com-prenant -hə, -χə et -ɕɣî signifiant 'dent' ainsi que le réflexe de notre étymon en question, -gɔ, -ɢɔ et -rqəm, ce qui montre que le sens d'origine de cet étymon n'était pas 'dent, molaire'. Il n'est pas rare pour le même sens soit dérivé du nom d'un instrument dans
	les langues du monde. Par exemple, le français molaire vient du latin molaris 'meule'. De même, le chinois 臼齒 jiù-chǐ est compsé de 臼 jiù 'mortier' et 齒 chǐ 'dent'. 5 La rime -ɔ en minyag correspond régulièrement à -Vm dans d'autres langues apparen-tées, comme le montre le tableau 5.
				3163	pʰowʶ¹
	"dans le sein" est suivi de 𗅁 2983	uʶ² 'LOC'. Le marqueur 𗅁 2983	uʶ² 'LOC' est attesté avec des formes

4787

ɴqoʶ² 'étau', avec le numéro de rime 54 (reconstruit comme -ō par partagent *s-ⁿqom, voir la section 2.1) 'étau, pinces', ce qui préserve le sens originel. Dans la section 2.1 précédente, on a vu que le réflexe tangoute était 𘖩 4787 ɴqoʶ² 'étau, pinces'. Ici, on propose 𘟗 0039 qowʶr² 'dent' qui est formellement distincte de ce réflexe-là. Néanmoins, il convient de noter qu'il y a une alternance de rime entre -o et -ow en tangoute, montrée dans le tableau 6. Forme à -o Sens Forme à -ow Sens 𗰖 0478 nɕʕoʶ¹ rassembler 𗒋 3163 pʰowʶ¹ (dans) (...) sein TAB. 6. Alternance -o∼-ow en tangoute Il se peut donc que 𘟗 0039 qowʶr² 'dent' soit un état lié de 𘖩 4787 ɴqoʶ² 'étau, pinces', originellement dans un composé avec le mot hérité pour 'dent' (comme dans le cas du minyag), et réanalysée comme une forme indépendante. La seule paire de mots dans le tableau 6 ayant des cognats, c'est celle de 'couvrir, dans le sein/le pli des vêtements'. Les cognats en wobzi sont mpʰəḿ 'couvrir' et pʰəm(gə) 'dans le sein, dans le pli des vêtements'. Notons qu'il est obligatoire que la forme cognate soit suivie de -gə, originellement un marqueur du locatif (=gə), pour donner le sens 'dans le sein'. En tangoute, il semble que ce soit le même cas, voir l'exemple (1a), où 𗫟 2396 dzūʶ² PST-s'asseoir 'Elle s'est assise dans son sein.' (Leilin 06.06B.3) b. Wobzi

4 'Vagin, vulve'

  lie cette forme avec 含 *Cə-m-kˤ[ə]m 'tenir dans la bouche'. Cette liaison morphologique étant vraisemblable, il faudrait expliquer l'initiale uvulaire dans 含 *Cə-m-kˤ[ə]m 'tenir dans la bouche' s'opposant à l'uvulaire dans 頷 *[ɢ]ˤ[ə]mʔ 'mâchoire', ce qui dépasse la portée de cet article. à consonne orale et une ancienne prénasalisation de l'initiale correspondant à mdans le rgyalronguique oriental. Le préfixe est peut-être apparenté au possessif inaliénable tV-, comme en japhug ta-mɢom, cogtse tə-mkɐm et somang tə-mkám.

	Choyul	Siyuewu	Japhug Sens
	pʰǒw kówków kʰəm spʰəḿ rów ɕɑrəḿ 'propre' rom	couvercle ferme sec
	TAB. 7. Correspondances de -ow en choyul
	Les initiales uvulaires ont aussi des correspondances régulières dans les autres langues
	apparantées. Voir le tableau 8.
	Choyul Siyuewu	Wobzi	Japhug Sens
	qʰóʶ χqó χqó qʰóʶ ɴɢwaeʶ	ʁû nʁîd nʁîd qʰóz 'poche'	ʁû nʁaeî nʁaeî ɴqɑtsəɣ	tɯ-ku ɴqa ɕɤkʰoz sac de cuir tête difficile fatigué

Le mot 頷 *[ɢ]ˤ[ə]mʔ 'mâchoire' semble être une forme dialectale. Sa description dans le Fangyan (Langues régionales, 2000 ans avant le présent) est la suivante : 頷 *[ɢ]ˤ[ə]mʔ, 頤 *[ɢ](r)ə, mâchoire. Cela s'appelle 頷 *[ɢ]ˤ[ə]mʔ au Chǔ du sud, et 頜 #gop 6 au Qín et au Jìn. 頤 *[ɢ](r)ə est le mot dans la langue commune. 7 Si la description ci-dessus est valide, nous pouvons proposer deux hypothèses. La première est qu'il est possible que 頷 *[ɢ]ˤ[ə]mʔ apparaisse relativement tard remplaçant 頤 *[ɢ](r)ə dans des parlers régionaux. La deuxième est que 頤 *[ɢ](r)ə, dont le sens originel était 'menton', évolue vers le sens de 'mâchoire dans la langue standard, et que 頷 *[ɢ]ˤ[ə]mʔ ne survive ausens de 'mâchoire' au Chǔ du sud. 8

2.

Le choyul de Pubarong présente un cas particulier qui semble très récent. Le réflexe, nommément nɢów, développe le sens 'vagin, vulve' dans cette variété.

La rime -ow en choyul correspond régulièrement à *-Vm dans les autres langues rgyalronguiques. Voir le tableau 7. boiteux TAB. 8. Correspondances des consonnes uvulaires du choyul La préinitiale ndans nɢów 'vagin, vulve' est sans doute une fusion d'un préfixe

3.2 'Étau, pinces' + 'mâchoire'

  Trois langues, nommément le yamphu, le pa'O et le deuri semble d'avoir préservé le sens d'origine, 'étau, pinces'. Voir le tableau 9. Sauf le dumi bak-, les parties accompagnant l'étymon en question dans toutes les langues dans le tableau signifient 'dent' : jiːen bokar, iːen galo, hàen moyon, tã en idu, haen lalung, waen garo, paen nocte, ho̥ uen khiamngan, paen konyak, waen jingpho et sʌen thakali. Ce qui implique que 'molaire' dans ces langues-là sont des composé de 'étau, pinces' et 'dent'. Réflexes au sens 'molaire' Trois langues kirtanti, le yakha, le limbou et le chamling, développent le sens verbal de 'mâcher', comme l'illustre le tableau 12.Le bantawa et le tshangla (motuo) permet les deux sens, 'étau, pinces' et 'mâchoire', de coexister avec le même étymon, voir le tableau 13.On atteste la colexification d''étau, pinces' et 'dent, molaire' en bahing, où kʰamta signifie 'étau, pinces', et ka.kam 'dent'.

	Forme	Branche Sens	Branche	Langue Langue Forme	Forme Forme Sens
	Kiranti Karenique Pa'O Yamphu kamba kīm Bodo-Garo Deuri gum.je TAB. 9. Réflexes au sens 'étau' Le tableau 10 montre les quatorze langues où l'étymon en question signifie 'mâchoire'. Branche Langue Forme Hrusais Miji gɯm.jaʔ Tani Apatani iː.gam Tani Bengni gam.mɯ Kuki-Chin Matupi kɑm⁴ Kuki-Chin Kom rem Naga Ao Naga Lotha Naga Naga Sangtam Angami-Pochuri Sema Mruique Mru Bodo-Garo Dimasa Manang Manang kʌnta² Tamanguique Tamang kam¹ Tamaguique Thakali (Syang) kʌm⁵⁵ Kiranti Yakha Kiranti Limbou Kirtanti Chamling kem-(u) kʰamt-kʰam.maː kʰampʰo 'barbe' Branche Langue Forme kam am.kʰum.hi gɯm¹rɯ² o.kan tə²-kəm² rʌ.kám Tani Bokar jiː.kam Tani Galo iː.gam Tani Hill Miri gam.ɲa Digaro-Mishmi Idu tã.ɟam.bo Naga Moyon hà.bʌ.kʌm Bodo-Garo Lalung ha.gam Bodo-Garo Garo wa-gam Konyak Nocte pa³kʌmʔ² Konyak Khiamngan ²ho̥ u²¹gɪ¹m Konyak Konyak kams-(u) put together kʰoms(-u) put together pakt-(u) plant, sow pʰot-(a) cultivate mamt-(u) close the lips momt-(u) wrap katt-(u) reach with hand kʰott-(u) point (finger) kʰapt-(u) thatch, to roof kʰopt-(u) close kʰatt-(u) satisfy s.o.'s hunger kʰott-(u) chew for s.o. latt-(u) take out for s.o. lott-(u) bring out, bring up ŋams-(u) chew, munch ŋoms-(u) put in mouth and chew ramt-(u) overheat romt-(u) roast pá-gɯm Kachinique Jingpho wā-gām Kiranti Dumi TAB. 14. Alternances vocalique / d'aspiration potentielles en bantawa bak.ham Tamanguique Thakali (Tukche) kʌm-sʌ² 3.3 'Étau, pinces' + 'dent, molaire'
	TAB. 12. Réflexes au sens 'mâcher'
	Branche Langue		Étau/pinces	Mâchoire
	Kiranti Bodique Tshangla (Motuo) kam⁵⁵pa⁵⁵ Bantawa kʰam.pʰak 'pinces en bambou' yon.kʰom kam¹³ti⁵⁵
	TAB. 13. Réflexes aux sens 'étau/pinces' + 'mâcher'
	Notons que le bantawa présente deux formes différentes, kʰam-et kʰom-, pour les deux sens. Il s'agit probablement une alternance vocalique montrée dans le tableau 14. La
	fonction de cette alternance est toutefois encore inconnue.

TAB. 10. Réflexes au sens 'mâchoire' Le bokar, le galo et onze autres langues attestent le sens 'molaire', comme le montre le tableau 11. TAB. 11.

  Chepanguique Chepang de l'est hluk.gam.səyk ma.gam.səyk Le khaling est la seule langue où l'on atteste trois sens pour cet étymon, kɛm.ŋa.lu 'molaire, mâchoire' et kɛm-u 'mâcher'.

	On trouve cinq langues colexifiant 'dent, molaire' et 'mâchoire', montrées dans le ta-
	bleau 15.
	Branche	Langue	'Dent, molaire' Mâchoire
	Himalayen Tamanguique Gurung Lepcha Kiranti Hayu	fo-gom kãːsa kʰoː.bum	fo-gom kãːq kʰoː.bum
	Nung	Trung	skam	skam
	TAB. 15. Réflexes aux sens 'dent/molaire' + 'mâchoir'
	3.6 'Dent, molaire' + 'mâcher'
	Le thulung, quant à celui-ci, colexifie kam-li 'dent/molaire' et 'kam-' 'mâcher'.
	3.7 'Dent, molaire' + 'mâchoire' + 'mâcher'

On peut citer les travaux suivant qui traitent de ce type de problèmes :[START_REF] Thurgood | The Sino-Tibetan copula *wəy[END_REF],[START_REF] Bauer | Sino-Tibetan *vulva[END_REF],[START_REF] Jacques | A note on the etymology of the Tangut name Ngwemi[END_REF]), Gong (2018b).

Les sources des formes citées sont : japhug(Jacques 2015), cogtse (notes de terrain de Lin You-Jing et Lin 2016), wobzi (mes notes de terrain), njorogs[START_REF] Yin | Yelong Lawurongyu Yanjiu (业隆拉坞戎语研究) [Study on the 'Jorogs Lavrung language[END_REF], 'brongrdzong[START_REF] Sun | Collection & analysis of oral texts in the rGyalrong language group III[END_REF], guanyinqiao[START_REF] Huang | Lawurongyu yanjiu (拉坞戎语研究) [Study on the Lavrung language[END_REF], siyuewu (mes notes de terrain), geshiza[START_REF] Honkasalo | A grammar of Eastern Geshiza[END_REF], tangoute[START_REF] Li | Xia-Han zidian 夏漢字典 [Tangut-Chinese dictionary[END_REF], minyag (Jiagenba)[START_REF] Bai | A grammar of munya[END_REF], minyag (Pengbuxi)[START_REF] Gao | Description de la langue menya : phonologie et syntaxe[END_REF], chinois archaïque[START_REF] Baxter | Old Chinese : A new reconstruction[END_REF], choyul de Pubarong (notes de terrain de Guan Xuan), choyul de Xinlong(Huang and Dai 1992), bantawa[START_REF] Winter | A Bantawa dictionary[END_REF]. Toutes les autres formes citées proviennent de la base de données STEDT[START_REF] Richard | The Sino-Tibetan Etymological Dictionary and Thesaurus : STEDT project data resources and protocols[END_REF].

L'état lié est une forme non-autonome en tant qu'élément non-final d'une forme composée(Jacques 2021 : 147). Dans les langues rgyalronguiques, l'état lié souvent présente un changement vocalique de sa contrepartie autonome.

Un relecteur propose qu'il ne soit pas impossible que le sens originel de l'étymon soit bien 'molaire', et que l'ajout du morphème signifiant 'dent' s'effectue après le nouveau sens 'mâchoire' apparaisse pour relever l'ambigüité. Bien que cette analyse alternative soit possible, il semble moins économique que la nôtre dans cet article.

Il s'agit de ma reconstruction provisoire, car ce mot n'a pas été reconstruit dans[START_REF] Baxter | Old Chinese : A new reconstruction[END_REF].

頷、頤,頜也。南楚謂之頷。秦晉謂之頷。頤,其通語也。

Je dois cette deuxième hypothèse à l'un des relecteurs anonymes.

Choyul (Xinlong) psue⁵⁵ 'brillant' :: siyuewu fsôɣ 'brillant (*fso(ˠ)k)

Voir aussi[START_REF] Bauer | Sino-Tibetan *vulva[END_REF] pour une discussion sur cet étymon.