# CODAG: Characterizing and Optimizing Decompression Algorithms for GPUs

Jeongmin Park  UIUC
<EMAIL_ADDRESS>Zaid Qureshi  UIUC
<EMAIL_ADDRESS>Vikram Sharma Mailthody  UIUC
<EMAIL_ADDRESS>Andrew Gacek  UIUC
<EMAIL_ADDRESS>Shunfan Shao  UIUC
<EMAIL_ADDRESS>Mohammad AlMasri  UIUC
<EMAIL_ADDRESS>Isaac Gelado  NVIDIA
<EMAIL_ADDRESS>Jinjun Xiong  University at Buffalo
<EMAIL_ADDRESS>Chris J Newburn  NVIDIA
<EMAIL_ADDRESS>I-hsin Chung  IBM Research
<EMAIL_ADDRESS>Michael Garland  NVIDIA
<EMAIL_ADDRESS>Nikolay Sakharnykh  NVIDIA
<EMAIL_ADDRESS>Wen-mei Hwu  UIUC / NVIDIA
<EMAIL_ADDRESS>
###### Abstract

Data compression and decompression have become vital components of big-data
applications to manage the exponential growth in the amount of data collected
and stored. Furthermore, big-data applications have increasingly adopted GPUs
due to their high compute throughput and memory bandwidth. Prior works presume
that decompression is memory-bound and have dedicated most of the GPU’s
threads to data movement and adopted complex software techniques to hide
memory latency for reading compressed data and writing uncompressed data. This
paper shows that these techniques lead to poor GPU resource utilization as
most threads end up waiting for the few decoding threads, exposing compute and
synchronization latencies.

Based on this observation, we propose CODAG, a novel and simple kernel
architecture for high throughput decompression on GPUs. CODAG eliminates the
use of specialized groups of threads, frees up compute resources to increase
the number of parallel decompression streams, and leverages the ample compute
activities and the GPU’s hardware scheduler to tolerate synchronization,
compute, and memory latencies. Furthermore, CODAG provides a framework for
users to easily incorporate new decompression algorithms without being
burdened with implementing complex optimizations to hide memory latency. We
validate our proposed architecture with three different encoding techniques,
RLE v1, RLE v2, and Deflate, and a wide range of large datasets from different
domains. We show that CODAG provides 13.46$\times$, 5.69$\times$, and
1.18$\times$ speed up for RLE v1, RLE v2, and Deflate, respectively, when
compared to the state-of-the-art decompressors from NVIDIA RAPIDS.

##  I Introduction

Compression has become a vital component of the data storage infrastructure
due to unprecedented data size growth in big-data, high-performance computing
(HPC), and data analytics applications [1, 2, 3, 4, 5, 6]. Compression aims to
reduce the number of bits used to represent the data through various encoding
techniques specialized for different data types and patterns. Compressed data
must be first decompressed before it can be processed by compute operations.

As GPUs are increasingly used to meet the compute throughput demands of
emerging big-data workloads [7, 8, 9], data analytics pipelines typically read
compressed data from the storage devices into the GPU memory and deploy one or
more GPU kernels for decompression before performing query computation on the
data. Prior works [10, 6] have shown this approach’s benefits over storing and
moving uncompressed data. Thus, high-throughput decompression on the GPU is
becoming a key requirement to avoid performance bottlenecks in these
pipelines. Our profiling of data-dependent queries with the state-of-the-art
GPU-enabled RAPIDS pipeline [1], e.g., “What is the average fare per trip for
ride starting from Williamsburg?” on the NYC taxicab dataset [11], a popular
dataset for evaluating data analytics pipelines, shows that approximately 91%
of the total GPU execution time is consumed by the decompression kernel.

Prior works [12, 13, 14] have proposed various data-layout transformations on
compressed datasets to reduce the GPU execution time spent on decompression.
However, applying data-layout transformations on existing massive datasets of
modern big-data applications [15, 2, 10, 16, 17] is an expensive proposition,
making decompression techniques that require such transformations impractical.
Therefore, it’s much more desirable for a decompression scheme to support the
standard compressed data formats, e.g., ORC [2], and not require data-layout
transformations, a design goal that we aim to achieve in this work.

To better understand why decompression accounts for so much of the GPU
execution time, we performed a detailed profiling study of the state-of-the-
art GPU decompressor in RAPIDS. As modern compression schemes support indexed
access to multiple chunks of compressed data that correspond to fixed-size
chunks of uncompressed input data, the state-of-the-art decompressor exploits
chunk-level parallelism by assigning compressed chunks to thread blocks.
However, data dependencies still serialize the decoding operations within each
chunk. Thus, the state-of-the-art decompressor dedicates a single thread, the
leader thread, in each thread block to decode the compressed data and
broadcast the decoded information to the remaining threads in the thread
block. All threads then collaborate to execute the memory writing or copying
operations based on the decoded information. The state-of-the-art decompressor
further enhances the memory bandwidth utilization and hides more of the memory
latency by dedicating one warp in each thread block to prefetch compressed
data.

By profiling the RAPIDS decompressor, we observe that (1) the memory-
bandwidth-driven resource provisioning strategy in the state-of-the-art
decompressor in fact fails to effectively hide compute operation, memory, and
synchronization latencies; and (2) the root cause of the exposed latencies is
the insufficient number of decoding (leader) threads that cannot take full
advantage of the GPU’s parallelism and provide the GPU warp schedulers with
sufficient numbers of concurrently available instructions to hide latencies.

Based on these key observations, we propose CODAG, a modular and flexible
framework for high-throughput decompression on GPUs. In CODAG, we assign
compressed chunks to warps rather than thread blocks to drastically increase
the exploitation of the GPU’s parallelism to process more chunks at any given
time, while maintaining the ability to exploit fine-grained parallelism within
each chunk. We refer to the assignment of each chunk to a warp as _warp-level
decompression_. We believe that the warp-level decompression offers a number
of benefits in CODAG. First, warp-level decompression results in a much higher
number of active decoding threads, which helps the hardware scheduler to
effectively hide memory, synchronization, and decoding compute operation
latencies. Second, since warp-level synchronization primitives have much lower
latency and higher throughput than block-level synchronization, warp-level
decompression reduces synchronization overhead.

On top of that, we further propose (1) a novel strategy for CODAG to ensure
that all on-demand reads or writes in each warp are performed on a cache-line
granularity and are thus coalesced, making efficient use of the available GPU
memory bandwidth; (2) a clean abstraction for common operations during
decompression with the associated Application Programming Interface (API)
functions as implemented in decompression kernels. This new abstraction
provides the flexibility for programmers to easily incorporate their own
encoding techniques into the decoding section of the CODAG kernels while
leveraging the rest of the kernel to achieve the desired level of latency
tolerance and memory bandwidth efficiency.

We make the following key contributions in this paper.

  * •

We conduct a systematic characterization study on the state-of-the-art GPU
decompressor, identify the performance bottlenecks, and determine that
decompression on modern GPUs is a compute-bound process.

  * •

We present CODAG, a novel and superior GPU resource provisioning strategy for
decompressors with various optimization techniques and abstraction. The
proposed APIs can support common decompression operations so that users can
easily integrate their compression algorithms into CODAG with low programming
effort.

  * •

We demonstrate CODAG’s effectiveness and flexibility by implementing three
encoding techniques, i.e., RLE v1, RLE v2, and Deflate. Results based on the
NVIDIA A100 GPUs show that CODAG achieves 13.46$\times$, 5.69$\times$, and
1.18$\times$ speedups over a state-of-the-art decompressor, respectively.

##  II Background

In this section, we provide an overview of common compression encoding
techniques and how their decompression algorithms have been traditionally
parallelized. We then briefly explain why the resource provisioning strategies
used in the state-of-the-art GPU decompressors limit the performance and
efficiency of decompression on GPUs.

###  II-A Traditional Encoding Techniques

Generally, compressors take a contiguous sequence of bytes and attempt to
encode that sequence in as few symbols as possible. The primary measure of
merit of a compressor is the ratio of the size of the uncompressed data over
that of the compressed data, referred to as the compression ratio. To maximize
compression ratio, compressors use an arsenal of encoding techniques each of
which targets different data types and patterns. Here, we describe some of the
most common encoding techniques.

The run-length encoding (RLE) technique takes advantage of the same byte or
value being repeated in a contiguous sequence, called a  _run_. The RLE v1
encoding scheme encodes each _run_ with two values: a literal and the length
of the run.

Delta encoding provides a high compression ratio when there is a contiguous
sequence of values and the differences (deltas) between adjacent values are
small. Delta encoding stores the starting value and the number of values in
the sequence followed by the deltas to create a compressed data
representation. The RLE v2 encoding scheme leverages delta encoding in
addition to RLE to capture more data patterns and further improve the
compression ratio over RLE v1.

Dictionary-based encoding generates compact encoding when sequences of values
are repeated between different sections of the uncompressed data. Given a
sequence of bytes, the decompressor can use a dictionary to look up the
sequence to see if it has been encountered in the past. If so, the repeated
instance can be encoded with 2 values, an offset into the dictionary, and the
length of the repeated sequence.

In dictionary encoding, some sequences are repeated more than others. Thus,
some references (offset/length pairs) appear more frequently than others [18]
in the compressed stream. The Deflate encoding scheme sorts the literals and
the length-offset pairs according to their frequency into a Huffman tree [19]
so that the number of bits used for each is inversely proportional to their
frequency, thus minimizing the total number of bits used.

(a) State-of-the-art

(b) CODAG

Figure 1: Comparison of CODAG with the state-of-the art decompressor NVIDIA
RAPIDS implementation. In CODAG design, the conditional statement A checks
whether there is enough space in the input buffer to store full cacheline
data, and B checks whether threads need to write uncompressed data based on
the decoded information.

###  II-B Parallelizing Decompression

Designing high-throughput parallelizable decompression schemes is challenging
due to the nature of compressed data. First, compressors output a sequence of
variable-length blocks of compressed data, often without any alignment
guarantees (See $\S$ II-A), resulting in size irregularities. Second, the
encoding of a bit or byte in a compressed stream depends on the actions taken
for the bits or bytes before it in the stream, creating data-dependencies in
the compressed stream. Thus, decompression requires starting at the first byte
of the compressed stream and decoding the stream sequentially, limiting
parallelism during decompression.

To expose data parallelism during decompression, modern compression data
formats divide the uncompressed input into fixed-size chunks, compress each
chunk independently, and keep the bytes of the compressed chunk sequential in
memory [2, 16, 1]. During decompression, each compressed chunk is assigned to
a processing unit so that it can be decompressed independently of other
chunks. The data format further specifies metadata for offsets in the
uncompressed and/or compressed streams.

###  II-C State of the Art in Exploiting GPUs for Decompression

As GPUs become increasingly adopted by big-data applications due to their high
compute throughput and memory bandwidth, there has been growing interest in
high-performance decompression on GPUs. Generally, high-performance
decompressors on GPUs assume modern data formats that support chunking. Let’s
take a look at how the decompression kernel is mapped to GPU resources in the
state-of-the-art NVIDIA RAPIDS implementation [1, 4, 5].

As shown in Figure 1(a), after finding the offset for each compressed chunk,
RAPIDS assigns each compressed chunk to a thread block, referred to as a
decompression unit. Each decompression unit consists of a prefetching warp,
batch buffers in the fast shared memory, and decoding/writing warps. RAPIDS
opts to use a specialized prefetching warp that fetches bytes of the
compressed stream from the slow global memory in a coalesced manner to the
batch buffers. Having a specialized prefetching warp in the same thread block
as the decoding/writing warps has the benefit of allowing the batch buffers to
reside in the fast shared memory and the prefetch operations to overlap with
the decoding operations for the same chunk. Both features help to reduce the
latency of decompressing each chunk.

In the decoding warp, a leader thread is responsible for decoding the bytes
from the compressed stream, as this is a sequential process. During decoding,
the leader thread reads the data from the batch buffers. When the information
has been decoded, the leader thread synchronizes with the other threads in the
thread block and directs them to collectively perform a writing operation to
the output decompressed data stream in a coalesced manner.

##  III Characterization Study of State-of-the-art Decompressor

Understanding and optimizing decompression on GPUs is essential as the
decompression can consume approximately 91% of total GPU time on a modern data
analytics pipeline. However, none of the prior works have performed a detailed
analysis of existing state-of-the-art decompression schemes on GPUs. To this
end, we conduct an in-depth characterization study of the state-of-the-art
decompressor in NVIDIA RAPIDS [1], and identify its key performance
bottlenecks.

Setup: We use a variety of datasets across different domains as shown in Table
IV and evaluate RAPIDS’s decompressor with RLE v1 and Deflate schemes. Due to
space limitations, we only present the results from two datasets (MC0 and TPC)
but similar trends are observed for the rest. To characterize decompressor
performance, we use the system configuration described in Table III and use
the NVIDIA Nsight profiling tools to capture important metrics.

Figure 2: The peak throughput percentages and the stalled instruction
distribution of RAPIDS decompressor for RLE v1.

Key insights: From the left side of Figure 2, we know that for the RLE v1
encoding scheme, both compute and memory resources are far from being fully
utilized. As we discussed in $\S$II-C, these results point to latency exposed
due to insufficient parallelism during the decoding process as all threads in
a thread block wait on a barrier for the leader thread to complete its
decoding operation. To understand the impact of the barrier, we measured the
distribution of conditions causing an SM scheduler stall cycle, which is shown
on the right side of Figure 2. This metric provides the percentage of threads
waiting for each type of stalling condition during the cycles where an SM
scheduler cannot find a warp ready for execution.

From Figure 2, during the SM scheduler stall cycles, most (up to 80%) of the
threads are not ready because they are waiting on a barrier for the leader
thread to complete its current decoding operation. Apart from barriers, up to
20% of the instructions are waiting for the resolution of branch conditions
(Branch Resolve) when there is a lack of ready threads for the SM scheduler to
dispatch for execution. The Wait metric indicates that a warp is stalled
waiting on a fixed latency (arithmetic/logic unit) execution dependency.
Moreover, some threads do not receive any work during the writing stage as run
lengths can be smaller than the size of the thread block, further degrading
the achievable performance.

For Deflate, unlike the traditional view that it is memory bandwidth bound [1,
13], our profiling results in Figure 3 show that its memory bandwidth
utilization is less than 25%, which indicates that it is not limited by memory
bandwidth but is compute limited. In contrast to RLE v1, Deflate has a much
higher utilization of compute throughput (up to 65%) than RLE v1 because of
its high computational complexity in dictionary look-up and Huffman tree
traversal. We need more detailed metrics and analysis in order to discern
whether the performance is limited by a particular type of compute pipeline,
memory latency, or compute latency.

We use the profiler to analyze the utilization of the different GPU compute
pipelines for Deflate decompression. The right side of Figure 3 shows that the
arithmetic logic units (ALU) are up to 55% utilized while the fused
multiply/add and load-store units are up to 35% and 20% utilized. This
modestly high degree of ALU and FMA utilization mainly comes from the fact
that during decoding, the leader thread executes a large number of arithmetic
instructions for every byte to identify the next action for generating the
output for the byte.

Nevertheless, all these results point to memory latency and/or compute latency
as the main limiter of Deflate performance. Due to the high computational
complexity of Deflate decoding process, it takes a long time for the leader
thread to complete the decoding of each byte while most threads in the thread
block are waiting for the decoding thread. During this time, there are not
enough decoding threads for the SM scheduler to tolerate the compute pipeline
latency. As a result, the execution speed of Deflate is limited by compute
pipeline latency in addition to barrier stalling.

Figure 3:  The compute/memory peak throughput percentages and compute unit
utilization of RAPIDS decompressor for Deflate.

Summary and take-aways: From the profiling results, we see that the
decompression speed of RLE v1 is limited by barrier latency and that Deflate
is limited by barrier stall latency as well as the fixed latency of compute
pipelines. Neither the large number of writing threads nor the specialized
prefetching warps help address this limitation. Any further improvement in
decompression speed for both RLE v1 and Deflate would require more
simultaneously running decoding threads (more chunks being decoded in
parallel), which is an important design goal of the CODAG architecture.

##  IV  CODAG Design

The goal of CODAG is to provide a flexible framework for high-throughput
decompression on GPUs without requiring developers to implement sophisticated,
GPU-specific optimizations. CODAG’s architecture, resource provisioning
strategy, and kernel APIs enable developers to easily incorporate their
encoding techniques of choice into CODAG and rely on the rest of the framework
to achieve high decompression speed. On the one hand, as discussed in Section
III, the key to improved decompression speed is to increase the number of
simultaneously running decompression threads so that the SM scheduler can
better hide the synchronization/memory/compute-pipeline overheads.

On the other hand, it is also critical to maintain the coalesced access
patterns in reading uncompressed data, reading dictionary value sequences, and
writing decompressed data so that the memory access cost does not become
excessive and creates new bottlenecks. Thus, the CODAG architecture must
address three key challenges in providing an efficient and effective solution
that can be used with a wide variety of encoding techniques:

  * •

CODAG must provision GPU execution resources so that many threads can
simultaneously decode more chunks in each SM to effectively tolerate latencies
and efficiently utilize hardware resources.

  * •

As the naïve on-demand reading and writing operations can make poor use of
memory bandwidth and exhibit excessive latency, CODAG must optimally coalesce
these accesses to minimize memory access overhead in spite of the irregular
consumption of compressed data and sporadic generation of decompressed data by
the decoder.

  * •

As the sophisticated optimizations for efficiently using memory bandwidth and
tolerating latencies require significant programming skills and efforts, the
optimizations in the CODAG framework must be generalized into reusable
primitives for a wide range of encoding techniques.

###  IV-A CODAG System Overview and Abstraction

Figure 1(b) shows an overview of the CODAG framework. CODAG uses warps as
decompression units instead of the larger thread blocks used by the state-of-
the-art decompressor, enabling CODAG to simultaneously decompress many
compressed chunks. This enables the hardware scheduler to overlap instruction
execution from several independent warps and efficiently hide compute pipeline
and memory latencies. In CODAG, all threads in each warp perform the decoding
operations as well as amalgamate to read uncompressed data and write
decompressed data. Thus, CODAG avoids the high barrier-synchronization
overhead that the state-of-the-art decompressor incurs.

In CODAG, each warp executes a decompression GPU kernel that implements a
control flow of three main operations: sequential decoding, coalesced reading,
and coalesced writing, as shown in Figure 1(b). All threads in a warp follow
the same control flow to avoid using any broadcast synchronization operations.
The control flow consists of a main loop for the decoding process and two
conditional statements for on-demand reading and writing operations. During
decoding, all threads sequentially execute the same decoding process in-
compliance with the serial nature of the compressed data stream. In our
current framework implementation, the sequential decoding code for different
combinations of pertinent encoding techniques can be easily incorporated into
the kernel as a CUDA device function.

For coalesced on-demand reading operations, all threads in the warp
collaborate to fetch the next cache-line section (e.g., 128 bytes in the A100
GPUs) of the compressed input data stream in a coalesced manner into an input
buffer. Having a buffer for the fetched bytes minimizes the number of fetches
in spite of the potentially irregular patterns of compressed data consumption
during decoding, which helps to reduce the cost of accessing the compressed
data. In comparison with the RAPIDS decompressor, CODAG does not have a
specialized prefetch warp in each thread block to asynchronously prefetch
input data while the thread block’s leader thread is decoding. However, with
many independent warps being simultaneously executed, CODAG can still tolerate
the memory latency.

When a writing operation is required during decoding, all threads in the warp
collaborate to execute the coalesced on-demand writing of a cache-line of the
uncompressed output data to the output buffer. Compared to the RAPIDS
decompressor, this design at first glance may not seem to offer as much memory
parallelism when writing into each output chunk, as RAPIDS uses thread blocks
with up to 1024 threads available for writing into each output chunk.

However, the benefit with CODAG is that warp-level barriers are significantly
less expensive than their block-level counterparts. Furthermore, with the same
thread block size, CODAG supports decompressing and writing up to $32\times$
more independent chunks in parallel. The increased data parallelism can easily
compensate for the loss of memory-level parallelism within a chunk.

As the reading pattern is the same across all encoding techniques, the on-
demand reading operation code is portable to all (combinations of) encoding
techniques. However, the on-demand writing operations can vary depending on
the encoding techniques and need to be specialized by decompressor developers.
To shield the decompressor developers from the complexity of GPU algorithmic
optimizations that are critical to high-performance writing operations, CODAG
provides several major types of optimized writing primitives to support a wide
range of popular encoding techniques.

Figure 4:  Illustration of the issued instructions between the baseline and
CODAG in a steady state.

###  IV-B CODAG Meta-data and Performance APIs

CODAG provides an infrastructure for realizing the warp-level decompression as
well as the optimized reading and writing operations to lower the programming
effort when using CODAG for different encoding techniques. To isolate the
complexity of GPU optimizations from decompressor developers, CODAG provides
input_stream and output_stream abstractions that transparently manage the
synchronization, coalescing and metadata for the on-demand reading and writing
operations. Developers can simply call the member functions, shown in Table I
and Table II, as needed to read and write data efficiently.

To support the highly efficient on-demand reading and writing operations,
input_stream and output_stream internally manage meta-data such as offsets for
the reading and writing pointers, information of input buffer, or fetched
bites in the top entry of the input buffer. With these provided reading and
writing primitives, CODAG only requires small changes to the sequential
decoding code as the memory access optimizations are hidden via the internal
functions. In the next few sections, we will discuss the CODAG decompression
scheme and its optimizations more in detail.

TABLE I: CODAG APIs for input_stream.  Member Functions |  Description
---|---
fetch_bits(int n) |  Fetch the next n bits in the compressed stream.
peek_bits(int n) |  Peek at the next n bits in the compressed stream.
TABLE II: CODAG APIs for output_stream.  Member Functions |  Description
---|---
write_byte(int8_t b) |  Write a literal (byte) to the output buffer
write_run(T init, size_t len, T delta) |  Write a run based on an initial value of type T, a length, and a delta.
memcpy(int off, int len) |  Memcpy operation for dictionary-based encoding with the offset and length.

###  IV-C Warp Level Decompression

The large granularity of the decompression units in the state-of-the-art
limits the scheduler’s ability to hide latency, because the vast majority of
threads end up spending most of their execution time waiting for a few leader
threads to finish sequential decoding. The left side of Figure 4 is a cartoon
illustration of how the instructions are issued during the decompression for
the state-of-the-art decompressor. We assume an SM with two schedulers and an
occupancy capacity of 4 warps for the example. In this tiny example, assume
that the state-of-the-art has a decompression unit of a thread block that
consists of two warps. Thus, only two independent compressed chunks (green and
orange) can be assigned to a single SM in the state-of-the-art where one of
the two warps is specialized to prefetch data from global memory into the
buffers, and both warps work together to perform writes to global memory after
the leader has decoded a symbol.

We assume that during the first part of the timeframe, each scheduler is
issuing instructions of the warp containing the leader thread of one of the
decompressor units and, for simplicity, that the batch buffer is full so there
is no prefetch activity. Since there is effectively only one thread active in
each scheduler, the latency of the function units is fully exposed, shown as
pipeline bubbles between the decoding operations. After the decoding
operations are complete, all threads synchronize, incurring the
synchronization pipeline bubbles before the write operations. The SM execution
resources are underutilized due to the pipeline bubbles.

The right side of Figure 4 shows the execution timing of CODAG. CODAG uses
warps as a decompression unit to allow more chunks (green, orange, purple,
blue) to be decompressed simultaneously. This gives the hardware scheduler
more independent warps to work with. Assume that the green and purple warps
are assigned to scheduler 1 whereas the orange and blue are assigned to
scheduler 2. During the first part of the time frame, the decoding operations
of the green and blue warps can be issued alternatively to tolerate each
other’s latency. Thus, the SM execution resources are much better utilized. In
reality, the reading, decoding, and writing operations of all the warps will
be mixed together to more fully utilize the memory bandwidth than what is
shown in Figure 4. That is, since there are no inter-warp dependencies in
CODAG, the SM’s scheduler has a wider pool of useful instructions to schedule
from and thus can better reduce pipeline bubbles and hide latencies.
Furthermore, there is no more synchronization overhead before the writing
operations.

###  IV-D CODAG Decoding Scheme

As the prefetch warp is removed in CODAG’s warp-level decompression scheme,
all threads in the warp collectively execute on-demand reading operations
during decoding whenever there are not enough fetched bytes. To do this, the
decoding state needs to be saved before the reading operation and restored
after the reading operation to enable the leader thread to continue decoding
after the reading operations.

To avoid the above-mentioned overhead during the transitions between the
decoding operations and the reading/writing operations, CODAG uses an all-
thread decoding technique by activating all threads in the warp to execute the
same serial decoding operations on the same chunk. In this decoding scheme,
all threads have a local copy of decoded information as they execute the same
decoding process thus broadcasts are no longer needed. As for the reading
operations, all threads in the warp only need to perform a low-cost warp-level
synchronization to ensure a consistent view of the input buffer before reading
as newer GPU architectures [20, 21] do not guarantee lock-step execution of
the threads in a warp. After the reading operation, all threads in the warp
can return naturally to continue the decoding operation, eliminating the need
for explicitly saving and restoring the decoding state.

Since the decoding operation is redundantly performed by all threads in a
warp, one might take a pessimistic view of the all-thread decoding technique.
Thus, we use a micro-benchmark to compare the ALU compute throughput of the
single-thread decoding and the all-thread decoding. For this micro-benchmark,
the number of arithmetic operations executed per global memory access is
varied from 1 to 100,000 to capture the compute throughput at different
compute intensity levels. The result shows that the difference in compute
throughput of the two decoding techniques never exceeds 0.1%. Thus, the
additional arithmetic operations executed due to the all-thread decoding
technique do not degrade the achieved compute throughput.

Algorithm 1 On-demand Reading With Shared Input Buffer

1: if len - counter $<$ 128 then

2: barrier()

3: read_idx = offset \+ lane * 4

4: enq_idx = (head + counter + lane) % len

5: Buff[enq_idx] = Input[read_idx]

6: counter += 128

7: barrier()

8: end if

9: req_val = Buff[$head$]

10: head = (head + 1) % len

11: counter–

12: return req_val

###  IV-E Coalesced On-Demand Reading Operation

To maximize the benefits of the warp-level decompression scheme, CODAG also
provides the APIs for highly efficient on-demand reading operations since it
is also critical to maintain the coalesced access patterns in reading
compressed data to optimize GPU memory bandwidth utilization. Thus, CODAG
implements coalesced on-demand reading operations behind the input_stream
abstraction.

CODAG’s on-demand reading operation with its default configuration is
illustrated with the pseudocode in Algorithm 1. By default, CODAG deploys an
input buffer in shared memory, to temporarily store the fetched bytes,
minimizing the number of global memory accesses. The size of the input buffer
is at least twice the cacheline size so that entire cachelines can be easily
enqueued in the input buffer. An input buffer also has a head pointer for
threads to determine where to read from in the input buffer, and a length to
keep track of the number of bytes left in the buffer.

For a given on-demand reading operation, all threads in a warp first
collectively refill the shared input buffer with bytes fetched from a
cacheline. Each thread in the warp fetches 4-bytes after a warp-
synchronization barrier to ensure coalesced memory accesses. Once the input
buffer is refilled, threads update the offset into the input in the
input_stream instance. Then, the threads read the requested bits from the
input buffer using the head pointer and update the internal state of its
input_stream for the next on-demand reading operation.

Using Registers: The shared memory usage for the input buffer might limit the
applicability of CODAG in applications with high shared memory utilization. To
address this, CODAG provides a configuration where the input buffer is
implemented with the threads’ local registers instead of the shared memory. In
this configuration, each thread uses two 32-bits registers to implement the
input buffer.

Similar to the default configuration, the input buffer is first refilled when
there is enough space to store a complete cacheline. This is done by
leveraging double buffering where two sets of 32 local registers across the
warp operate as a double buffer. Since the data is stored in a local register,
CODAG exploits warp-wise broadcast when fetching requested bytes from the
input buffer. For this operation, the head pointer is first used to identify
which thread’s local register is holding the requested bytes so that the
requested bytes can be broadcast by the appropriate thread to all other
threads in a warp.

###  IV-F Coalesced On-demand Writing Operation

Algorithm 2 Memcpy(offset, len)

1: if !is_4B_aligned(write_ptr) then

2: pad_bytes = (4 - (write_ptr % 4)) % 4

3: byte_memcpy(offset, pad_bytes)

4: end if

5: barrier()

6: num_itr = ceil((len - pad_bytes) / 128)

7: for (int i = 0; i ¡ num_itr; i++) do

8: read_idx = ceil((read_ptr + lane*4) / 4)

9: r1 = out[read_idx-1]

10: r2 = out[read_idx]

11: rdata = funnel_shift(r1, r2, read_ptr % 4)

12: out[(write_ptr + 4*tid) / 4] = rdata

13: read_ptr, write_ptr += 32

14: barrier()

15: end for

Similar to the on-demand reading operations, CODAG supports highly efficient
coalesced on-demand writing operations to reduce the cost of memory accesses.
Ideally, all threads in a warp should collectively execute an on-demand
writing operation of a full cacheline and then synchronize. However, as the
required writing operations can vary based on the encoding techniques, CODAG
provides optimized writing primitives to support the most common encoding
techniques: 1) writing a literal value, 2) writing a run with delta (delta is
0 for a simple run), and 3) memory copy (memcpy).

Writing a literal value: As this operation only requires a few bytes to write
to the output buffer, only one thread executes this operation.

Writing a run: This operation is commonly used for RLE-based encoding
techniques. For this operation, each thread independently computes its
decompressed value based on the initial value, delta, and thread index and
then writes it to the output buffer at the appropriate offset to maximize
memory utilization.

Memcpy: Common dictionary-based encoding techniques compress data with a
length specifying the number of bytes to copy from the dictionary to the
output buffer and an offset into the output buffer from the last byte written
before the current memcpy starts. Unlike a regular memcpy, using a warp to
perform a memcpy operation in a dictionary-based compression algorithm has a
couple of challenges.

First, as the dictionary is implicit in the previously written output, the
reading and the writing pointer both point to the same output buffer. Second,
common algorithms only provide byte alignment for both the input and the
output. These characteristics make it difficult to exploit the warp’s
parallelism in memcpy, especially when we would like each thread to copy more
than 1-byte (e.g. 4-bytes) at a time for better coalescing, reduced executed
memory instructions, and improved GPU cacheline and memory bandwidth
utilization.

To address these challenges, we design CODAG’s fast memcpy operation, shown in
Algorithm 2. For a given memcpy operation, the threads first fetch the writing
pointer from output_stream and check whether it is 4-byte aligned (line 1 in
Algorithm 2). If it is not aligned, the threads collectively copy up to 3
bytes to guarantee 4-byte alignment (line 3). Once the writing pointer is
aligned, the warp synchronizes to ensure that the copied bytes are visible to
the whole warp. Then, the warp executes copy operation in iterations to cover
the entire memcpy length, with each thread writing one 4-byte value to the
output each iteration. Each thread independently calculates its reading
pointer for a given iteration based on its warp lane index (line 8). Since the
reading pointer is not guaranteed to be 4-byte aligned, each thread reads two
consecutive 4-byte elements (line 9-10), then bit-shifts and concatenates them
to generate the correct 4-byte output element (line 11) to write to its
assigned offset (line 12).

Note, in each iteration, each thread in the warp independently reads two
4-byte elements written in previous iterations and writes one complete 4-byte
element. We found that using 4-byte granularity for the memcpy loop body
provided the best trade-off between overhead to align the writing pointer, the
amount of parallelism exposed, and cacheline utilization. All threads in a
warp synchronize at the end of each loop iteration, preventing race conditions
between iterations.

Moreover, with dictionary compression, it is possible that the length of a
memcpy is larger than the offset. In this special case, the byte sequence from
the offset till the end of the output (before the memcpy executes) is to be
appended to the output, possibly multiple times. Just like in Algorithm 2, we
first align the writing pointer to 4-bytes and then assign each thread one
4-byte output element to generate per memcpy loop iteration. However, in this
special case, each 4-byte output element is created by combining two
consecutive 4-byte elements in a fixed circular window of 4-byte elements,
starting from the 4-byte element pointed to by the offset until the last
4-byte element of the output buffer. Thus, each thread can independently use
modulo-arithmetic to find the two consecutive 4-byte elements in the specified
window, and bit-shift and concatenate them to generate its output element each
iteration. We don’t show the algorithm listing for this case due to space
constraints.

##  V Evaluation

Our evaluation demonstrates that: (1) CODAG framework improves the performance
of decompression by efficiently utilizing the hardware resources provided by
the GPU. CODAG provides 13.46$\times$, 5.69$\times$, 1.18$\times$
decompression throughput speed up (geomean) on ORC RLE v1, ORC RLE v2, and
Deflate compared to the state-of-the art implementation in NVIDIA RAPIDS [1].
(2) A decompressor implemented in the CODAG framework can rival the state-of-
the-art GPU implementation without application-specific optimizations.

###  V-A Experiment Setup

We evaluate CODAG using three different encoding techniques, namely RLE v1 and
RLE v2 from the Apache ORC file format, and Deflate. These encoding schemes
are chosen due to wide applicability but other encodings can also be enabled.
We compare CODAG and the state-of-the-art NVIDIA RAPIDS baseline on an
environment shown in Table III. We lock the GPU’s clock frequency to peak
before running the experiments.

###  V-B Evaluation Datasets

For each algorithm, we evaluate CODAG and the baseline on seven real-world
sample datasets curated from four industrial application domains as shown in
Table IV. Each dataset has unique properties enabling us to evaluate the
proposed scheme on not only different types of encoding techniques but also
different data types and patterns. For example, MC0 and MC3 have long run-
lengths while TPC and TPT contains many repeated patterns. CD2, and TC2 follow
power law distributions for the data. HRG is human genome dataset and has
repeated text pattern containing A,T,C,G and N characters.

For the baseline, we use the official data management tools provided by ORC
[2] to compress data with RLE v1 and RLE v2 algorithms. To compress the data
with deflate, we use the zlib library [22] with a compression level of 9, the
highest compression level. The chunk size for the original data is fixed to be
128KB for both CODAG and the baseline.

TABLE III: Configuration used to evaluate CODAG. Configuration |  Specification
---|---
CPU |  AMD EPYC 7702 64-Core Processor
Memory |  1TB DDR4
GPU 1 |  NVIDIA Tesla V100 HBM2 32GB
GPU 2 |  NVIDIA A100 HBM2 40GB
TABLE IV: Real-world dataset used for evaluating CODAG.  Dataset |  Category |  DType |  Size (GB)
---|---|---|---
Mortgage Col 0 (MC0) [23] |  Analytics |  uint_64 |  4.86
Mortgage Col 3 (MC3) [23] |  Analytics |  fp32 |  2.43
NYC Taxi Passenger Count (TPC) [11] |  Analytics |  int_8 |  3.07
NYC Taxi Payment Type (TPT) [11] |  Analytics |  char |  7.41
Criteo Dense 2 (CD2) [17] |  Recommenders |  uint_32 |  0.73
Twitter COO Col 1 (TC2) [24] |  Graph |  uint_64 |  5.47
Human Reference Genome (HRG) [25] |  Genomics |  char |  3.1
TABLE V: Compression Ratios and average compressed symbol length when the datasets are compressed by RLE v1, RLE v2, and Deflate. The compressed symbol length for RLE v2 is excluded due to various compressed types for RLE v2.  Dataset | Compression Ratio | Avg Comp Sym Len
---|---|---
RLE v1 |  RLE v2 |  Deflate |  RLE v1 |  Deflate
MC0 |  0.023 |  0.022 |  0.017 |  29.7 |  81.3
MC3 |  0.038 |  0.039 |  0.015 |  40.5 |  156.6
TPC |  0.867 |  0.637 |  0.119 |  1.00 |  10.1
TPT |  1.41 |  0.99 |  0.042 |  1.00 |  32.6
CD2 |  0.286 |  0.308 |  0.625 |  20.9 |  5.89
TC2 |  0.087 |  0.075 |  0.0172 |  34.3 |  67.3
HRG |  0.975 |  0.972 |  0.305 |  1.00 |  7.76
Figure 5: Stalled instruction distribution comparison between CODAG and RAPIDS
baseline. Figure 6: Comparing the compute/memory peak throughput of CODAG and
RAPIDS. Comp/Mem stand for compute/memory throughput percentage.

###  V-C Latency Analysis

To analyze how effectively CODAG hides latency and utilizes resources, we
measured the stalled instruction distribution and compute/memory throughput
percentages of CODAG and the baseline. For the stalled instruction
distribution, we measured the number of instructions stalled due to
synchronization (SB) and Math Pipeline Throttle (MPT). Measuring SB helps to
understand thread occupancy during decompression, and measuring MPT helps us
to understand the pressure on computing units as MPT means the instructions
are stalled to wait on a specific oversubscribed math pipeline.

In this experiment, MC0 and TPC datasets are used for the measurement as they
have drastically different compression ratios. As shown in Figure 5, we
observe a dramatic difference in SB percentages between the baseline and
CODAG. The lower SB percentage shows that there are more active threads during
decompression. Moreover, we observe substantially more instructions are
stalled due to MPT in CODAG, showing CODAG utilizes more compute resources and
shifts decompression from latency-bound to compute-bound on GPUs.

Figure 6 shows that CODAG achieves a higher memory throughput percentage
compared to the baseline. This is because CODAG optimally utilizes memory
bandwidth by leveraging coalesced on-demand reading/writing operations and
enabling more memory instructions to be executed in parallel by higher ILP.

Figure 7: Comparing decompression throughput of CODAG and RAPIDS baseline on
A100 GPU. Figure 8: CODAG speedup without and with a prefetching warp when
compared to RAPIDS in the latest two generations of NVIDIA GPUs. First two
bars represent the speedup compared to RAPIDS on A100, and the baseline for
the third bar is RAPIDS on V100.

###  V-D Decompression Throughput

We now evaluate the decompression throughput (output bytes per second) of
CODAG and RAPIDS on three encoding techniques. Figure 7 shows the
decompression throughput of each encoding technique for both CODAG and the
baseline, and Figure 8 shows the CODAG speedup over the baseline decompressor.
As shown in Figure 7, CODAG achieves 38.07 GBps, 26.87 GBps, and 51.96 GBps
geo-mean decompression throughput while the state-of-the-art decompressor
achieves 2.83 GBps, 4.72 GBps, and 44.18 GBps. Thus, CODAG achieves
13.46$\times$, 5.69$\times$, and 1.18$\times$ geo-mean speed ups compared to
the baseline as shown in Figure 8.

CODAG sees significantly higher decompression throughput in MC0 and MC3
datasets. This is because highly compressable datasets amplify the effective
memory bandwidth, requiring fewer bytes to decompress the chunk and enabling
more active threads during writing operations. In our datasets, the mortgage
datasets (MC0 and MC3) are highly compressable datasets, thus showing the
highest decompression throughput.

We note that CODAG achieves substantially higher speedups over the baseline
for RLE-based encoding techniques than Deflate. This is because with the
state-of-the-art the average run length needs to be longer for RLE v1 for best
thread occupancy; else, most threads become inactive during writing
operations.

###  V-E Single-Thread Decoding vs All-Thread Decoding

Next, we present the benefit of performing all-thread decoding over single-
thread decoding for RLE v1 (RLE v2 has similar behavior) and Deflate
decompression. When we compare the decompression throughput of the two
decoding techniques, the result shows the CODAG with the all-thread decoding
technique achieves 1.17$\times$ and 1.19$\times$ higher geo-mean decompression
throughput than the single-thread decoding technique. As discussed in $\S$
IV-D, this shows CODAG successfully reduces the number of broadcast operations
while maintaining the achieved compute throughput of ALUs.

###  V-F Impact of Finer Decompression Unit

This section presents an ablation study to show that (1) reducing the number
of threads responsible for decoding a chunk is key to CODAG’s performance
gains, and (2) CODAG’s warp-level decompression, i.e., without a prefetching
warp, further improves performance.

For this evaluation, we add a prefetch warp to CODAG, such that two warps are
scheduled per chunk - one warp to prefetch compressed input data and one warp
to perform the decompression. This makes the CODAG architecture more similar
to the baseline except that the baseline allocates more threads to decompress
each chunk - 1024 threads for RLE v1 and RLE v2 and 128 threads for Deflate.

Figure 8 shows CODAG with a prefetching warp is able to achieve 7.10$\times$,
4.33$\times$, and 1.02$\times$ higher geo-mean decompression throughputs
compared to the baseline for the three respective encoding techniques. Based
on the result, finer decompression units enable the GPU hardware to tolerate
compute operation and memory latency more effectively.

However, the performance gains for each of the three encoding techniques is
lower than the performance gains achieved by CODAG’s single-warp decompression
unit, i.e. without a prefetching warp. This is because the decompression
process is not bounded by memory bandwidth, thus manually allocated hardware
resources for prefetching to tolerate memory latency limits hardware resource
utilization. Additionally, CODAG minimizes the cost of memory access latency
by coalesced on-demand memory operations. Thus, the benefits from reallocating
the resources used for prefetching to decoding out-scales the benefits from
prefetching units.

###  V-G Impact of Hardware Resources

In the previous sections, we establish that CODAG offers superior
decompression bandwidth on A100 GPU compared to the baseline. In this section,
we constrain the CODAG and baseline design by downgrading the GPU to NVIDIA
Volta V100. From this experiment, we analyze the impact of hardware resources
on CODAG and the scalability of CODAG and the baseline design.

As shown in Figure 8, CODAG achieves 11.19$\times$, 4.39$\times$, and
1.10$\times$ geo-mean speed up compared to the baseline on V100 while it
achieves 13.46$\times$, 5.69$\times$, and 1.18$\times$ geo-mean speed up
compared to the baseline on A100. The impact of compute under-utilization
becomes more critical as newer GPU provides higher compute resources and
memory bandwidth. Thus, we observe CODAG scales better than the baseline with
the available hardware resources as it more efficiently exploits parallelism
provided by the GPUs.

##  VI Related Work

Parallel compression/decompression with CPUs: Multiple efforts to parallelize
compression and decompression algorithms on CPUs have been proposed in the
past [26, 27, 28]. pigz decompression is inherently hard to parallelize due to
data dependency between the variable-length compressed blocks [29, 26]. Zstd
[28] and pbzip2 [27] implementations add decompressed size as part of their
headers to exploit parallelism. However, such coarse-grain parallelism across
the compressed blocks is insufficient to fully utilize GPUs [13].

GPU Decompression: Several prior works have proposed implementing efficient
compression or encoding algorithms for GPUs [30, 4, 31, 32, 13, 33, 34, 35,
12, 14]. The prior works [13, 12, 14, 30, 33], define their own file formats
to achieve high decompression throughput for Deflate algorithm on GPU and
propose to split the input data into smaller data blocks. However, compared to
CODAG, the datasets should go through expensive pre-processing steps to
recompress the entire datasets into their own data format. Other prior works
[34, 35, 32, 31] leverage algorithm-specific optimizations, limiting
flexibility.

##  VII Conclusion

In this work, we introduced CODAG, a flexible framework for high throughput
decompression on GPUs without requiring decompressor developers to implement
sophisticated, GPU-specific optimizations. CODAG exploits warp-level
decompression and highly efficient on-demand reading and writing operations to
effectively hide latencies. Moreover, CODAG completely removes the broadcast
operations by dedicating all threads to participate in decoding. Our
evaluations show that CODAG can achieve 13.46$\times$, 5.69$\times$, and
1.18$\times$ speed up compared to the state-of-the-art decompressor for RLE
v1, RLE v2, and Deflate encoding techniques.

## References

  * [1] “CUDA RAPIDS: GPU-Accelerated Data Analytics and Machine Learning,” https://developer.nvidia.com/rapids, 2020.
  * [2] A. S. Foundation, “Orc file format specification,” https://orc.apache.org/specification/, 2021.
  * [3] W. McKinney and the Pandas Development Team, “pandas: powerful python data analysis toolkit,” https://pandas.pydata.org/docs/, 2022.
  * [4] “A library for fast lossless compression/decompression on the gpu,” https://github.com/NVIDIA/nvcomp, 2020.
  * [5] “Accelerating Lossless GPU Compression with New Flexible Interfaces in NVIDIA nvCOMP,” https://developer.nvidia.com/blog/accelerating-lossless-gpu-compression-with-new-flexible-interfaces-in-nvidia-nvcomp/, 2022\.
  * [6] “Boosting Data Ingest Throughput with GPUDirect Storage and RAPIDS cuDFP,” https://developer.nvidia.com/blog/boosting-data-ingest-throughput-with-gpudirect-storage-and-rapids-cudf/, 2022\.
  * [7] D. Wang, F. Zhang, W. Wan, H. Li, and X. Du, “Finequery: Fine-grained query processing on cpu-gpu integrated architectures,” in _2021 IEEE International Conference on Cluster Computing (CLUSTER)_ , 2021, pp. 355–365.
  * [8] S. Floratos, M. Xiao, H. Wang, C. Guo, Y. Yuan, R. Lee, and X. Zhang, “Nestgpu: Nested query processing on gpu,” in _2021 IEEE 37th International Conference on Data Engineering (ICDE)_ , 2021, pp. 1008–1019.
  * [9] Y. Liu, J. Wang, and S. Swanson, “Griffin: Uniting cpu and gpu in information retrieval systems for intra-query parallelism,” _SIGPLAN Not._ , vol. 53, no. 1, p. 327–337, feb 2018. [Online]. Available: https://doi.org/10.1145/3200691.3178512
  * [10] J. Li, H.-W. Tseng, C. Lin, Y. Papakonstantinou, and S. Swanson, “HippogriffDB: Balancing I/O and GPU Bandwidth in Big Data Analytics,” _Proceedings of the VLDB Endowment_ , vol. 9, no. 14, p. 1647–1658, Oct. 2016.
  * [11] The City of New York, “TLC Trip Record Data,” https://www1.nyc.gov/site/tlc/about/tlc-trip-record-data.page.
  * [12] S. Funasaka, K. Nakano, and Y. Ito, “Adaptive loss-less data compression method optimized for gpu decompression,” _Concurrency and Computation: Practice and Experience_ , vol. 29, no. 24, p. e4283, 2017.
  * [13] E. Sitaridi, R. Mueller, T. Kaldewey, G. Lohman, and K. A. Ross, “Massively-parallel lossless data decompression,” in _2016 45th International Conference on Parallel Processing (ICPP)_ , 2016, pp. 242–247.
  * [14] A. Shanbhag and B. W. Yogatama, “Tile-based lightweight integer compression in gpu,” _Proceedings of the 2022 International Conference on Management of Data_ , pp. 1390–1403, 2022.
  * [15] Z. Qureshi, V. S. Mailthody, I. Gelado, S. W. Min, A. Masood, J. Park, J. Xiong, C. Newburn, D. Vainbrand, I.-H. Chung, M. Garland, W. Dally, and W. mei Hwu, “Bam: A case for enabling fine-grain high throughput gpu-orchestrated access to storage,” 2022.
  * [16] A. S. Foundation, “Parquet file format specification,” https://parquet.apache.org/docs/file-format/, 2022.
  * [17] Criteo AI Lab, “Criteo 1TB Click Logs dataset ,” https://ailab.criteo.com/download-criteo-1tb-click-logs-dataset/.
  * [18] P. Deutsch, “Rfc8478: Zstandard compression and the application/zstd media type,” https://tools.ietf.org/html/rfc1951, 1996.
  * [19] D. A. Huffman, “A method for the construction of minimum-redundancy codes,” _Proceedings of the IRE_ , vol. 40, no. 9, pp. 1098–1101, 1952.
  * [20] “NVIDIA Tesla V100 GPU Architecture Whitepaper,” https://images.nvidia.com/content/volta-architecture/pdf/volta-architecture-whitepaper.pdf, 2017\.
  * [21] “NVIDIA Tesla A100 Tensor Core GPU Architecture,” https://www.nvidia.com/content/dam/en-zz/Solutions/Data-Center/nvidia-ampere-architecture-whitepaper.pdf, 2020\.
  * [22] J. G. G. Roelofs and M. Adler, “A massively spiffy yet delicately unobtrusive compression library,” https://www.zlib.net/, 2022.
  * [23] F. Mae, “Fannie mae single-family loan performance data,” https://docs.rapids.ai/datasets/mortgage-data.
  * [24] H. Kwak, C. Lee, H. Park, and S. Moon, “What is Twitter, a Social Network or a News Media?” in _Proceedings of the 19th International Conference on World Wide Web_ , ser. WWW ’10, Raleigh, NC, 2010.
  * [25] N. C. for Biotechnology Information, “Human reference genome (grch38-fasta),” https://www.ncbi.nlm.nih.gov/genome/guide/human/.
  * [26] M. Adler, “Parallel gzip,” https://zlib.net/pigz/, 2010.
  * [27] J. Gilchrist, “Parallel bzip2 file compressor,” https://linux.die.net/man/1/pbzip2, 2015.
  * [28] Y. Collet and M. Kucherawy, “Rfc8478: Zstandard compression and the application/zstd media type,” https://tools.ietf.org/html/rfc8478 and https://facebook.github.io/zstd/, 2015.
  * [29] H. Jang, C. Kim, and J. W. Lee, “Practical speculative parallelization of variable-length decompression algorithms,” in _Proceedings of the 14th ACM SIGPLAN/SIGBED Conference on Languages, Compilers and Tools for Embedded Systems_ , ser. LCTES ’13. New York, NY, USA: Association for Computing Machinery, 2013, p. 55–64.
  * [30] T. Lloyd, K. Barton, E. Tiotto, and J. N. Amaral, “Run-length base-delta encoding for high-speed compression,” in _Proceedings of the 47th International Conference on Parallel Processing Companion_ , ser. ICPP ’18. New York, NY, USA: Association for Computing Machinery, 2018.
  * [31] R. A. Patel, Y. Zhang, J. Mak, A. Davidson, and J. D. Owens, “Parallel lossless data compression on the gpu,” in _2012 Innovative Parallel Computing (InPar)_ , 2012, pp. 1–9.
  * [32] A. Ozsoy, M. Swany, and A. Chauhan, “Optimizing lzss compression on gpgpus,” _Future Generation Computer Systems_ , vol. 30, pp. 170 – 178, 2014.
  * [33] M. A. O’Neil and M. Burtscher, “Floating-point data compression at 75 gb/s on a gpu,” in _Proceedings of the Fourth Workshop on General Purpose Processing on Graphics Processing Units_ , ser. GPGPU-4. New York, NY, USA: Association for Computing Machinery, 2011\.
  * [34] A. Ozsoy, M. Swany, and A. Chauhan, “Pipelined parallel lzss for streaming data compression on gpgpus,” in _2012 IEEE 18th International Conference on Parallel and Distributed Systems_ , 2012, pp. 37–44.
  * [35] M. Olano, D. Baker, W. Griffin, and J. Barczak, “Variable bit rate gpu texture decompression,” in _Proceedings of the Twenty-Second Eurographics Conference on Rendering_ , ser. EGSR ’11. Goslar, DEU: Eurographics Association, 2011, p. 1299–1308.


# Modulated phases in a spin model with Dzyaloshinskii-Moriya interactions

William de Castilho* and S. R. Salinas
Instituto de Física, Universidade de São Paulo,
São Paulo, SP, Brazil

###### Abstract

We analyze the phase diagram of an elementary statistical lattice model of
classical, discrete, spin variables, with nearest-neighbor ferromagnetic
isotropic interactions in competition with chiral interactions along an axis.
At the mean-field level, we show the existence of paramagnetic lines of
transition to a region of modulated (helimagnetic) structures. We then turn to
the analysis of the analogous problem on a Cayley tree. Taking into account
the simplicity introduced by the infinite-coordination limit of the tree, we
explore several details of the phase diagrams in terms of temperature and a
parameter of competition. In particular, we characterize sequences of
modulated (helical) structures associated with devil’s staircases of a fractal
character.

††footnotetext<EMAIL_ADDRESS>
Keywords: helimagnetism; Dzyaloshinkii-Moriya interactions; competing
interactions

##  1 Introduction

Competing interactions are known to lead to very rich phase diagrams, with
critical and multicritical points and sequences of spatially modulated
structures [1][2]. The ANNNI model, which stands for axial-next-nearest
neighbor Ising model, is perhaps the best investigated statistical lattice
model with competing interactions [3][4][5]. In the standard ANNNI model,
there are ferrromagnetic interactions between pairs of nearest-neighbor Ising
spin variables on the sites of a cubic lattice, with the addition of competing
antiferromagnetic interactions between second-neighbor spin pairs along an
axial direction. In terms of temperature and a parameter to gauge the
competition, the ANNNI model exhibits perhaps the richest phase diagram of the
literature, with a wealth of modulated phases and the onset of fractal
structures, which have been called devil’s staircases [6][7].

In the present work we were motivated by recent interest in different forms of
helimagnetism [8][9], which can be explained in terms of a chiral mechanism
proposed by Dzyaloshinskii and Moriya (DM) a long time ago [10]. According to
this DM mechanism, interactions are restricted to first neighbors, but with
vector spin variables on each lattice site. Modulation effects come from the
competition between the usual Heisenberg exchange, of scalar nature, and an
additional exchange associated with the vector products of the DM
interactions. Different formulations of spin models including competing
exchange and DM interactions have indeed been used to account for the
monoaxial helical patterns of magnetic compounds [11][12]. We were also
motivated by early work for special clock models with chiral interactions,
which have been investigated by different techniques [13][14][15]. In contrast
to the ANNNI model, elementary lattice model systems associated with the DM
mechanism have not been subjected to more detailed statistical mechanics
calculations. We then believe that there is still room to analyze some simple
and schematic version of these spin models with DM interactions, and to show
that they are capable of displaying a phase diagram with many sequences of
modulated structures.

We consider a ferromagnetic version of a classical Heisenberg spin model with
DM interactions, in zero external field, and restrict the competition to a
single lattice axis. The spin Hamiltonian of this simple monoaxial lattice
system is written as

$\mathcal{H}_{DM}=-J_{0}\sum_{\left(\overrightarrow{r},\overrightarrow{r}^{\prime}\right)}^{\left(x,y\right)}\overrightarrow{S}_{\overrightarrow{r}}\cdot\overrightarrow{S}_{\overrightarrow{r}^{\prime}}-J\sum_{\left(\overrightarrow{r},\overrightarrow{r}^{\prime}\right)}^{\left(z\right)}S_{\overrightarrow{r}}^{z}\,S_{\overrightarrow{r}^{\prime}}^{z}-D\sum_{\left(\overrightarrow{r},\overrightarrow{r}^{\prime}\right)}^{\left(z\right)}\left(\overrightarrow{S}_{\overrightarrow{r}}\times\overrightarrow{S}_{\overrightarrow{r}^{\prime}}\right)\cdot\widehat{z},$
(1)

where $\overrightarrow{S}_{\overrightarrow{r}}$ is a vector spin on site
$\overrightarrow{r}$ of a simple cubic crystal lattice, we assume positive
exchange parameters, $J_{0}$, $J>0$, the first sum is over nearest-neighbor
sites on the $x-y$ planes, and the second and third sums are over nearest-
neighbor sites along the $z$ direction. The last term in this spin Hamiltonian
provides the competition with the usual exchange term, which is the essential
ingredient to mimic a monoaxial DM model system.

We now make another simplifying assumption. Instead of considering continuous
variables, we assume discrete vector spin variables, with just six
possibilities, along the directions of the Cartesian axes. We then write the
six spin components,

$\overrightarrow{S}_{\overrightarrow{r}}=\begin{pmatrix}1\\\ 0\\\
0\end{pmatrix},\,\begin{pmatrix}0\\\ 1\\\ 0\end{pmatrix},\,\begin{pmatrix}0\\\
0\\\ 1\end{pmatrix},\,\begin{pmatrix}-1\\\ 0\\\
0\end{pmatrix},\,\begin{pmatrix}0\\\ -1\\\
0\end{pmatrix},\,\begin{pmatrix}0\\\ 0\\\ -1\end{pmatrix}.$ (2)

Although this is a drastic simplification of the original spin model system,
it keeps the essential features of the DM interactions. This formulation,
however, is still not amenable to exact calculations. We then resort to well-
known mean-field approximate schemes, and analyze the main features of the
phase diagrams in terms of temperature and a parameter that gauges the
competition.

In Section 2, we formulate a solution of this problem on the basis of a layer-
by-layer mean-field calculation, which has been successfully used to
investigate the phase diagram of the ANNNI model [16]. Although it is simple
to formulate and provides very sensible results, this mean-field calculation
demands a considerable numerical effort to unveil the fine details of the
phase diagram. We then restrict the analysis to the overall aspects of the
phase diagrams. In particular, we locate the paramagnetic transition lines,
which already indicate the presence of modulated structures.

In Section 3, we consider an analog of this model system on a Cayley tree. We
formulate the statistical problem as a non linear discrete map along the
generations of the tree. Attractors of this map, in the deep interior of the
tree, are known to correspond to physically realistic solutions, which come
from a standard pair approximation for the analogous model system on a Bravais
lattice. In the limit of infinite coordination of the tree, the problem
becomes considerably simpler, with solutions that should be close to the usual
mean-field results. According to earlier calculation for the analog of the
ANNNI model on a Cayley tree [17], we take advantage of this infinite-
coordination limit to analyze a number of fine details of the phase diagram,
including the onset of sequences of modulated structures and the existence of
devil’s staircases of a fractal character.

##  2 Layer-by-layer mean-field calculation

We consider a cubic lattice of $N^{3}$ sites, and split the Hamiltonian (1)
into two terms,

$\mathcal{H}=\mathcal{H}_{1}+\mathcal{H}_{2},$ (3)

where the first term includes the nearest-neighbor pair interactions on the
$x-y$ planes of the lattice,

$\mathcal{H}_{1}=-J_{0}\sum_{z=1}^{N}\sum_{\left(x,y;\,\,x^{\prime}y^{\prime}\right)}\left[S_{x,y,z}^{x}S_{x^{\prime},y\prime,z}^{x}+S_{x,y,z}^{y}S_{x^{\prime},y\prime,z}^{y}+S_{x,y,z}^{z}S_{x^{\prime},y\prime,z}^{z}\right],$
(4)

and the second term refers to the interactions along the axial $z$ direction,

$\mathcal{H}_{2}=-J\sum_{x,y}\sum_{z=1}^{N}\left[S_{x,y,z}^{x}S_{x,y,z+1}^{x}+S_{x,y,z}^{y}S_{x,y,z+1}^{y}+S_{x,y,z}^{z}S_{x,y,z+1}^{z}\right]-$
$-D\sum_{x,y}\sum_{z=1}^{N}\left[S_{x,y,z}^{x}S_{x,y,z+1}^{y}-S_{x,y,z}^{y}S_{x,y,z+1}^{x}\right].$
(5)

We now use the Bogoliubov inequality,

$G\left(\mathcal{H}\right)\leq
G_{0}\left(\mathcal{H}_{0}\right)+\left\langle\mathcal{H}-\mathcal{H}_{0}\right\rangle_{0}=\Phi,$
(6)

where $G\left(\mathcal{H}\right)$ is the free energy associated with
$\mathcal{H}$, $G_{0}\left(\mathcal{H}_{0}\right)$ is the free energy
associated with a trial Hamiltonian $\mathcal{H}_{0}$, and
$\left\langle...\right\rangle_{0}$ is an expected value in a canonical
ensemble defined by the trial Hamiltonian. The mean-field approximation
consists in assuming that the free energy of this system comes from the
minimization of the upper bound $\Phi$ with respect to the parameters of the
trial Hamiltonian.

The simplest form of a trial Hamiltonian, which still preserves the
possibility of modulated structures along the $z$ axis, is given by

$\mathcal{H}_{0}=-\sum_{x,y}\sum_{z}\overrightarrow{\eta}_{z}\cdot\overrightarrow{S}_{x,y,z}=-\sum_{x,y}\sum_{z}\left[\eta_{z}^{x}S_{x,y,z}^{x}+\eta_{z}^{y}S_{x,y,z}^{y}+\eta_{z}^{z}S_{x,y,z}^{z}\right],$
(7)

where
$\overrightarrow{\eta}_{z}=\left(\eta_{z}^{x},\eta_{z}^{y},\eta_{z}^{z}\right)$
is a three-component trial field. We then have

$G_{0}=-\frac{N^{2}}{\beta}\sum_{z}\ln\left[2\cosh\left(\beta\eta_{z}^{x}\right)+2\cosh\left(\beta\eta_{z}^{y}\right)+2\cosh\left(\beta\eta_{z}^{z}\right)\right]$
(8)

and

$\frac{1}{N^{2}}\left\langle\mathcal{H}-\mathcal{H}_{0}\right\rangle_{0}=-2J_{0}\sum_{z}\left[\left(m_{z}^{x}\right)^{2}+\left(m_{z}^{y}\right)^{2}+\left(m_{z}^{z}\right)^{2}\right]-$
$-J\sum_{z}\left[m_{z}^{x}m_{z+1}^{x}+m_{z}^{y}m_{z+1}^{y}+m_{z}^{z}m_{z+1}^{z}\right]-$
$-D\sum_{z}\left[m_{z}^{x}m_{z+1}^{y}-m_{z}^{y}m_{z+1}^{x}\right]+\sum_{z}\left[\eta_{z}^{x}m_{z}^{x}+\eta_{z}^{y}m_{z}^{y}+\eta_{z}^{z}m_{z}^{z}\right],$
(9)

in which we have introduced the definitions

$m_{z}^{\nu}=\frac{\sinh\left(\beta\eta_{z}^{\nu}\right)}{\cosh\left(\beta\eta_{z}^{x}\right)+\cosh\left(\beta\eta_{z}^{y}\right)+\cosh\left(\beta\eta_{z}^{z}\right)},$
(10)

where $\nu=x$, $y$, $z$, and $\beta=1/\left(k_{B}T\right)$ is the inverse of
temperature.

The minimization of $\Phi$ with respect to the trial fields
$\left\\{\overrightarrow{\eta}_{z}\right\\}$ leads to the equations

$\eta_{z}^{x}=4J_{0}m_{z}^{x}+J\left(m_{z-1}^{x}+m_{z+1}^{x}\right)+D\left(m_{z-1}^{y}-m_{z+1}^{y}\right),$
(11)
$\eta_{z}^{y}=4J_{0}m_{z}^{y}+J\left(m_{z-1}^{y}+m_{z+1}^{y}\right)-D\left(m_{z-1}^{x}-m_{z+1}^{x}\right),$
(12)

and

$\eta_{z}^{z}=4J_{0}m_{z}^{z}+J\left(m_{z-1}^{z}+m_{z+1}^{z}\right).$ (13)

We now insert these expressions for the trial fields into the definitions of
the local magnetizations, given by (10), and write an infinite system of
nonlinear coupled mean-field equations of state for the local magnetizations.

As in the mean-field calculations for the ANNNI model [16], given the
temperature and the parameters of the system, the problem consists in finding
a multiple set of numerical solutions, and choosing the solution that
minimizes the free energy. Except in the immediate neighborhood of a critical
line, this search for the physically acceptable solutions becomes a purely
numerical problem. Usually, we find a numerical solution with the assumption
of a periodicity of $L$ lattice spacings along the $z$ direction, and take
note of the associated value of the free energy. We then repeat this procedure
for a sequence of values of $L$. Physically acceptable solutions are shown to
minimize the free energy with respect to the length $L$. Instead of carrying
out this cumbersome numerical calculation, we now turn to an expansion of the
free energy.

It is easy to write the mean-field free energy as a power series in the spin
magnetizations, and make contact with a Landau-Ginzburg expansion. Keeping
terms up to order $4$, we have

$\frac{1}{N^{3}}\Phi=-\frac{1}{\beta}\ln
6+\frac{1}{N}\sum_{z}\left(-2J_{0}+\frac{3}{2\beta}\right)\left[\left(m_{z}^{x}\right)^{2}+\left(m_{z}^{y}\right)^{2}+\left(m_{z}^{z}\right)^{2}\right]-$
$-\frac{J}{N}\sum_{z}\left[m_{z}^{x}m_{z+1}^{x}+m_{z}^{y}m_{z+1}^{y}+m_{z}^{z}m_{z+1}^{z}\right]-\frac{D}{N}\sum_{z}\left[m_{z}^{x}m_{z+1}^{y}-m_{z}^{y}m_{z+1}^{x}\right]+$
$+\frac{9}{4\beta
N}\sum_{z}\left[\left(m_{z}^{x}\right)^{2}\left(m_{z}^{y}\right)^{2}+\left(m_{z}^{x}\right)^{2}\left(m_{z}^{z}\right)^{2}+\left(m_{z}^{y}\right)^{2}\left(m_{z}^{z}\right)^{2}\right]+...$
(14)

We then assume periodic boundary conditions, and use a Fourier representation,

$m_{z}^{\nu}=\frac{1}{\sqrt{N}}\sum_{q}m_{q}^{\nu}\exp\left(iqz\right),$ (15)

where $\nu=x,y,z$, and the sum is restricted to the first, and symmetric,
Brillouin zone. Also, it is convenient to write

$m_{q}^{\nu}=\frac{1}{\sqrt{2}}(R_{q}^{\nu}+iI_{q}^{\nu}),$ (16)

with

$R_{q}^{\nu}=R_{-q}^{\nu},\qquad I_{q}^{\nu}=-I_{-q}^{\nu}.$ (17)

Keeping terms up to second order, we have

$\frac{1}{N^{3}}\Phi-\frac{1}{\beta}\ln 6=$
$=\frac{1}{2}\sum_{q}\left(\frac{3}{2\beta}-2J_{0}-J\cos
q\right)\left[\left(R_{q}^{x}\right)^{2}+\left(I_{q}^{x}\right)^{2}+\left(R_{q}^{y}\right)^{2}+\left(I_{q}^{y}\right)^{2}+\left(R_{q}^{z}\right)^{2}+\left(I_{q}^{z}\right)^{2}\right]-$
$-D\sum_{q}\left(\operatorname{sen}q\right)\left[I_{q}^{x}R_{q}^{y}-I_{q}^{y}R_{q}^{x}\right]+\cdot\cdot\cdot,$
(18)

which is a quadratic form in terms of the real normal modes in Fourier space.
Note the symmetry along the $q$ axis. Also, note that the coupling is
restricted to the $x-y$ real modes only. It is then convenient to define

$A=\frac{3}{4\beta}-J_{0}-\frac{1}{2}J\cos q;\qquad B=D\operatorname{sen}q,$
(19)

and to introduce the simplified notation

$x_{1}=R_{q}^{x};\quad x_{2}=I_{q}^{x};\quad y_{1}=R_{q}^{y};\quad
y_{2}=I_{q}^{y}\quad z_{1}=R_{q}^{z};\quad z_{2}=I_{q}^{z}.$ (20)

We then have to analyze the quadratic form

$Q=A\left(x_{1}^{2}+x_{2}^{2}+y_{1}^{2}+y_{2}^{2}+z_{1}^{2}+z_{2}^{2}\right)+B\left(y_{2}x_{1}-x_{2}y_{1}\right)=\left(\overrightarrow{v}\right)^{t}\mathbf{M}\left(\overrightarrow{v}\right),$
(21)

in which $\overrightarrow{v}$ is a column vector,
$\widetilde{\mathbf{v}}=\left(x_{1},x_{2},y_{1},y_{2},z_{1},z_{2}\right)$, and
the matrix $\mathbf{M}$ is given by

$\mathbf{M}=\left(\begin{array}[c]{cccccc}A&0&0&B/2&0&0\\\ 0&A&-B/2&0&0&0\\\
0&-B/2&A&0&0&0\\\ B/2&0&0&A&0&0\\\ 0&0&0&0&A&0\\\
0&0&0&0&0&A\end{array}\right).$ (22)

This matrix has two trivial real and degenerate eigenvalues,

$\lambda_{1}=\lambda_{2}=A=\frac{3}{4\beta}-J_{0}-\frac{1}{2}J\cos q.$ (23)

The remaining eigenvalues are associated with a $4\times 4$ supersymmetric
matrix. We then have

$\lambda_{3,4}=A\pm\frac{1}{2}B=\frac{3}{4\beta}-J_{0}-\frac{1}{2}J\cos
q\pm\frac{1}{2}D\sin q.$ (24)

Taking into account the form of these eigenvalues, it is immediate to obtain
an expression for the border of the disordered region,

$\frac{1}{\beta_{c}}=k_{B}T_{c}=\max\limits_{q}\left[\frac{4}{3}J_{0}+\frac{2}{3}J\,\cos
q\pm\frac{2}{3}D\,\sin q\right].$ (25)

For $D=0$, we recover the (mean-field) critical temperature of a simple
ferromagnet,

$\frac{1}{\beta_{c}^{F}}=k_{B}T_{c}^{F}=\max\limits_{q}\left(\frac{4}{3}J_{0}+\frac{2}{3}J\cos
q\right)=\frac{4}{3}J_{0}+\frac{2}{3}J,$ (26)

For $D\neq 0$, we have

$\tan q=\pm\frac{D}{J},$ (27)

which corresponds to propagating waves along the two directions of the $z$
axis, with the critical temperature

$\frac{1}{\beta_{c}^{M}}=k_{B}T_{c}^{M}=\frac{4}{3}J_{0}+\frac{2}{3}\left(J^{2}+D^{2}\right)^{1/2}.$
(28)

From equations (26) and (28) we have the paramagnetic transition lines,

$k_{B}T_{c}^{F}-\frac{4}{3}J_{0}=\frac{2}{3}J\text{ },\qquad
k_{B}T_{c}^{M}-\frac{4}{3}J_{0}=\frac{2}{3}\left(J^{2}+D^{2}\right)^{1/2}.$
(29)

These equations indicate that there is always a transition to a spatially
modulated structure, except at a trivial, ferromagnetic, multicritical point,
at $D=0$.

##  3 Analysis of the DM model on a Cayley tree

We now turn to calculations for the analogous DM model on a Cayley tree, which
is perhaps the simplest way to unveil the rich structure of the ordered
phases. We then consider the monoaxial Hamiltonian,

$\mathcal{H}=-J\sum_{\left(\overrightarrow{r},\overrightarrow{r}^{\prime}\right)}\overrightarrow{S}_{\overrightarrow{r}}\cdot\overrightarrow{S}_{\overrightarrow{r}^{\prime}}-D\sum_{\left(\overrightarrow{r},\overrightarrow{r}^{\prime}\right)}\left(\overrightarrow{S}_{\overrightarrow{r}}\times\overrightarrow{S}_{\overrightarrow{r}^{\prime}}\right)\cdot\widehat{z},$
(30)

where $J>0$, and $\overrightarrow{r}$ and $\overrightarrow{r}^{\prime}$ are
nearest-neighbor sites along the successive generations of a Cayley tree (and
we discard ferromagnetic interactions between sites belonging to the same
generation of the tree). Although it is not entirely similar to the previously
analyzed mean-field model, the formulation on a Cayley tree is known to lead
to essentially the same qualitative results [18]. Also, the problem is written
as a nonlinear dissipative discrete map, along the generations of a Cayley
tree, which is much easier to analyze than the corresponding area-preserving
map associated with the mean-field solutions. Attractors of this dissipative
mapping problem are known to mimic the phase structure of the ordered region
of the mean-field phase diagrams.

Figure 1: Sketch of three generations of a Cayley tree with coordination
three.

In Figure 1, we sketch a Cayley tree of coordination $z=3$ (and ramification
$r=z-1=2$). We then consider a tree of arbitrary ramification $r$, and use
well-known procedures to write six recursion relations,

$Z_{1}^{\prime}=\left(e^{\beta J}Z_{1}+e^{\beta D}Z_{2}+Z_{3}+e^{-\beta
J}Z_{4}+e^{-\beta D}Z_{5}+Z_{6}\right)^{r},$ (31)
$Z_{2}^{\prime}=\left(e^{-\beta D}Z_{1}+e^{\beta J}Z_{2}+Z_{3}+e^{\beta
D}Z_{4}+e^{-\beta J}Z_{5}+Z_{6}\right)^{r},$ (32)
$Z_{3}^{\prime}=\left(Z_{1}+Z_{2}+e^{\beta J}Z_{3}+Z_{4}+Z_{5}+e^{-\beta
J}Z_{6}\right)^{r},$ (33) $Z_{4}^{\prime}=\left(e^{-\beta J}Z_{1}+e^{-\beta
D}Z_{2}+Z_{3}+e^{\beta J}Z_{4}+e^{\beta D}Z_{5}+Z_{6}\right)^{r},$ (34)
$Z_{5}^{\prime}=\left(e^{\beta D}Z_{1}+e^{-\beta J}Z_{2}+Z_{3}+e^{-\beta
D}Z_{4}+e^{\beta J}Z_{5}+Z_{6}\right)^{r},$ (35)
$Z_{6}^{\prime}=\left(Z_{1}+Z_{2}+e^{-\beta J}Z_{3}+Z_{4}+Z_{5}+e^{\beta
J}Z_{6}\right)^{r},$ (36)

where $Z_{i}^{\prime}$ is a partial partition function associated with a
central site in the spin state $i$, given by equation (2), which is related to
the partial partition functions, $Z_{1}$ to $Z_{6}$, associated with the
nearest-neighbor sites belonging to the next generation of the tree. The
cycle-free structure of the Cayley tree enormously simplifies the form of
these recursion relations. It is worth to remark that, for $r=1$, we regain
the transfer matrix associated with the solution of the corresponding model on
a chain.

Although it is not difficult to analyze these recursion relations, we further
simplify the problem, and at the same time emphasize the analogy with the
mean-field results, by resorting to the infinite coordination limit of the
tree. We then take $r\rightarrow\infty$, and $J$, $D\rightarrow 0$, but keep
fixed the products $Jr$ and $Dr$.

In this infinite-coordination limit, the problem is reduced to the analysis of
just three recursion relations, which involve the vector components of the
average of an “effective magnetization”,

$m_{x}^{\prime}=\frac{1}{M}\left[\sinh(\beta Jrm_{x}+\beta Drm_{y})\right],$
(37) $m_{y}^{\prime}=\frac{1}{M}\left[\sinh(-\beta Drm_{x}+\beta
Jrm_{y})\right],$ (38) $m_{z}^{\prime}=\frac{1}{M}\left[\sinh(\beta
Jrm_{z})\right],$ (39)

with the denominator

$M=\cosh(\beta Jrm_{x}+\beta Drm_{y})+\cosh(-\beta Drm_{x}+\beta
Jrm_{y})+\cosh(\beta Jrm_{z}).$ (40)

We now use analytical and numerical techniques to analyze the attractors of
this map.

###  3.1 Attractors of the map

There is always a trivial solution, $m_{x}^{\ast}=m_{y}^{\ast}=m_{z}^{\ast}=0$
of equations (37)-(39). The linearization about this trivial fixed point leads
to the matrix equation

$\overrightarrow{m}^{\prime}=\mathbf{M}\,\overrightarrow{m},$ (41)

where

$M=\frac{1}{3}\begin{pmatrix}\beta Jr&\beta Dr&0\\\ -\beta Dr&\beta Jr&0\\\
0&0&\beta Jr\end{pmatrix}.$ (42)

Therefore, we have the eigenvalues

$\lambda_{1}=\frac{1}{3}\beta Jr,$ (43) $\lambda_{2}=\frac{1}{3}(\beta
Jr+i\beta Dr),$ (44)

and

$\lambda_{3}=\frac{1}{3}(\beta Jr-i\beta Dr).$ (45)

The conditions of linear stability of this trivial (disordered) fixed point
are given by

$T>\frac{1}{3},$ (46)

and

$T>\frac{\left(1+p^{2}\right)^{1/2}}{3},$ (47)

where we have introduced a simplified notation, $T=1/\left(\beta Jr\right)$
and $p=D/J$. We then have a similar picture as in the previous mean-field
calculation. There is always a paramagnetic transition, which is associated
with the complex eigenvalues of the trivial fixed point, except at $p=0$,
which leads to a simple ferromagnetic transition.

In order to go beyond the linear analysis, and investigate the stability of
additional attractors, we have to resort to numerical calculations. We then
characterized two types of behavior: (i) in one of these possibilities, the
magnetization settles in a specific value (one of the components settles in
the value $1$, and the other components vanish); (ii) the second possibility
is a cyclic pattern in the representation of $m_{x}$ versus $m_{y}$, as we
indicate in Figure 2 (the component $m_{z}$ of the magnetization flows to a
fixed value, as in the first case). The first possibility leads to a quite
trivial ordered structure. The second possibility, however, leads to a
spatially modulated structure.

Figure 2: Graph of $200$ iterations of $m_{y}$ versus $m_{x}$, for $T=0.2$ and
$p=0.9$. We assume the initial values $m_{x}=m_{y}=m_{z}=1$, and discard the
initial $10000$ iterations. This modulated structure is characterized by the
wave number of the main harmonic component, $q/2\pi=0.13$ (assuming a unit
lattice paramater).

We have used this numerical procedure to sketch the phase diagram of Figure 3.
Attractors of the map correspond to disordered (D), ferromagnetically ordered
(O), and spatially modulated (M) structures. In general, given the values of
reduced temperature $T$ and the parameter $p=D/J$, we choose arbitrary initial
conditions ($0<m_{x},m_{y},m_{z}<1$) and iterate equations (37) to (39) about
$10^{4}$ times. Final values of the magnetizations $m_{x}$ and $m_{y}$ are
determined with a precision of about $10^{-5}$ (the component $m_{z}$ does not
display any modulation).

Figure 3: Phase diagram in terms of temperature ($T$) and the parameter of
chirality ($p=D/J$). We indicate disordered (D), ferromagnetically ordered (O)
and some spatially modulated phases (in the white region). We use the
dimensionless wave number, $q/2\pi$, which is written in terms of $2\pi$, to
indicate some of the main modulated phases. In the groud state, for $p>1$, we
have $q/2\pi=0.25$.

##  4 The devil’s staircases

The devil’s staircase is a fractal structure that has been associated with the
sequences of modulated structures in the ordered region of the phase diagrams
of lattice statistical models, in particular the ANNNI model, with competing
interactions. We show that the same kind of behavior comes from the analysis
of the iterates of equations (37) to (39).

Given the parameters $T$ and $p$, we plot a graph of $m_{y}$ versus $m_{x}$,
as shown in Figure 2, and note that the attractor is a periodic function,
which gives rise to the definition of a period $q$ in units of $2\pi$. In a
more precise calculation, we perform a Fourier analysis of the attractor and
find the main Fourier component, which sets the period $q$. In Figure 4 we
plot the wave number $q$ versus the chiral parameter $p$, for a fixed value of
temperature, $T=0.2$, which already displays the characteristic shape of a
devil’s staircase.

Figure 4: Sketch of a devil’s staircase at $T=0.2$.

We now turn to the calculation of the Hausdorff dimension associated with
these plots of $q$ versus $p$, which leads to the characterization of the
fractal nature of the devil’s staircase. According to a box-counting algorithm
[7], we choose a value $\epsilon>0$, and then calculate the sum of the step
sizes, along the $p$ axis, of all steps with width (size along the $p$ axis)
larger than than $\epsilon$. We then subtract the total width of the $p$
interval and obtain the function $X(\epsilon)$. The slope of a plot of
$\log(X(\epsilon)/\epsilon)$ versus $\log(1/\epsilon)$, in the limit of small
$\epsilon$, leads to an estimate of the Hausdorff dimension $D_{H}$ of the
staircase. In the modulated region, we found several examples of plots of $q$
versus $p$, at fixed temperature, with $D_{H}<1$, which characterizes a
fractal object. Also, we observed that $D_{H}$ increases with temperature,
which indicates a route to an incommensurate structure. In figure 5, we show
one of these plots of the Hausdorff dimension versus normalized temperature.

Figure 5: Plot of the Hausdorff dimension as a function of normalized
temperature $T$, calculated in the interval $0.65<p<1.15.$

The increase of the Hausdorff dimension $D_{H}$ with $T$ ($D_{H}\rightarrow 1$
as $T\rightarrow T_{c}$) indicates the typical route to an incommensurate
structure.

##  5 Conclusions

We used mean-field approximations to analyze the phase diagram of an
elementary system of discrete classical spin variables on a crystal lattice,
with ferromagnetic nearest-neighbor exchange interactions and the addition of
competing chiral interactions along an axis. This is perhaps the simplest
three-dimensional lattice system with the inclusion of Dzyaloshinskii-Moriya
interactions. We show that this monoaxial system displays the characteristic
modulated phases of helimagnetic compounds. Besides writing the equations of a
conventional mean-field layer-by-layer calculation, which already indicates
the existence of transitions to modulated structures, we performed some
detailed calculations for the analogous model on a Cayley tree, in the limit
of infinite coordination. Calculations on the tree confirm the general
features of the mean-field results, and provide a simpler way to check fine
details of the phase diagrams, in terms of temperature and a parameter of
chirality. In particular, we show the existence of many sequences of modulated
structures. Also, we confirm the existence of devil’s staircases, which are
associated with a non-integer Hausdorff dimension.

Acknowledgment

We acknowledge many fruitful discussions with Doctor Eduardo S. Nascimento.
Also, we acknowledge the financial support of the Brazilian agencies CNPq and
CAPES.

## References

  * [1] M. Seul and D. Andelman, Science, 267, 476–483 (1995).
  * [2] D. Andelman and R. E. Rosensweig, J. Phys. Chem. B 113, 3785–3798 (2009).
  * [3] P. Bak, Repts. Progr. Phys. 45, 587, 1982.
  * [4] W. Selke, Phys. Rep. 170, 213, 1988.
  * [5] J. Yeomans, Sol. State Phys. 41, 151, 1988.
  * [6] W. Selke, “Spatially modulated structures in systems with competing interactions” in Phase Transitions and Critical Phenomena, Vol. 15, p.1-72, ed. by C. Domb and J.L. Lebowitz (Academic Press, 1992).
  * [7] E. S. Nascimento, J. P. de Lima and S. R. Salinas, Physica A 409, 78–86 (2014).
  * [8] U. K. Rossler, A. N. Leonov, A. N. Bogdanov, J. Phys. Conf. Ser. 200, 022029 (2010).
  * [9] N. Nagaosa, Y. Tokura, Nature Nanotechnology 8, 899 (2010).
  * [10] Yu. A. Izyumov, Sov. Phys. Usp. 27, 845, 1984.
  * [11] J. Kishine, K. Inoue, Y. Yoshida, Progr. Theor. Phys. Suppl. 159, 82, 2005.
  * [12] M. Shinozaki, S. Hoshino, Y. Masaki, J. Kishine, Y. Kato, J. Phys. Soc. Jpn. 85, 074710, 2016; T. Toretsume, T. Kikuchi, R. Anita, J. Phys. Soc. Jpn. 87, 041011, 2018.
  * [13] D. A. Huse, Phys. Rev. B24, 5180, 1981.
  * [14] S. Ostlund, Phys. Rev. B24, 398, 1981.
  * [15] H. C. Öttinger, J. Phys. C16, L257, 1983; J. Phys. C16, L597, 1983.
  * [16] C. S. O. Yokoi, M. D. Coutinho-Filho, S. R. Salinas, Phys. Rev. B24, 4047, 1981.
  * [17] C. S. O. Yokoi, M. J. de Oliveira, S. R. Salinas, Phys. Rev. Lett. 54, 163, 1985.
  * [18] M. J. de Oliveira and S. R. Salinas, J. Phys. A18, L1157, 1985.