Junjie Lai
email: junjie.lai@inria.fr

André Seznec
email: andre.seznec@inria.fr

Performance Upper Bound Analysis and Optimization of SGEMM on Fermi and Kepler GPUs

Keywords: C.1.4 [Processor Architectures]: Parallel Architectures, C.4 [Performance of Systems]: Modeling techniques, F.2.1 [Numerical Algorithms and Problems]: Computations on matrices Measurement, Performance Kepler GPU, Fermi GPU, SGEMM, CUDA, Performance Upper Bound Analysis

In this paper, we present an approach to estimate GPU applications' performance upper bound based on algorithm analysis and assembly code level benchmarking. As an example, we analyze the potential peak performance of SGEMM (Single-precision General Matrix Multiply) on Fermi (GF110) and Kepler (GK104) GPUs. We try to answer the question of how much optimization space is left for SGEMM and why. According to our analysis, the nature of Fermi (Kepler) instruction set and the limited issue throughput of the schedulers are the main limitation factors for SGEMM to approach the theoretical peak performance. The estimated upper-bound peak performance of SGEMM is around 82.5% of the theoretical peak performance on GTX580 Fermi GPU and 57.6% on GTX680 Kepler GPU. Guided by this analysis and using the native assembly language, on average, our SGEMM implementations achieve about 5% better performance than CUBLAS in CUDA 4.1 SDK for large matrices on GTX580. The achieved performance is around 90% of the estimated upper-bound performance of SGEMM on GTX580. On GTX680, the best performance we achieve is around 77.3% of the estimated performance upper bound. We also describe how to use native assembly language directly in the CUDA runtime source code.

Introduction

Nowadays, vendors are putting more and more computing cores per die and the GPU architecture is evolving rapidly. Fermi and Kepler GPUs are the present generations of NVIDIA GPUs [START_REF]Fermi Whitepaper[END_REF][START_REF]GTX680 Whitepaper[END_REF]. There are many studies about optimizing specific kernels on Fermi GPU. However, since the architecture is changing with each generation, we may need to repeat the optimization work again very soon. Unfortunately, no practical performance upper bound evaluation is available to the developers. In practice, developers apply several optimization techniques based on the analysis to the algorithm or serial code, and their expert experience. Then developers may modify the optimizations with feedback provided by tools like NVIDIA Visual Profiler [START_REF] Nvidia | Visual Profiler[END_REF]. However, they can not be sure how far the obtained performance is from the best achievable performance. In this paper, we present an approach to project performance upper bound using algorithm analysis and assembly code level benchmarking.

To understand GPU performance result, there exist many works about how to project/predict CUDA applications' performance using analytical or simulation methods. Meng et al. [START_REF] Meng | Grophecy: Gpu performance projection from cpu code skeletons[END_REF] proposed a GPU performance projection framework based on annotated code skeletons. Hong and Kim [START_REF] Hong | An analytical model for a gpu architecture with memory-level and thread-level parallelism awareness[END_REF] introduced the MWP-CWP model to predict CUDA application performance using PTX code. Bakhoda et al. [START_REF] Bakhoda | Analyzing cuda workloads using a detailed gpu simulator[END_REF] developed a detailed GPU simulator and the simulator also uses the PTX code as input. Recently, Sim et al. [START_REF] Sim | A performance analysis framework for identifying potential benefits in gpgpu applications[END_REF] extended the MWP-CWP model and utilize the assembly code of CUDA kernel to predict performance. The quantitative GPU performance model proposed by Zhang and Owens [START_REF] Zhang | A quantitative performance analysis model for gpu architectures[END_REF] is also based on the native assembly code. Since very little information about the underlying GPU architecture is disclosed, it becomes very unlikely to build accurate simulators for each new GPU generation. Luckily, the results [START_REF] Hong | An analytical model for a gpu architecture with memory-level and thread-level parallelism awareness[END_REF][START_REF] Meng | Grophecy: Gpu performance projection from cpu code skeletons[END_REF][START_REF] Sim | A performance analysis framework for identifying potential benefits in gpgpu applications[END_REF][START_REF] Zhang | A quantitative performance analysis model for gpu architectures[END_REF] show that we can have very good approximation of GPU performance using analytical approaches. However existing GPU performance models all rely on certain level of an application's implementation (C++ code, PTX code, assembly code. . . ) and do not answer the question of how good the current optimized version is and whether further optimiza-tion effort is worthwhile or not. Different from existing GPU performance models, our approach does not project the possible performance from certain implementations, but the performance upper bound that an application cannot exceed.

Researchers are also interested in the outcome of different optimization combinations on GPUs. The roofline model [START_REF] Williams | Roofline: an insightful visual performance model for multicore architectures[END_REF] is well known for estimating the optimization effects. The recent work by Sim et al. [START_REF] Sim | A performance analysis framework for identifying potential benefits in gpgpu applications[END_REF] studied the effects of different optimization techniques on GPUs using the similar approach as the roofline model. However, the chosen optimizations normally rely on the initial code version and different optimizations are likely to have complex impacts on each other. Our approach tackles the problem from the opposite angle as the roofline method. We first assume an optimistic situation on GPUs (no shared memory bank conflict, global memory accesses are all coalesced, all the auxiliary operations like address calculations are neglected, etc.). Then we try to predict a performance upper bound when mapping an application on the GPU based on the constraints introduced by the architecture, the instruction set and the application itself, or the constraints that we are not able to eliminate using optimization techniques. With a tight performance upper bound of an application, we have an evaluation on how much optimization space is left and can decide the optimization effort. Also, with the analysis, we can understand which parameters are critical to the performance and have more insights into the performance result. Hence, with these knowledge, it would be easier for the community to move to the new architecture.

As an example, we analyze the potential peak performance of SGEMM (Single-precision General Matrix Multiply) on Fermi (GF110) and Kepler (GK104) GPUs. GEMM1 operation is essential for Level 3 BLAS (Basic Linear Algebra Subprograms) [2] routines and generally represents the practical best performance of a computer system. If we compare the performance of SGEMM from CUBLAS with the theoretical peak performance, on Fermi, it achieves around 70% and on Kepler, only around 42% of the theoretical peak performance. There are already some articles about optimizing GEMM kernels on Fermi GPU [START_REF] Nath | An improved magma gemm for fermi gpus[END_REF] [START_REF] Tan | Fast implementation of dgemm on fermi gpu[END_REF], and an auto-tuning framework has also been presented [START_REF] Kurzak | Autotuning gemm for the fermi gpu[END_REF]. In this research, the focus is to answer the question of how much optimization space is left for SGEMM and why. We also show that the analysis can help optimization efforts since it uncovers critical parameters. Only single precision SGEMM is evaluated, since we could only access the GTX580 Fermi and the GTX680 Kepler Geforce cards, which have much poorer double precision performance than Tesla products. It is not really worth the effort to study the DGEMM performance on Geforce GPU. In November 2012, NVIDIA has announced the new Tesla K20X Kepler GPU (GK110) and the documented SGEMM efficiency is around 73% of the theoretical peak performance [START_REF]NVIDIA Tesla K20/K20X GPU Accelerators Application Performance Technical Brief[END_REF]. The K20X Kepler GPU (GK110) architecture is different from the GTX680 (GK104) and uses a different instruction set (each thread can utilize maximum 255 registers on the new architecture while the limit is 63 on GTX680 GPU). With a Tesla GPU card, it should not be difficult to extend the analysis to SGEMM and DGEMM on the Tesla GPU using our approach.

This paper is organized as follows: Section 2 introduces our target GPU architectures. Section 3 introduces our assembly level benchmarking approach. Section 4 presents our analysis for performance upper bound of SGEMM on Fermi and Kepler GPUs. In Section 5 assembly code level optimization methods and performance result of SGEMM are presented. Section 6 concludes this study.

GPU Architecture Characteristics

NVIDIA GPUs are composed of a cluster of independent SMs (Streaming Multiprocessor). Each SM includes several SPs (Streaming Processor), which is the basic computing component, scheduler, SFUs (Special Functional Unit), LD/ST (Load/Store) Units, register file and a unified shared memory/L1 cache. The SPs inside one SM are similar to a lane of SIMD engines and they share the memory resource of the SM like the registers and shared memory.

A comparison of the three generations of NVIDIA GPUs is illustrated in Table 1. From GT200 to Kepler GPU, the number of SPs increases dramatically, from 240 (GTX280, 65nm) to 1536 (GTX680, 28nm) [START_REF]Fermi Whitepaper[END_REF][START_REF]GTX680 Whitepaper[END_REF]. However, if we consider the memory resource (registers and shared memory) per SP, the on-die storage per SP actually decreases. Previous generations have two clock domains in the SM, the core clock for the scheduler and the shader clock for the SPs. The shader clock is roughly twice the speed of the core clock. On Kepler (GK104) GPU, shader clock no longer exists, the functional units with SMs run at the same core clock. However, to compare the different generations more easily, we still use the term shader clock on Kepler GPU. In the rest of this paper, all throughput data is calculated with the shader clock.

A typical CUDA [START_REF]Nvidia cuda c programming guide 4[END_REF] program normally creates thousands of threads to hide memory access latency or math pipeline latency. The threads are grouped into 1D to 3D blocks, and further into 1D or 2D grids. Each block is assigned to one SM at execution time. The warp is the basic execution and scheduling unit of a SM, and is composed of 32 threads. We define a warp instruction as the same instruction shared by all threads in the same warp, and a thread instruction as the instruction executed by one thread. So a warp instruction launches 32 operations or consists of 32 thread instructions. On the SM, only a limited set of threads can run concurrently (active threads). On one hand, the increased SPs require more active threads to hide latency. On the other hand, the register and the shared memory resource For Fermi (and Kepler GK104) instruction set, there is a hard limit of maximum 63 registers per thread (for GT200 generation the limit is 127 registers per thread) since in the instruction encoding, only 6 bits are left for one register.

CUDA Programming with Native Assembly Code

For performing this study, we have to develop some software components and reverse engineer many characteristics of the hardware. We use the GPU assembly code directly with an assembly tool Asfermi [1]. Asfermi was first developed to work on Fermi GPU. We patch Asfermi to support Kepler GPU (GK104) and manage to use native assembly language directly in the CUDA runtime source code.

Using Native Assembly Code in CUDA Runtime API Source Code

Programming in assembly code on NVIDIA GPUs is not publicly supported by the company. However, our analysis is requiring such programming. With an assembly tool for Fermi GPU called Asfermi [1] and a little hacking into the CUDA programming compiling stages, we manage to use hand-tuned GPU assembly code in CUDA projects using CUDA runtime APIs .

There are several advantages of using assembly code or native machine code directly instead of using high level languages like C++. First, we can carefully control the register allocation since the register resource per thread is very limited and sometimes the compiler may spill many registers for programs utilizing much register resource per thread like SGEMM. Second, the instruction order can be carefully designed to better prefetch data from global memory and mix different instruction types to get better throughput. Third, SIMD-like instructions (LDS.64 or LDS.128) could be used intentionally to reduce the instruction number. Also, we can control the exact behavior of the machine code. For example, the compiler might choose to use wider load instructions (LDS.64 or LDS.128) based on the data alignment in shared memory. However, using wide load instructions does not always benefit the performance (Section 4.1).

Kepler GPU Binary File Format

Asfermi was first developed to work on Fermi GPU. We patched Asfermi to support CUDA sm 30 (GK104 Kepler GPU). However, although the CUDA program using assembly code kernel can still run correctly on Kepler GPU, the performance is very poor. The reason is that new control information is embedded into the CUDA binary file to help processor scheduling. According to the GTX680 white paper [START_REF]GTX680 Whitepaper[END_REF], the compiler places the scheduling information along with the actual instructions in the CUDA binary file.

The scheduling information (we call it control notation) is placed before each group of 7 instructions and has the format of 0xXXXXXXX7 0x2XXXXXXX. 0x7 and 0x2 are identifiers and the rest of the notation is separated into 7 fields and associated with each following instruction. It is similar to the explicit-dependence lookahead used in Tera computer system [START_REF] Alverson | The tera computer system[END_REF]. Unfortunately, NVIDIA does not disclose the encoding of the control notation. We do not know how to generate the control notation exactly as the nvcc compiler. In our implementation of SGEMM on Kepler GPU, as a compromise, we use the same control notation for same kind of instructions and try to find the best combination of those notations for major instruction types. However, our decryption of the notations is still not enough.

Math Instruction Throughput on Kepler GPU

Understanding and modeling the behavior of math instructions on Kepler GPU is a major difficulty. We use two approaches to test the throughput of math instructions. First, a kernel is written in C++ code and compiled into binary with control notations embedded by nvcc. Second, a kernel is written in assembly code directly and the controlling Instruction FFMA performs single precision fused multiplyadd operation (FFMA RA, RB, RC, RD performs the operation RA := RB * RC + RD). With first approach, the instruction throughput of FFMA R9, R8, R9, R5 is measured as 129.2 operations per shader cycle 2 . With the second approach and the control notation of 0x25, the throughput is 132.0 operations per shader cycle (The actual shader clock cannot be obtained during execution. All throughput data is calculated by boost clock of 1058MHz [START_REF]GTX680 Whitepaper[END_REF]). Some math instructions' throughput is illustrated in Table 2 measured with the second approach. In these cases, the scheduling function units on one SM can only issue about maximum 132 thread instructions per shader cycle, which is much lower than the SP's processing throughput (192 thread instructions per shader cycle). If some of the three source registers are the same (like FFMA RA, RB, RB, RA), with some carefully designed code structures, the FFMA throughput can approach around 178 thread instructions per shader cycle. However, considering 'useful' FFMA's throughput, that is (FFMA RA, RB, RC, RA), the maximum single precision performance for many applications like SGEMM on GTX680 GPU (GK104) cannot exceed around 68.75% (132/192) of the claimed performance (3090GFlops) by NVIDIA.

Our benchmark result also shows that the instruction throughput is related to register indices. According to some other experiments, we speculate that the registers reside on four banks. Take the instruction FFMA RA, RB, RC, RD for instance, if there are two different source registers on the same bank, the throughput drops by 50%, and if all three source registers RB, RC, RD are different registers on the same bank, the throughput is around 33.3% of the best 2 The actual implementation is not 8192 FFMA R9, R8, R9, R5 instructions per thread but 4 independent FFMA instructions like FFMA R9, R8, R9, R5 unrolled by 2048 times.

(R index %8 < 4 && R index %2 == 0), even 1 (R index %8 ≥ 4 && R index %2 == 0), odd 0 (R index %8 < 4 && R index %2 == 1)
, and odd 1(R index %8 ≥ 4 && R index %2 == 1). Since we implement SGEMM with assembly code directly, the register indices have to be carefully chosen to make sure there is no bank conflict. The detailed optimization is illustrated in section 5.4.

Analysis of Potential Peak Performance of SGEMM

The general analysis approach can be similar for all applications while the detailed analysis process may differ from application to application. Our method is applicable for applications which use a few major instruction types and a simple execution path. Many high-performance computing kernels have this characteristic, especially linear algebra routines. Our analysis requires characteristics of the architecture such as register file size, maximum number of registers per thread, shared memory size, instruction throughput for different instruction mix, etc. Those characteristics need to be collected on the real hardware and are independent of the effective application. First, we should analyze the instruction types and percentage of a routine. Second, we should find the critical parameters which affect the different instructions' mixing percentage. Third, we analyze how the instruction throughput changes when we vary these critical parameters. Fourth, we can use the instruction throughput with critical parameters' optimal combination to estimate the performance upper bound. With this approach, not only we can have the performance upper bound estimation, know how much performance gap is left and decide the optimization effort, but we can also understand what parameters are essential to the performance and how to distribute our optimization effort.

For SGEMM, all well-implemented SGEMM kernels actually utilize shared memory on the GPU to reduce the global memory pressure as illustrated in Figure 1. First, data is loaded from global memory to shared memory and then threads within one block can share the loaded data in the shared memory. For Fermi (GF110) and Kepler (GK104) GPUs, arithmetic instructions like FFMA cannot take operands from the shared memory. Since LDS instructions are needed to load data first from shared memory into registers, most of the instructions executed in SGEMM are FFMA and LDS instructions. For instance, in our SGEMM implementation with 1024x1024 matrix size, 80.5% of instructions executed are FFMA instructions and 13.4% are LDS.64 instructions. So essentially, in our analysis, we define a few key parameters and study the instruction throughput mixing FFMA and LDS.X instructions while varying these parameters.

The rest of this section is our analysis of SGEMM's performance upper bound. We show that the analysis can give good insights about how to optimize a specific kernel (SGEMM) and help us to understand the performance result.

Using Wider Load Instructions

To achieve better performance, it is essential to minimize auxiliary instructions' percentage. By auxiliary instructions, we mean non-math instructions, especially LDS instruction. The assembly code for CUDA sm 20 (GF110 Fermi GPU) and sm 30 (GK104 Kepler GPU) provides SIMD-like LDS.64 and LDS.128 instructions to load 64bit and 128bit data from the shared memory. Using wider load instructions can reduce the total number of LDS instructions.

According to our benchmarks, on Fermi GPU, the peak throughput for LDS instruction is 16 32bit-operations per shader clock per SM. Using LDS.64 instructions does not increase the data throughput and the LDS.128 instruction normally leads to 2-way shared memory bank conflict on Fermi GPU. LDS.128 has the throughput of only 2 thread instructions per shader cycle on one SM. On Kepler GPU, the throughput for LDS operation is measured as 33.1 64bit operations per shader clock per SM. Using the 32bit LDS operation actually decreases the data throughput in half comparing with using LDS.64 instructions and properly used LDS.128 instruction does not introduce penalty. Figure 2 illustrates the instruction throughput of mixing FFMA and LDS.X instructions. While gradually increasing the ratio of FFMA instructions to LDS instructions, the overall instruction throughput approaches the FFMA's peak processing throughput. The instruction ratio of FFMA to LDS.X depends on the algorithm parameters such as register blocking size. Apparently, the overall performance does not always benefit from using wider load instructions. However, the compiler might choose to use the wider load instructions based on the data alignment in the shared memory. With the native assembly language, it is possible for us to carefully design the data layout and use the best instruction type.

Register Blocking

As in Table 1, the scheduler of GT200 GPU can issue one warp instruction like FFMA per core cycle and since there are 8 SPs per SM, SPs need 4 shader cycles to process one warp instruction. Apparently, as the issue throughput is higher than the SP's processing throughput, math instructions executed in SPs cannot saturate the scheduler's issue throughput. So the scheduler has some 'free cycles' to issue instructions to other functional units.

On Fermi GPUs, SM are redesigned with 2 warp schedulers and 32 SPs. Each warp scheduler, equipped with one dispatch unit, issues instructions to 16 SPs. With an issue rate of one warp instruction per shader cycle, the schedulers' ability could be fully utilized by 32 SPs. The theoretical peak performance for math instructions comes from the SPs' performance. The percentage of other instructions becomes an issue when there are many auxiliary instructions: there are fewer cycles left for schedulers to issue useful instructions like FFMA. In the worst case, without any register reuse, 2 LDS instructions are needed to fetch data for 1 FFMA instruction in the SGEMM main loop. In that case, only 1/3 of the instructions are floating point operations. Blocking is a well-known technique to better utilize memory hierarchy for scientific programs [START_REF] Lam | The cache performance and optimizations of blocked algorithms[END_REF][START_REF] Mckellar | Organizing matrices and matrix operations for paged memory systems[END_REF]. To increase the percentage of math instructions, register blocking is needed. We illustrate the percent-age of FFMA instructions varying register blocking factors in Figure 3.

If 6-register blocking is used (which is the case of our SGEMM implementation on Fermi GPU), the FF-MA/LDS.X ratios are 3:1, 6:1, and 12:1 if shared memory accesses are implemented with LDS, LDS.64 and LDS.128 respectively. The percentage of FFMA instructions is 75%, 85.7% and 92.3%. On Fermi GPU, the overall instruction throughputs for one SM in these cases are 31.3, 30.4 and 24.5 thread instructions per shader clock. Because using LDS.128 instruction may lead to extra penalties, even if all the accesses to shared memory are implemented with LDS.128, in the best case we can only achieve around 71% ( 24.5 32 * 92.3%) of SMs' single precision floating point performance. Also, in many cases, a lot of padding in shared memory has to be used to get proper data alignment. Apparently, it is not worth the programming effort to mix FFMA with LDS.128 for SGEMM on the Fermi GPU.

Active Threads on SM

Normally, the more active threads one SM executes, the higher performance the GPU can achieve. Since register and shared memory resource is limited per SM, only a limited set of warps can be executed concurrently (T SM ).

TSM * RT ≤ RSM (1)
The registers that each thread can utilize (R T ) is less than or equal to 63 on Fermi and Kepler GPUs (R M ax ). Furthermore, the register budget of the active warps cannot exceed the SM's register amount (R SM ) (Equation 1). Figure 4 illustrates the instruction throughput mixing FFMA and LDS.64 instructions with ratio 6:1 under different number of active threads on one SM. We tested two cases. In the first case ( 'independent' in Figure 4), 6 FFMA and 1 LDS.64 instructions are all independent. In the second case ( 'dependent' in Figure 4), 6 FFMA instructions are dependent on one LDS.64 instruction. The second case is closer to the actual implementation of SGEMM. On Fermi GPU, with 512 active threads, the instruction throughput of the second case is already close to the best situation. On Kepler GPU, however, with fewer than 1024 active threads, the Kepler GPU is very sensitive to the dependences between instructions.

Register and Shared Memory Blocking Factors

Larger register blocking size can introduce more register reuse within one thread and higher percentage of FFMA instructions. However, the register blocking size is limited by the register resource on the SM and the instruction set constraint. With a register blocking factor B R , if we only consider the registers needed for blocking, we can describe the resource constraint as Equation 2.

B 2 R + BR + 1 < RT ≤ RMax (2)
This loose condition for register blocking factor B R can be used to roughly estimate B R . B 2 R is the register set needed to hold C sub-matrix per thread, B R is one column/row of A or B sub-matrix. For instance, with maximum 63 registers per thread, B R ≤ 7.

As depicted in Figure 1, T B * B 2 R is the size of the C sub-matrix per block (each block has T B threads) and T B * B 2 R * L is the size of a sub-matrix for A or B (L is the stride). To overlap the data transfer and the computation, extra registers are needed to fetch data from global memory to shared memory since no direct data transfer is provided between the two memory space. The stride L needs to be chosen such that each thread loads the same amount of data (Equation 3).

( √ TB * BR * L)%TB = 0 (3)
Considering data prefetching from global memory and a few registers to store the addresses of matrices in global memory and shared memory (R addr ), the overall strict constraint for register blocking factor can be described as Equation 4.

B 2 R + 2 * √ TB * BR * L TB +BR +1+R addr ≤ RT ≤ RMax (4)
Since shared memory is allocated in block granularity, for Blk active blocks, Blk * 2 * √ T B * B R * L is needed to store prefetched global memory data (Equation 5). The shared memory blocking factor can be defined as

B Sh = T B * B 2 R .
With the shared memory blocking factor B Sh , the performance bounded by global memory bandwidth can be roughly estimated using Equation 6.

Blk * 2 * √ TB * BR * L ≤ ShSM (5)
P M emBound #GlobalM em bandwidth = 2 * B Sh 2 2 * B Sh * 4 (6)

Potential Peak Performance of SGEMM

The instruction factor F I is the ratio of FFMA instructions in the SGEMM main loop (We only consider FFMA and LDS.X instructions here). It depends on the choice of LDS.X instruction and register blocking factor B R (Figure 3). For instance, if LDS.64 is used with register blocking factor 6, F I = 0.5.

The throughput factor F T is a function of register blocking factor (B R ), number of active threads (T SM ), throughput of SPs (#SP T P ), LD/ST units (#LDS T P ) and dispatch units (#Issue T P )) (Equation 7). The function f for Fermi and Kepler GPUs is illustrated in Figure 2 and in Figure 4 (only shows LDS.64) and obtained through benchmarks varying these parameters.

FT = f (BR, #Issue T P, #SP T P, #LDS T P, TSM ) [START_REF] Kurzak | Autotuning gemm for the fermi gpu[END_REF] With the register blocking factor B R , the instruction factor F I and the throughput factor F T , the performance bounded by SMs' processing throughput is estimated as Equation 8and the overall performance is as Equation 9.

P SM Bound = B 2 R B 2 R + BR * 2 * FI * FT * P theoretical (8)
P potential = min(P M emBound , P SM Bound ) (9)
With the previous analysis, we can estimate the performance upper bound of SGEMM on Fermi and Kepler GPUs. On the Fermi GPU for instance, because of the hard limit of 63 registers (R M ax ) per thread, considering prefetching and using the strict condition of Equation 4, the maximum blocking factor is only 6. The detailed register allocation is illustrated in Section 5.2. With the register blocking factor of 6, the register resource per SM can support up to 512 threads. Using Equation 3, we choose 256 threads per block.

To easily program the data prefetching, according to Equation 3, L could be 8, 16, 24, . . . . Considering the condition in Equation 4, we choose L as 16. With a 6-register blocking factor, mixing LDS or LDS.64 with FFMA instructions, the throughput can achieve close to 32 thread instructions per shader clock per SM. Using a LDS.64 instruction can increase the FFMA instruction percentage to 85.7% from 75% (using LDS). Though LDS.128 instruction can provide higher percentage of FFMA instructions, the instruction processing throughput is too low.

According to Equations 6, 8 and 9, the performance is bounded by SMs' processing throughput, and the potential peak is about 82.5% (

6 2
6 2 +6 * 2 * 0.5 * 30.8 32 ) of the theoretical peak performance for SGEMM. The main limitation comes from the nature of the Fermi instruction set and the limited issue throughput of schedulers.

It is similar to estimate the performance upper bound of SGEMM on Kepler GPU as Fermi GPU. The Kepler GPU (GK104) instruction set is very close to that of Fermi GPU. It means that the limit of 63 registers per thread still exists. Thus, 6-register blocking is also applicable. And the register resource can support 1024 active threads per SM (64K 32bit registers per SM). We can choose either 256 or 1024 threads per block. Similarly, if we use LDS.64 instructions, the FFMA instruction percentage is 85.7%. If we use LDS.128 instructions (need padding or data layout transform), the FFMA instruction percentage is 92.3%.

Similarly, according to Equations 6, 8 and 9, the performance is bounded by SM's processing throughput, and the potential peak is about 54.6% (

6 2 6 2 +6 * 2 * 0.5 * 122.4
192 ) of the theoretical peak performance for SGEMM using LDS.64 instructions. Using LDS.128 instructions, the potential peak is about 57.6% (

6 2
6 2 +6 * 2 * 0.25 * 119.9 192 ) of the theoretical peak. The main limitation factors are still the nature of instruction set and the limited issue throughput of schedulers.

Assembly Code Level Optimization

The estimated performance upper bound is a limit that an actual implementation cannot exceed. It can be a little optimistic since we only consider the major performance degradation factors. Besides the considered parameters, there might be other aspects which can limit the performance. The 'real' upper bound or the best possible performance is between the estimated upper bound and the achieved performance.

Depending on whether to apply transpose operation on input matrix A or B, there are 4 variations for GEMM kernel. Figure 5 illustrates the performance of four SGEMM variations from CUBLAS and our implementation (ASM) with 2400x2400 and 4800x4800 matrices. On GTX580 GPU, we achieve around 74.2% of the theoretical peak performance, i.e., about 90% of the estimated performance upper bound, which we think is good enough. In our analysis, we only study the two main instruction types. There are other auxiliary instructions which do not devote to the GFLOPS. And also, we do not consider the effect of barriers which will harm the performance too. We show that the 'real' upper bound is within this 10% and future optimization is unlikely to achieve a lot of speedup. On Kepler GPU, although we cannot provide the optimal controlling information as discussed in section 3.2, we achieve around 77.3% of the estimated upper bound. Similar to Fermi GPU, there are some factors we do not consider in our analysis. The larger gap between our achieved performance might be due to our very limited knowledge of the undisclosed scheduling information of Kepler GPU, which is critical to performance or to some hidden characteristics that we are not able to discover due to limited documentation. Figure 6 illustrates the performance comparison on Fermi GPU between our implementation (assembly), CUBLAS from CUDA 4.1 and MAGMA library [START_REF] Nath | An improved magma gemm for fermi gpus[END_REF]. Figure 7 is the performance comparison on Kepler GPU between our implementation (assembly), CUBLAS from CUDA 4.2 and MAGMA library.

The rest of the section briefly describes our optimizations on assembly code level of SGEMM.

Optimization of Memory Accesses

Assembly code level optimization of memory accesses is similar to high level language optimization. Global memory requests from the threads within a warp could be grouped (coalesced) into one or more memory transactions depending on the compute capability of the device and the memory accessing pattern. To access global memory efficiently, generally it is better to let threads in a warp access continuous data elements in global memory to get coalescing. Considering the majority of instructions in the SGEMM main loop are FFMA and LDS, and it is critical to reduce the number of LDS instructions (using LDS.64 or LDS.128), sub-matrices in shared memory should be grouped such that each thread accesses continuous B R data elements. Also, proper padding needs to be applied to reduce shared memory access conflicts and satisfy the alignment restriction of the LDS instruction.

Register Spilling Elimination

The register resource is 32K 32-bit registers per SM for the Fermi GPU and each thread can use a maximum of 63 registers. The register R1 is normally occupied as stack pointer. According to our analysis, the number of per-thread registers with prefetching is at least

B 2 R + 2 * √ T B * B R * L T B + B R + 1 + R index .
With the register blocking factor of 6 for Fermi GPU, the register allocation of our implementation is as the following. Note that we use 32bit addressing to save address registers. In all, 63 registers are used. Since we do not need thread stack, R1 is used to store the loop end condition in our code. Therefore, we are able to fully eliminate the register spilling.

Instruction Reordering

Generally, we try to interleave different instruction types to get better balance between functional units within one SM and better instruction throughput. We apply the following simple reordering optimizations:

1. In the main loop, between the 2 barriers are all shared memory accesses. By moving address calculation from start of the loop to mix with the shared memory accesses, we can achieve better performance.

2. Interleaving prefetching from global memory with FFMA and LDS instructions can benefit performance.

Register

. Register Allocation

As in Figure 8, around 30% of the FFMA instructions in the MAGMA [START_REF] Nath | An improved magma gemm for fermi gpus[END_REF] SGEMM binary for the Kepler GPU (nvcc generated) have the 2-way register bank conflict and 1% of the FFMA instructions have the 3-way register bank conflict. In our first version of SGEMM NN on GTX680, which achieves around 1100GFLOPS, 68.8% of the FFMA instructions have the 2-way register bank conflict, and 10.6% of the 3-way conflict. After applying the optimization, the modified version, which achieve around 1300GFLOPS, has only 1.2% of the 2-way FFMA register bank conflict and the 3-way conflict is fully removed.

Our optimization is depicted in Figure 9. In the SGEMM main loop, at each stage, one column from matrix A and one row from matrix B are processed. To use the register blocking and the LDS.64 instructions, at least 6 and 2 different registers are needed for column A and row B. Of course, there are many possible implementations, here we describe one possibility. We select registers from E0 and O0 for column A. Row B uses registers from E1 and O1. Then we use the first table in Figure 9 as the constraints of register allocation. In the final mapping stage, we make sure that 36 registers of C sub-matrix have 9 registers on each bank and had our register allocation as the second table, which does not have any register bank conflict to compute the 36 elements from the C sub-matrix.

Opportunity for Automatic Tools

Our study emphasizes that for Fermi and Kepler GPUs, it is essential to study the impact of algorithm parameters on instruction throughput to get insight into the performance result. The main optimization opportunity comes from the allocation of registers. For example, the four SGEMM variations of MAGMA library compiled with nvcc spill at least 10 registers (40 Bytes) on the Kepler GPU. When the active thread number is 512, at least 20KB L1 cache is needed to make sure that the spilled data stays in the L1 cache. However, since normally the unified 64KB shared memory/L1 cache is configured as 48KB shared memory and 16KB L1 cache, some data will be spilled out of L1 cache. As the active threads increase, more data is spilled out of L1 cache and the performance will be harmed. We already show that with careful design, register spilling could be eliminated. We also show that around 30% of FFMA instructions in the nvccgenerated SGEMM binary from MAGMA library have register bank conflict. We propos a simple solution in Section 5.4. It is possible for optimizers to detect the loop structure and remove the conflicts with proper register allocation.

An automatic tuning tool normally needs to explore a large design space and evaluate the performance of many configurations [START_REF] Kurzak | Autotuning gemm for the fermi gpu[END_REF][START_REF] Meng | Grophecy: Gpu performance projection from cpu code skeletons[END_REF][START_REF] Ryoo | Program optimization space pruning for a multithreaded gpu[END_REF]. It may take a significant amount of time. Normally, the automatic tuning tool is applicationdependent and each includes several efficient optimizations for the specific application. To build the tool relies on the developers' understanding of the application and optimization experience. With the proposed analysis approach, we can better understand which parameters are critical to the performance. The estimated upper bound actually corresponds to a set of parameters and optimization options. This knowledge can help an automatic tool to explore the design space in a relatively small region. And of course by comparing the performance of an automatic tool's output code and the estimated performance upper bound, we can judge whether the optimized version is good enough.

In our analysis, to study the instruction throughput mixing FFMA and LDS.X instructions, we manually write some benchmarks varying several key parameters such as instruction type choice (LDS.X), the mixing ratio, the blocking factor, the instructions' dependence, active threads and study these parameters' impact on the instruction throughput. For many applications with few major instruction types, a similar approach can be used to estimate the performance upper bound. The difference would be the chosen instruction types and their mixing pattern (mixing ratio, dependence, etc.). Systematic and automatic development of a set of microbenchmarks to help to estimate the performance upper bound of other applications is possible. A family of assembly level microbenchmarks could be defined and evaluated in order to provide a small database of performance references that could be used by the auto-tuning tool, and also the developer to transform the code for performance. Generally, the assembly level microbenchmarks can also help to un-derstand the difference between different GPU architectures. For example, the benchmarks illustrated in Figure 4 show the increasing need for active threads on Kepler GPU. Assembly level benchmarking requires an assembly tool chain which is missing from the official support. We manage to make it work on Fermi GPU. But on Kepler, there are some issues like the hidden scheduling information, which we cannot fully decrypt.

Conclusion

In this work, we have proposed an approach to analyze GPU applications' performance upper bound. Different from existing works on GPU performance models, our approach relies on application analysis and assembly level benchmarking. As an example, we analyze the potential peak performance of SGEMM on Fermi and Kepler GPUs. We show that the nature of the Fermi (Kepler) instruction set and the limited issue throughput of schedulers are the main limitation factors for SGEMM to approach the theoretical peak performance. The general guideline is to reduce the auxiliary instructions and increase the FFMA instruction's percentage. Proper register allocation, shared memory data layout and memory access pattern need to be carefully designed to minimize the impact of memory accesses on performance. We also show that our analysis can help to decide some critical algorithm parameters and show how much optimization space exists. Guided by the analysis, we further optimize the four SGEMM kernel variations and achieve better performance on Fermi GPU (around 5% on average for large matrices) than highly optimized routine provided by NVIDIA.

Figure 1 .

 1 Figure 1. SGEMM Implementation case. We name the four banks as even 0(R index %8 < 4 && R index %2 == 0), even 1 (R index %8 ≥ 4 && R index %2 == 0), odd 0 (R index %8 < 4 && R index %2 == 1), and odd 1(R index %8 ≥ 4 && R index %2 == 1). Since we implement SGEMM with assembly code directly, the register indices have to be carefully chosen to make sure there is no bank conflict. The detailed optimization is illustrated in section 5.4.

Figure 2 .

 2 Figure 2. Thread Instruction Throughput Mixing FFMA and LDS.X

Figure 3 .

 3 Figure 3. FFMA Instruction Percentage in SGEMM Mainloop with Different Register Blocking Factors

Figure 4 .

 4 Figure 4. Instruction Throughput Mixing FFMA and LDS.64 with Ratio of 6:1

Figure 5 .Figure 6 .

 56 Figure 5. SGEMM Performance of CUBLAS and Our Implementation on Fermi and Kepler GPUs

Figure 7 .

 7 Figure 7. SGEMM NN Performance on GTX680

Table 2 .

 2 Examples of Math Instruction Throughput on Ke-

pler GPU with Various Operand Register Indices notations are embedded with our parsing tool. Each thread executes the same 8192 math instructions. Each block has 1024 threads without synchronization and 40960 blocks are spawned to keep the GPU busy.

  Allocation for Kepler GPU Register Conflict of FFMA InstructionAs we describe in Section 3.3, to get the best throughput, the 3 source registers of FFMA instructions should reside on 3 different banks if they are different. In our current implementation, 6-register blocking is used. 6 registers are used to load A from the shared memory and 2 registers to load B from the shared memory in the main loop. 36 different registers (R26∼R61) hold the C sub-matrix. In this implementation, register spilling is eliminated.

	A	C B	Col A	E0 O0 E0 O0 E0 O0	E1 O O1 /E0 O O1 /E0 O O1 /E0	O1 E1/ O0 E E1/ O0 E E1/ O0 E	0% 20% 40% 60% 80% 100% Percentage O1 E1/ O0 E E1/ O0 E E1/ O0 E Figure 8. R6 ma gm a_ NN ma gm a_ NT ma gm a_ TN ma gm a_ TT as m_ NN mo d_ as m_ NN No Conflict 2-way Conflict 3-way Conflict R7 R27 R28 R33 R30 R35 R36 R29 R26 R31 R32 R37 R34 R41 R38 R43 R44 R49 R46 R39 R40 R45 R42 R47 R48 R51 R52 R57 R54 R59 R60 R53 R50 R55 R56 R61 R58 R8 R9 R10 R11 R18 R19 E1 O O1 /E0 O O1 /E0 O O1 /E0 Row B O1 E1/ O0 E E1/ O0 E E1/ O0 E Sub-matrix C E1 O O1 /E0 O O1 /E0 O O1 /E0

GEMM performs the matrix-matrix operation C := alpha * op(A) * op(B) + beta * C. alpha and beta are scalars, and A, B and C are matrices. op(X) is op(X) = X or op(X) = X T .

Acknowledgments

Many thanks to CJ Newburn, our shepherd and the anonymous reviewers for their insightful comments on improving this paper. This work is supported by French National Research Agency (ANR) through COSINUS program (project PETAQCD N o ANR-08-COSI-010).


Gregory Kucherov
email: gregory.kucherov@univ-mlv.fr

Yakov Nekrich
email: yakov.nekrich@googlemail.com

Tatiana Starikovskaya
email: tat.starikovskaya@gmail.com

Cross-Document Pattern Matching

We study a new variant of the string matching problem called cross-document string matching, which is the problem of indexing a collection of documents to support an efficient search for a pattern in a selected document, where the pattern itself is a substring of another document. Several variants of this problem are considered, and efficient linear-space solutions are proposed with query time bounds that either do not depend at all on the pattern size or depend on it in a very limited way (doubly logarithmic). As a side result, we propose an improved solution to the weighted level ancestor problem.

Introduction

In this paper we study the following variant of the string matching problem that we call cross-document string matching: given a collection of strings (documents) stored in a "database", we want to be able to efficiently search for a pattern in a given document, where the pattern itself is a substring of another document. More formally, assuming we have a set of documents T 1 , . . . , T m , we want to answer queries about the occurrences of a substring T k [i..j] in a document T ℓ .

This scenario may occur in various situations when we have to search for a pattern in a text stored in a database, and the pattern is itself drawn from a string from the same database. In bioinformatics, for example, a typical project deals with a selection of genomic sequences, such as a family of genomes of evolutionary related species. A common repetitive task consists then in looking for genomic elements belonging to one of the sequences in some other sequences. These elements may correspond to genes, exons, mobile elements of any kind, regulatory patterns, etc., and their location (i.e. start and end positions) in the sequence of origin is usually known from a genome annotation provided by a sequence data repository (such as GenBank or any other). A similar scenario may occur in other application fields, such as the bibliographic search for example.

In this paper, we study different versions of the cross-document string matching problem. First, we distinguish between counting and reporting queries, asking respectively about the number of occurrences of T k [i..j] in T ℓ or about the occurrences themselves. The two query types lead to slightly different solutions. In particular, the counting problem uses the weighted level ancestor problem [START_REF] Farach | Perfect hashing for strings: Formalization and algorithms[END_REF][START_REF] Amir | Dynamic text and static pattern matching[END_REF] to which we propose a new solution with an improved complexity bound.

We further consider different variants of the two problems. The first one is the dynamic version where new documents can be added to the database. In another variant, called document counting and reporting, we only need to respectively count or report the documents containing the pattern, rather than counting or reporting pattern occurrences within a given document. This version is very close to the document retrieval problem previously studied (see [START_REF] Muthukrishnan | Efficient algorithms for document retrieval problems[END_REF] and later papers referring to it), with the difference that in our case the pattern is itself selected from the documents stored in the database. Finally, we also consider succinct data structures for the above problems, where we keep involved index data structure in compressed form.

Let m be the number of stored strings and n the total length of all strings. Our results are summarized below: (i) for the counting problem, we propose a solution with query time O(t + log log m), where t = min( log occ/ log log occ, log log |P |), P = T k [i..j] is the searched substring and occ is the number of its occurrences in T ℓ , (ii) for the reporting problem, our solution outputs all the occurrences in time O(log log m + occ), (iii) in the dynamic case, when new documents can be dynamically added to the database, we are able to answer counting queries in time O(log n) and reporting queries in time O(log n + occ), whereas the updates take time O(log n) per character, (iv) for the document counting and document reporting problems, our algorithms run in time O(log n) and O(t + ndocs) respectively, where t is as above and ndocs is the number of reported documents, (v) finally, we also present succinct data structures that support counting, reporting, and document reporting queries in cross-document scenario (see Theorems 6 and 7 in Section 4.3).

For problems (i)-(iv), the involved data structures occupy O(n) space. Interestingly, in the cross-document scenario, the query times either do not depend at all on the pattern length or depend on it in a very limited (doubly logarithmic) way.

Throughout the paper positions in strings are numbered from 1. Notation T [i..j] stands for the substrings T [i]T [i + 1] . . . T [j] of T , and T [i..] denotes the suffix of T starting at position i.

Preliminaries

Basic Data Structures

We assume a basic knowledge of suffix trees and suffix arrays.

Besides using suffix trees for individual strings T i , we will also be using the generalized suffix tree for a set of strings T 1 , T 2 , . . . , T m that can be viewed as the suffix tree for the string T 1 $ 1 T 2 $ 2 . . . T m $ m . A leaf in a suffix tree for T i is associated with a distinct suffix of T i , and a leaf in the generalized suffix tree is associated with a suffix of some document T i together with the index i of this document. We assume that for each node v of a suffix tree, the number n v of leaves in the subtree rooted at v, as well as its string depth d(v) can be recovered in constant time. Recall that the string depth d(v) is the total length of strings labelling the edges along the path from the root to v.

We will also use the suffix arrays for individual documents as well as the generalized suffix array for strings T 1 , T 2 , . . . , T m . Each entry of the suffix array for T i is associated with a distinct suffix of T i and each entry of the generalized suffix array for T 1 , . . . , T m is associated with a suffix of some document T i and the index i of the document the suffix comes from. We store these document indices in a separate array D, called document array, such that D[i] = k if the i-th entry of the generalized suffix array for T 1 , . . . , T m corresponds to a suffix coming from T k .

For each considered suffix array, we assume available, when needed, two auxiliary arrays: an inverted suffix array and another array, called the LCP-array, of longest common prefixes between each suffix and the preceding one in the lexicographic order.

Weighted Level Ancestor Problem

The weighted level ancestor problem, defined in [START_REF] Farach | Perfect hashing for strings: Formalization and algorithms[END_REF], is a generalization of the level ancestor problem [START_REF] Berkman | Finding level-ancestors in trees[END_REF][START_REF] Bender | The level ancestor problem simplified[END_REF] for the case when tree edges are assigned positive weights.

Consider a rooted tree T whose edges are assigned positive integer weights. For a node w, let weight(w) denote the total weight of the edges on the path from the root to w; depth(w) denotes the usual tree depth of w.

A weighted level ancestor query wla(v, q) asks, given a node v and a positive integer q, for the ancestor w of v of minimal depth such that weight(w) ≥ q (wla(v, q) is undefined if there is no such node w).

Two previously known solutions [START_REF] Farach | Perfect hashing for strings: Formalization and algorithms[END_REF][START_REF] Amir | Dynamic text and static pattern matching[END_REF] for weighted level ancestors problem achieve O(log log W ) query time using linear space, where W is the total weight of all tree edges. Our data structure also uses O(n) space, but achieves a faster query time in many special cases. We prove the following result.

Theorem 1. There exists an O(n) space data structure that answers weighted ancestor query wla(v, q) in O(min( log g/ log log g, log log q)) time, where g = min(depth(wla(v, q)), depth(v) -depth(wla(v, q))).

If every internal node is a branching node, we obtain the following corollary.

Corollary 1. Suppose that every internal node in T has at least two children. There exists an O(n) space data structure that finds w = wla(v, q) in O( log n w / log log n w ) time, where n w is the number of leaves in the subtree of w.

Our approach combines the heavy path decomposition technique of [START_REF] Amir | Dynamic text and static pattern matching[END_REF] with efficient data structures for finger searching in a set of integers. The proof is given in the Appendix.

3 Cross-document Pattern Counting and Reporting

Counting

In this section we consider the problem of counting occurrences of a pattern T k [i..j] in a document T ℓ .

Our data structure consists of the generalized suffix array GSA for documents T 1 , . . . , T m and individual suffix trees T i for every document T i .

For every suffix tree T ℓ we store a data structure of Theorem 1 supporting weighted level ancestor queries on T ℓ . We also augment the document array D with an O(n)-space data structure that answers queries rank(k, i) (number of entries storing k before position i in D) and select(k, i) (i-th entry from the left storing k). Using the result of [START_REF] Golynski | Rank/select operations on large alphabets: a tool for text indexing[END_REF], we can support such rank and select queries in O(log log m) and O(1) time respectively. Moreover, we construct a data structure that answers range minima queries (RMQ) on the LCP array: for any 1 ≤ r 1 ≤ r 2 ≤ n, find the minimum among LCP [r 1 ], . . . LCP [r 2 ]. There exists a linear space RMQ data structure that supports queries in constant time, see e.g., [START_REF] Bender | The lca problem revisited[END_REF]. An RMQ query on the LCP array computes the length of the longest common prefix of two suffixes GSA[r 1 ] and GSA[r 2 ], denoted LCP (r 1 , r 2 ).

Our counting algorithm consists of two stages. First, using GSA, we identify a position p of T ℓ at which the query pattern T k [i..j] occurs, or determine that no such p exists. Then we find the locus of T k [i..j] in the suffix tree T ℓ using a weighted ancestor query.

Let r be the position of T k [i..] in the GSA. We find indexes r 1 = select(ℓ, rank(r, ℓ)) and

r 2 = select(ℓ, rank(r, ℓ)+1) in O(log log m) time. GSA[r 1 ] (resp. GSA[r 2 ]
) is the closest suffix from document T ℓ that precedes (resp. follows) T k [i..] in the lexicographic order of suffixes. Observe now that T k [i..j] occurs in T ℓ if and only if either LCP (r 1 , r) or LCP (r, r 2 ) (or both) is no less than j -i + 1. If this holds, then the starting position

p of GSA[r 1 ] (respectively, starting position of GSA[r 2 ]) is the position of T k [i..j] in T ℓ . Once such a position p is found, we jump to the leaf v of T ℓ that contains the suffix T ℓ [p..].
The weighted level ancestor u = wla(v, (j

-i + 1)) is the locus of T k [i..j] in T ℓ . This is because T ℓ [p..p + j -i] = T k [i..j].
By Corollary 1, we can find node u in O( log n u / log log n u ) time, where n u is the number of leaf descendants of u. Since u is the locus node of T k [i..j], n u is the number of occurrences of T k [i..j] in T ℓ . By Theorem 1, we can find u in O(log log(j -i + 1)) time.

Summing up, we obtain the following Theorem.

Theorem 2. For any 1 ≤ k, ℓ ≤ m and 1 ≤ i ≤ j ≤ |T k |, we can count the number of occurrences of T k [i..j] in T ℓ in O(t + log log m) time, where t = min( log occ/ log log occ, log log(j -i + 1)) and occ is the number of occurrences. The underlying indexing structure takes O(n) space and can be constructed in O(n) time.

Reporting

A reporting query asks for all occurrences of a substring T k [i..j] in T ℓ .

Compared to counting queries, we make a slight change in the data structures: instead of using suffix trees for individual documents T i , we use suffix arrays. The rest of the data structures is unchanged.

We first find an occurrence of T k [i..j] in T ℓ (if there is one) with the method described in Section 3.1. Let p be the position of this occurrence in T ℓ . We then jump to the corresponding entry r of the suffix array SA ℓ for the document T ℓ . Let LCP ℓ be the LCP-array of SA ℓ . Starting with entry r, we visit adjacent entries t of SA ℓ moving both to the left and to the right as long as LCP ℓ [t] ≥ j -i + 1. While this holds, we report SA ℓ [t] as an occurrence of T k [i..j]. It is easy to observe that the procedure is correct and that no occurrence is missing. As a result, we obtain the following theorem.

Variants of the Problem

Dynamic Counting and Reporting

In this section we focus on a dynamic version of counting and reporting problems, where the only dynamic operation consists in adding a document to the database 4 .

Recall that in the static case, counting occurrences of T k [i..j] in T ℓ is done through the following two steps (Section 3.1):

1. compute position p of some occurrence of T k [i..j] in T ℓ , 2. in the suffix tree of T ℓ , find the locus of string T ℓ [p..p + j -i], and retrieve the number of leaves in the subtree rooted at u.

For reporting queries (Section 3.2), Step 1 is basically the same, while Step 2 is different and uses an individual suffix array for T ℓ .

In the dynamic framework, we follow the same general two-step scenario. Note first that since Step 2, for both counting and reporting, uses data structures for individual documents only, it trivially applies to the dynamic case without changes. However, Step 1 requires serious modifications that we describe below.

Since the suffix array is not well-suited for dynamic updates, at Step 1 we will use the generalized suffix tree for T 1 , T 2 , . . . , T m hereafter denoted GST . For each suffix of T 1 , T 2 , . . . , T m we store a pointer to the leaf of GST corresponding to this suffix.

We maintain a dynamic doubly-linked list EL corresponding to the Euler tour of the current GST . Each internal node of GST is stored in two copies in EL, corresponding respectively to the first and last visits of the node during the Euler tour. Leaves of GST are kept in one copy. Observe that the leaves of GST appear in EL in the "left-to-right" order, although not consecutively.

On EL, we maintain the data structure of [START_REF] Bender | Two simplified algorithms for maintaining order in a list[END_REF] which allows, given two list elements, to determine their order in the list in O(1) time (see also [START_REF] Dietz | Two algorithms for maintaining order in a list[END_REF]). Insertions of elements in the list are supported in O(1) time too.

Furthermore, we maintain a balanced tree, denoted BT , whose leaves are elements of EL. Note that the size of EL is bounded by 2n (n is the size of GST ) and the height of BT is O(log n). Since the leaves of GST are a subset of the leaves of BT , we call them suffix leaves to avoid the ambiguity.

Each internal node u of BT stores two kinds of information: (i) the rightmost and leftmost suffix leaves in the subtree of BT rooted at u, (ii) minimal LCP value among all suffix leaves in the subtree of BT rooted at u.

Finally, we will also need an individual suffix array for each inserted document T i .

We are now in position to describe the algorithm of Step 1. Like in the static case, we first retrieve the leaf of GST corresponding to suffix T k [i..]. To identify a position of an occurrence of T k [i..j] in T ℓ , we have to examine the two closest elements in the list of leaves of GST , one from right and from left, corresponding to suffixes of T ℓ . To find these two suffixes, we perform a binary search on the suffix array for T ℓ using order queries of [START_REF] Bender | Two simplified algorithms for maintaining order in a list[END_REF] on EL. This step takes O(log |T ℓ |) time.

We then check if at least one of these two suffixes corresponds to an occurrence of T k [i..j] in T ℓ . In a similar way to Section 3, we have to compute the longest common prefix between each of these two suffixes and T k [i..], and compare this value with (j -i+1). This amounts to computing the minimal LCP value among all the suffixes of the corresponding range.

This can be done in O(log n) time by using a standard range trees approach [?]: for any sublist of EL we can retrieve O(log n) nodes v i that cover it. The least among all minimal LCP values stored in nodes v i is the minimal LCP value for the specified range of suffixes.

The query time bounds are summarized in the following lemma. When a new node v is added to a suffix tree, the following updates should be done (in order):

(i) insert v at the right place of the list EL (in two copies if v is an internal node), (ii) rebalance the tree BT if needed, (iii) if v is a leaf of GST (i.e. a suffix leaf of BT ), update LCP values and rightmost/leftmost suffix leaf information in BT , To see how update (i) works, we have to recall how suffix tree is updated when a new document is inserted. Two possible updates are creation of a new internal node v by splitting an edge into two (edge subdivision) and creating a new leaf u as a child of an existing node. In the first case, we insert the first copy of v right after the first copy of its parent, and the second copy right before the second copy of its parent. In the second case, the parent of u has already at least one child, and we insert u either right after the second (or the only) copy of its left sibling, or right before the first (or the only) copy of its right sibling.

Rebalancing the tree BT (update (ii)) is done using standard methods. Observe that during the rebalancing we may have to adjust the LCP and rightmost/leftmost suffix leaf information for internal nodes, but this is easy to do as only a constant number of local modifications is done at each level.

Update (iii) is triggered when a new leaf u is created in GST and added to EL. First of all, we have to compute the LCP value for u and possibly to update the LCP value of the next suffix leaf u ′ to the right of u in EL. This is done in O(1) time as follows. At the moment when u is created, we memorize the string depth of its parent D = d(parent(u)). Recall that the parent of u already has at least one child before u is created. If u is neither the leftmost nor the rightmost child of its parent, then we set LCP (u) = D and LCP (u ′ ) remains unchanged (actually it also equals D). If u is the leftmost child of its parent, then we set LCP (u) = LCP (u ′ ) and then LCP (u ′ ) = D. Finally, if u is the rightmost child, then LCP (u) = D and LCP (u ′ ) remains unchanged.

We then have to follow the path in BT from the new leaf u to the root and possibly update the LCP and rightmost/leftmost suffix leaf information for all nodes on this path. These updates are straightforward. Furthermore, during this traversal we also identify suffix leaf u ′ (as the leftmost child of the first right sibling encountered during the traversal), update its LCP value and, if necessary, the LCP values on the path from u ′ to the root of BT . All these steps take time O(log n).

Thus, updates of all involved data structures take O(log n) time per symbol. The following theorem summarizes the results of this section. Theorem 4. In the case when documents can be added dynamically, the number of occurrences of T k [i..j] in T ℓ can be computed in time O(log n) and reporting these occurrences can be done in time O(log n + occ), where occ is their number. The underlying data structure occupies O(n) space and an update takes O(log n) time per character.

Document Counting and Reporting

Consider a static collection of documents T 1 , . . . , T m . In this section we focus on document reporting and counting queries: report or count the documents which contain at least one occurrence of T k [i..j], for some 1 ≤ k ≤ m and i ≤ j.

For both counting and reporting, we use the generalized suffix tree, generalized suffix array and the document array D for T 1 , T 2 , . . . , T m . We first retrieve the leaf of the generalized suffix tree labelled by T k [i..] and compute its highest ancestor u of string depth at least j -i + 1, using the weighted level ancestor technique of Section 2.2. The suffixes of T 1 , T 2 , . . . , T m starting with T k [i..j] (i.e. occurrences of T k [i..j]) correspond then to the leaves of the subtree rooted at u, and vice versa. As shown in Section 3.1, this step takes O(t) time, where t = min( log occ/ log log occ, log log(j -i + 1)) and occ is the number of occurrences of T k [i..j] (this time in all documents).

Once u has been computed, we retrieve the interval [lef t(u)..right(u)] of ranks of all the leaves under interest. We are then left with the problem of counting/reporting distinct values in D[lef t(u)..right(u)]. This problem is exactly the same as the color counting/ color reporting problem that has been studied extensively (see e.g., [START_REF] Gagie | Colored range queries and document retrieval[END_REF] and references therein).

For color reporting queries, we can use the solution of [START_REF] Muthukrishnan | Efficient algorithms for document retrieval problems[END_REF] based on an O(n)-space data structure for RMQ, applied to (a transform of) the document array D. The pre-processing time is O(n). Each document is then reported in O(1) time, i.e. all relevant documents are reported in O(ndocs) time, where ndocs is their number. The whole reporting query then takes time O(t + ndocs) for t defined above.

For counting, we use the solution described in [START_REF] Bozanis | New upper bounds for generalized intersection searching problems[END_REF]. The data structure requires O(n) space and a color counting query takes O(log n) time. The following theorem presents a summary. Theorem 5. We can store a collection of documents T 1 , . . . , T m in a linear space data structure, so that for any pattern P = T k [i..j] all documents that contain P can be reported and counted in O(t + ndocs) and O(log n) time respectively. Here t = min( log occ/ log log occ, log log |P |), ndocs is the number of documents that contain P and occ is the number of occurrences of P in all documents.

Compact Counting, Reporting and Document Reporting

In this section, we show how our reporting and counting problems can be solved on succinct data structures [START_REF] Navarro | Compressed full-text indexes[END_REF].

Reporting and Counting. Our compact solution is based on compressed suffix arrays [START_REF] Grossi | Compressed suffix arrays and suffix trees with applications to text indexing and string matching (extended abstract)[END_REF]. A compressed suffix array for a text T uses |CSA| bits of space and enables us to retrieve the position of the suffix of rank r, the rank of a suffix T [i..], and the character T [i] in time Lookup(n). Different trade-offs between space usage and query time can be achieved (see [START_REF] Navarro | Compressed full-text indexes[END_REF] for a survey).

Our data structure consists of a compressed generalized suffix array CSA for T 1 , . . . , T m and compressed suffix arrays CSA i for each document T i . In [START_REF] Sadakane | Compressed suffix trees with full functionality[END_REF] In [START_REF] Sadakane | Succinct data structures for flexible text retrieval systems[END_REF], it was shown how to report, for any 1 ≤ r 1 ≤ r 2 ≤ n, all distinct documents T f such that at least one suffix of T f occurs at position r, r 1 ≤ r ≤ r 2 , of the generalized suffix array. The construction uses O(n + m log n m ) additional bits, and all relevant documents are reported in O(Lookup(n) • ndocs) time, where ndocs is the number of documents that contain T k [i..j]. Summing up, we obtain the following result.

Theorem 3 .

 3 All the occurrences of T k [i..j] in T ℓ can be reported in O(log log m + occ) time, where occ is the number of occurrences. The underlying indexing structure takes O(n) space and can be constructed in O(n) time.

Lemma 1 .

 1 Using the above data structures, counting and reporting all occurrences of T k [i..j] in T ℓ can be done respectively in time O(log n) and time O(log n + occ), where occ is the number of reported occurrences. We now explain how the involved data structures are updated. Suppose that we add a new document T m+1 . Extending the generalized suffix tree by T m+1 is done in time O(|T m+1 |) by McCreight's or Ukkonen's algorithm, i.e. in O(1) amortized time per symbol.

Theorem 6 .

 6 it was shown that using O(n) extra bits, the length of the longest common prefix of any two suffixes can be computed in O(Lookup(n)) time. Besides, the ranks of any two suffixes T k [s..] and T ℓ [p..] can be compared in O(Lookup(n)) time: it suffices to compareT ℓ [p + f ] with T k [s + f ] for f = LCP (T k [s..], T ℓ [p..]).Note that ranks of the suffixes of T ℓ starting with T k [i..j] form an interval [r 1 ..r 2 ]. We use a binary search on the compressed suffix array of T ℓ to find r 1 and r 2 . At each step of the binary search we compare a suffix ofT ℓ with T k [i..]. Therefore [r 1 ..r 2 ] can be found in O(Lookup(n) • log n) time. Obviously, the number of occurrences of T k [i..j] in T ℓ is r 2 -r 1 .To report the occurrences, we compute the suffixes of T ℓ with ranks in interval[r 1 ..r 2 ]. All occurrences of T k [i..j] in T ℓ can be counted in O(Lookup(n)•log n) timeand reported in O((log n+occ)Lookup(n)) time, where occ is the number of those. The underlying indexing structure takes 2|CSA| + O(n + m log n m ) bits of memory. Document Reporting Again, we use a binary search on the generalized suffix array to find the rank interval [r 1 ..r 2 ] of suffixes that start with T k [i..j]. This can be done in O(Lookup(n) • log n) time.

Theorem 7 .

 7 All documents containing T k [i..j] can be reported in O((log n + ndocs)Lookup(n)) time, where ndocs is the number of those. The underlying indexing structure takes 2|CSA| + O(n + m log n m ) bits of space.

document deletions are also possible to support but require some additional constructions that are left to the extended version of this paper

Acknowledgments: T.Starikovskaya has been supported by the mobility grant funded by the French Ministry of Foreign Affairs through the EGIDE agency and by a grant 10-01-93109-CNRS-a of the Russian Foundation for Basic Research. Part of this work has been done during a one-month stay of Y.Nekrich at the Marne-la-Vallée University supported by the BEZOUT grant of the French government. The authors also would like to thank Tsvi Kopelowitz for fruitful discussions and valuable comments.

Appendix: Proof of Theorem 1

Here we prove Theorem 1. We use the heavy path decomposition technique of [START_REF] Amir | Dynamic text and static pattern matching[END_REF].

A path π in T is heavy if every node u on π has at most twice as many nodes in its subtree as its child v on π. A tree T can be decomposed into paths using the following procedure: we find the longest heavy path π r that starts at the root of T and remove all edges of π r from T . All remaining vertices of T belong to a forest; we recursively repeat the same procedure in every tree of that forest.

We can represent the decomposition into heavy paths using a tree T. Each node v j in T corresponds to a heavy path π j in T . A node v j is a child of a node v i in T if the head of π j (i.e., the highest node in π j ) is a child of some node u ∈ π i . Some node in π i has at least twice as many descendants as each node in π j ; hence, T has height O(log n). O(n log n)-Space Solution. Let p j denote a root-to-leaf path in T. For a node v in T let weight(v) denote the weight of the head of π, where π is the heavy path represented by v in T. We store a data structure D(p j ) that contains the values of weight(v) for all nodes v ∈ p j . D(p j ) contains O(log n) elements; hence, we can find the highest node v ∈ p j such that weight(v) ≥ q in O(1) time. This can be achieved by storing the weights of all nodes from p j in the q-heap [START_REF] Fredman | Trans-dichotomous algorithms for minimum spanning trees and shortest paths[END_REF].

For every heavy path π j , we store the data structure E(π j ) from [START_REF] Andersson | Dynamic ordered sets with exponential search trees[END_REF] that contains the weights of all nodes u ∈ π j and supports the following queries: for an integer q, find the lightest node u ∈ π j such that weight(u) ≥ q. Using Theorem 1.5 in [START_REF] Andersson | Dynamic ordered sets with exponential search trees[END_REF], we can find such a node u ∈ π j in O( log n ′ / log log n ′ ) time where

Moreover, we can also find the node u in O(log log q) time; we will show how this can be done in the full version of this paper. Thus E(π j ) can be modified to support queries in O(min( log n ′ / log log n ′ , log log q)) time.

For each node u ∈ T we store a pointer to the heavy path π that contains u and to the corresponding node v ∈ T.

A query wla(v, q) can be answered as follows. Let v denote the node in T that corresponds to the heavy path containing v. Let p j be an arbitrary root-to-leaf path in T that also contains v. Using D(p j ) we can find the highest node u ∈ p j , such that weight(u) ≥ q in O(1) time. Let π t denote the heavy path in T that corresponds to the parent of u, and π s denote the path that corresponds to u. If the weighted ancestor wla(v, q) is not the head of π s , then wla(v, q) belongs to the path π t . Using E(π t ), we can find u = wla(v, q) in O(min( log n ′ / log log n ′ , log log q)) time where

All data structures E(π i ) use linear space. Since there are O(n) leaves in T and each path

Lemma 2. There exists an O(n log n) space data structure that finds the weighted level ancestor u in O(min( log n ′ / log log n ′ , log log q)) time.

O(n)-Space Solution. We can reduce the space from O(n log n) to O(n) using a micro-macro tree decomposition. Let T 0 be a tree induced by the nodes of T that have at least log n/8 descendants. The tree T 0 has at most O(n/ log n) leaves. We construct the data structure described above for T 0 ; since T 0 has O(n/ log n) leaves, its heavy-path tree T 0 also has O(n/ log n) leaves. Therefore all structures D(p j ) use O(n) words of space. All E(π i ) also use O(n) words of space. If we remove all nodes of T 0 from T , the remaining forest F consists of O(n) nodes. Every tree T i , i ≥ 1, in F consists of O(log n) nodes. Nodes of T i are stored in a data structure that uses linear space and answers weighted ancestor queries in O( 1) time. This data structure will be described later in this section. Suppose that a weighted ancestor wla(v, q) should be found. If v ∈ T 0 , we answer the query using the data structure for T 0 . If v belongs to some T i for i ≥ 1, we check the weight w r of root(T i ). If w r ≤ q, we search for wla(v, q) in T i . Otherwise we identify the parent v 1 of root(T i ) and find wla(v 1 , q) in T 0 . If wla(v 1 , q) in T 0 is undefined, then wla(v, q) = root(T i ).

Data Structure for a Small Tree. It remains to describe the data structure for a tree T i , i ≥ 1. Since T i contains a small number of nodes, we can answer weighted level ancestor queries on T i using a look-up table V . V contains information about any tree with up to log n/8 nodes, such that node weights are positive integers bounded by log n/8. For any such tree T , for any node v of T , and for any integer q ∈ [1, log n/8], we store the pointer to wla(v, q) in T . There are O(2 log n/4 ) different trees T (see e.g., [START_REF] Bender | The level ancestor problem simplified[END_REF] for a simple proof); for any T , we can assign weights to nodes in less than (log n/8)! ways. For any weighted tree T there are at most (log n) 2 /64 different pairs v, q. Hence, the table V contains O(2 log n/4 (log n) 2 (log n/8)!) = o(n) entries. We need only one look-up table V for all mini-trees T i .

We can now answer a weighted level ancestor query on T i using reduction to rank space. The rank of a node u in a tree T is defined as rank(u, T ) = |{ v ∈ T | weight(v) ≤ weight(u) }|. The successor of an integer q in a tree T is the lightest node u ∈ T such that weight(u) ≥ q. The rank rank(q, T ) of an integer q is defined as the rank of its successor. Let rank(T ) denote the tree T in which the weight of every node is replaced with its rank. The weight of a node u ∈ T is not smaller than q if an only if rank(u, T ) ≥ rank(q, T ). Therefore we can find wla(v, q) in a small tree T i , i ≥ 1, as follows. For every T i we store a pointer to T i = rank(T i ). Given a query wla(v, q), we find rank(q, T i ) in O(1) time using a q-heap [START_REF] Fredman | Trans-dichotomous algorithms for minimum spanning trees and shortest paths[END_REF]. Let v ′ be the node in T i that corresponds to the node v. We find u ′ = wla(v ′ , rank(q, T i )) in T i using the table V . Then the node u in T i that corresponds to u ′ is the weighted level ancestor of v.