Xinzhe Wu
email: xinzhe.wu.etu@univ-lille.fr.

Serge G Petiton
email: serge.petiton@univ-lille.fr.

A Distributed and Parallel Asynchronous Unite and Conquer Method to Solve Large Scale Non-Hermitian Linear Systems with Multiple Right-hand Sides

Keywords: Linear systems, Krylov subspace methods, unite and conquer, multiple right-hand sides, asynchronous communication 2010 MSC: 00A69, 65F10 BGMRES Krylov Subspace size

Many problems in the field of science and engineering often require to solve simultaneously large-scale non-Hermitian sparse linear systems with multiple right-hand sides (RHSs). Efficiently solving such problems on extreme-scale platforms requires the minimization of global communications, reduction of synchronization points and promotion of asynchronous communications. We develop an extension of the Unite and Conquer GMRES/LS-ERAM (UCGLE) method [1] by combining it with Block GMRES method to solve non-Hermitian linear systems with multiple RHSs. UCGLE is a hybrid method consisting of three computing algorithms with asynchronous communication that allow the use of approximate eigenvalues to accelerate to solve of linear systems and to improve their fault tolerance. In this paper, the variant of UCGLE with novel components and manager engine implementations is introduced. This engine is capable of allocating multiple Block GMRES at the same time, each Block GMRES solving the linear systems with a subset of RHSs and accelerating the convergence using the eigenvalues approximated by other eigensolvers. Dividing the entire linear system with multiple RHSs into subsets and solving them simultaneously with different allocated linear solvers allow localizing calculations, reducing global communication, and improving parallel performance. Meanwhile, the asynchronous preconditioning using eigenvalues is able to speed up the convergence and improve the fault tolerance and reusability. Numerical experiments using different test matrices on supercomputer ROMEO indicate that the proposed method achieves a substantial decrease in both computation time and iterative steps with good scaling performance.

Introduction

In this paper, we consider solving the system

AX = B, (1)
where A ∈ C n×n is a large, sparse and non-Hermitian matrix of order n, X

= [x 1 , • • • , x s ] ∈ C n×s and B = [b 1 , • • • , b s ] ∈
C n×s are rectangular matrices of dimension n × s with s ≤ n.

In this paper, the rectangular matrices such as B is also called multi-vector, which can be seen as the combination of s vectors b i ∀i ∈ 1, 2, • • • , s. This kind of linear systems with multiple RHSs arise from a variety of applications in different scientific and engineering fields, such as the Lattice Quantum Chromodynamics (QCD) [START_REF] Sakurai | Application of block krylov subspace algorithms to the wilson-dirac equation with multiple right-hand sides in lattice qcd[END_REF][START_REF] Nakamura | Modified block bicgstab for lattice qcd[END_REF][START_REF] Fiebach | Variants of the block-qmr method 740 and applications in quantum chromodynamics[END_REF], the wave scattering and propagation simulation [START_REF] Malhotra | Iterative solution of multiple radiation and scattering problems in structural acoustics using a block 745 quasi-minimal residual algorithm[END_REF], dynamics of structures [START_REF] Barbella | Block krylov subspace methods for the computation of structural response to turbulent wind[END_REF][START_REF] Ferraz | Block-krylov component synthesis and minimum rank perturbation theory for damage detection in complex structures[END_REF][START_REF] Nour-Omid | Short communication block lanczos method for dynamic analysis of structures[END_REF], etc. The block Krylov methods are good candidates if we want to solve these large linear systems at the same time because the block methods can expand the search space associated with each RHS and may accelerate the convergence. Another feature of block Krylov methods is that they can be implemented up not only the number of compute nodes and central processing unit (CPU) cores but also the heterogeneity of components by introducing graphics processing units (GPUs) and manycore processors. This results in the tendency of transition to multi-and many cores within computing nodes, which commu-30 nicate explicitly through faster interconnection networks. These hierarchical supercomputers can be seen as the intersection of distributed and parallel computing. Indeed, for a large number of cores, the communication of overall reduction operations and global synchronization of applications are the bot-35 tleneck. When solving linear systems by block Krylov methods on large-scale distributed memory platforms, the cost of using BLAS3 operations to enlarge search space and reduce the memory requirement is apparent: the communication bound of SpGEMM in each step of Arnoldi projection damages heavily 40 their performance, which cannot be compensated by the advantages of the block methods. Even using classic Krylov methods, such as GMRES (Gen-eralized Minimal Residual method), to solve a large-scale problem on parallel clusters, the cost per iteration of them becomes the most significant concern, typically caused by communication and synchronization overheads. Consequently, large scalar products, overall synchronization, and other operations involving communication among all cores have to be avoided. The numerical applications should be optimized for more local communication, less global communication and synchronization. That is the reason for the recent tendance to study the communication avoiding techniques for linear algebra operations [START_REF] Demmel | Avoiding communication in sparse matrix computations[END_REF][START_REF] Hoemmen | Communication-avoiding krylov subspace methods[END_REF][START_REF] Carson | Communication-avoiding krylov subspace methods in theory and practice[END_REF] and different pipelined strategies for Krylov methods [START_REF] Ghysels | Hiding global 770 communication latency in the gmres algorithm on massively parallel machines[END_REF][START_REF] Morgan | A stochastic performance model for pipelined krylov methods[END_REF][START_REF] Cools | The communication-hiding pipelined bicgstab method for the parallel solution of large unsymmetric linear systems[END_REF]. For benefiting the full computational power of such hierarchical systems, it is central to explore novel parallel methods and models for the solving of linear systems. These methods should not only accelerate the convergence but also have the abilities to adapt to multi-grain, multi-level memory, to improve the fault tolerance, reduce synchronization and promote asynchronization.

We propose to combine the Block GMRES (BGMRES) [START_REF] Vital | Etude de quelques méthodes de résolution de problemes linéaires de grande taille sur multiprocesseur[END_REF] with UCGLE [START_REF] Wu | A distributed and parallel asynchronous unite and conquer method to solve large scale non-hermitian linear systems[END_REF] to solve Equation [START_REF] Wu | A distributed and parallel asynchronous unite and conquer method to solve large scale non-hermitian linear systems[END_REF] in parallel on modern computer architectures. UCGLE is a hybrid method which composes three computing components with asynchronous communication: ERAM (Explicitly Restarted Arnoldi Method), GM-RES and LS (Least Squares) components. GMRES Component is used to solve the systems, and LS Component performs as a preconditioner which uses the eigenvalues approximated by ERAM Component to speed up the convergence. LS Component is implemented based on the Least squares polynomial method proposed by Saad [START_REF] Saad | Least squares polynomials in the complex plane and their use for solving nonsymmetric linear systems[END_REF] in 1987. Compared with the traditional hybrid methods, the key features of UCGLE are its distributed and parallel asynchronous communications and the implementation of a manager engine among three components, which target at the extreme-scale supercomputing platforms. The asynchronous communication overlaps the overhead of overall synchronization and improves the fault tolerance and reusability. In this paper, firstly, we develop a block version of Least Squares Polynomial (B-LSP) method based on [START_REF] Saad | Least squares polynomials in the complex plane and their use for solving nonsymmetric linear systems[END_REF], then replace the three computing components of UCGLE respectively by the BGMRES, Shifted Krylov-Schur (s-KS), and B-LSP. Additionally, in order to solve linear systems with multiple RHSs and reduce the global communication, we design and implement a new manager engine to replace the former one in UC-GLE. This novel engine allows to allocate and deploy multiple BGMRES and/or s-KS Components at the same time and support their asynchronous communications. Each allocated BGMRES is assigned to solve the linear systems with a subset of RHSs. This extension is denoted as multiple-UCGLE or m-UCGLE even though the ERAM Component is replaced by s-KS method.

The paper is organized as follows. In Section 2, we present the related work, including the basic Krylov methods, BGM-RES and UCGLE. The theoretical work and algorithms inside m-UCGLE are presented in Section 3, especially the mathematical definitions fo B-LSP extended from the Least Squares polynomial method. In Section 4, we present the implementation of a newly designed manager engine for m-UCGLE. In Section 5, we give the information of experimental hardware/software and test matrices. We evaluate the convergence performance of our m-UCGLE implementation 1 compared with conventional restarted BGMRES in Section 6 and we show that our approach has up to 9.5× speedup against the classic method. The parallel performance of m-UCGLE presented in Section 7 show its 105 better strong scalability against the classic restarted BGMRES, and the time to the solution can be significantly decreased. We give the conclusions and perspectives in Section 8.

Related Work

In this section, we give a glance at the basic Krylov Sub-110 space methods, BGMRES, and UCGLE.

Krylov Subspace Methods

In linear algebra, the m-order Krylov subspace generated by a n × n matrix A and a vector b of dimension n is the linear subspace spanned by the images of b under the first m powers

115 of A, that is K m (A, b) = span(b, Ab, A 2 b, • • • , A m-1 b),
The Krylov subspace methods are often used to solve largescale linear systems and eigenvalue problems.

GMRES is an efficient Krylov subspace method to solve 120 non-Hermitian linear systems. It approximates the solution x m of Ax = b starting from an initial guessed solution x 0 , with the minimal residual in a subspace K n (A, r 0 ), where r 0 = b -Ax 0 is the initial residual and K n is the n th Krylov subspace spanned by A from r 0 . This method was developed by Saad and Schultz 125 in 1986 [START_REF] Saad | Gmres: A generalized minimal residual algorithm for solving nonsymmetric linear systems[END_REF]. If GMRES method is restarted after a number, say m g , of iterations with x m g as a new initial guess, to avoid enormous memory and computation requirements with the increase of Krylov subspace projection number, it is called the restarted GMRES. A well-known difficulty of restarted GMRES is that it 130 can suffer from stagnation in convergence when the matrix A is not positive definite as the restarted subspace is often closer to the earlier one. The typical methods are to use preconditioning techniques, to deflate the eigenvalues/eigenvectors, to recycle the Krylov subspaces, etc.

135

Arnoldi algorithm is a widely used Krylov subspace method to approximate the eigenvalues of large sparse matrices under the Krylov subspace K n (A, v), where v is an initial given vector. The numerical accuracy of the computed eigenpairs by Arnoldi method depends highly on the size of K m (A, v) and the orthog-140 onality of selected basis. Generally, the larger the subspace is, the better approximation of the eigenpairs is. The problem is that firstly the orthogonality of the computed basis tends to degrade with each basis extension. Also, the larger the subspace size is, the larger memory is required. Hence available memory 145 may also limit the subspace size, and the achievable accuracy of the Arnoldi process. In order to overcome this, it is necessary to restart the Arnoldi algorithm, with explicit [START_REF] Saad | Numerical methods for large eigenvalue problems: revised edition[END_REF] or implicit [START_REF] Sorensen | Implicitly restarted arnoldi/lanczos methods for large scale eigenvalue calculations[END_REF] deflation of the unwanted eigenvalues. Krylov-Schur method proposed by Stewart [START_REF] Stewart | A krylov-schur algorithm for large eigenproblems[END_REF] is another variant of Arnoldi U j = AV j 6:

for i = 1, 2, • • • , j do 7:
H i, j = V T i U j 8:

U j = U j -V i H i, j 9:
end for 10:

U j = V j+1 H j+1, j 11:
end for 12:

W m = [V 1 , V 2 , • • • , V m ], H m = {H i, j } 1≤i≤ j+1;1≤ j≤m 13: Find Y m , s.t. ||β -H m Y m || 2 is minimized 14: X = X 0 + W m Y m 15: end function
BGMRES is a variant of GMRES, which aims to solve sparse non-Hermitian linear systems with multiple RHSs at the same time. As shown in Algorithm 1, the procedure of BGMRES is similar to classic GMRES except replacing SpMV in the Arnoldi projection operation with SpGEMM. Block implementation of GMRES also introduces extra scalar works comparing with running multiple times of GMRES in sequence. Regarding numerical performance, the main attraction of BGMRES is that it can enlarge the search space and potentially speed up the convergence. When solving Equation (1) by BGMRES with Krylov subspace size fixed as m, the search space would be:

K m = span(b 1 , Ab 1 , • • • , A m-1 b 1 , • • • , b s , Ab s , • • • , A m-1 b s )
As discussed in Section 1, the first limitation of BGMRES is caused by the instability and bad parallel performance of its SpGEMM on large-scale distributed memory platforms with the increase of RHS number s. Secondly, for special cases, the convergence by BGMRES cannot always be guaranteed even by enlarging the search spaces. Thus the techniques to deflate of eigenpairs should be applied to accelerate the convergence. That is the motivation for us to propose m-UCGLE which can reduce the global communication and accelerate the convergence by the deflation of eigenvalues.

Unite and Conquer GMRES/LS-ERAM Method

UCGLE is a hybrid method implemented with asynchronous communication. The design of UCGLE is inspired by the Unite and Conquer (UC) approach.

Unite and Conquer Approach

UC approach is proposed by Emad [START_REF] Emad | Unite and conquer approach for high scale numerical computing[END_REF] to explore the novel methods for modern computer architectures. It is a model for the design of numerical methods by combining different computing components using asynchronous communication, and all these components working for the same objective. Unite im-185 plies the combination of components, and conquer represents they work together to solve one problem. In the UC methods, different parallel computing components work independently, which can be deployed on various platforms such as P2P, cloud and the supercomputer systems, or different computing units on 190 the same platform.

UCGLE Implementation

UCGLE is a kind of UC approach. As presented in Section 1, it composes three computing components: ERAM, GM-RES, and LS. GMRES Component is used to solve the sys- vector for solving non-Hermitian linear systems. Better convergence acceleration and parallel performance of UCGLE for several test matrices compared with preconditioned GMRES are given in [START_REF] Wu | A distributed and parallel asynchronous unite and conquer method to solve large scale non-hermitian linear systems[END_REF].

The convergence acceleration of UCGLE is similar to a de- mial should be extended to support the solving of linear systems with multiple RHSs. Moreover, as shown in Figure 1, the manager engine in [START_REF] Wu | A distributed and parallel asynchronous unite and conquer method to solve large scale non-hermitian linear systems[END_REF] is static. In detail, four communicators of MPI COMM WORLD are created for the computing components through MPI Split using the different colors and keys.

235

The asynchronous communication among them is supported by the MPI non-blocking sending and receiving operations of inter-communicators. This former implementation of manager engine cannot fulfill the requirements of m-UCGLE. Indeed, in order to reduce the global communications and localize the computation, m-UCGLE should be able to allocate dynamically different numbers of asynchronous computing components, and each linear solver solves a subset of linear systems.

m-UCGLE for Multiple Right-hand Sides

The m-UCGLE proposed in this paper is an extension of UCGLE and is designed to solve linear systems with multiple RHSs. The three types of components inside UCGLE are respectively replaced by restarted BGMRES that we presented in Section 2.2, s-KS, and B-LSP methods. In this section, firstly, we introduce the s-KS used in m-UCGLE, then the mathematical definitions of B-LSP to support for linear systems with multiple RHSs is given. Finally, we present the numerical analysis and parameter definition of m-UCGLE.

Shifted Krylov-Schur Algorithm

UCGLE uses the dominant eigenvalues to accelerate the convergence of GMRES, and theoretically, the more eigenvalues are applied, the acceleration of Least Squares Polynomial will be more significant [START_REF] Wu | A distributed and parallel asynchronous unite and conquer method to solve large scale non-hermitian linear systems[END_REF]. In order to approximate more eigenvalues by ERAM Component, the easiest way is to enlarge the size of related Krylov Subspace. In m-UCGLE, we replace ERAM Component by s-KS method which is another variant of the Arnoldi algorithm with an effective and robust restarting scheme and numerical stability [START_REF] Stewart | A krylov-schur algorithm for large eigenproblems[END_REF]. The Krylov subspace of s-KS cannot be too large. Otherwise BGMRES Component is not able to receive the eigenvalues in time to perform the B-LSP acceleration. With the novel developed manager engine of m-UCGLE in this paper, several different s-KS components can be allocated at the same time to approximate efficiently the different part of dominant eigenvalues of matrix A, by the shift with different values and thickly restarting with smaller Krylov subspace sizes. The algorithm of s-KS is given in Algorithm 2.

Least Squares Polynomial for Multiple Right-hand sides

The Least Squares polynomial method is an iterative method proposed by Saad [START_REF] Saad | Least squares polynomials in the complex plane and their use for solving nonsymmetric linear systems[END_REF] to solve linear systems. It is applied to calculate a new preconditioned residual for restarted GM- If not satisfied, go to step 3 9: end function method, which is a block extension of Least Squares polynomial method to solve linear systems with multiple RHSs at the same time. The iterates of B-LSP method can be written as X n = X 0 + P d (A)R 0 , where X 0 ∈ C N×s is a selected initial guess 280 to the solution, R 0 ∈ C N×s the corresponding residual multivector, and P d a polynomial of degree d -1. We set a polynomial of n degree R n such that

R d (λ) = 1 -λR d (λ)
.

285

The residual of n th steps iteration R n can be expressed as equation

R n = R d (A)R 0 , with the constraint R d (0) = 1. We want to find a kind of polynomial which can minimize all ||R d (A)R (p) 0 || 2 , with p ∈ 0, 1, • • • , s -1, R (p)
0 the p th vector in the multi-vector R 0 and ||.|| 2 the Euclidean norm.

290

Suppose A is a diagonalizable matrix with its spectrum denoted as σ(A) = λ 1 , • • • , λ n , and the associated eigenvectors u 1 , • • • , u n . Expanding the each component of R n in the basis of these eigenvectors as as R (p) n = n i=1 R d (λ i )ρ i u i , which allows to get the upper limit of

||R (p) n || 2 with p ∈ 0, 1, • • • , s -1 as: 295 ||R (p) 0 || 2 max λ∈σ(A) |R d (λ)| (2)
In order to minimize the norm of R (p) n , it is possible to find a polynomial P d which can minimize the Equation ( 2

) ∀p ∈ 0, 1, • • • , s -1.
As presented in [START_REF] Wu | A distributed and parallel asynchronous unite and conquer method to solve large scale non-hermitian linear systems[END_REF], P d can be expanded with a basis of

300 Chebyshev polynomial t j (λ) = T j λ-c b T j c b
, where t i is constructed by an ellipse englobing the convex hull formulated by the computed eigenvalues, with c the centre of ellipse, and b the focal distance of this ellipse. P d is under form that P d = d-1 i=0 η i t i . The selected Chebyshev polynomials t i meet still the three terms 305 recurrence relation [START_REF] Nakamura | Modified block bicgstab for lattice qcd[END_REF].

t i+1 (λ) = 1 β i+1 [λt i (λ) -α i t i (λ) -δ i t i-1 ] (3)
For the computation of parameters H = (η 0 , η

H d = (η 0 , • • • , η d-1
) is presented in Algorithm 3, where A is a n × n matrix, B represents the multi-vector of RHSs, d is the degree of Least Squares polynomial, Λ r the collection of approximate eigenvalues, a, c, b the required parameters to fix an ellipse in the plan, with a the distance between the vertex and centre, c the centre position and b the focal distance. The iterative reccurence implementation of Equation ( 5) and ( 6) using the parameters gotten from the pre-treatement procedure to construct the restarted residual for BGMRES by B-LSP is given in Algorithm 4. In this algorithm, X 0 is the temporary solution in BGMRES before performing the restart. Compared with Least Squares polynomial method for single RHS, the difference in B-LSP is to replace the SpMV in each iteration step with SpGEMM, as shown in Equation [START_REF] Barbella | Block krylov subspace methods for the computation of structural response to turbulent wind[END_REF].

R n = m i=1 s j=1 ρ((R ( j) d )(λ i ) ι )u i + n i=m+1 s j=1 ρ((R ( j) d )(λ i ) ι )u i (7)
The first part is constructed with the m known eigenvalues used to compute the convex hull in B-LSP Component, and the second part represents the residual with unknown eigenpairs. Get X 0 , which is temporary solution in BGMRES 3:

R 0 = B -AX 0 , Ω 1 = R 0 4: for k = 1, 2, • • • , l do 5: for i = 1, 2, • • • , d -1 do 6: Ω i+1 = 1 β i+1 [AΩ i -α i Ω i -δ i Ω i-1 ]
7:

X i+1 = X i + η i+1 Ω i+1 8:
end for 9:

end for 10:

Update GMRES restarted residual by X d 11: end function In the practical implementation, for each time preconditioning by the B-LSP method, it is often repeated for several times to 345 improve its acceleration of convergence, that is the meaning of parameter l in Equation [START_REF] Ferraz | Block-krylov component synthesis and minimum rank perturbation theory for damage detection in complex structures[END_REF]. The B-LSP preconditioning applies R d as a deflation vector for each time restart of BGMRES. The first part in Equation ( 7) is small since the B-LSP finds R d minimizing |R d (λ)| in the convex hull, but not with the second 350 part, where the residual will be rich in the eigenvectors associated with the eigenvalues outside H k . As the number of approximated eigenvalues k increasing, the first part will be much closer to zero, but the second part keeps still large. This results in an enormous increase of restarted BGMRES precondi-355 tioned vector norm. Meanwhile, when BGMRES restarts with the combination of some eigenvectors, the convergence will be faster even if the residual is enormous, and the convergence of BGMRES can still be significantly accelerated.

m-UCGLE is a combination of different methods. Thus it 360 has a large number of parameters, which have impacts on the convergence. These parameters are listed and classified according to their relations with different components as follows:

requirement of m-UCGLE. Thus, in order to extend UCGLE method to solve non-Hermitian linear systems and to reduce 1) can be decomposed as:

B = [B 1 , B 2 , • • • , B k , • • • , B n g ] (8)
Each BGMRES[k] will solve the linear systems with multiple RHSs B k , which is a subgroup of B: Table 1 gives the comparison of memory and communication complexity of SpGEMM operation inside m-UCGLE and BGMRES with the same number of RHSs. The factors n, nnz, s, P g and n g represent respectively the matrix size, the number of non-zero entries of matrix, the number of multiple RHSs, the to-405 tal number of computing units for BGMRES Components, and the number of BGMRES Components allocated by the manager engine of m-UCGLE. The average memory requirement for BGMRES on each computing unit is O( nnz(1+s)

AX k = B k (9

P g

). For m-UCGLE, the matrix is duplicated n g times into different allo-410 cated linear solvers. Thus the required memory to store the matrix should be scaled with the factor n g comparing with BGM-RES. Due to the localization of computation inside m-UCGLE, the total number global communication can be reduced with a factor 1 n g comparing with BGMRES. In practice, the selection of the number to allocate the BGMRES Components should end while 10: end function make a balance between the increase of memory requirement and the reduction of global communication.

Here we present in detail the workflow of this new engine. In the beginning, the manager will simultaneously allocate the 420 required number of three kinds of computing components. For each BGMRES[k], it will load a full matrix A and its related subgroup B k , and then start to solve Equation ( 9) separately. Meanwhile, each s-KS[q] load a full matrix A from local, and start to find the required part of eigenvalues of A, through the 425 s-KS method, using different parameters such as the shift value σ q , the Krylov subspace size (m a ) q , etc. If the eigenvalues of the required number are approximated on s-KS[q], these values will be asynchronously sent to the manager process. The manager process will always check if new eigenvalues are available 430 from different s-KS Components, if yes, it will collect and update the new coming eigenvalues together and send them to B-LSP Component. B-LSP Component will use all the eigenvalues received from manager process to do the pre-treatment of the B-LSP, the parameters gotten will be sent back to the man-435 ager process. Immediately, these parameters will be distributed to BGMRES [k]. BGMRES[k] can use the B-LSP residual constructed by these parameters to speed up the convergence. If the exit signals from all BGMRES Components are received by manager process, it will send a signal to all other components 440 to terminate their executions.

The allocation of a different number of computing components is implemented with MPI SPAWN, and their asynchronous communication is assured by the MPI non-blocking sending and receiving operations between the manager process and each if Recv(Λ r ) then

count = 0 3: BGMRES(input: A, m, X 0 , B, output: X m ) 4: if ||B -AX m || < g then 5:
return X m 6:

Send (exit == T RUE) to manager process

Hardware/Software Settings and Test Sparse Matrices

After the implementation of m-UCGLE, we test it on the supercomputer with selected test matrices. The purpose of this section is to give the details about the hardware/software settings and test sparse matrices.

Hardware and Software Settings

Experiments were obtained on the supercomputer ROMEO 2 , a system located at University of Reims Champagne-Ardenne, France. Made by Atos, this cluster relies on totally 115 the Bull Sequana X1125 hybrid servers, powered by the Xeon Gold 6132 (products formerly Skylake) and NVidia P100 cards. Each dense Bull Sequana X1125 server accommodates 2 Xeon Scalable Processors Gold bi-socket nodes, and 4 NVidia P100 cards connected with NVLink. In total, this supercomputer includes 3,220 Xeon cores and 280 Nvidia P100 accelerators.

The MPI (Message Passing Interface) used is the OpenMPI 3.1.2, all the shared libraries and binaries were compiled by gcc (version 6.3.0). The released scientific computational libraries Trilinos (version 12.12.1) and LAPACK (version 3.8.0) were also compiled and used for the implementation of m-UCGLE.

Test Sparse Matrices

In order to test m-UCGLE with matrices of high dimensions, we use the Scalable Matrix Generator with Given Spectra (SMG2S) [START_REF] Wu | A parallel generator of non-hermitian matrices computed from given spectra[END_REF][START_REF] Wu | SMG2S Manual v1.0[END_REF] to generate different test matrices. SMG2S 3is an open source package implemented and optimized using MPI and C++ templates, which allows generating efficiently large-scale sparse non-Hermitian test matrices with customized eigenvalues or spectral distribution by users to evaluate the impacts of spectra on the convergence of linear solvers. Figure 3 gives an example of SMG2S to generate test matrix. SMG2S has good scalability and acceptable accuracy to keep the given spectra. One more important benefit of SMG2S is that the matrices are generated in parallel, the data are already allocated to different processes. These distributed data can be used directly by the users to efficiently evaluate the numerical methods without concerning the I/O operation. Various interfaces provided by SMG2S to different programming languages and computational libraries on various computer architectures make it easier to do the tests and benchmarks.

In order to generate the test matrices with given dimensions and user-defined spectral distribution, three parameters p, d and h should be specified. The definition in detail of these parameters can be referred to [START_REF] Wu | A parallel generator of non-hermitian matrices computed from given spectra[END_REF]. As in shown Fig. 3, h determine the lower band's bandwidth of generated matrices as the same value, and their upper band's bandwidth is bounded by 2pd. In this paper, denote SMG2S(p, d, h, spec) a collection of matrix generation suite determined by p, d, h and the prescribed spectral distribution function spec, which can create sparse non-505 Hermitian test matrices with different sizes.

Convergence Evaluation

In this section, we evaluate the numerical performance of m-UCGLE for solving non-Hermitian linear systems compared with conventional BGMRES. In the experiments, the total number of RHSs of linear systems to be solved for each test is fixed as 64. As shown in Figure 4, we propose three strategies to divide these systems into various subgroup: Moreover, for m-UCGLE with 4 or 16 allocated components, they can be applied either with or without the preconditioning of B-LSP using approximate eigenvalues. Denote the special variant of m-UCGLE without B-LSP preconditioning 525 as m-BGMRES. m-BGMRES is also able to reduce the global communications through allocating multiple BGMRES components by the manager engine. Table 2 gives the naming of the five alternatives to solve linear systems with 64 RHSs and 3 are generated by SMG2S(1, 3, 4, spec), the definition of spec functions for different tests are shown in Table 3. For example, the spec of TEST-1 is given as (rand(21.0, 66.0), rand(-21.0,24.0)), the first part rand(21.0, 66.0) defines that the real parts of given eigenvalues for TEST-1 are the floating numbers generated randomly in the fixed interval [21.0, 66.0], similarly its imaginary parts are randomly generated in the fixed interval [-21.0, 24.0]. The relative tolerance for the convergence test is fixed as 1.0 × 10 -8 , the Krylov subspace size m g is given as 40 for all tests, the number of times that B-LSP applied in m-UCGLE l, the degree of Least Squares polynomial in B-LSP, the number of BGMRES restarts between two times of B-LSP preconditioning are respectively set as 10, 15 and 1. For m-BGMRES(16)×4, m-BGMRES(4)×16, m-UCGLE(16)×4 and m-UCGLE(16)×4, their iteration steps in Table 3 are defined as the maximal ones among their allocated components to solve the systems.

Result Analysis

The iteration steps of different methods for convergence are shown in Table 3 (dnc in this table significates do not converge in 5000 iteration steps). In this table, the blue and red cells represent respectively the worst and the best cases for each test.

From Table 3, firstly we can conclude that the enlargement of the Krylov subspace by more RHSs in BGMRES is effective to accelerate the convergence for TEST-1, TEST-2, TEST-3, and TEST-6. However, this acceleration cannot always be guaranteed. Referring to TEST-4 and TEST-5, m-BGMRES( 16))×4 converge faster than BGMRES(64). Secondly, m-UCGLE components converge much faster than their related m-BGMRES with the same number of RHSs, except in TEST-6. In the TEST-2, TEST-3, TEST-4, and TEST-5 of this table, m-UCGLE with 16 RHSs works even much better than BGMRES with 64

RHSs. An extreme special case in TEST-4 shows that convergence of m-UCGLE(4)×16 and m-UCGLE(16)×4 have the speedup respectively 3× and 9.5× over m-BGMRES [START_REF] Vital | Etude de quelques méthodes de résolution de problemes linéaires de grande taille sur multiprocesseur[END_REF].

In conclusion, for most tests, the method that converges the fastest is m-UCGLE( 16

Parallel Performance Evaluation 580

The iteration step for convergence is not the only concern about the iterative methods. After the convergence comparison of m-UCGLE and BGMRES, in this section, we evaluate its strong scalability and performance on supercomputer ROMEO. One important concern of BGMRES is the time cost per iteration, due to the communication bound of SpGEMM. In order to evaluate the parallel performance of m-UCGLE on large clusters, we compare the strong scaling of m-BGMRES(4)×16, m-UCGLE(4)×16, m-BGMRES( 16)×4, m-UCGLE(16)×4 and BGMRES(64) by their average time cost per iteration.

For the strong scaling evaluation on CPUs, the methods tested are set up the same as in Section 6. The test matrix is generated by SMG2S(1, 3, 4, spec) with size fixed as 1.792 × 10 6 . No larger matrices are tested, due to the memory limitation during the Arnoldi projection of BGMRES. The Krylov subspace size m g for all methods is set as 40. The average time cost for these methods is computed by 100 iterations. Time per iteration is suitable for demonstrating scaling behavior. The total CPU core number for B-GMRES( 64 In Figure 5, we can conclude that the strong scaling of m-BGMRES(4)×16 and m-UCGLE(4)×16 perform very well, the strong scalability of the rest are bad, especially BGMRES(64). In the beginning, the scalability of m-BGMRES(16)×4 and m-UCGLE(16)×4 is good, but it turns bad quickly with the increase of CPU number. It is demonstrated that the properties of m-UCGLE to promote the asynchronous communication and localize of computation can improve significantly the parallel performance of BGMRES for solving systems with multiple RHSs. Additionally, for the m-BGMRES and m-UCGLE with the same number of RHSs, the time per iteration of former is a little less the latter, since m-UCGLE introduces the iterative operations (SpGEMM) by B-LSP preconditioning.

Since m-UCGLE uses additional computing units for other components especially s-KS Component, it is unfair only to compare the scaling performance that total CPU number of all BGMRES Components in m-UCGLE equals to these numbers of m-GMRES and BGMRES(64). Thus we plot two more curves of m-UCGLE(4)×16 and m-UCGLE( 16)×14 with all their CPU numbers (including the CPU of s-KS Component) in Figure 5.

The two additional curves are respectively the blue and green ones with the marker set as the cross. It is shown that m-630 UCGLE(4)×16 and m-UCGLE( 16)×4 can still have respectively up to 35× and 4× speedup per iteration against BGMRES(64).

Time Consumption Evaluation

After the evaluation of parallel performance, we compare the time consumption of all methods to solve large-scale lin-635 ear systems with multiple RHSs on CPUs. The test matrices are generated by SMG2S(1, 3, 4, spec), the same as the convergence evaluation in Section 6. The size of matrices are all 1.792 × 10 6 for the experiments on CPUs, the total numbers of CPU cores for different methods (either BGMRES Compo-640 nents in m-UCGLE and m-GMRES or conventional BGMRES) are respectively fixed as 1792, and the Krylov subspace size m g is set as 40. The results on CPUs is given in Table 4, where the blue and red cells represent respectively the worst and the best cases for each test.

645

We can find that for the TEST-2, TEST-3, TEST-4, TEST-4, m-UCGLE(4)×16 take the least time to converge. For TEST-1 and TEST-6, m-BGMRES(4)×16 takes a little less time than m-UCGLE(4)×16 to get the convergence. For an extremely special case TEST-4, m-UCGLE(4)×16 has about 100× speedup in 650 time consumption over BGMRES(64) to solve the linear systems with 64 RHSs.

Performance Analysis

In conclusion, m-UCGLE(4)×16 and m-GMRES(4)×16 with the most decrease of global communications have the best strong 655 scaling performance. m-UCGLE cost a little more time per iteration compared with m-BGMRES with the same number of RHSs, but this might be made up by its decrease of iterative steps with B-LSP preconditioning. The experiments in this section demonstrate the benefits of m-UCGLE to reduce global 660 communication and promote the asynchronization. Two more important points that cannot be concluded from the experiments in this paper are:

(1) The increase of memory requirement should be considered when dividing the whole RHSs into subsets;

665

(2) The number of RHSs per component of m-UCGLE to enlarge the search space and the computation time per iteration should be balanced to achieve the best performance.

Requirement of Workflow Environment for UC Approach

After the validation of the numerical and parallel perfor-670 mance of m-UCGLE, a major difficulty to profit from UC methods including m-UCGLE is to implement the manager engine which can well handle their fault tolerance, load balance, asynchronous communication of signals, arrays and vectors, the management of different computing units such as GPUs, etc. In this paper, we tried to give a naive implementation of the engine to support the management tasks on the homogeneous platforms based on MPI SPAWN and MPI non-blocking sending and receiving functionalities. Further development of this engine should be done to support the manager of the computing algorithms and their asynchronous communication on more different computer architectures and platforms, such as the multi-GPU connected with NVLink on ROMEO. The stability of this implementation of the engine cannot always be guaranteed. Thus we are also thinking about to select the suitable workflow/task based environments to manage all these aspects in the UC approach. YML 4 is a good candidate, which is a workflow environment to provide the definition of the parallel application independently from the underlying middleware used. The special middleware and workflow scheduler provided by YML al-690 lows defining the dependencies of tasks and data on the supercomputers [START_REF] Delannoy | Yml: un workflow scientifique pour le calcul haute performance[END_REF]. YML, including its interfaces and compiler to various programming languages and libraries, will facilitate the implementation of UC based methods with different numerical components. We plan to replace the manager engine of 695 m-UCGLE by YML in the future.

Conclusion and Perspectives

This paper presents m-UCGLE, an extension of distributed and parallel method UCGLE to solve large-scale non-Hermitian sparse linear systems with multiple RHSs on modern supercom-700 puters. m-UCGLE is implemented with three kinds of computational components which communicate by the asynchronous communication. A special engine is proposed to manger the communication and allocate multiple different components at the same time. m-UCGLE is able to accelerate the convergence, minimize the global communication, cover the synchronous points for solving linear systems with multiple RHSs on large-scale platforms. The experiments on supercomputers prove the good numerical and parallel performance of this method.

Various parameters have impacts on the convergence. Thus, the auto-tuning scheme is required in the next step, where the systems can select the dimensions of Krylov subspace, the numbers of eigenvalues to be computed, the degrees of Least Squares polynomial according to different linear systems and cluster architectures. The fault tolerance and reusability of m-UCGLE and the influence of multiple s-KS components with different shifted values need to be evaluated. Parallel performance on heterogeneous clusters (e.g., multi-GPU) should be tested with further development of the manager engine, and workflow based model can be applied to replace the proposed manager engine

Algorithm 1 for j = 1 , 2 ,

 112 Non-Restarted Block GMRES Algorithm 1: function BGMRES(input:A, m, B, X 0 ∈ C N×s of full rank, output: X) • • • , m do 5:

  195 tems, LS and ERAM Components work as the preconditioning part. The asynchronous communication of this hybrid method among three components reduces the number of overall synchronization points and minimizes global communication. The workflow of UCGLE with three computing components is given 200 by Figure 1: ERAM Component computes the desired number of dominant eigenvalues, and then sends them to LS Component; LS Component uses these received eigenvalues to generate a new residual vector, and sends it to GMRES Component; GMRES Component uses this residual as a new restarted initial 205

210

  flated preconditioner. The difference is that the asynchronous communication among the different computing components of the latter reduces the global communication and synchronization, improves the fault tolerance and the reusability of this method. The three computing components work independently 215 from each other, when errors occur inside of ERAM, GMRES or LS Component, UCGLE can continue to work as a standard restarted GMRES method. The materials for accelerating the convergence are eigenvalues. These values can be improved on the fly, introducing a continuous improvement of convergence 220 while solving the linear systems. With the help of asynchronous communications, approximated eigenvalues by ERAM Component can be stored into a local file and reused for the other solving procedures with the same matrix. UCGLE is a distributed and parallel method which can profit both shared memory and 225 distributed memory of computational architectures with multiple levels of parallelism.UCGLE is used to solve non-Hermitian with single RHS, thus its GMRES Component should be replaced by a BGMRES algorithm, and a deflation technique by Least Squares polyno-

  230

Figure 1 :

 1 Figure 1: Workflow of UCGLE method [1].

5 :

 5 275RES in UCGLE. In this section, we will present the B-LSP Algorithm 2 Shifted Krylov-Schur Method 1: function s-KS(input: A, x 1 , m, σ, output: Λ k with k ≤ p)2:A ← A -σI 3: Build an initial Krylov decompostion of order m 4: Apply orthogonal transformations to get a Krylov-Schur decompostion Reorder the diagonal blocks of the Krylov-Schur decompostion 6: Truncate to a Krylov-Schur decompostion of order p 7:Extend to a Krylov decomposition of order m 8:

Algorithm 3

 3 Least Square Polynomial Pre-treatement 1: function LS(input: A, B, d, Λ r , output: A d , B d , ∆ d , H) 2: construct the convex hull C by Λ r 3: construct ellispe(a, c, b) by the convex hull C 4: compute parameters A d , B d , ∆ d by ellispe(a, c, b) 5: construct matrix T (d + 1) × d matrix by A d , B d , ∆ d 6: construct Gram matrix M d by Chebyshev polynomials basis 7: Cholesky factorization M d = LL T 8: F d = L T T 9: H d satisfies min l 11 e 1 -F d H 10: end function 3.3. Analysis Suppose that the computed convex hull by B-LSP contains eigenvalues λ 1 , • • • , λ m , the restarted residual for BGMRES generated by B-LSP for solving Equation (1) can be divided into two parts:

Algorithm 4

 4 Update BGMRES residual by LS Polynomial 1: function LSUpdateResidual(input:A, B, A d , B d , ∆ d , H d ) 2:

385

  the global communication and synchronization points, we design and implement a new manager engine for m-UCGLE. As shown in Figure2, the new engine allows creating a number of different computing components at the same time. Suppose that we have allocated n g BGMRES Components, n k s-KS Components and 1 B-LSP Component. The exact implementation for s-KS, B-LSP, BGMRES Components and manager process are respectively given in Algorithms 5, 6, 7 and 8. Denote the BGMRES Components as BGMRES[k] with k ∈ 1, 2, • • • , n g , and the s-KS Components as s-KS[q] with q ∈ 1, 2, • • • , n a . The matrix B in Equation (

# 1 s-KS # 2 Figure 2 :

 122 Figure 2: Manager Engine Implementation for m-UCGLE.

445

  computing components. Same as UCGLE, m-UCGLE has multiple levels of parallelism for distributed memory systems: 1. Coarse Grain/Component level: it allows the distribution of different numerical components, including the precon-450 ditioning part (B-LSP and s-KS) and the solving part Algorithm 6 B-LSP Component 1: function B-LSP-EXEC(input: A, b, d) 2:

  2 https://romeo.univ-reims.fr Algorithm 8 Manger of m-UCGLE with MPI Spawn 1: function Master(Input : n g , n a )

Figure 4 :

 4 Figure 4: Different strategies to divide the linear systems with 64 RHSs into subsets: (a) divide the 64 RHSs into to 16 different components of m-UCGLE, each holds 4 RHSs; (b) divide the 64 RHSs into to 4 different components of m-UCGLE, each holds 16 RHSs; (c) One classic BGMRES to solve the linear systems with 64 RHSs simultaneously.

510 6 . 1 .

 61 Specific Experimental Setup

515( 1 )

 1 1 group with all 64 RHSs solved by classic BGMRES (shown by Figure 4(c)); (2) 4 allocated BGMRES Components in m-UCGLE with each holding 16 Rhs (shown by Figure 4(a)); (3) 16 allocated BGMRES Components in m-UCGLE with each 520 holding 4 Rhs (shown in Figure 4(b)).

Table 3 :

 3 Iteration steps of convergence comparison (SMG2S generation suite SMG2S(1, 3, 4, spec), relative tolerance for convergence test = 1.0 × 10 -8 ), Krylov subspace size m g = 40, s use = 10, d = 15, L = 1, dnc = do not converge in 5000 iteration steps). Method spec m-BGMRES(4)×16 m-UCGLE(4)×16 m-BGMRES(16)×4 m-UCGLE(16)×4 BGMRES(64) TEST-1 (rand(21.0, 66.0), rand(-21.0,24.their allocated components and the numbers of RHSs per component. These 64 RHSs for the tests are all generated in random with different given seed states. All the test matrices from TEST-1 to TEST-6 in Table

  )×4. The combination of enlarging the search space by enough number of RHSs and preconditioning by B-LSP makes substantial acceleration on the convergence of solving linear systems. Moreover, compared with classic BGMRES, due to the preconditioning of B-LSP, m-UCGLE with less RHSs and smaller search space can still have a better 575 acceleration of the convergence. The potential damages caused by the localization of computation and reduction of global communications with less RHSs of each component in m-UCGLE can be covered by the B-LSP preconditioning.

  ) and all BGMRES Components in m-UCGLE and m-BGMRES ranges from 224 to 1792. Thus for each BGMRES Component of m-BGMRES(4)×16 and m-UCGLE(4)×16, the number ranges from 14 to 112. Similarly, for each BGMRES Component of m-BGMRES(16)×4 and m-UCGLE(16)×4, this number ranges from 56 to 448. All the tests allocate only 1 s-KS Component which always has the same number of CPU cores with each BGMRES Component inside m-UCGLE.

  LL T by the Cholesky factorization. The parameters H can be computed by a least squares problem of the formula:min l 11 e 1 -F d H(4)With the definition of Ω i ∈ R N×s by Ω i = t i (A)R 0 , we can and

, • • • , η d-1 ), we construct a modified gram matrix M d with dimension d × d, and matrix T d with dimension (d + 1) × d by the three terms

Table 1 :

 1 Memory and Communication Complexity Comparison between m-UCGLE and BGMRES.

			m-UCGLE	BGMRES	ratio
		Memory	O(	nnz(ng+s) Pg	)	O( nnz(1+s) Pg	)	ng+s 1+s
		Communication	O(	nnzsPg ng )		O(nnzsP g )	1 ng
	Algorithm 5 s-KS Component
	1: function s-KS-EXEC(input: A, m a , ν, r, a )
	2:	while exit==False do
	3:	s-KS(A, r, m a , ν, a , output: Λ r )
	4:	Send (Λ r ) to LS
	5:	if Recv (exit == T RUE) then
	6:	Send (exit) to LS Component
	7:	stop
	8:	end if
	9:

3 :

 3 LS(input:A, b, d, Λ r , output: A d , B d , ∆ d , H d )

	4:	Send (A d , B d , ∆ d , H d ) to GMRES Component
	5:	end if
	6:	if Recv (exit == T RUE) then
	7:	stop
	8:	end if
	9: end function
	Algorithm 7 BGMRES Component
	1: function BGMRES-EXEC(input: A, m g , X 0 , B, g , L, l,
		output: X m )
	2:

Table 2 :

 2 Alternative methods for experiments, and the related number of allocated component, Rhs number per component and preconditioners.

	Method	Component nb	RHS nb per component	Preconditioner
	BGMRES(64)	1	64	None
	m-BGMRES(16)×4	4	16	None
	m-BGMRES(4)×16	16	4	None
	m-UCGLE(16)×4	4	16	B-LSP
	m-UCGLE(4)×16	16	4	B-LSP

  BGMRES(64); test matrix size is 1.792 × 10 6 ; X-axis refers to the total number of CPUs from 224 to 1792; Y-axis refers to the average execution time per iteration. A base 2 logarithmic scale is used for all X-axis, and a base 10 logarithmic scale is used for all Y-axis.

	585	7.1. Strong Scalability Evaluation
			m-UCGLE(4)x16 (total CPUs)	m-BGMRES(4)x16
			m-UCGLE(4)x16	m-BGMRES(16)x4
			m-UCGLE(16)x4 (total CPUs)	BGMRES(64)
			m-UCGLE(16)x4
			10 1
		Time (s)	10 0
			10 -1
			224	448	896	1792
				CPU core number
		Figure 5: Strong scalability test on CPUs of solving time per iteration for m-
		BGMRES(4)×16, m-UCGLE(4)×16, m-BGMRES(16)×4, m-UCGLE(16)×4,

Table 4 :

 4 Consumption time (s) comparison on CPUs (SMG2S generation suite SMG2S(1, 3, 4, spec), the size of matrices = 1.792 × 10 6 , relative tolerance for convergence test = 1.0 × 10 -8 ), Krylov subspace size m g = 40, l = 10, d = 15, L = 1, dnc = do not converge in 5000 iteration steps).

	Method	spec	m-BGMRES(4)×16	m-UCGLE(4)×16 m-BGMRES(16)×4 m-UCGLE(16)×4 BGMRES(64)
	TEST-1 (rand(21.0, 66.0), rand(-21.0,24.0))	34.2	35.3	133.9	98.9	362.8
	TEST-2	(rand(0.5, 3.0), rand(-0.5,2.0))	dnc	40.9	231.3	111.5	580.6
	TEST-3	(rand(0.2, 5.2), rand(-2.5,2.5))	dnc	66.0	dnc	145.8	522.5
	TEST-4	(rand(-5.2, -0.2), rand(-2.5,2.5))	dnc	68.2	768.3	178.2	6829.3
	TEST-5	(rand(-0.23, -0.03), rand(-0.2,0.2))	132.3	50.1	209.4	120.5	1959.5
	TEST-6	(rand(-9.3, -3.2), rand(-2.1,2.1))	11.4	34.1	87.8	91.7	275.8

algorithm with an effective and robust restarting scheme and numerical stability. 1 https://github.com/brunowu/m-UCGLE

recurrence of the basis t i . M d can be factorized to be M d =

https://smg2s.github.io

Acknowledgment

The authors would like to thank ROMEO HPC Center Champagne Ardenne for their support in providing the use of cluster ROMEO. This work is funded by the MYX project of French 725 National Research Agency (ANR) (Grant No. ANR-15-SPPE-003) under the SPPEXA framework.

obtain the Equation [START_REF] Malhotra | Iterative solution of multiple radiation and scattering problems in structural acoustics using a block 745 quasi-minimal residual algorithm[END_REF], and in the end iteration [START_REF] Barbella | Block krylov subspace methods for the computation of structural response to turbulent wind[END_REF].

The pre-treatement of this method to obtain the parameters in the future. 4 http://yml.prism.uvsq.fr/


est destinée au dépôt


Anne-Sophie Michaud

Michel Pichavant

The H/F ratio as an indicator of contrasted wolframite deposition mechanisms

Keywords: wolframite, hübnerite, W deposition, magmatic control, W ore deposit, veins

Understanding wolframite deposition mechanisms is a key to develop reliable exploration guides for W. In quartz veins from the Variscan belt of Europe and elsewhere, wolframites have a wide range of compositions, from hübnerite-(MnWO 4 ) to ferberite-rich (FeWO 4 ). Deposition style, source of Mn and Fe, distance from the heat/fluid source and temperature have been proposed to govern the wolframite H/F (hübnerite/ferberite ratio) defined as 100 at. Mn / (Fe + Mn). The Argemela mineralized district, located near the world-class Panasqueira W mine in Portugal, exposes a quartz-wolframite vein system in close spatial and genetic association with a rare-metal granite.

Wolframite is absent as a magmatic phase, but W-rich whole-rock chemical data suggest that the granite magma is the source of W. Wolframite occurs as large homogeneous hübnerites (H/F = 64-75%) coexisting with montebrasite, K-feldspar and cassiterite in the latest generation of intragranitic veins corresponding to magmatic fluids exsolved from the granite. Locally, early hübnerites evolve to late more Fe-rich compositions (H/F = 45-55%). In a country rock vein, an early generation of Fe-rich hübnerites (H/F = 50-63%) is followed by late ferberites (H/F = 6-23%). Most Argemela wolframites have H/F ratios lower than at Panasqueira and other Variscan quartz-vein deposits which dominantly host ferberites. In greisens or pegmatitic veins, wolframites generally have intermediate H/F ratios. In those deposits, fluid-rock interactions, either involving country rocks (quartz-veins) or granite (greisens) control W deposition. In contrast, at Argemela, wolframite from intragranitic veins was deposited from a magmatic fluid. Differentiation of highly evolved peraluminous crustal magmas led to high Mn/Fe in the fluid which promoted the deposition of hübnerite. Therefore, the H/F ratio can be used to distinguish between contrasted deposition environments in perigranitic W ore-forming systems. Hübnerite is a simple mineralogical indicator for a strong magmatic control on W deposition.

Introduction

W is a critical metal, essential for many industrial applications. As for other metals and commodities, exploration for W now focuses on sub-surface to deep contexts in complex geological settings (e.g., Poitrenaud, 2018, www.almonty.com, www.wolfminerals.com.au). Targeting of mineralized systems requires the integration of geological, geophysical and geochemical data sets and their interpretation in the light of an ore deposition model [START_REF] Mccuaig | Translating the mineral systems approach into an effective exploration targeting system[END_REF][START_REF] Mccuaig | The mineral system concept: the key to exploration targeting[END_REF]. With the development of analytical tools, mineralogical compositions of both gangue and ore minerals are increasingly used to footprint the presence of proximal orebodies [START_REF] Belousova | Apatite as an indicator mineral for mineral exploration: trace-element compositions and their relationship to host rock type[END_REF][START_REF] Codeço | Chemical and boron isotopic composition of hydrothermal tourmaline from the Panasqueira W-Sn-Cu deposit[END_REF]Harlaux et al., 2018a;[START_REF] Neiva | Geochemistry of cassiterite and wolframite from tin and tungsten quartz veins in Portugal[END_REF].

W mineral deposits include quartz veins, stockworks, greisens, skarns and breccia pipes spatially associated with granitic bodies, and are thought to result from magmatic and hydrothermal processes, combined [START_REF] Bobos | Fe-, Fe,Mn-and Fe,Mg-chlorite: A genetic linkage to W, (Cu,Mo)mineralization in the magmatic-hydrothermal system of Borralha, Northern Portugal[END_REF][START_REF] Harlaux | Nb-ti-y-hree-w-u Oxide Minerals with Uncommon Compositions Associated with the Tungsten Mineralization in the Puy-lesvignes Deposit (Massif Central, France): Evidence For Rare-metal Mobilization By Late Hydrothermal Fluids with a Peralkaline Signature[END_REF][START_REF] Kelly | Geologic, fluid inclusion, and stable isotope studies of the tin-tungsten deposits of Panasqueira, Portugal[END_REF][START_REF] Ramboz | The P-V-T-X-f O2 evolution of H2O-CO2-CH4-bearing fluid in a wolframite vein: Reconstruction from fluid inclusion studies[END_REF][START_REF] Smith | Fluid inclusion and stable isotope constraints the genesis of the Cligga Head Sn-W deposit, SW England[END_REF]. Although the main parameters controlling the transport of W have been identified [START_REF] Che | Tungsten solubility in evolved granitic melts: An evaluation of magmatic wolframite[END_REF][START_REF] Heinrich | The chemistry of hydrothermal tin (-tungsten) ore deposition[END_REF][START_REF] Wood | The Hydrothermal Geochemistry of Tungsten in Granitoid Environments: I. Relative Solubilities of Ferberite and Scheelite as a Function of T, P, pH, and mNaCl[END_REF], there is as yet no comprehensive model for W deposition.

Quartz veins such as found in Southeast China and in the Variscan belt of Europe form an important class of W deposits. W, in the form of the dominant ore mineral wolframite ([Fe,Mn]WO4), is generally associated with other metals such as Sn, Cu, Ag. The group of quartz-wolframite deposits is diverse. It includes intragranitic (e.g., [START_REF] Smith | Fluid inclusion and stable isotope constraints the genesis of the Cligga Head Sn-W deposit, SW England[END_REF] and country rock (e.g., [START_REF] Kelly | Geologic, fluid inclusion, and stable isotope studies of the tin-tungsten deposits of Panasqueira, Portugal[END_REF] veins with a wide range of geometries, widths and extensions (e.g., [START_REF] Giuliani | Le gisement de tungstène de Xihuashan (Sud-Jiangxi, Chine): Relations granites, altérations deutériques-hydrothermales, minéralisations[END_REF]. Quartz and wolframite can be accompanied by various gangue (e.g., tourmaline, [START_REF] Codeço | Chemical and boron isotopic composition of hydrothermal tourmaline from the Panasqueira W-Sn-Cu deposit[END_REF] and ore (e.g., cassiterite, [START_REF] Noronha | Characterization and timing of the different types of fluids present in the barren and ore-veins of the W-Sn deposit of Panasqueira, Central Portugal[END_REF] minerals. Wolframites have compositions ranging from hübnerite (MnWO 4 ) to ferberite (FeWO 4 , e.g., Harlaux et al., 2018a;[START_REF] Lecumberri-Sanchez | Fluid-rock interaction is decisive for the formation of tungsten deposits[END_REF]. Factors of controls of wolframite composition are the subject of an intense debate. The deposition style, source of Mn and Fe, distance from the heat/fluid source and temperature have been proposed to govern the wolframite H/F (i.e., hübnerite/ferberite) ratio, hereafter calculated as 100 at. Mn / (Fe + Mn). [START_REF] Hosking | The significance of the variable iron content of the "Wolframites" of the west of England[END_REF] pointed out that wolframites in Cornish pegmatites are relatively Mn-rich while wolframites from quartz veins have high Fe contents. Greisens generally have wolframites with intermediate Fe/Mn (e.g., [START_REF] Breiter | Assessment of magmatic vs. metasomatic processes in rare-metal granites: A case study of the Cínovec/Zinnwald Sn-W-Li deposit[END_REF][START_REF] Hosking | The significance of the variable iron content of the "Wolframites" of the west of England[END_REF]. Several studies have suggested that Mn-rich wolframites precipitate at higher temperatures than Fe-rich wolframites (e.g., [START_REF] Oelsner | Über arzebirge Wolframite: Ber. Freib[END_REF][START_REF] Leutwein | Die Wolframit Grappe[END_REF][START_REF] Taylor | Manganese-iron ratios in wolframite, South Crofty Mine, Cornwall[END_REF]. However, attempts to use the H/F ratio as a geothermometer [START_REF] Moore | On the application of the hübnerite:ferberite ratio as a geothermometer[END_REF] have failed as temperature is not the only factor controlling wolframite composition [START_REF] Amossé | Physicochemical study of the hubnerite-ferberite (MnWO4-FeWO4) zonal distribution in wolframite (Mn X Fe (1-X) WO 4 ) deposits[END_REF].

Recent studies in the Central Iberian Zone (which includes the world class Panasqueira W deposit) have identified a quartz-wolframite vein system (Michaud et al., 2018b) spatially and genetically associated with a Rare-Metal Granite (RMG), the Argemela granite [START_REF] Charoy | Multistage growth of a rare-element, volatile-rich microgranite at Argemela (Portugal)[END_REF]. Field, structural, mineralogical and geochemical evidence suggest that the intragranitic wolframite-bearing veins are the expression of magmatic fluids exsolved from the granite at the magmatic-hydrothermal transition. Below, the main characteristics of the Argemela veins are summarized and compared with other W vein systems. Focus is placed on wolframite compositions, essentially Mn-rich at Argemela and Fe-rich in the neighboring Panasqueira and most other Variscan country rock vein systems. This first-order compositional contrast between wolframites suggests major differences in W deposition mechanisms. It is proposed that the H/F ratio can be used as an exploration tool to help distinguishing between contrasted types of W deposition environments.

Analytical techniques

Samples from the Argemela W veins were collected and thin as well as polished sections prepared.

Preliminary petrographic and textural observations were conducted with a petrographic microscope. Backscattered electron (BSE) images and EDS element mapping were obtained using a Merlin compact ZEISS Microscope equipped with a Bruker EDS detector and working under an acceleration voltage of 15 kV. Wolframites were analyzed by electron microprobe at ISTO, Orléans (France). Punctual analyses, traverses and elements distribution maps for major and minor elements were obtained with a CAMECA SX Five microprobe operated under an acceleration voltage and a beam current of 20 kV and 40 nA respectively. Standards included SnO 2 (SnLα), Fe 2 O 3 (FeKα), MnTiO 3 (MnKα and TiKα) and pure Nb (NbLα), Ta (TaLα) and W (WLα) metals.

The Argemela W vein system

The Variscan Argemela Mineralized District (AMD) occurs within the Central Iberian Zone (CIZ, Central Portugal) and is divided in two subsystems at < 0.5 km from each other but evolving separately. The first (the Argemela tin mine, ATM) consists of a vein swarm emplaced in low grade metamorphic rocks from the Neoproterozoic Schist and Greywacke Complex (SGC, [START_REF] Meireles | New data on the lithostratigraphy of Beiras Group (Schist Greywackes Complex) in the region of Góis-Arganil-Pampilhosa da Serra (Central Portugal)[END_REF]. Veins can be up to 2 m wide, are anastomosed and contain Sn and Li mineralization in the form of cassiterite and montebrasite. No W concentration has been recognized. The second (the Argemela Hill Top, AHT) comprises a RMG with disseminated Li, Sn, Nb, Ta mineralization in the form of magmatic montebrasite, cassiterite and columbo-tantalite. The granitic intrusion is associated with a network of intragranitic veins, and a few quartz veins occur in country rocks near the granite contact. Although wolframite is absent as a magmatic phase, W enrichment is marked by granite whole-rock concentrations ranging up to ~20 ppm W [START_REF] Charoy | Multistage growth of a rare-element, volatile-rich microgranite at Argemela (Portugal)[END_REF]. Three generations of intragranitic veins occur within the granite. Field evidence (confinement of veins inside the granite), structural (granite and veins are emplaced in the same structural context), mineralogical (veins lack hydrothermal alteration aureoles and the earliest have a "pegmatitic" texture and a granitic mineral assemblage) and stable isotope data indicate that the intragranitic veins are the expression of fluids exsolved from the granite at the magmatic-hydrothermal transition [START_REF] Michaud | From magmatic to hydrothermal Sn-Li-Nb-Ta-W mineralization: the Argemela District (Portugal)[END_REF]. Wolframite occurs in variable amounts in the third and latest generation of intragranitic veins which are numerous but of restricted size, never exceeding a 10 cm width. It is mostly associated with quartz, K-feldspar and montebrasite, with quartz and montebrasite only in some other veins, and cassiterite is occasionally present. Wolframites form mostly large homogeneous crystals with hübnerite compositions (early WII, H/F between 64 and 75%, Table 1 and Fig. 1). On the eastern border, the granite is affected by a dextral shear zone and the intragranitic veins host wolframites with a patchy texture made of early hübnerites evolving to more Fe-rich compositions (late WII, H/F from 45 to 55%, Table 1,Fig.1,e.g. Baumann 1964). This late wolframite ore deposition stage is followed by sulfides (arsenopyrite) and other phases (vivianite); cassiterite is transformed into stannite.

Wolframite is also present in one country rock vein crosscut by the granite. It appears as large crystals and in variable amounts associated with milky geodic quartz typical of hydrothermal veins in the AMD. In contrast to intragranitic veins, K-feldspar is absent. Wolframite crystals can be homogeneous but most show textural evidence for recrystallization (Fig. 2). Early wolframites are Fe-rich hübnerites (early WI, H/F from 50 to 63%, Table 1 and Fig. 1). The second generation crystallizes in fractures or near the margins of early wolframites and is a ferberite with high Fe contents (late WI, H/F = 6-23%, Table 1, Fig. 1 and2). This late stage is followed by the deposition of iron oxides.

Comparison with other wolframite vein systems

Wolframite occurs as the dominant ore mineral in several types of Variscan W deposits and is characterized by highly variable compositions, from hübnerite to ferberite. The Panasqueira W world class mine, located at only 13 km west from Argemela, provides an excellent comparison point with the AMD. The deposit consists of > 1000 nearly horizontal veins emplaced in SGC metamorphic rocks on top of a granite cupola [START_REF] Neiva | Geochemistry of cassiterite and wolframite from tin and tungsten quartz veins in Portugal[END_REF]Foxford et al., 2000). The hydrothermal veins, which can exceed 1 m widths, host wolframite, quartz, muscovite, tourmaline, topaz, arsenopyrite and early cassiterite [START_REF] Kelly | Geologic, fluid inclusion, and stable isotope studies of the tin-tungsten deposits of Panasqueira, Portugal[END_REF]Rye, 1979, Polya et al., 2000). Wolframites form large chemically homogeneous crystals with H/F ratios generally < 23% (e.g. ferberite), in strong contrast with most AMD wolframites (Table 2 and Fig. 1). Same as at Panasqueira, in the French Massif Central (FMC), most Variscan W deposits show wolframites with ferberite compositions [START_REF] Harlaux | Tracing metal sources in peribatholitic hydrothermal W deposits based on the chemical composition of wolframite: The example of the Variscan French Massif Central[END_REF] and references therein, Table 2 and Fig. 1). For example, the stockwork of La Bosse (north FMC) consists in a network of interconnected subhorizontal quartz veins emplaced in metamorphic rocks; wolframite has H/F ratios between 21 and 32% [START_REF] Aubert | Les coupoles granitiques de Montebras et d'Échassières:(Massif Central français) et la genèse de leurs minéralisations en étain, lithium, tungstène et béryllium[END_REF]Aïssa et al., 1987a, Table 2 and Fig. 1). Other representative examples of ferberite-bearing vein-type W deposits in the FMC include Enguialès, Leucamp and St-Goussaud (e.g., Harlaux et al., 2018 and references therein). Both the Panasqueira and the FMC deposits are representative of perigranitic hydrothermal systems developed in metamorphic country rocks. In the case of Panasqueira, ferberite deposition is controlled by fluid-rock interactions involving SGC country rocks [START_REF] Lecumberri-Sanchez | Fluid-rock interaction is decisive for the formation of tungsten deposits[END_REF] and, so, there is a significant non-magmatic contribution to the ore-forming system (e.g., [START_REF] Polya | Evolution and paragenetic context of low δD hydrothermal fluids from the Panasqueira W-Sn deposit, Portugal: new evidence from microthermometric, stable isotope, noble gas and halogen analyses of primary fluid inclusions[END_REF].

Other Variscan W deposits show a closer spatial association with granitic rocks. This is the case of Cligga Head (Cornwall, UK) where wolframite is deposited together with cassiterite, arsenopyrite and other sulfides in a dense network of intragranitic veins [START_REF] Charoy | Greisenisation, minéralisation et fluides associés à Cligga Head, Cornwall (sud-ouest de l'Angleterre)[END_REF]. Wolframite is ferberite-rich (H/F = 35-40%, Table 2, [START_REF] Charoy | Greisenisation, minéralisation et fluides associés à Cligga Head, Cornwall (sud-ouest de l'Angleterre)[END_REF]Michaud, unpublished) and, although the mineralizing fluids have been interpreted as magmatic [START_REF] Smith | Fluid inclusion and stable isotope constraints the genesis of the Cligga Head Sn-W deposit, SW England[END_REF], W deposition at Cligga Head is associated with hydrothermal alteration around the veins indicating fluid-granite equilibration at subsolidus temperatures (300-400°C, [START_REF] Charoy | Greisenisation, minéralisation et fluides associés à Cligga Head, Cornwall (sud-ouest de l'Angleterre)[END_REF]. At Cínovec (Bohemia, DE and CZ), wolframite occurs disseminated both in greisens and in quartz veins within the host RMG. H/F ratios are intermediate (between 45 and 63%, Table 2 and Fig. 1, [START_REF] Breiter | Assessment of magmatic vs. metasomatic processes in rare-metal granites: A case study of the Cínovec/Zinnwald Sn-W-Li deposit[END_REF] consistent with the presence of subequal concentrations of Fe and Mn in fluid inclusions [START_REF] Korges | Depressurization and boiling of a single magmatic fluid as a mechanism for tin-tungsten deposit formation[END_REF]. Sn-W deposition is associated either with greisenisation of the granite (i.e., with fluid-rock reaction at subsolidus temperatures) or with boiling of the magmatic fluid [START_REF] Korges | Depressurization and boiling of a single magmatic fluid as a mechanism for tin-tungsten deposit formation[END_REF]. Last, deposits from Cornwall (UK) demonstrate a systematic increase of H/F ratios from pegmatites and early feldspathic veins (H/F = 47-58), greisens (H/F = 40-47) to quartz lodes (H/F = 26-40, Table 2 and Fig. 1, [START_REF] Hosking | The significance of the variable iron content of the "Wolframites" of the west of England[END_REF]. In the South Crofty Mine, the "pegmatitic vein swarms", mainly composed of quartz, feldspar and other accessory minerals, contain wolframites of intermediate compositions (H/F = 37-60, [START_REF] Taylor | Manganese-iron ratios in wolframite, South Crofty Mine, Cornwall[END_REF].

In some deposits, several wolframite crystallization stages can be identified, indicating distinct deposition mechanisms. For example, the La Bosse stockwork contains a late generation of hübneritic wolframite (H/F = 44-79%, Table 2 and Fig. 1) developed in minor amounts on early ferberites (H/F = 21-32%, see above) and attributed to the influence of the late Beauvoir RMG [START_REF] Aïssa | Géologie profonde de la Franche Echassières : le forage scientifique. Une clé pour la compréhension des mécanismes magmatiques et hydrothermaux associées aux granites à métaux rares[END_REF]. Isolated hübnerite-bearing veins occur north of the La Bosse deposit (e.g., Le Mazet and Gilbert veins, Table 2) and they have been interpreted as footprints of another Beauvoir-type RMG at depth [START_REF] Aubert | Les coupoles granitiques de Montebras et d'Échassières:(Massif Central français) et la genèse de leurs minéralisations en étain, lithium, tungstène et béryllium[END_REF]. In the Portuguese CIZ, the Vale Das Gatas deposit consists of quartz veins cutting the contact between SGC country rocks and a muscovite>biotite granite [START_REF] Neiva | Geochemistry of cassiterite and wolframite from tin and tungsten quartz veins in Portugal[END_REF]. The veins are characterized by at least two ore deposition stages, the first involving cassiterite and wolframite and the second only wolframite. Wolframites are mostly hübnerites (H/F = 89-96%) except rims that range from Fe-rich hübnerite to ferberite (H/F = 18-54%, Table 2 and Fig. 1).

Wolframite deposition mechanism

None of the wolframite deposition mechanisms mentioned above is directly applicable to the AMD since none explains the precipitation of hübnerite directly from a magmatic fluid. At Argemela, the W-bearing intragranitic veins show no hydrothermal interaction aureole at their margins. The earliest veins exhibit "pegmatitic" textures and crystallize a granitic mineral assemblage [START_REF] Michaud | From magmatic to hydrothermal Sn-Li-Nb-Ta-W mineralization: the Argemela District (Portugal)[END_REF]. Quartz in the granite and veins both have high δO 18 (respectively 13.7-14.5 and 14.5-15.9 ‰). Therefore, the intragranitic veins formed from fluids exsolved from the crystallizing granite and fluid-granite equilibration was originally established at relatively elevated temperatures. Using solidus temperatures for RMG (Pichavant et al., 1987a) and oxygen isotope thermometry on quartz and muscovite from the granite, magmatic fluids were originally at equilibrium with the granite at temperatures of 500-550°C. Their subsequent evolution within the intragranitic veins most probably followed a fluid-buffered instead of a rock-buffered path (e.g., [START_REF] Heinrich | The chemistry of hydrothermal tin (-tungsten) ore deposition[END_REF]. There is no evidence for postmagmatic hydrothermal interaction processes such as greisenisation around the veins. Quartz from the relatively late wolframite-bearing veins have the same stable isotope signature (δO 18 : 15.5-15.9 ‰ ; δD: -25 to -31 ‰) than quartz from the early veins (δO 18 : 14.5-15 ‰ ; δD: -21 to -35 ‰). Furthermore, the surrounding Beira schists are W-depleted and cannot be the source of W [START_REF] Charoy | Multistage growth of a rare-element, volatile-rich microgranite at Argemela (Portugal)[END_REF]. Although the intragranitic vein system was not closed to external non-magmatic contributions (as suggested by the mineralogical data summarized above), wolframite deposition in the Argemela intragranitic veins resulted from cooling of the igneous fluid in isolation from the enclosing rock, in contrast with mechanisms proposed for W deposition in other contexts (e.g., [START_REF] Lecumberri-Sanchez | Fluid-rock interaction is decisive for the formation of tungsten deposits[END_REF].

Wolframite deposition necessitates a source of Fe and Mn, in addition to W. [START_REF] Amossé | Physicochemical study of the hubnerite-ferberite (MnWO4-FeWO4) zonal distribution in wolframite (Mn X Fe (1-X) WO 4 ) deposits[END_REF] demonstrated that the composition of the deposited wolframite depends essentially on two factors, temperature and aMn/aFe of the ore fluid. The Mn/Fe ratio is always higher (and so the H/F ratio higher) in the deposited wolframite than in the equilibrium fluid. Because lowering temperature slightly increases Mn/Fe in wolframite, [START_REF] Amossé | Physicochemical study of the hubnerite-ferberite (MnWO4-FeWO4) zonal distribution in wolframite (Mn X Fe (1-X) WO 4 ) deposits[END_REF] emphasized the importance of the Mn/Fe ratio of the mineralizing fluid as the main factor controlling wolframite composition. At Argemela, the Mn/Fe of the ore fluid is defined by fluid/melt partitioning for Mn and Fe. Both elements partition similarly between melt and fluid as shown by their nearly identical ratios in coexisting silicate melt and fluid inclusions from Ehrenfriedersdorf [START_REF] Zajacz | Determination of fluid/melt partition coefficients by LA-ICPMS analysis of co-existing fluid and silicate melt inclusions: Controls on element partitioning[END_REF]. Thus, the Mn/Fe of the magmatic fluid that deposited hübnerite at Argemela was essentially that of the crystallizing granitic melt. Although this Mn/Fe ratio is unknown, Variscan RMG (e.g., Beauvoir, [START_REF] Raimbault | Geochemical evidence for a multistage magmatic genesis of Ta-Sn-Li mineralization in the granite at Beauvoir, French Massif Central[END_REF] and pegmatites (Chèdeville, [START_REF] Raimbault | Composition of complex lepidolite-type granitic pegmatites and of constituent columbitetantalite, Chedeville, Massif Central, France[END_REF], and equivalent volcanic products (Macusani, Pichavant et al., 1987b), demonstrate a general increase of Mn/Fe during magmatic differentiation. The whole-rock data for the Argemela granite show such an increase [START_REF] Charoy | Multistage growth of a rare-element, volatile-rich microgranite at Argemela (Portugal)[END_REF], also marked in the enrichment of disseminated columbo-tantalites in MnTa 2 O 6 [START_REF] Michaud | From magmatic to hydrothermal Sn-Li-Nb-Ta-W mineralization: the Argemela District (Portugal)[END_REF], see also [START_REF] Raimbault | Composition of complex lepidolite-type granitic pegmatites and of constituent columbitetantalite, Chedeville, Massif Central, France[END_REF][START_REF] Linnen | Granite-related rare-element deposits and experimental constraints on Ta-Nb-W-Sn-Zr-Hf mineralization[END_REF]. Therefore, the magmatic evolution of the Argemela RMG led to the exsolution of relatively high Mn/Fe igneous fluids. We suggest that the deposition of hübnerite in perigranitic magmatic-hydrothermal systems is the signature of metals and fluids originating from highly evolved peraluminous crustal magmas similar to the Argemela granite.

W is known to be transported in hydrothermal fluids as an oxyanion, WO 4 2-, complexed with Na+, K+, or H+ [START_REF] Wood | The Hydrothermal Geochemistry of Tungsten in Granitoid Environments: I. Relative Solubilities of Ferberite and Scheelite as a Function of T, P, pH, and mNaCl[END_REF][START_REF] Heinrich | The chemistry of hydrothermal tin (-tungsten) ore deposition[END_REF]. Substitution of Na for Li does not change wolframite solubilities in granitic melts (all the other parameters being equal, [START_REF] Linnen | The solubility of Nb-Ta-Zr-Hf-W in granitic melts with Li and Li + F; constraints for mineralization in rare metal granites and pegmatites[END_REF] and, so, Li is a W complexing agent with an efficiency comparable to Na. Thus, destabilization of Li-W complexes as a result of precipitation of a Li-bearing phase will cause wolframite deposition. Montebrasite (LiAlPO 4 (OH, F)) is a systematic phase in the Argemela intragranitic veins where it is accompanied by K-feldspar (whose crystallization removes K from the fluid).

Muscovite never occurs as a mineral in the veins suggesting that the pH of the fluid was never strongly acid.

Compositions of montebrasites (Michaud et al., 2018b) in K-feldspar-bearing veins show no variations in F content [START_REF] Loh | Synthesis and fluorine-hydroxyl exchange in the amblygonite series[END_REF] and so deposition was probably not driven by a drastic increase of the fluid pH. However, montebrasites from montebrasite-wolframite veins can have F contents nearly twice as elevated as montebrasites in K-feldspar-bearing veins, which indicates a significant pH decrease (thus counteracting wolframite precipitation). Therefore, a decrease in temperature is probably the main factor controlling wolframite deposition [START_REF] Heinrich | The chemistry of hydrothermal tin (-tungsten) ore deposition[END_REF].

One element specific of the Argemela magmatic-hydrothermal system is phosphorus. The Argemela granite, with whole-rock concentrations of up to ~2 wt% P2O5 [START_REF] Charoy | Multistage growth of a rare-element, volatile-rich microgranite at Argemela (Portugal)[END_REF], is a good example of a peraluminous high phosphorus (PHP) rare-metal granite [START_REF] Linnen | Granite-related rare-element deposits and experimental constraints on Ta-Nb-W-Sn-Zr-Hf mineralization[END_REF]. Montebrasite crystallizes as a magmatic phase in the Argemela granite and is also a key phase in intragranitic veins [START_REF] Charoy | Multistage growth of a rare-element, volatile-rich microgranite at Argemela (Portugal)[END_REF].

Phosphorus increases wolframite solubility in peraluminous granitic melts [START_REF] Che | Tungsten solubility in evolved granitic melts: An evaluation of magmatic wolframite[END_REF]. Its influence is however indirect, through the formation of AlPO4 species which increases the availability of Na, K and Li to complex with W. Therefore, removal of such species by phosphate crystallization will cause wolframite deposition. The overall mechanism of wolframite precipitation in the Argemela intragranitic vein system is summarized with the following reaction: An equivalent reaction can be written for the fluorine Li phosphate end-member (amblygonite). Equation ( 1) is probably driven to the right by a decrease in temperature, i.e., by cooling of the intragranitic vein fluids.

Discussion and implications for exploration

The model proposed for Argemela emphasizes the differences with other W deposition environments, and has direct implications for exploration. The model applies to hübnerite deposition from a magmatic fluid. Wolframite composition reflects enrichment in Mn relative to Fe in the crystallizing granite which is the source of both W and fluids. Wolframite deposition calls upon a specific chemical environment (presence of both Li and P) in a magmatic-hydrothermal system that evolves under fluid-buffered conditions. Li, together with F and P, is systematically enriched in Variscan granitoïds from the CIZ (Roda-Robles et al., 2018 and references therein).

Such anomalies reflect the specific chemistry (Al-, P-, F-and Li-rich, Ca-poor) of the metasedimentary rock sources of granitic magmas [START_REF] Villaseca | Petrogenetic relationships between Variscan granitoids and Li-(F-P)-rich aplite-pegmatites in the Central Iberian Zone: Geological and geochemical constraints and implications for other regions from the European Variscides[END_REF]. The Li-and P-rich geochemical signature exhibited by the Argemela granite is a reflection of such specific crustal sources. Therefore, the W deposition mechanism proposed at Argemela also probably applies to other similar W occurrences in the CIZ.

The model for the AMD differs from others where a magmatic source for W is assumed but fluid interactions with country rocks are considered necessary for W deposition [START_REF] Lecumberri-Sanchez | Fluid-rock interaction is decisive for the formation of tungsten deposits[END_REF]. One first-order difference between the two types of models is that, in the former case, veins are intragranitic (Argemela) and, in the latter, within country rocks (Panasqueira). Therefore, each corresponds to a specific type of W deposition environment, as reflected by the composition of wolframite, hübnerite in the former and ferberite in the latter (Fig. 3). In the case of intragranitic vein systems, a magmatic origin for the ore fluid has been generally proposed (e.g., [START_REF] Korges | Depressurization and boiling of a single magmatic fluid as a mechanism for tin-tungsten deposit formation[END_REF][START_REF] Smith | Fluid inclusion and stable isotope constraints the genesis of the Cligga Head Sn-W deposit, SW England[END_REF]. However, one major difference between Argemela and these systems (Cligga Head and Cínovec respectively) is the temperature and mechanism of fluid-rock equilibration. At Argemela, the composition of the ore fluid is defined by fluid/melt equilibrium at near-solidus temperatures; the magmatic fluid evolves along a fluid-buffered path. At Cligga Head and Cínovec, the fluid is buffered by fluid-rock interaction at subsolidus temperatures (greisenisation). Therefore, Argemela can be viewed as a high temperature variant of the Cligga Head and Cínovec cases. Preservation of the high temperature intragranitic vein system at Argemela was probably favored by the small volume of the granite intrusion which prevented the development of a large long-lived hydrothermal system involving adjacent country rocks and fluids.

Our results suggest that the H/F ratio can be used as an indicator of contrasted W deposition environments in perigranitic ore-forming systems (Fig. 3). Hübnerite crystallization is associated with the following key features for exploration: (1) source of W resides in nearby igneous rocks, (2) highly evolved granitic magmas are present, (3) ore fluids have a magmatic origin and (4) a "high temperature" magmatic-hydrothermal stage is preserved.

Overall, hübnerite can be viewed as a relatively simple indicator of a strong magmatic control on W deposition.

The fact that, in W deposits, ferberite compositions are generally more common than hübnerite could be the consequence of early hübnerite recrystallizing to lower H/F ratios as observed for example in Fig. 2. Although, at Argemela hübnerite deposition is the consequence of a high Mn/Fe in the ore fluid, hübnerite crystallization progressively depletes the fluid in Mn [START_REF] Amossé | Physicochemical study of the hubnerite-ferberite (MnWO4-FeWO4) zonal distribution in wolframite (Mn X Fe (1-X) WO 4 ) deposits[END_REF] and, therefore, precipitation of wolframites with progressively lower H/F ratios is expected (Fig. 3). However, large variations of the H/F ratio within a given deposit (for example at La Bosse, Vale Das Gatas, Argemela, Fig. 3) are the mark of a complex multi-phase deposition history involving ore fluids with contrasted Fe/Mn and different, superimposed, wolframite deposition environments. Wolframite can recrystallize to lower H/F ratios (WI in the Argemela country rock vein, Fig. 3) when the magmatic-hydrothermal system becomes open to fluids derived from Fe-, S-rich metamorphic country rocks [START_REF] Lecumberri-Sanchez | Fluid-rock interaction is decisive for the formation of tungsten deposits[END_REF]. Alternatively, overprinting of early low H/F ratios in country rock veins by magmatic fluids derived from a crystallizing RMG can be also observed in the FMC (La Bosse, Fig. 3).

Hübnerite is a tracer of proximal highly evolved igneous bodies such as RMG and/or pegmatite with potential for commodities such as quartz and feldspar, kaolin and disseminated Li-Sn-Nb-Ta ores (e.g., [START_REF] Aubert | Les coupoles granitiques de Montebras et d'Échassières:(Massif Central français) et la genèse de leurs minéralisations en étain, lithium, tungstène et béryllium[END_REF]Fig.3), in addition to W. Wolframite with intermediate H/F ratios is characteristic of medium temperature rock-buffered hydrothermal systems with a strong magmatic signature, both for the fluids and the metals (Fig. 3, [START_REF] Smith | Fluid inclusion and stable isotope constraints the genesis of the Cligga Head Sn-W deposit, SW England[END_REF][START_REF] Korges | Depressurization and boiling of a single magmatic fluid as a mechanism for tin-tungsten deposit formation[END_REF]. These compositions are typical of greisen-type (both vein and massive) deposits where mostly Li-Sn-Nb-Ta are economically interesting and of intragranitic vein type deposits that follow a rock-buffered path. Lastly, H/F ratios < 40 and down to 0 indicate W deposition environments where a significant non-magmatic contribution to the ore fluid and W deposition is necessary (Fig. 3, [START_REF] Lecumberri-Sanchez | Fluid-rock interaction is decisive for the formation of tungsten deposits[END_REF].

Ferberites showing a low H/F ratio are tracers of W-bearing veins/stockwork deposits of great economic importance.

The comparison can be extended to deposits worldwide (Fig. 3). Fractionated Li-bearing pegmatites from Mount Begbie (Bristish Columbia; [START_REF] Dixon | Mineralogy and geochemistry of pegmatites on Mount Begbie, British Columbia[END_REF] exhibit scarce hübnerites with a high H/F ratio, similar to Argemela (Fig. 3). Intragranitic veins of the Piaotang and Xihuashan W-deposits [START_REF] Zhang | In-situ LA-ICP-MS trace element analyses of scheelite and wolframite: Contraints on the genesis of veinlet-disseminated and vein-type tungsten depostis, South China[END_REF] are characterized by an alteration halo at vein margins which indicates a rock-buffered path (i.e., granite-fluid interaction). These vein systems behave similarly to Cligga Head and both show comparable intermediate H/F ratios. In comparison, Piaotang country rock veins show a lower H/F ratio [START_REF] Zhang | In-situ LA-ICP-MS trace element analyses of scheelite and wolframite: Contraints on the genesis of veinlet-disseminated and vein-type tungsten depostis, South China[END_REF], most probably a consequence of interaction with these rocks. The Bugarama, Nyakabingo and Gifurwe vein-type deposits from the "tungsten belt" of Rwanda are all hosted in country rocks and wolframite shows very low H/F ratio [START_REF] Goldmann | Mineralogy and Trace Element Chemistry of Ferberite/Reinite from Tungsten Deposits in Central Rwanda[END_REF], comparable to the Panasqueira compositions. Thus, the H/F ratio appears as a reliable exploration tool worldwide despite differences in ages and geodynamic contexts between deposits.

Conclusions

Intragranitic veins in the Argemela Mineralized District remarkably illustrate the case of wolframite precipitation directly from relatively "high temperature" magmatic fluids. High H/F ratios wolframites (hübnerites) are characteristic of this specific W deposition environment which contrasts with other mechanisms that lead to intermediate to low H/F ratios. It is proposed that the H/F ratio can be used as an indicator of wolframite deposition style for Variscan deposits (and maybe beyond) and thus as a relevant exploration tool. Lastly, it should be emphasized that the genesis of Variscan W deposits does not only result from locally favorable country rocks, as for Panasqueira (e.g., [START_REF] Lecumberri-Sanchez | Fluid-rock interaction is decisive for the formation of tungsten deposits[END_REF]. Chemically specific source rocks and granite magma types are also necessary.

Figure captions

Fig. 1: H/F ratio vs atomic W concentration for wolframites from representative Variscan deposits. Argemela data are from this study (Table 1). For other deposits see references in Table 2.

(1.5 fitting image) 2 for Variscan deposits and see [START_REF] Zhang | In-situ LA-ICP-MS trace element analyses of scheelite and wolframite: Contraints on the genesis of veinlet-disseminated and vein-type tungsten depostis, South China[END_REF], [START_REF] Dixon | Mineralogy and geochemistry of pegmatites on Mount Begbie, British Columbia[END_REF] and [START_REF] Goldmann | Mineralogy and Trace Element Chemistry of Ferberite/Reinite from Tungsten Deposits in Central Rwanda[END_REF] for worldwide deposits. Three main H/F domains corresponding to different wolframite deposition mechanisms are identified: (A) for H/F ratio > 60, wolframite precipitates from a Mn-rich magmatic fluid evolving under a fluid-buffered path [START_REF] Heinrich | The chemistry of hydrothermal tin (-tungsten) ore deposition[END_REF]; (B) for H/F ratio between 40 and 60, wolframite precipitates from a fluid buffered by granite/fluid interactions [START_REF] Smith | Fluid inclusion and stable isotope constraints the genesis of the Cligga Head Sn-W deposit, SW England[END_REF]; (C) for H/F ratio < 40, wolframite precipitates from fluid carrying a significant non magmatic signature derived from country rocks (Lecumberri-Sanchez et al., 2017). Three main mechanisms controlling variations of the H/F ratio in a single deposit are illustrated: (1) step increase of the H/F ratio resulting from addition of Fe through fluid-rock interactions [START_REF] Lecumberri-Sanchez | Fluid-rock interaction is decisive for the formation of tungsten deposits[END_REF]; (2) Progressive increase of the H/F ratio due to Mn removal during hübnerite crystallization [START_REF] Amossé | Physicochemical study of the hubnerite-ferberite (MnWO4-FeWO4) zonal distribution in wolframite (Mn X Fe (1-X) WO 4 ) deposits[END_REF] and (3) step decrease of the H/F ratio due to influx of Mn-rich magmatic fluid. Deposit type and commodity associated with the three domains are detailed. Eng-Leu: Enguialès-Leucamp; FMC: French Massif Central; "peg": pegmatitic veins [START_REF] Hosking | The significance of the variable iron content of the "Wolframites" of the west of England[END_REF].

(2 column fitting image) Highlights -Mn-rich wolframite precipitates from a Mn-rich magmatic fluid evolving through a fluid-buffered path.

-The H/F ratio can be used as an indicator of wolframite deposition mechanisms and as a tracer of the fluid source. -The H/F ratio is a simple guide tool for ore deposits exploration.

Li

  NaWO 4 + (Mn, Fe) 2+ + AlPO 4 + H 2 O = LiAlPO 4 OH + (Mn, Fe)WO 4 + Na + + H + (1) fluid fluid fluid fluid montebrasite wolframite fluid fluid

Fig. 2 :

 2 Fig. 2: Patchy texture of Argemela WI (schist-hosted vein) illustrating evolution of the H/F ratio between early and late crystallization stages. EDS map of (a) Fe and (b) Mn. (single column fitting image) Fig. 3: Environments and mechanisms of W deposition as revealed by the H/F ratio for several representative Variscan deposits and extension to worldwide deposits. Source of data in Table2for Variscan deposits

Acknowledgements

The authors acknowledge P. Moutela and P. Ferraz for field access and I. Di Carlo for assistance with the electron microprobe and SEM. This work was supported by the Labex VOLTAIRE (ANR-10-LABX-100-01), the ERAMIN project NewOres and the ANR project VARPEG (ANR-15-CE01-0001). We acknowledge P.

Lecumberri-Sanchez for her constructive remarks.