Kai Feng

Wentao Ma

Wei Huang

Qing Zhang

Yili Gong
email: yiligong@whu.edu.cn

Speeding up Galois Field Arithmetic on Intel MIC Architecture

Keywords: Galois Field Arithmetic, MIC Architecture, SIMD, Open-MP, Speedup

Galois Field arithmetic is the basis of LRC, RS and many other erasure coding approaches. Traditional implementations of Galois Field arithmetic use multiplication tables or discrete logarithms, which limit the speed of its computation. The Intel Many Integrated Core (MIC) Architecture provides 60 cores on chip and very wide 512-bit SIMD instructions, attractive for data intensive applications. This paper demonstrates how to leverage SIMD instructions and shared memory multiprocessing on MIC to perform Galois Field arithmetic. The experiments show that the performance of the computation is significantly enhanced.

Introduction

From disk arrays [START_REF] Blaum | EVENODD: An efficient scheme for tolerating double disk failures in RAID architectures[END_REF], cloud platforms [START_REF] Huang | Erasure coding in windows azure storage[END_REF] to archival systems [START_REF] Tansley | DSpace as an open archival information system: Current status and future directions[END_REF] storage systems must have fault tolerance to protect themselves from data loss. Erasure codes provide the basic technology for the fault tolerance of a storage system. The classic Reed-Solomon code [START_REF] Plank | A tutorial on Reed-Solomon coding for fault-tolerance in RAID-like systems[END_REF] organizes a storage system as a set of linear equations whose arithmetic is Galois Field arithmetic, termed GF(2 w ). W is the length of a word, the basic computing unit. Encoding and decoding of a storage system for fault tolerance are implemented by computing these linear equations by multiplying large regions of bytes by various w -bit constants in GF(2 w ) and combining the products using bitwise exclusive-or (XOR).

Traditional implementations of Galois Field arithmetic use multiplication tables or discrete logarithms, which limit the speed of its computation. The performance using multiplication is at least four times slower than using XOR [START_REF] Plank | A performance evaluation and examination of open-source coding libraries for storage[END_REF]. James S. Plank et al. fast Galois Field arithmetic using 128-bit SIMD instruction [START_REF] Plank | Screaming fast Galois Field arithmetic using Intel SIMD instructions[END_REF].

In late 2012, Intel released its commercial products based on the Many Integrated Core (MIC) architecture [START_REF]Intel? C++ Compiler XE 13.1 User and Reference Guides[END_REF], targeting to High Performance Computing field for the PetaFLOPS era. It is based on the streamlined x86 core and similar to the architecture of the existing CPUs. Since its architectural compatibility, it can utilize existing parallelization software tools, including OpenMP [START_REF]OpenMP Application Program Interface[END_REF], etc. and specialized versions of Intel's Fortran, C++ and math libraries [START_REF]Intel Math Kernel Library for Linux* OS User's Guide[END_REF]. Its SIMD instructions are further extended to very wide 512-bit and allow 512-bit numbers to be manipulated on a core simultaneously. MIC's 60 cores also greatly enhance its parallel computing capabilities.

To the best of our knowledge, how to use a computing unit as powerful as a MIC coprocessor for Galois Field arithmetic has not been discussed yet. When the operator size of SIMD instructions extends from 128 bits to 512 bits, though the number of elements keeps at 16, the size of each element changes from 8 bits to 32 bits. With smaller w, e.g. w = 4, the spatial utilization ratio is only 1/8 for the multiplication table. The obvious waste needs to be avoided to save memory usage. As to larger w, e.g. w = 32, the existed algorithm [START_REF] Plank | Screaming fast Galois Field arithmetic using Intel SIMD instructions[END_REF] maps a word into 4 8-bit parts since the element size of 128-bit SIMD instructions is 8-bit, which in-creases complexity and decreases performance. With 32-bit elements, the over-head should be reduced. This paper will detail how to leverage 512-bit SIMD instructions and shared memory multiprocessing to multiply regions of bytes by constants in GF(2 w ) for w ∈ {4, 8, 16, 32}. Each value of w has similar but still different implementation techniques. We will present these techniques and compare the performance of our algorithms on MIC with other approaches on other platforms.

The rest of this paper is organized as follows. The next section describes related work. Section 3 gives description about Erasure Codes and Galois Fields. Section 4 introduces 512-bit instructions used in our algorithms. Section 5 details our algorithms leveraging 512-bit SIMD instructions and OpenMP to multiply regions of bytes by constants in GF(2 w ) for w varying from 4 to 32. Section 6 compares and analyzes the performance of our algorithms and the others. Section 7 is the conclusion and future work.

Related Work

Erasure coding is an alternative to replication for fault tolerance as storage systems scale. Traditionally used in the communication field, erasure codes have gained their popularity due to lower spatial requirement under the same reliability.

Many erasure codes are based on Galois Field arithmetic, such as Pyramid codes [START_REF] Huang | Pyramid codes: Flexible schemes to trade space for access efficiently in reliable data storage systems[END_REF], LRC codes [START_REF] Huang | Erasure coding in windows azure storage[END_REF], RS codes [START_REF] Kalcher | Accelerating Galois Field arithmetic for Reed-Solomon erasure codes in storage applications[END_REF] and F-MSR codes [START_REF] Hu | NCCloud: Applying network cod-ing for the storage repair in a cloud-of-clouds[END_REF], among which the most common one is RS codes. RS codes are used in Bigtable [START_REF] Chang | Bigtable: A Dis-tributed Storage System for Structured Data[END_REF] from Google, Cassandra [START_REF] Lakshman | Cassandra: A Decentralized Structured Storage System[END_REF] from Facebook and Cleversafe [START_REF] Resch | AONTRS: blending security and performance in dispersed storage systems[END_REF]. Microsoft Azure uses LRC codes [START_REF] Huang | Erasure coding in windows azure storage[END_REF].

Traditional implementations of Galois Field arithmetic adopt multiplication tables or discrete logarithms. There are methods proposed to improve Galois Field arithmetic, such as Kevin M. Greenan et al. using split multiplication tables and composite fields [START_REF] Greenan | Optimizing Galois Field arithmetic for diverse processor architectures and applications[END_REF], Jianqiang Luo et al. using bit-grouping tables [START_REF] Luo | Efficient software implementations of large finite fields GF(2 n ) for secure storage applications[END_REF] and H. Peter Anvins approach based on fast multiplication by two [START_REF]OpenMP Application Program Interface[END_REF][START_REF] Anvin | The mathematics of RAID-6[END_REF] and so on.

Recently in [START_REF] Plank | Screaming fast Galois Field arithmetic using Intel SIMD instructions[END_REF] James S. Plank et al. present the algorithms of Galois Field arithmetic on CPUs using 128-bit SIMD instructions. As with [START_REF] Plank | Screaming fast Galois Field arithmetic using Intel SIMD instructions[END_REF], this paper focuses solely on multiplying regions of bytes by constants. We will exploit 512bit SIMD instructions as well as OpenMP on MIC coprocessors.

Erasure Codes and Galois Fields Arithmetic

Fault tolerance of a storage system is enabled by redundancy. For Galois Field Arithmetic based erasure codes, n disks are partitioned into k disks for original data and m disks for coding information, which is calculated from the original data. When no more than m disks fail, the lost data can be recovered through the remaining disks. For example, RAID-6 has two (m = 2) coding disks (C 0 and C 1 ), which are created from k data disks (D i , 0≤i <n) as shown in Fig. 1 (a). Content of every disk is composed of w -bit words, such as d ih and c ih (0≤i <k, 0≤j <2, 0≤h<l). Here l is the number of words in a disk. The coding disks are created by a set of linear equations on the right.

The arithmetic of redundant code generation mainly includes Galois Field multiplication and addition, which correspond to multiplication and XOR operations. Taking C 1 as an example, every word d ih is multiplied by a constant a i , shown in Fig. 1 (b). The products of d ih and a i (0≤i <k) are added (XOR-ed) and the sum is c ih (0≤h<l). Since the speed of XOR operations is very fast for modern computes, multiplication becomes the dominant concern with code calculating.

The selection of w decides the number of disks in the storage system for protection. For example, when using Reed-Solomon codes, w = 4 means the disk number cannot be larger than 16; w = 16 sets the limit to 65,536 disks. The value of w also greatly impacts the computation performance. Larger values of w perform much more slowly than smaller ones. Usually w is a power of 2 to match the size of machine words. Combining all the factors together, typically w is 4 or 8 for storage systems [START_REF] Huang | Erasure coding in windows azure storage[END_REF][START_REF] Resch | AONTRS: blending security and performance in dispersed storage systems[END_REF] and could be 32 and 64 for security and erasure coding purpose [START_REF] Luo | Efficient software implementations of large finite fields GF(2 n ) for secure storage applications[END_REF].

512-bit SIMD Instructions

The Intel Many Core not only has ordinary vector floating-point units, but also uses special registers that enable packed data of up to 512 bits in length for optimal vector graphic streaming SIMD processing. These 512-bit instructions [START_REF]Intel? C++ Compiler XE 13.1 User and Reference Guides[END_REF] can manipulate sixteen elements of 32 bits or eight elements of 64 bits at a time. In this paper, we use manipulation of 16 elements of 32 bits simultaneously. We leverage the following instructions in our implementations: mm512 setzero epi32(void): sets all the elements of the 512-bit vector to zero. Returns a 512-bit vector with all elements set to zero. mm512 set1 epi32(int a): sets all 16 elements of an int32 result vector to an equal integer value specified by a. Returns an int32 vector with 16 elements each equal to integer value specified by a. mm512 slli epi32( m512i v2, unsigned int count): performs an element-byelement logical left shift of int32 vector v2, shifting by the number of bits given by immediate count. If the shift value specified by this parameter is greater than 31 then the result of the shift is zero. mm512 srli epi32( m512i v2, unsigned int count): performs an element-byelement logical right shift. mm512 and epi32( m512i v2, m512i v3): performs a bitwise AND operation between int32 vectors v2 and v3. mm512 xor epi32( m512i v2, m512i v3): performs a bitwise XOR operation between int32 vectors v2 and v3. mm512 loadunpackhi epi32( m512i v1 old, void const* mt): the high 64byte-aligned portion of the double word stream starting at the elementaligned address mt is loaded. It usually works together with the intrinsic mm512 loadunpacklo epi32( m512i v1 old, void const* mt) to load 64 bytes in memory into a 512-bit variable. mm512 permutevar epi32( m512i v2, m512i v3): this is the real enabling SIMD instruction for GF(2 w ). It permutes 32-bit blocks of int32 vector v3 according to indices in the int32 vector v2. The i th element of the result is the j th element of v3, where j is the i th element of v2.

Galois Field Arithmetic on MIC

In this section, calculating yA in GF(2 4 ), GF(2 8 ), GF (2 16 ) and GF(2 32 ) on MIC are presented respectively.

Calculating y A in GF(2 4 )

When w = 4, each word is composed of four bits, and there are only 16 values that a word may be. All operations are based on a 16 16 multiplication table that is small enough to fit into main memory and can be calculated in advance. A table lookup is needed every four bits, i.e. 2K lookups for a region of 1K bytes.

The SIMD intrinsics operates on operators composed of 16 32-bit elements simultaneously. In the original table, each entry corresponds to the 16 4-bit results of a number y multiplied by 16 numbers from 0 to 15. Storing only 4-bit in a 32-bit element is obviously a waste. Thus we try to merge multiple entries into one in the multiplication table, which is showed in Fig. 2. The products of y and 0x0 to 0xf from 8 entries are placed in 16 elements from the lowest to highest, and in each element the product from entry 7 on the high end and the one from entry 0 at the low end. Compressing entries 8-15, 16-23 is similar.

Since the processing element of SIMD instructions is 32-bit while w = 4, every 32 bits in an element are split into 8 4-bit unit using mask[i], shown in Fig. 3 step [START_REF] Plank | Screaming fast Galois Field arithmetic using Intel SIMD instructions[END_REF]. Step ( 7)-( 9) calculated tmp[i] and should be executed for 0≤i<8. Finally, perform XOR operation on all tmp values and get yA. Thus 40 SIMD instructions fulfill 128 multiplication operations.

In general the amounts of data to be computed are huge. Dividing data into basic units of 512 bits and there are no data dependence among them. Thus it is natural to parallelize Galois Field Arithmetic by OpenMP exploiting 60 cores on MIC and opens up to 240 threads.

Calculating y A in GF(2 8 )

When w = 8, each word is 8-bit and there are 256 values that a word may have. In principal the method used in GF(2 4 ) is applicable to the one in GF(2 8 ). The difference is that the instruction mm512 permutevar epi32() only works on 16element tables (each element is 32-bit), 256 values are too large to fit into a Fig. 2. Merge eight entries into one in the multiplication table when w = 4. Four entries in the original table are merged into one to fit 512-bit registers and variables on MIC. The upper line is high-order 256-bits and the lower line is low-order 256-bits. All variables are presented in hex. 16-element variable. Let a be an 8-bit word and a h and a l be the high-order 4 bits and low-order 4 bits of a respectively, and we have:

a = (a h 4) ⊕ a l . (1)
Thus ya = y(a h 4) ⊕ ya l .

(2)

Based on the above analysis, the multiplication table is divided into two, table high which stores the result of y(a h 4) and table low which storage the result of ya l . As with GF(2 4 ), multiplication tables are compressed and occupy 8KB memory. Fig. 4 shows the steps to extract the corresponding content from the compressed lookup tables for mm512 permutevar epi32() to permute. Since the lookup content for y = 7 is at 24-31 bit of each element in the compressed table entry, both table high and table low , it is extracted by right-shifting 24 bits and masked by 0xff. After acquiring the lookup tables, the remaining steps are similar to the ones with w = 4 in Fig. 3, except for step [START_REF]OpenMP Application Program Interface[END_REF] and [START_REF]Intel Math Kernel Library for Linux* OS User's Guide[END_REF]. For w = 8, eight 4-bits in an element is indexed by i (0≤i <8). When i is odd, it means that these 4 bits are high-order of a word; when it is even, these 4 bits are low-order of a word. High-order 4 bits and low-order 4 bits are subject to looking up different tables, table high and table low , as well as left-shifting different bits. The revisions are as follows: (8) for the high-order 4 bits i.e. i is odd tmp[i] = mm512 permutevar epi32(tmp[i], th). for the low-order 4 bits i.e. i is even tmp[i] = mm512 permutevar epi32(tmp[i], tl). ( 9) When i is odd:

tmp[i] = mm512 slli epi32(tmp[i], (i-1)
2). When i is even:

tmp[i] = mm512 slli epi32(tmp[i], i
2).

Calculating y A in GF(2 16 )

For GF(2 16 ) each 16-bit word may have 2 16 = 64K values. Since the instruction mm512 permutevar epi32() only works on 16-element tables, word a is divided into 4-bit sub-words, named a 3 through a 0 :

a = (a 3 12) ⊕ (a 2 8) ⊕ (a 1 4) ⊕ a 0 . (3)
Then ya = y(a 3 12) ⊕ y(a 2 8) ⊕ y(a 1 4) ⊕ ya 0 . (4)
Thus, we need perform 4 table lookup operations for a 16-bits word. We use compressed tables for data storage. The entries from four tables for a constant y take up 256 bytes and the total memory usage is 8 MB.

Calculating y A in GF(2 32 )

For w = 32, the processing is similar. We split each word a (32 bits) into 4-bit sub-words, named a 7 through a 0 :

a = (a 7 28) ⊕ (a 6 24) ⊕ (a 5 20) ⊕ (a 4 16) ⊕ (a 3 12) ⊕(a 2 8) ⊕ (a 1 4) ⊕ a 0 . (5)
Then

ya = y(a 7 28) ⊕ y(a 6 24) ⊕ y(a 5 20) ⊕ y(a 4 16) ⊕ y(a 3 12) ⊕y(a 2 8) ⊕ y(a 1 4) ⊕ ya 0 . (6)
Thus we need perform 8 table lookup operations for a 32-bit word. Since the element size is 32-bit and the same as the size of Galois Field arithmetic word, w, there is no need for compression. The entries from eight tables for a constant y take up 512 bytes and the total size is 2 TB, which is too large to fit into main memory.

Performance Evaluation

The performance of our proposed algorithms on a MIC coprocessor is evaluated and for comparison the Multiplication Table algorithms [START_REF] Plank | A performance evaluation and examination of open-source coding libraries for storage[END_REF] and the 128-bit SIMD algorithms from [START_REF] Plank | Screaming fast Galois Field arithmetic using Intel SIMD instructions[END_REF] are run on a CPU machine.

The MIC machine used in the experiments is Intel Xeon Phi coprocessor 5110p, 60 cores, core frequency 1.053 GHz, 8 GB GDRR5 memory, 32 KB L1 Instruction Cache, 32 KB L1 Data Cache, 512 KB unified L2 Cache. When the cores do not share data or code, the effective L2 Cache is 30 MB. The comparing machine is Intel Xeon CPU E5620 * 2, 2.4 GHz, 32 KB L1 Instruction Cache, 32 KB L1 Data Cache, 256 KB L2 Cache, 12 MB L3 Cache, 32 GB memory.

The multiplication table algorithms and 128-bit SIMD algorithms are tested on CPU and MIC machines. Our proposed 512-SIMD algorithms are run on MIC with native mode. In all algorithms, regions of random values are multiplied by constants in GF(2 w ). For OpenMP accelerated algorithms the region size varies from 1 MB to 1 GB, while for Multiplication Table and SIMD only algorithms the size range is 1 KB to 1 GB. The results are shown in Fig. 5 -Fig. 9.

From Fig. 5 (MulTa is the abbreviation for multiplication table) it can be seen that the SIMD algorithms (128-bit SIMD on CPU and 512-bit SIMD on MIC) greatly outperform the multiplication table algorithms. When w = 4, the performance using SIMD on MIC is 13 times more than that of using multiplication table, and 10.6 times on CPU. We can also conclude that the performance of both algorithms on CPU is better than that on MIC, mainly because the core on CPU is more powerful than the one on MIC (2.4 GHz over 1.053 GHz). For example the multiplication table algorithm on CPU is about 1.8 times faster than on MIC and the SIMD is 1.3 times faster. With w = 8, 16 and 32 we have similar results and the details are omitted. Fig. 6 presents the performance under different w values. We can see that the performance does not change much as w grows which is quite different from the conclusion from [START_REF] Plank | Screaming fast Galois Field arithmetic using Intel SIMD instructions[END_REF]. In [START_REF] Plank | Screaming fast Galois Field arithmetic using Intel SIMD instructions[END_REF] w = 4 and w = 8 perform roughly the same, w = 16 slightly slower and w = 32 slower still. This is because MIC SIMD instructions can operate on more bits (512 bits over 128 bits) simultaneously thus fewer operations needed for a word processing, which benefits larger w. For a certain w, when the region size reaches a point between 256 KB and 512 KB, the performance peaks and then drops dramatically. This is because L2 cache saturation impacts the performance greatly. The results of OpenMP-based acceleration on the algorithms are shown in Fig. 7 -Fig. 9. For the multiplication table algorithm, it is always CPU-intensive thus changing the region size does little impact on performance as given in Fig. 7. For the 128-bit SIMD algorithms, before L3 cache saturates 8 threads are better than 4 threads; after the saturation they are of the same since it is I/O bound now. In the best case, the 128-bit SIMD outperforms the multiplication table by 9.5. Fig. 8 -Fig. 9 compare the performance of the multiplication table algorithm (w = 4) with the 512-bit SIMD (w = 4, 8, 16 and 32) on MIC. Though each core on MIC is capable of 4-way hardware multi-threading, 240 threads do not have the best performance while generally speaking 180 threads are the best. The 512-bit SIMD + OpenMP algorithm is better than the multiplication table + OpenMP on MIC by 6.8 times and better than the 128-bit SIMD + OpenMP on CPU by 7.2 times.

The peak speedups for all algorithms and conditions are summarized in Table 1 with w = 4. Here we take the performance of the single-threaded multiplication table algorithm on CPU as the base 1.

From Fig. 9 (a) -(d) right before the combined 32 MB L2 cache saturates the computing peak can be about 220 GB/s. MIC works as a coprocessor and is connected to the host by standard PCIe x16 which has one-way bandwidth 8 GB/s theoretically. In practice, we have tested that the peak bandwidth from MIC to CPU is 7.0 GB/s and that is 6.7 GB/s from CPU to MIC. Obviously I/O is the bottleneck of Galois Field arithmetic.

Fig. 1 .

 1 Fig. 1. The RAID-6 data disks and coding disks (k = 4). (a) The composition of the RAID-6. (b) How to create code disk C1.

Fig. 3 .

 3 Fig. 3. Multiplying a 512-bit region A by y = 7 in GF(2 4 ).

Fig. 4 .

 4 Fig. 4. Multiplying a 512-bit region A by y = 7 in GF(2 8 ).

  (a) MulTa algorithms. (b) SIMD algorithms.

Fig. 5 .

 5 Fig. 5. The performance of MulTa algorithms and SIMD algorithms on CPU and MIC with w = 4.

Fig. 6 .

 6 Fig. 6. The performance of SIMD algorithms on MIC with w = 4, 8, 16 and 32.

Fig. 7 .

 7 Fig. 7. The performance of OpenMP accelerated MulTa and SIMD algorithms on CPU with 4 and 8 threads, when w = 4.

Fig. 8 .

 8 Fig. 8. The performance of OpenMP accelerated MulTa algorithm on MIC, when w = 4.

Fig. 9 .

 9 Fig. 9. The performance OpenMP accelerated 512-bit SIMD algorithms with different thread numbers on MIC.

Table 1 .

 1 The speedups with w = 4 (taking the performance of the single-threaded multiplication table on CPU as the base 1; MulTa is the abbreviation for multiplication table).In this paper, we detail how to apply 512-bit SIMD instructions with OpenMP on MIC to Galois Field arithmetic. The algorithms are evaluated with different w from 4 to 32. The performance of our algorithms is about 7.2 to 35.2 times faster than the implementations using 128-bit SIMD with OpenMP on CPU.With 512-bit SIMD and OpenMP, cache, main memory and I/O to host become bottlenecks. In future we focus on improving the I/O performance and coordination between computation and data transfer.

	Peak Speedup	MulTa MulTa+OpenMP SIMD SIMD+OpenMP
	CPU(128-bit SIMD)	1	5.5	10.6	52.2
	MIC(512-bit SIMD) 0.56	55.4	7.2	373.4
	7 Conclusion and Future Work

Acknowledgments. We would like to thank Inspur (Beijing) Electronic Information Industry Co., Ltd offering a server with Intel Xeon Phi coprocessor for experiments.

This work is supported by the National Natural Science Foundation of China under Grant No. 61100020.


Baolei Cheng
email: chengbaolei@suda.edu.cn

Jianxi Fan
email: jxfan@suda.edu.cn

Jiwen Yang
email: jwyang@suda.edu.cn

Xi Wang

Circular Dimensional-Permutations and Reliable Broadcasting for Hypercubes and Möbius cubes

Keywords: dimensional-permutation, reliable broadcasting, hypercube, Möbius cube, independent spanning tree

Reliable broadcasting for interconnection networks can be achieved by constructing multiple independent spanning trees(ISTs) rooted at the same node. In this paper, we prove that there exists (n -1)! sets of ISTs rooted at arbitrary node for Qn and Mn based on circular dimensional-permutations of 0, 1, . . . , n -1 and n ≥ 1. At the same time, we give an parallel algorithm, called BCIST, which is the further study of IST problem for Qn and Mn in literature. Furthermore, simulation experiments of ISTs based on JUNG framework and different sets of disjoint paths between node 1 and any node v ∈ V (0-M4)\{1} for 0-M4 are also presented.

Introduction

It is well known that hypercubes are widely used in parallel computing systems, which have many advantageous properties such as lower node degree and diameter, higher connectivity, symmetry, and etc [START_REF] Lee | Logic and computer design in nanospace[END_REF], [START_REF] Schlosser | HyperCuP-hypercubes, ontologies, and efficient search on peer-to-peer networks[END_REF]. Furthermore, by changing their links between some nodes, the variants of hypercubes, such as Möbius cubes [START_REF] Cull | The Möbius cubes[END_REF], crossed cubes [START_REF] Kulasinghe | Multiply-twisted hypercube with five or more dimensions is not vertex-transitive[END_REF], and twisted cubes [START_REF] Abraham | The twisted cube topology for multiprocessors: a study in network asymmetry[END_REF] were proposed, which have better properties [START_REF] Fan | Hamilton-connectivity and cycle-embedding of the Möbius cubes[END_REF], [START_REF] Fan | Embedding meshes into crossed cubes[END_REF], [START_REF] Fan | Optimal embeddings of paths with various lengths in twisted cubes[END_REF], [START_REF] Hsieh | Pacyclicity on Möbius cubes with maximal edge faults[END_REF], [START_REF] Xu | Paths in Möbius cubes and crossed cubes[END_REF].

Independent spanning trees(ISTs for short) have been used in reliable broadcasting, secure message distribution [START_REF] Bao | Reliable broadcasting in product networks[END_REF], reliable communication protocols [START_REF] Itai | The multi-tree approach to reliability in distributed networks[END_REF], one-to-all broadcasting [START_REF] Tseng | Efficient broadcasting in wormholerouted multicomputers: A network-partitioning approach[END_REF], the multi-node broadcasting [START_REF] Chen | Multi-node broadcasting in all-ported 3-D wormhole-routed torus using an aggregation-then-distribution strategy[END_REF], and diagnosis [START_REF] Cheng | Parallel construction of independent spanning trees and an application in diagnosis on Möbius Cubes[END_REF]. Therefore, the problem to construct ISTs for a given network is becoming an important issue.

However, there is a well-known conjecture on the existence of ISTs for any network [START_REF] Itai | The multi-tree approach to reliability in distributed networks[END_REF][31]:

Corresponding author.

Conjecture 1. Let G be an n-node-connected network with n ≥ 1. Then, there exist n node-independent spanning trees rooted at any node for G.

In what follows, we use independent to represent node-independent. For n ≤ 4, Conjecture 1 was solved [START_REF] Cheriyan | Finding nonseparating induced cycles and independent spanning trees in 3-connected graphs[END_REF], [START_REF] Curran | Finding four independent trees[END_REF], [START_REF] Itai | The multi-tree approach to reliability in distributed networks[END_REF], [START_REF] Zehavi | Three tree-paths[END_REF], but when n ≥ 5, it has remained open. Consequently, researchers are interested in the study of ISTs for various special networks. Conjecture 1 has been solved for some restricted classes of networks, such as planar networks [START_REF] Huck | Independent trees in planar graphs[END_REF], product networks [START_REF] Obokata | Independent spanning trees of product graphs and their construction[END_REF], hypercubes [START_REF] Tang | Optimal independent spanning trees on hypercubes[END_REF], [START_REF] Werapun | An efficient parallel construction of optimal independent spanning trees on hypercubes[END_REF], [START_REF] Yang | Parallel construction of optimal independent spanning trees on hypercubes[END_REF], Möbius cubes [START_REF] Cheng | Constructive algorithm of independent spanning trees on Möbius cubes[END_REF], [START_REF] Cheng | Parallel construction of independent spanning trees and an application in diagnosis on Möbius Cubes[END_REF], locally twisted cubes [START_REF] Liu | Constructing independent spanning trees for locally twisted cubes[END_REF], crossed cubes [START_REF] Cheng | Independent spanning trees in crossed cubes[END_REF], [START_REF] Cheng | Dimension-adjacent trees and parallel construction of independent spanning trees on crossed cubes[END_REF], twisted cubes [START_REF] Wang | Independent spanning trees on twisted cubes[END_REF], even networks [START_REF] Kim | Independent spanning trees on even networks[END_REF], odd networks [START_REF] Kim | Optimal independent spanning trees on odd graphs[END_REF], and etc.

We say that a sequence of n integers is a permutation if it contains all integers from 0 to n -1 exactly once. Considering the results for hypercubes and Möbius cubes, each paper in literature only considers a set of ISTs for the special network and lacks the discussion of the relation between the permutations and ISTs. Question 1. Can all permutations of 0, 1, . . . , n-1 be used to construct spanning tree and ISTs for the n-dimensional hypercube Q n and the n-dimensional Möbius cube M n ?

To solve this question, we adopt the definition of circular dimensional-permutation and prove that any circular dimensional-permutation of 0, 1, . . . , n -1 can be used to construct n ISTs for Q n and M n in this paper, which is the further discussion of spanning trees and ISTs for Q n and M n comparing with the results in literature.

The rest of this paper is organized as follows. Section 2 presents some definitions and graph terminologies and notations. Section 3 discusses the IST problem for Q n and M n rooted at any node. We draw the conclusion of the paper in the last section.

Preliminaries

Definition of Hypercubes, Möbius Cubes, and ISTs

We use a unique binary string of length n to denote the address of each node in the n-dimensional hypercube Q n and the n-dimensional Möbius cube M n . In what follows, nodes and their addresses will be used alternatively. Q n is a network consists of 2 n nodes. Any two nodes of Q n are adjacent whenever their corresponding addresses differ in exactly one place.

M n is a variant of the Q n , which has two types, 0-type n-dimensional Möbius cube and 1-type n-dimensional Möbius cube. We adopt the following definition of M n in [START_REF]Diagnosability of the Möbius Cubes[END_REF]. Definition 1. [START_REF]Diagnosability of the Möbius Cubes[END_REF] 0-M 1 and 1-M 1 are both the complete graph on two nodes whose addresses are 0 and 1. For any integer n with n ≥ 2, both 0-M n and 1-M n contain one 0-type (n -1)-dimensional sub-Möbius cube M 0 n-1 and one

1-type (n -1)-dimensional sub-Möbius cube M 1 n-1 . The nodes in M 0 n-1 have a common prefix 0; the nodes in M 1 n-1 have a common prefix 1. For two nodes x = x n-1 x n-2 . . . x 0 ∈ V (M 0 n-1 ) and y = y n-1 y n-2 . . . y 0 ∈ V (M 1 n-1 ), where x n-1 = y n-1 = 0, (1) (x, y) ∈ E(0-M n ) if and only if x i = y i , i = 0, 1, . . . , n -2; (2) (x, y) ∈ E(1-M n ) if and only if x i = y i , i = 0, 1, . . . , n -2.
A binary string x of length n is denoted by

x n-1 x n-2 . . . x 0 . Suppose that u = u n-1 u n-2 . . . u 0 and v = u n-1 u n-2 . . . u l u l-1 v l-2 v l-3 . . . v 0 are two nodes in X n ∈ {Q n , M n }.
We say that u and v have a leftmost differing bit at position l -1. We use LDF(u, v) to denote the leftmost differing bit of two nodes u and v. Given two adjacent nodes u and v, if LDF(u, v)=d, we say that v is the d-neighbor of u or that the edge (u, v) is an edge of dimension d. For this purpose, let N d (u) denote the d-neighbor of u. We follow the definitions of path and ancestor in [START_REF] Cheng | Dimension-adjacent trees and parallel construction of independent spanning trees on crossed cubes[END_REF].

Two paths P and P starting from a node u and ending with another node v are said to be internally disjoint if E(P )∩E(P ) = ∅ and V (P )∩V (P ) = {u, v}. Two spanning trees for a network G are independent if they are rooted at the same node, said u, and for each node v ∈ V (G)\{u}, the two paths starting at u and ending with v are internally disjoint. A set of spanning trees of G rooted at v are called independent spanning trees if they are pairwisely independent.

Definition of Dimensional-Permutation

The sequence of n integers is called a dimensional-permutation if it contains all integers from 0 to n -1 exactly one (Noting that each node in Q n or M n has n neighbors, which are 0-neighbor, 1-neighbor, . . ., (n -1)-neighbor). A circular dimensional-permutation (CDP for short) is a type of permutation to put all integers from 0 to n -1 along a closed circle in the clockwise order [START_REF] Cheng | Parallel construction of independent spanning trees and an application in diagnosis on Möbius Cubes[END_REF]. Suppose that {a 0 , a 1 , . . . a n-1 }= {0, 1, . . . , n -1}. All cyclic permutations of integers are equivalent in the circle,

a 0 → a 1 → . . . → a n-1 , a 1 → a 2 → . . . → a n-1 → a 0 , . . . , a n-1 → a 0 → a 1 → . . . → a n-2
belong to the same CDP. The total number of CDPs is n!/n = (n -1)!. For example, the six CDPs of 0, 1, 2, 3 are

3 → 0 → 1 → 2, 3 → 0 → 2 → 1, 3 → 1 → 0 → 2, 3 → 1 → 2 → 0, 3 → 2 → 0 → 1, 3 → 2 → 1 → 0.

An Reliable Broadcasting Algorithm Based on ISTs for Hypercubes and Möbius Cubes

In this section, we point out that every CDP of n integers 0, 1, 2, . . . , n -1 can be used to construct ISTs for X n ∈ {Q n , M n }. We now present the following observation.

Observation 1. For the n-dimensional Möbius cube, we proved the correctness of ISTs rooted at any node based on the descending CDP n -1, n -2, . . . , 0 [START_REF] Cheng | Parallel construction of independent spanning trees and an application in diagnosis on Möbius Cubes[END_REF].

In essence, the set of optimal ISTs for Q n in [START_REF] Tang | Optimal independent spanning trees on hypercubes[END_REF] can be obtained by the ascending CDP 0, 1, . . . , n -1; the set of ISTs in [START_REF] Yang | Parallel construction of optimal independent spanning trees on hypercubes[END_REF] is similar to that in [START_REF] Werapun | An efficient parallel construction of optimal independent spanning trees on hypercubes[END_REF] for Q n , which can be constructed by the descending CDP n -1, n -2, . . . , 0. Thus, the result in this section is the further study of spanning trees and ISTs for Q n and M n .

ISTs for Q n and M n with any Circular Dimensional-Permutation

In what follows, we always let u denote any node in Q n or M n . Now we present an algorithm, called BCIST, to construct n ISTs rooted at an arbitrary node u for X n ∈ {Q n , M n }. Fig. 1 demonstrates the construction procedures of n spanning trees

T 0 , T 1 , . . . , T n-1 rooted at node u for X n ∈ {Q n , M n } in radial style.
Algorithm BCIST Input: An array S = {a 0 ,a 1 ,..., a n-1 }, where the permutation a 0 , a 1 , ...,

a n-1 is a CDP of 0, 1, ..., n -1; an arbitrary node u in X n ∈ {Q n , M n }; Output: T 0 ,T 1 ,..., T n-1 rooted at u for X n ; Begin 1: for i = 0 to n -1 do in parallel 2: V (T i ) = N ai (u); /*N ai (u) denotes the a i -neighbor of u.*/ 3: E(T i ) = ∅; 4:
for l = 0 to n -1 do 5:

for each node v ∈ V (T i ) do in parallel 6:

d = S[(i + l + 1) mod n]; /*Indexing of S is counted from 0 to n -1,
where n is length of the array S.*/ 7:

E(T i ) = E(T i ) ∪ {(v, N d (v)}; 8: V (T i ) = V (T i ) ∪ {N d (v)}; 9:
end for 10: end for 11: end for end

The construction procedures of T 0 , T 1 , . . . , T n-1 are similar. Take T 0 for example, the construction procedures are described as follows (See Fig. 1).

At first, there is only one node N a0 (u) in tree T 0 ; during the 1st iteration (l =0), node

N a0 (u) in T 0 is connected to its a 1 -neighbor node N a1 (N a0 (u)). There- fore, V (T 0 ) = {N a0 (u), N a1 (N a0 (u))} and E(T 0 ) = {(N a0 (u), N a1 (N a0 (u)))};
during the 2nd iteration (l = 1), each node v in T 0 is connected to its a 2neighbor node N a2 (v). Thus, the edges (N a1 (N a0 (u)), N a2 (N a1 (N a0 (u)))) and (N a0 (u), N a2 (N a0 (u))) are appended to T 0 and V (T 0 ) = {N a0 (u), N a2 (N a0 (u)), N a1 (N a0 (u)), N a2 (N a1 (N a0 (u)))}; during the 3rd iteration, each node v in T 0 is connected to its a 3 -neighbor node N a3 (v). As a result, it has doubled the number of nodes in T 0 ; During the l-th iteration, each node v in T 0 is connected to its

S[l + 1]-neighbor node N S[l+1] (v) with 4 ≤ l ≤ n -2; in the last iteration, each node v in T 0 is connected to its a 0 -neighbor node N a0 (v).
Consequently, T 0 is a spanning tree for

X n ∈ {Q n , M n }.
More examples will be shown in the next subsection. Now we give the following lemma about the relation of adjacent nodes in M n and Q n . Lemma 1. For any node x n-1 x n-2 . . . x 0 and its k-neighbor node

y n-1 y n-2 . . . y 0 with 0 ≤ k ≤ n-1 in X n ∈ {M n , Q n }, we have x n-1 x n-2 . . . x k+1 = y n-1 y n-2 . . . y k+1 and x k = y k .
Based on Definition 1, the definition of Q n , and Lemma 1, we have the following lemma. Lemma 2. For any two nodes x, y in V (X

xn-1 n-1 ), if LDF(x, y) = k with 0 ≤ k ≤ n -2 and X xn-1 n-1 ∈ {M xn-1 n-1 , Q xn-1 n-1 }, then LDF(N n-1 (x), N n-1 (y))=k and N n-1 (x), N n-1 (y) ∈ X xn-1 n-1 . Lemma 3.
[5] Given a walk W : u (0) , u (1) = N m1 (u (0) ), u (2) = N m2 (u (1) ), . . . ,

u (k) = N m k (u (k-1) ) in M n for any integer k with 1 ≤ k ≤ n, if m 1 , m 2 , . . . , m k differ from one another and 0 ≤ m i ≤ n -1 for i = 1, 2, . . . , k, then W is a path.
For the convenience of proof, we define a vector < β 1 , β 2 , . . . , β n > such that the set {β 1 , β 2 , . . . , β n } equals to the set {0, 1, . . . , n -1}. Lemma 4. Suppose that P : u (0) , u (1) = N a1 (u (0) ), u (2) = N a2 (u (1) ), . . . ,

u (k) = N a k (u (k-1) ) and P : u (0) , v (1) = N a1 (u (0) ), v (2) = N a2 (v (1) ), . . . , v (t) = N at (v (t-1) ) are two paths in X n ∈ {Q n , M n } for any two integers k, t with 1 ≤ k, t ≤ n and a 1 = a 1 . If the following conditions hold: (1) < β i1 , β i2 , . . . , β i k >=< a 1 , a 2 , . . . , a k >, where 1 ≤ i 1 < i 2 < . . . < i k ≤ n; (2) < β j1 , β j2 , . . . , β jt >=< a 1 , a 2 , .
. . , a t >, where 1 ≤ j 1 < j 2 < . . . < j t ≤ n, then V (path(P, u (1) , u (k) )) ∩ V (path(P , v (1) , v (t) )) = ∅.

Proof. Suppose that there exists a node v, such that v ∈ V (P ) ∩ V (P ). We denote v as v = u (i) and v = v (j) with 1 ≤ i ≤ k and 1 ≤ j ≤ t. Let A = {a 1 , a 2 , . . . , a i } and B = {a 1 , a 2 , . . . , a j }. We have M 1 = max((A∪B)\(A∩B)). By Lemma 2, we can verify that the M 1 -bit of v in P is different from that of v in P , which is a contradiction.

Based on Lemma 4, we have the following corollary. Corollary 1. Suppose that P : u (0) , u (1) = N a1 (u (0) ), u (2) = N a2 (u (1) ), . . . , u (k) = N a k (u (k-1) ) = v and P : v (0) , v (1) = N a1 (u (0) ), v (2) = N a2 (v (1) ), . . . , v (t) = N at (v (t-1) ) = v are two paths in X n ∈ {Q n , M n } for any two integers k, t with 1 ≤ k, t ≤ n and a k = a t . If the following conditions hold:

(1) < β i1 , β i2 , . . . , β i k >=< a 1 , a 2 , . . . , a k >, where 1 ≤ i 1 < i 2 < . . . < i k ≤ n; (2) < β j1 , β j2 , . . . , β jt >=< a 1 , a 2 , . . . , a t >, where 1 ≤ j 1 < j 2 < . . . < j t ≤ n, then V (path(P, u (0) , u (k) )) ∩ V (path(P , v (0) , v (t) )) = {v}.
Lemma 5. Let the input a 0 , a 1 , . . . , a n-1 of algorithm BCIST be any dimensionalpermutation of integers 0, 1,. . . ,n -1 and X n ∈ {Q n , M n }. T i obtained by Algorithm BCIST is a spanning tree for X n with integer i = 0, 1, . . . , n -1.

Proof. Without loss of generality, we consider the tree T 0 obtained by algorithm BCIST. After the n iterations, we have 1 + 2 0 + 2 1 + ... + 2 n-1 = 2 n nodes in tree T 0 . Choosing arbitrary two nodes v (1) and v (2) from T 0 , the < N a0 (u), v (1) >path can be denoted by N a0 (u), x (1) = N ai 1 (N a0 (u)), x (2) = N ai 2 (x (1) ), . . . ,

x (k) = N ai k (x (k-1)
) and the < N a0 (u), v (2) >-path can be denoted by N a0 (u), x (1) = N aj 1 (N a0 (u)), x (2) = N aj 2 (x (1) ), . . . , x (m) = N aj m (x (m-1) ), where 0 ≤ i 1 < i 2 < . . . < i k ≤ n -1 and 0 ≤ j 1 < j 2 < . . . < j m ≤ n -1. By algorithm BCIST, it is easy to verify < N a0 (u), v (1) >-path and < N a0 (u), v (2) >-path satisfy the conditions in Lemma 4, which implies that v (1) = v (2) . Thus, we can say that T 0 is a spanning tree rooted at N a0 (u) for X n ∈ {Q n , M n }. Noting that u is the child of node N a0 (u) and the leaf node in T 0 , that is, T 0 is a spanning tree rooted at u. Lemma 6. [START_REF] Cheng | Independent spanning trees in crossed cubes[END_REF] Let T and T be two spanning trees rooted at node u for a network G. T and T are independent if and only if for every node v ∈ V (G)\{u}, ancestor (v, T ) ∩ ancestor (v, T ) = {u} and ancestor (v, T ) ∪ ancestor (v, T ) ⊃ {u}.

Lemma 7. T 0 , T 1 , . . . , T n-1 obtained by Algorithm BCIST are n ISTs for X n ∈ {Q n , M n }.

Proof. We have the following two cases.

Case 1. X n is M n . By Lemma 5, T i obtained by Algorithm BCIST is a spanning tree for V (M n ) for integer i = 0, 1, . . . , n -1. We only need to prove that T i and T j are independent for 0

≤ i ≤ j ≤ n -1.
The trivial cases is n = 1 and n = 2. Now we consider n with n ≥ 3. The longest path P 1 in T i and the longest path P 2 in T j can be denoted by

P 1 : u, x 0 = N ai (u), x 1 = N ai+1 (x 0 ), . . . , x n-i-1 = N an-1 (x n-i-2 ), x n-i = N a0 (x n-i-1 ), x n-i+1 = N a1 (x n-i ), . . . , x n = N ai (x n-1 ) and P 2 : u, y 0 = N aj (u), y 1 = N aj+1 (y 0 ), . . . , y n-j-1 = N an-1 (x n-j-2
), y n-j = N a0 (y n-j-1 ), y n-j+1 = N a1 (y n-j ), . . . , y n = N aj (y n-1 ), respectively, where 0 ≤ i < j ≤ n -1.

Let a c0 = a i and a d0 = a j . By lemma 6, we only need to prove that for any v ∈ V (T i ) ∩ V (T j ), ancestor (v, T i ) ∩ ancestor (v, T j ) = {u} and ancestor (v, T i ) ∪ ancestor (v, T j ) ⊃ {u}. Any path in T i and any path in T j can be denoted by P 3 and P 4 , respectively, as follows.

P 3 : u, x 0 = N ac 0 =ai (u), x 1 = N ac 1 (xi) , . . . , x k = N ac k (x k-1 ), x k+1 = N ac 0 (x k ) and P 4 : u, y 0 = N a d 0 =aj (u), y 1 = N a d 1 (y i ), . . . , y l = N a d l (y l-1
), y l+1 = N a d 0 (y l ), Without loss of generality, suppose that a j > a i and v ∈ V (P 3 ) ∩ V (P 4 ). Let v = a cw = a dz where 1 ≤ u ≤ k and 1 ≤ z ≤ m. Based on P 1 , P 2 , P 3 , and P 4 , we define walks W 1 , W 2 , W 3 , and W 4 as follows.

W 1 : a i , a i+1 , . . . , a n-1 , a 0 , a 1 , . . . , a i , W 2 : a j , a j+1 , . . . , a n-1 , a 0 , a 1 , . . . , a j , W 3 : a c0 = a i , a c1 , . . . , a cw , and W 4 : a d0 = a j , a d1 , . . . , a dz . Since a c0 = a d0 , we have x 0 = y 0 , which implies that ancestor (v, P 3 ) ∪ ancestor (v, P 4 ) ⊃ {u}. We only need to prove that ancestor (v, P 3 ) ∩ ancestor (v, P 4 ) = {u}. We have the following Cases.

Case 1.1. max(V (W 3 )) = max(V (W 4 )). Since max(V (W 4 )) ≥ a j , then we have the following subcases.

Case 1.1.1. max(V (W 3 )) < a j and max(V (W 4 )) = a j . Then, each node in V (P 4 )\{y m+1 } and each node in V (P 3 ) have a leftmost different bit at position a j . Then we have v = y m+1 and ancestor (v,

P 3 ) ∩ ancestor (v, P 4 ) = {u}. Case 1.1.2. max(V (W 3 ∪ W 4 )) > a j . Then, the max(V (W 3 ∪ W 4 ))-bit of v in P 3 is different from that of v in P 4 . It is a contradiction. Case 1.2. max(V (W 3 )) = max(V (W 4 ))
. Then, we have the following cases. Case 1.2.1. max(V (W 3 )) = max(V (W 4 )) = a j . We can verify that v = y m+1 and a j ∈ {a d1 , a d2 , . . . , a dm }. We can divide path ancestor (v, P 3 ) into P 31 and P 32 as follows.

P 31 : u, x 0 = N ac 0 (u),

x 1 = N ac 1 (x 0 ), . . . , ,x f -1 = N ac f -1 (x f -2
) and P 32 :

x f = N ac f =aj (x f -1 ), x f +1 = N ac f +1 (x f ), . . . , x w = N ac w (x u-1 )
. By Lemma 2, each node in V (P 31 ) and each node in V (P 4 ) have a leftmost different bit at position a j . Then, V (P 4 ) ∩ V (P 31 ) = ∅. By Lemma 2, LDF(x f , y 0 )= max({a c0 = a i , a c1 , . . . , a c(f -1) }). Let A = {a c f +1 , a c f +2 , . . . , a cw } and B = {a d1 , a d2 , . . . , a dz }. Furthermore, we have the following cases.

Case

1.2.1.1. max((A ∪ B)\(A ∩ B)) > LDF(x f , y 0 ). Then, by Lemma 2, the max((A ∪ B)\(A ∩ B))-bit of v in P 3 is different from that of v in P 4 , which is a contradiction. Case 1.2.1.2. max((A∪B)\(A∩B)) = LDF(x f , y 0 )
. By Algorithm BCIST, T i and T j are constructed based on the same CDP. Then, we can verify that a cw = a dz . By Corollary 1,

V (P 4 ) ∩ V (P 32 ) = {v}. Case 1.2.1.3. max((A ∪ B)\(A ∩ B)) < LDF(x f , y 0 ). Then, by Lemma 2, the LDF(x f , y 0 ) bit of v in P 3 is different from that of v in P 4 , which is a contradiction. Case 1.2.2. max(V (W 3 )) = max(V (W 4 )) = M 1 > a j .
We can divide path sub-path < u, v >-path of P 3 into P 31 and P 32 as follows.

P 31 : u, x 0 = N ac 0 (u), x 1 = N ac 1 (x 0 ), . . . , ,x f -1 = N ac f -1 (x f -2 ) and P 32 :

x f = N a c f =M 1 (x f -1 ), x f +1 = N ac f +1 (x f ), . . . , x w = N ac w (x u-1 ).
We can divide path sub-path < u, v >-path of P 4 into P 41 and P 42 as follows.

P 41 : u, y j = N a d 0 =aj (u), y 1 = N a d 1 (y i ), . . . , y h-1 = N a d h-1 (y h-2 )
, and P 42 : y h = N a d h (y h-1 ), y h+1 = N a d h+1 (y h ), . . . , y z = N az (y z-1 ), By Lemma 2, V (P 31 ) ∩ V (P 41 ) = {u}. Clearly, the M 1 -bit of each node in V (P 42 ) is different from each node in V (P 31 ), thus V (P 42 ) ∩ V (P 31 ) = ∅. Similarly, we have V (P 32 ) ∩ V (P 41 ) = ∅. We only need to prove that V (P 32 ) ∩ V (P 42 ) = ∅. Similarly to Case 1.2.1, we can verify that a cw = b dz . By Corollary 1, V (P 32 ) ∩ V (P 42 ) = {v}.

As a consequence, we have ancestor (v, P 3 ) ∩ ancestor (v, P 4 ) = {u}. Case 2. X n is Q n . The proof is similar to that of Case 1. Based on the above discussion, by Lemma 6, the lemma holds.

Since there are (n -1)! CDPs of 0, 1, . . . , n -1, based on Lemma 7, we have the following theorem.

Theorem 1. Based on Algorithm BCIST, there are (n -1)! sets of ISTs for

X n ∈ {Q n , M n }.
Comparing with the result in [START_REF] Tang | Optimal independent spanning trees on hypercubes[END_REF], all the (n-1)! sets of ISTs can provide optimal reliable broadcasting for Q n . As far as the symmetry is concerned, all the (n -1)! sets of ISTs can also provide optimal reliable broadcasting for M n .

Simulation Experiments of ISTs and Disjoint Paths for 0-M 4

As it is well-known that JUNG framework is a software library which provides a common and extendible language for the modeling, analysis, and visualization of data that can be represented as a network. Using Java multi-thread and JUNG framework, we can easily construct ISTs for Q n and M n . For example, Fig. 2 illustrates some construction procedures of four ISTs rooted at 1 for 0-M 4 . The four trees can be obtained by four steps with parallel fashion.

Based on different CDPs, we can obtain multiple sets of ISTs. Thus, there are different sets of disjoint paths between arbitrary two nodes. Fig. 3(a)-(f) show six sets of ISTs rooted at 1 for 0-M 4 based on the six CDPs mentioned in Section 2.2. Clearly, there are six sets of n disjoint paths between node 1 and any node v ∈ V (0-M 4 )\{1}, among which the total length of four disjoint paths between node 1 and node v may be different. For example, the four disjoint paths between nodes 1 and 13 in Fig. 3(a 1→ 0→ 2→ 10→ 13, 1→ 3→ 11→ 12→ 13, 1→ 5→ 13, and 1→ 9→ 14→ 13. The total length of four disjoint paths between nodes 1 and 13 in Fig. 3(c) is shorter than that in Fig. 3(a). Thus, basing on CDPs, we may easily obtain a set of n disjoint paths with better performance between some pairs of nodes.

Conclusions

In this paper, we study the circular dimensional-permutations and reliable broadcasting based on ISTs for the n-dimensional hypercube Q n and the n-dimensional Möbius cube M n . We prove that any circular dimensional-permutation of 0, 1, . . . , n -1 can be used to construct n ISTs, which is the further discussion of IST problem for Q n and M n comparing with the results in literature. Moreover, we also conduct simulation experiments of ISTs based on JUNG framework and list some disjoint paths for 0-M 4 .

Acknowledgment

This

Fig. 1 .

 1 Fig. 1. Construction procedures based on algorithm BCIST.

Fig. 2 .

 2 Fig. 2. Simulation experiments of ISTs based on JUNG framework. (a) Initialization; (b) The 2nd iteration; (c) The 3rd iteration; (d) The 4th iteration.

  Fig. 3. Example of ISTs rooted at node 1 and disjoint paths between nodes 1 and 13 in 0-M4.

		0				0	1					3		1					5		1				9		1
		3		8		4	2			2	11		7			6		4		13			14	11		8
	(a)		15	10	9		7	5	3	10	6	0		12		9	4	7	14	2	12	0	10	13	15	6	10		3	0
		2	12	14	11		6		13	8	5		15		15	3	9			11		12	5	7		2
		1	13						14								8							4
		0				0		1					3		1					5		1				9		1
		3		8		2		4			7	2	11			4	13		6			11	14		8
	(b)		10	15	9		6	3	5	6	15	4		10		0	9	12	7	0		14	10	2	12	10	3	15	6	0
		1	13	11	14		7		14	5	12		8			15	11	3			9		13	4	2		7
		2	12						13								8							5
		0				0		1					3		1					5		1				9		1
	(c)	1	10	2 6	3	8	15	4 9	5	12	11 10	9	7	6	2	4	0	7	4 12	0	6	14	13 2	10	15	14 13	11 6	10	8	0	3
		3	13	11	7		14		13	15	8		5			15	3	11		9			12	7	5		2
		2	12						14								8							4
		0				0		1					3		1					5		1				9		1
		1		2		4		8			7	11		2			13		4		6			8	11	14
	(d)		6	10	3		12	5	9	15	6	4		10		9	0	12	14	10	7	0	2	10	15	0	12		3	6
		2	14	7	11		13		14	12	5		8			15	11	9			3			13	2	7		4
		3	15						13								8							5
		0				0		1					3		1					5		1				9		1
		2		4		8		2			2		7	11		13		6		4			11	8	14
	(e)		12	7	5		10	9	3	6	10	0		15		4	9	14	12	10	7	2	0	10	12	3	15		0	6
		3	15	13	6		11		14	5		8		12		15	9	11		3			13	2	4		7
		1	14						13								8							5
		0				0	1					3		1					5		1				9		1
		2		4		2		8			11		2		7			6	13		4			8	14	11
	(f)		7	12	5		10	3	9	10	12	9		6		0	4	14	7	2	12	10	0	15	10	0	13		6	3
		1	15	6	13		11		13	8	15		5			15	9	3			11		12	7	2		5
		3	14						14								8							4

  work is supported by National Natural Science Foundation of China (No. 61170021), Specialized Research Fund for the Doctoral Program of Higher Education (No. 20103201110018), Application Foundation Research of Suzhou of China (SYG201240), the 2011 Program for Postgraduates Research Innovation in University of Jiangsu Province (No. CXZZ11 0100), the 2012 Science and technology innovation team building program of Soochow University (SDT2012B02) and sponsored by Qing Lan Project.


Yamin Li

Wanming Chu

Total Exchange Routing on Hierarchical Dual-Nets

Keywords:

ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

Recently, because of the advances in computer and networking technologies, supercomputers containing hundreds of thousands of nodes have been built [START_REF]Supercomputer Sites[END_REF]. It was predicted that the parallel systems of the next decade will contain 10 to 100 millions of nodes [START_REF] Beckman | Looking toward exascale computing, keynote speaker[END_REF]. The interconnection network plays an important role for achieving high-performance in such ultra-scale parallel systems. The performance of an ultra-scale parallel computers depends largely on the time complexities of communication schemes, and in turn depends on the diameter of the network. An interconnection network consists of switches with multiple communication ports and cables connecting ports by following certain topologies. For an ultrascale parallel computer, the traditional interconnection networks may no longer satisfy the requirements for the high-performance computations or efficient communications. For such an ultra-scale parallel computer, the node degree and the diameter will be the critical measures for the effectiveness of the interconnection networks. The node degree is limited by the hardware technologies and the diameter affects all kinds of communication schemes directly. The number of communication ports (node degree) in the network-on-chip (NoC) is typically 4 to 8 in current implementations. The off-chip interconnect switches can have tens of ports, but the cost becomes expensive as the number of ports increases. Other important measures for the effectiveness of the interconnection networks include symmetricity, scalability, and efficient routing algorithms.

The following two categories of interconnection networks have attracted a great research attention and been used in many supercomputers' implementations. One is the hypercube-like family that has the advantage of short diameters for high-performance computing and efficient communications. The other is the 2D/3D mesh or torus family that has the advantage of small and fixed node degrees and easy implementations [START_REF] Adiga | Blue gene/l torus interconnection network[END_REF]. Traditionally, most supercomputers including those built by CRAY, IBM, SGI, and Intel use 3D tori or hypercubes. However, the node degree of the hypercube increases logarithmically as the number of nodes in the systems increases; the diameter of the 2D/3D torus becomes large in an ultra-scale parallel system. To solve these problems, the hierarchical (cluster-based) architectures are proposed in literature [START_REF] Ghose | Hierarchical cubic networks[END_REF][START_REF] Li | Efficient collective communications in dual-cube[END_REF]. The supercomputer Roadrunner built by IBM adopts a new approach for the interconnection network [START_REF]IBM: Roadrunner: Hardware and Software Overview[END_REF]. It is a cluster-based architecture: the connection among clusters is fully connected, and the fat-tree is used for the connection inside a cluster.

In this paper, we first present a flexible interconnection network, called Hierarchical Dual-Net (HDN) [START_REF] Li | Hierarchical dual-net: A flexible interconnection network and its routing algorithm[END_REF]. The HDN is symmetric and can connect a large number of nodes with a small node degree, meanwhile keeping the diameter short. The HDN was motivated by recursive dual-net (RDN) [START_REF] Li | Recursive dual-net: A new universal network for supercomputers of the next generation[END_REF]. The RDN can be viewed as a special case of HDN. The RDN has merits of low node degree and short diameter. The problem of the RDN is that it grows too fast in size, and there is no mechanism to control the rate of its growth. Different from the RDN, the scale of the HDN can be controlled by setting a set of suitable parameters while generating an expanded network through dual-construction. The HDN also adapts the cluster-based architecture. Compared to the Roadrunner, the HDN is symmetric, uses small number of links, and meanwhile keeps the diameter short. The HDN structure is also better than other popular existing networks such as hypercube and 2D/3D torus with respect to the degree and diameter. We investigate the topological properties of the HDN and show some examples of HDNs with simple base networks of small size. Then we compare them to other networks such as three-dimensional torus used in IBM Blue Gene/L [START_REF] Adiga | Blue gene/l torus interconnection network[END_REF], and hypercube. The total exchange, or all-to-all personalized communication, is one of the most dense communication patterns and is at the heart of numerous applications and programming models in parallel computing. In this paper, we present an efficient total exchange routing algorithm on a hierarchical dual-net. The time complexity T k (m) of the algorithm for an HDN(B, k, S) is

T k (m) = (2 k+1 -2)(t s + t w mn k /2) + 2 k T 0 (mn k /2)
, where n k is the total number of nodes, t s is startup latency, m is the message length in words, and t w is the per-word transfer time.

The rest of this paper is organized as follows. Section 2 introduces the hierarchical dual-net in details. Section 3 describes the routing algorithm. Section 4 gives the total exchange routing algorithm on a hierarchical dual-net. Section 5 concludes the paper.

The Hierarchical Dual-Net

We begin with a brief introduction to the recursive dual-net (RDN). The RDN is constructed recursively by a dual-construction. The dual-construction is a way to expand a given symmetric graph G of size n to a new symmetric graph G * of size 2n 2 . It generates 2n copies of G as subgraphs (denoted as clusters) of G * . Half of them, n clusters, are of class 0 and the others are of class 1. The connection method is described below.

If G is symmetric then the expanded graph G * is unique and symmetric. Therefore, the dual-construction can be applied recursively from a symmetric network (the base network). RDN(m, k) denotes an RDN generated from a base network of size m by applying dual-construction k times. The problem about an RDN is that its growth rate is super-exponential ((2m) 2 k ). There is very little space for selection of the size of an RDN. For example, let the base network be a 3-cube, then the sizes of RDN(8, k) will be 2 7 , 2 15 , and 2 31 for k = 1, 2, and 3, respectively. In HDN, we provide a mechanism to control the growth rate through its expansion from a base network. This new interconnection network has a very flexible way for adjusting its size.

The hierarchical dual-net, HDN(B, k, S), contains three sets of parameters: B is a symmetric product graph, we call it base network; k is an integer that indicates the level of the HDN (the number of dual-constructions applied); and

S = {G 1 , G 2 , . . . , G k }, where G i is a sub-graph of HDN(B, k -1, S) and s i = |G i |
is the number of nodes in a super-node at the level i for 1 ≤ i ≤ k. All these terminologies will be defined in the following paragraphs.

Given r graphs

G i = (V i , E i ), 1 ≤ i ≤ r, their product graph G = G 1 × G 2 × . . . × G r is defined as the graph G = (V, E), where V = {(v 1 , v 2 , . . . , v r )|v i ∈ V i , 1 ≤ i ≤ r} and E = {[(u 1 , u 2 , .
. . , u r ), (v 1 , v 2 , . . . , v r )]| for some j, (u j , v j ) ∈ E j and for i = j, u i = v i }.

In other words, the nodes of the product graph G are labeled with r-tuples, where the ith element of the r-tuples is chosen from the node set of the ith component graph. The edges of the product graph connect pairs of nodes whose labels are identical in all but the jth element, and the two nodes corresponding to the jth elements in the jth component graph are connected by an edge.

Meshes/tori or hypercubes are typical examples of product graphs. Given a product graph G = G 1 ×G 2 ×. . .×G r , we define a quotient graph

Q as Q = G/G where G is a sub-product graph of G such that G = G ×Q. A node in a product graph G = G 1 ×. . .×G i ×. . .×G r can be represented by (a 1 , . . . , a i , . . . , a r ) with 0 ≤ a i ≤ |G i |-1. We define a sub-graph G as G = G 1 ×. . .×G j ×. . .×G q with G j = G i for 1 ≤ j ≤ q ≤ r and 1 ≤ i ≤ r, G j = G k if j = k for 1 ≤ j, k ≤ q.
Then a node in the sub-graph G can be represented by (b 1 , . . . , b i , . . . , b q ) with 0 ≤ b i ≤ |G i | -1. We can consider a quotient graph Q as a reduced graph of G with G being mapped into a single node (a super-node).

A graph G is symmetric (node-symmetric) if all its nodes looks alike. A product graph is symmetric if all its component graphs are symmetric. We use the symmetric product graph as the base network for generating a hierarchical dual-net through dual-constructions. We denote the base network as B = B 1 × B 2 × . . . × B r where all the B i , 1 ≤ i ≤ r, are symmetric. We define a super-node of B, denoted as SN as a sub-product graph of B. That is, SN = B i1 × B i2 × . . . × B iq , where i j , 1 ≤ j ≤ q, are distinct and q ≤ r.

Let |B i | = b i be the number of nodes in B i for 1 ≤ i ≤ r. The HDN(B, 0, S) = B is the base network. For i > 0, the HDN(B, i, S) is generated from HDN(B, i-1, S) by a construction to be explained below. Note that S = {G 1 , G 2 , . . . , G k }, where G i is a sub-graph of HDN(B, k -1, S) and s i = |G i | is the number of nodes in a super-node at the level i for 1 ≤ i ≤ k. First, we define a super-node of level i, denoted as SN i , to be a sub-product graph G i of size s i in B. Then, we define graph Q i as the quotient graph HDN(B, i -1, S)/SN i . Suppose that there are N i-1 nodes in the HDN(B, i -1, S), then the number of nodes

n i in Q i is N i-1 /s i . The s i can be 1 or q j=1 |B ij |, where 1 ≤ i j ≤ r and q ≤ r. That is, s i can be a product of any number of integers in {b 1 , b 2 . . . . , b r }. For example, if r = 3, b 1 = 2, b 2 = 3, and b 3 = 5, the possible s i can be 1, 2, 3, 5, 2 × 3, 2 × 5, 3 × 5, or 2 × 3 × 5.
The construction of HDN(B, i, S), 1 ≤ i ≤ k, can be defined by a twostep process: First, we perform a dual-construction on the quotient graph

Q i-1 = HDN(B, i -1, S)/SN i (HDN(B, 0, S)= B).
Let the graph generated by the dual-construction be Q i , and the subgraph of two nodes that is connected by a cross-edge of level i be K 2 . Second, to get the HDN(B, i, S), we replace every That is, the number of clusters in each class is equal to the number of supernodes in a cluster. At level i, each super-node in a cluster has s i new links to a super-node in a distinct cluster of the other class. Because there are s i nodes in a super-node, one node contributes a new link. The dual-construction of an RDN is a special case of the construction of an HDN with

K 2 in Q i by a product graph K 2 × SN . We call HDN(B, i -1, S) cluster of HDN(B, i, S). HDN(B, i -1, S) 0 1 ni -1 0 1 ni -1
s i = 1 for 1 ≤ i ≤ k.
The indexes of the nodes in HDN(B, k, S) can be defined as follows. Let SN k id be a super-node id in a cluster of HDN(B, k, S) and N k id be a node id in a supernode, then a node in the HDN(B, k, S) can be represented by

(C k , U k id , SN k id , N k id )
where C k is the class id (0 or 1) and U k id is the cluster id. A cross-edge at level

k connects node (C k , U k id , SN k id , N k id ) and node (C k , SN k id , U k id , N k id ).
If we use a 2 × 3 × 5 torus as the base network, Table 1 lists the number of nodes in HDN(B, 1, S) and HDN(B, 2, S) under the different configurations of S. The node degrees are 7 and 8 for HDN(B, 1, S) and HDN(B, 2, S), respectively, because the node degree of B is 6. From the table, we can see that the HDN covers the nodes range from several hundreds to several millions. Suppose that the node degree of the base network B is d 0 , the node degree of the HDN(B, k, S) is d 0 + k. Let N i-1 be the number of nodes in the HDN(B, i -1, S). There are

k = 2 s2 = 1 s2 = 2 s2 = 3 s2 = 5 s2 = 6 s2 = 10 s2 = 15 s2 = 30 s1 = 1 6,
N i = 2(N i-1 /s i )N i-1 = 2N 2
i-1 /s i nodes in the HDN(B, i, S) for 1 ≤ i ≤ k, where N i-1 /s i is the number of clusters in each class. That is, the number of nodes in the HDN(B, k, S) is (2n 0 ) 2 k /(2 k i=1 s i ), where n 0 is the number of nodes in the base network.

Let the diameter of the HDN(B, i -1, S) be D i-1 and the diameter of the super-node (SN) be D(SN i ). Then, if we map a super-node into a single node, the diameter of the quotient graph

Q i-1 is D(Q i-1 ) = D i-1 -D(SN i ).
To route a node u in a cluster of class 0 (or 1) to a node v in a different cluster of the same class, we can route u along with a direct link of level i to a node u in a cluster of class 1 (or 0). This takes one step. Then, we route u inside the cluster to a node w that can reach a node w in the same cluster of node v along with direct link of level i. The longest distance between nodes u and w is D(Q i-1 ).

Similarly, we can route node w to a node w (by one step) and then to a node v which is in the same super-node of v (by D(Q i-1 ) steps). Finally, we route v to node v, this takes D(SN i ) steps. Therefore, we have the following recurrence:

D i = 2(1 + D(Q i-1 )) + D(SN i ) = 2D i-1 -D(SN i ) + 2
Solving the above recurrence, we get the diameter D k of HDN(B, k, S) as below:

D k = 2 k D(B) - k-1 j=0 2 j D(SN k-j ) + 2 k+1 -2
where D(B) and D(SN i ), 1 ≤ i ≤ k, are the diameters of the base network and the super-nodes, respectively. The results of the analysis in this section are summarized in the following theorem.

Theorem 1 Assume that the base network B is a symmetric, product graph and SN i , 1 ≤ i ≤ k, are sub-product graphs of B with |SN i | = s i . Let the number of nodes, the node-degree, and the diameter of B be n 0 , d 0 , and D 0 , respectively. Let the diameters of

SN i , 1 ≤ i ≤ k, be D(SN i ). Let S = {G 1 , G 2 , . . . , G k }, where G i is a sub-graph of HDN(B, k -1, S) and s i = |G i | is the number of nodes in a super-node at the level i for 1 ≤ i ≤ k. Then, the number of nodes of HDN(B, k, S) is (2n 0 ) 2 k /(2 k i=1 s i ), the node-degree is d 0 + k, and the diameter is D k = 2 k D(B) - k-1 j=0 2 j D(SN k-j ) + 2 k+1 -2,
where N is the number of nodes in HDN(B, k, S).

Table 2 lists the topological properties of the torus, n-cube, CCC [START_REF] Preparata | The cube-connected cycles: a versatile network for parallel computation[END_REF], Dual-Cube [START_REF] Li | Efficient collective communications in dual-cube[END_REF], RDN, and HDN. The CCC (cube-connected cycles) is obtained by replacing a node in an n-cube with an n-node cycle. The Dual-Cube is a special case of RDN with k = 1 and a base network of an n-cube.

(x + y + z)/2 n-cube 2 n n n CCC(n) n * 2 n 3 2n + n/2 -2 Dual-Cube(n) 2 2n-1 n 2n RDN(m, k) (2m) 2 k /2 d0 + k 2 k * D0 + 2 k+1 -2 HDN(B, k, S) (2|B|) 2 k /(2 k i=1 si ) d0 + k 2 k (D(B) -k-1 j=0 2 j (D(SN k-j )) + 2 k+1 -2
In [START_REF] Li | Recursive dual-net: A new universal network for supercomputers of the next generation[END_REF], we introduced the CR (cost ratio) for measuring the combined effects of the hardware cost (node degree) and the software efficiency (diameter) of an interconnection network. Instead of CR, this paper uses a more general measure, namely weighted cost ratio CR w (G), for the evaluation. The CR w (G) is defined as below. Let |(G)|, d(G), and D(G) be the number of nodes, the node degree, and the diameter of G, respectively. We define CR w (G) as

CR w (G) = w 1 d(G) + w 2 D(G) log 2 |(G)|
where w 1 and w 2 are weights for node degree and diameter, respectively. We have w 1 + w 2 = 100%.

The weighted cost ratio CR w of an n-cube is always 1 regardless of its size and weights. The CR w for some HDN(B, k, S) is shown in Table 3 where B is a 2×3×5 torus and we assume w 1 = w 2 = 50%. For simplicity, we use the number of nodes in super-nodes to represent S, instead of sub-graphs. From the table, we can see that the HDNs are more effective than hypercubes and tori measured by the weighted cost ratio although as the s i increases, the CR w becomes larger. The minimum CR w shown in the list is 0.69. Unfortunately, we do not know the theoretical or experimental optimal value of CR w up to the date we wrote this paper and it can be an open question for the future.

Routing on HDN

Given two nodes u and v in HDN(B, k, S), we first present a simple routing algorithm that finds a shortest path from u to v. In Section 2, we defined the product and quotient graphs. Now, we define the difference graph as follows.

Let SN 1 and SN 2 are two super-nodes in base network B, the difference graph SN 1 -SN 2 is the sub-product graph of B such that B i ,

1 ≤ i ≤ r, is in SN 1 -SN 2 if and only if B i ⊂ SN 1 and B i ⊂ SN 2 . For example, if B = C 2 × C 3 × C 5 , SN 1 = C 2 × C 3 , and SN 2 = C 3 × C 5 then SN 1 -SN 2 = C 2 .
We also need a re-indexing process of nodes in the cluster, which is an HDN(B, i -1, S), for routing via cross-edges of level i since the indexes of nodes in HDN(B, i -1, S) is based on SN i-1 and the cross-edge of level i is defined based on SN i . The index of a node in HDN(B, i -1, S) contains four parts (C i-1 , U i-1 id , SN i-1 id , N i-1 id ) as explained in the previous section. At the construction of the ith level, HDN(B, i -1, S) becomes a cluster containing only two parts, SN i id and N i id , of the node index in HDN(B, i, S). The other two parts, C i and U i id , are generated from the construction at the ith level. The re-indexing process that generates a 1-to-1 mapping between (C i-1 , U i-1 id , SN i-1 id , N i-1 id ) and (SN i id , N i id ) on an HDN(B, i -1, S) is necessary for the proposed routing algorithm.

Since the number of super-nodes

SN i in HDN(B, i -1, S) equals to N i-1 /s i , the range of SN i id is 2|U i-1 /(SN i -SN i-1 )| × |(SN i-1 -SN i )|. If s i-1
= s i then the re-indexing is simple: 1-1 mapping between SN i id and the 3-tuple (C i-1 id , U i-1 id , SN i-1 id ). However, when s i-1 = s i , the re-indexing is a little complicated and is explained below.

Let the q-tuple, (b i1 , . . . , b iq ) be the index of a node in a super-node SN , where b i1 ×. . .×b iq = |SN |. Then the re-indexing from id , SN 1 id , N 1 id ), where C 1 = 0 or 1, 0 ≤ U 1 id < 5, 0 ≤ SN 1 id < 5, and 0 ≤ N 1 id < 6. For the indexes of the nodes in HDN(B, 2, S), we perform re-indexing of nodes in HDN(B, 1, S), which is a cluster of HDN(B, 2, S), to get (SN 2 id , N 2 id ), where 0 ≤ SN 2 id < 2 × 5 × 2 = 20, and 0 ≤ N 2 id < 3 × 5 = 15, obtained by swapping

(C i-1 , U i-1 id , SN i-1 id , N i-1 id ) to (SN i id , N i id ) moves the indexes of those B j ⊂ SN i -SN i-1 into N i id and the indexes of those B j ⊂ SN i-1 -SN i into SN i id . For example, let B = C 2 × C 3 × C 5 , s 1 = |C 2 | × |C 3 | = 6,
|B 1 | and |B 3 |. That is, |SN 2 | = |C 1 | × |U 1 | × |B 1 | = 2 × 5 × 2 = 20, and |N 2 | = |B 2 | × |B 3 | = 15.
Table 4 shows four examples of re-indexing in detail for a cluster in the HDN(B, 2, S) with B = C 2 ×C 3 ×C 5 , s 1 = 2×3 = 6, and s 2 = 3×5 = 15. In the HDN(B, 1, S), the node representation Assume that the point-to-point routing algorithm in the base network is available. The proposed algorithm for routing node u to node v in HDN(B, k, S) works as follows. We first perform re-indexing of u and v if k > 1. Then, there are three cases: the two nodes are in the same cluster (Case 1), in the distinct clusters of the same class (Case 2), and in the distinct clusters of distinct classes (Case 3). Case 1 is trivial. Case 3 can be reduced to Case 2 by routing u via a cross-edge of level k. Therefore, we explain only the Case 2: The two nodes are in the distinct clusters with the same class. We first identify the super-nodes, denoted as SN k u and SN k v , in the two Q k-1 s containing u and v, respectively, such that SN k u and SN k v are connected by a unique cross-edge of level k in Q k from the dual-construction. Then, we route node u to node u , and node v to node v inside the clusters of level k, respectively. Notice that, u and v are not unique although SN k u and SN k v are unique. The algorithm finds the u and v that leave u k 3 and v k 3 unchanged if possible. And then, the routing from u to v is done by routing u to u ∈ SN k v via a cross-edge of level k in HDN(B, k, S) and routing from u to v inside SN k v . The algorithm is formally presented as Algorithm 1. The correctness of the algorithm and its time complexity are given in Theorem 2.

(C 1 , U 1 id , SN 1 id , N 1 id ) can be converted to a serial number i by i = C 1 ×(|B|/s 1 ) 2 ×s 1 +U 1 id ×(|B|/s 1 ) 1 ×s 1 +SN 1 id ×s 1 +N 1 id = C 1 × 150 + U 1 id × 30 + SN 1 id × 6 + N 1 id . Similarly, the (SN 2 id , N 2 id ) can be converted to a number SN 2 id × s 2 + N 2 id = SN 2 id × 15 + N 2 id .
(C 1 , U 1 id , SN 1 id , N 1 id ) Serial number (SN 2 id , N 2 id ) Serial number (0, 0, 0, 0) 0 × 150 + 0 × 30 + 0 × 6 + 0 = 0 (0, 0) 0 × 15 + 0 = 0 (1,
Theorem 2 Assume that the routing algorithms in the base network B is available. In HDN(B, k, S) for k > 0, routing between any two nodes can be done in at most

2 k R(B) - k-1 j=0 2 j R(SN k-j ) + 2 k+1 -2 steps, where R(B) and R(SN i ), 1 ≤ i ≤ k,
are the time complexities of the routing in B and SN i , respectively.

Proof: We show the correctness of Algorithm 1 by induction on k. Assume that the algorithm is correct for k -1 ≥ 0. From the algorithm, it is clear that we need to consider only Case 2. In Case 2, nodes u and u are in the same cluster by the definition of u . They can be connected by the induction hypothesis. Similarly, nodes v and v can be connected. The node u that is connected to u by a cross-edge of level k and node v are in the same super-node as can be seen from their IDs. Therefore, they can be connected by Base routing algorithm. Next, we derive the time complexity R k of the algorithm. In Case 2, there are two recursive calls to connect u to u and v to v , respectively. Since the nodeIDs of u and u are the same (so are v and v ), a recursive call takes only R k-1 -R(SN k ) time. Since the SupernodeIDs of u and v are the same, the last call to Base route to connect u to v takes only R(SN k ) time. In Case 3, there is an additional routing step via a cross-edge. Therefore, the time complexity R k of HDN Routing(HDN(B, k, S), u, v) satisfies the recurrence

R k = 2(R k-1 -R(SN k ))+R(SN k )+2 for k > 0. Solving this recurrence, we have R k = 2 R(B)- k-1 j=0 2 j R(SN k-j )+2 k+1 -2 where R(B) and R(SN i ), 1 ≤ i ≤ k,
are the time complexities of the routing in B and SN i , respectively.

Total Exchange Routing on HDN

Design of efficient routing algorithms for communications is the key issue in parallel computers or networks. Collective communications are required in load balancing, event synchronization, and data exchange. Based on the number of sending and receiving processors, these communications can be classified into one-to-many, one-to-all, many-to-many and all-to-all. The nature of the messages to be sent can be classified as personalized or non-personalized (multicast or broadcast). The all-to-all personalized communication (total exchange) is at the heart of numerical applications.

An important metric used to evaluate efficiency of communication is transmission latency, or communication time. The communication time depends on many factors such as contentions, switching techniques, network topologies etc. Therefore, we first define the communication model used in this paper.

We assume that the communication links are bidirectional, that is, two directly-connected processors can send messages to each other simultaneously. We also assume the processor-bounded model (one-port model) in which each node can access the network through a single input port and a single output port at a time. The port model of a network system refers to the number of internal channels at each node. In order to reduce the complexity of communication hardware, many systems support one-port communication architecture. We also assume the linear cost model in which the transfer time for a message is linearly proportional to the length of the message.

There are many switching methods. In this paper, we assume the packet switching model [START_REF] Kumar | Introduction to parallel computing: design and analysis of algorithms[END_REF]. In this model, each packet is maintained as an entity that is passed from node to node as it moves through the network. The long message can be partitioned and transmitted as fixed-length word w. The first few bytes of a packet contains routing and control information and are referred as packet header. A packet is completely buffered at each intermediate node before it is forwarded to the next node (for this reason, the model is also called store-andforward switching). In this paper, we allow packages that are headed for the same destination to be combined into a single message. The time to pack and unpack messages is included in the startup latency. The packet switching model is suitable for collective communication in MPP since it is safer than other switching models such as virtual cut-through switching. With packet switching model, the communication time for a message of length m (number of fixedlength words) to be sent to a node of distance d is d(t s + mt w ), where t s is startup latency, the time required for the system to handle the message at the sending node, t w is the per-word transfer time (1/t w is the bandwidth of the communication links). Through this paper, we will use the formula above for estimating the communication times of the proposed algorithms.

In total exchange, each node sends a distinct message of size m to every other node. The total number of messages is p 2 (a node also has a message for itself). Referring to Figures 2 and3, the algorithm for total exchange in HDN can be described in four stages. Note that all nodes in the figures do the same operations but only the cases of node 0 are shown for clarity. Also note that the commas of node address in Figure 2 are omitted for saving spaces. 1. In the first stage, we first divide M my id into two parts, M 1 my id and M 2 my id ,

where M 1 my id contains all messages to be sent to the nodes in the clusters of my type, and M 2 my id contains the rest message. Then, the first part of personalized messages M 1 is exchanged between my id and partner, the neighbor via cross-edge of level k. The time this stage takes is t s + t w mn k /2. 2. In the second stage, we first pack all messages that are to be sent to the nodes in the cluster of level k with clusterID = q into a single message msg q . Then, we perform total exchange inside each cluster, where msg q is to be sent to node with nodeID = q. The time this stage takes is denoted as T k-1 (n k /n k-1 ), the time for total exchange inside the cluster. 3. In the third stage, each node packs the received messages into a single message of length n k m and sends the packed message to its neighbor along the cross-edge of level k. After receive the message, each node unpacks the message received from its neighbor into n k-1 messages, msg q , where msg q is the collection of all messages destinated to node with nodeID = q in the cluster. The time this stage takes is also t s + t w mn k /2. 4. In the last stage, we perform total exchange again within each cluster of level k. This can be done since the packed messages sent through the level-k cross-edge are all destinated to the nodes inside the cluster. The time this stage takes is also T k-1 (n k /n k-1 ), which is the time for total exchange inside the cluster. The algorithm is showed in Algorithm 2. All nodes execute the algorithm concurrently. In Algorithm 4, my id is the id of the node. The initial message to be sent is M my id which contains p messages of length m. At the end of the algorithm, each node stores the collection of all p messages in result. Pack all messages in M my id that are to be sent to the nodes in the cluster of level k with cluster ID = q into a single message msgq to be sent to the node with node ID = q; T my id ← TOTAL EXCHANGE (HDN(B, k -1, S), M my id ); send message T my id to partner; receive message Tpartner from partner; T my id ← Tpartner; Unpack T my id into n k-1 messages, msg q , such that msg q is the collection of messages destinated to the node with node ID = q ; result ← TOTAL EXCHANGE (HDN(B, k -1, S), T my id ); endif end

The time to complete the total exchange on an HDN(B, 1, S) is T 1 (m) = (t s + mt w n 1 /2) + T 0 (mn 1 /2) + (t s + mt w n 1 /2) + T 0 (mn 1 /2) = 2(t s + mt w n 1 /2) + 2T 0 (mn 1 /2).

Generally, on an HDN(B, k, S), the time to complete the total exchange is T k (m) = (t s +t w mn k /2) +T k-1 (n k /n k-1 )+(t s +t w mn k /2)+T k-1 (n k /n k-1 ) = 2(t s + t w mn k /2) + 2T k-1 (n k /n k-1 ). That is, T k (m) = (2 k+1 -2)(t s + t w mn k /2) + 2 k T 0 (mn k /2)

where n k is the total number of nodes and T 0 (m) is the time complexity for total exchange in B.

In the examples of Figures 2 and3 where B is a 2-cube, T 0 (m) = 2(t s + t w m). If B is an n-cube, then T 0 (m) = n(t s + t w m2 n /2).

For the HDN(B, 1, S) shown in Figures 2, T 1 = 2(t s +8t w m)+4(t s +8t w m) = 6(t s + 8t w m). In contrast, T = 4(t s + 8t w m) for a 4-cube of same size. For the HDN(B, 2, S) shown in Figures 3, T 2 = 6(t s + 64t w m) + 4 × 2(t s + 64t w m) = 14(t s + 64t w m). In contrast, T = 7(t s + 64t w m) for a 7-cube of same size. The times of total exchange for HDNs are longer than that for hypercubes but an HDN has much less links than a hypercube of the same size. We summarize this result in the following theorem. Theorem 3 Assume that the time complexity T 0 (m) for total exchange in the base network B is known, where m is the length of each message. The time complexity T k (m) for total exchange on an HDN(B, k, S), k > 0, is T k (m) = (2 k+1 -2)(t s + t w mn k /2) + 2 k T 0 (mn k /2), where n k = (2n 0 ) 2 k /(2 k i=1 s i ).

Concluding Remarks

The hierarchical dual-net can connect a large number of nodes with a small nodedegree and a short diameter. It is a potential candidate for the interconnection network of the supercomputers of the next generation that may have more than one million of nodes. We can select a popular network of small size that is a product graph as the base network and then connect multiple base modules with cross links (cables) to construct a very large-scale hierarchical dual-net. We can also select a suitable set of integers based on the base network to control the number of nodes in the supercomputer. The base networks can be implemented in a NoC VLSI and high-speed line cables may be used as the cross links to connect PCB modules in cabinets. We presented an efficient algorithm for total exchange on recursive dual-net. There are many problems, such as disjoint path and fault-tolerant routing, on recursive dual-net that are worth further research.

Fig. 1 .

 1 Fig. 1. Build an HDN(B, i, S) from HDN(B, i -1, S)Referring to Figure1, an HDN(B, i, S) consists of 2n i clusters which are divided into two classes: class 0 and class 1 with each class containing n i clusters. That is, the number of clusters in each class is equal to the number of supernodes in a cluster. At level i, each super-node in a cluster has s i new links to a super-node in a distinct cluster of the other class. Because there are s i nodes in a super-node, one node contributes a new link. The dual-construction of an RDN is a special case of the construction of an HDN with s i = 1 for 1 ≤ i ≤ k.The indexes of the nodes in HDN(B, k, S) can be defined as follows. Let SN k

Table 1 .

 1 Number of nodes in HDN(B, k, S) where B is a 2 × 3 × 5 torus k = 1 s1 = 1 s1 = 2 s1 = 3 s1 = 5 s1 = 6 s1 = 10 s1 = 15 s1 =

and s 2 =

 2 |C 3 | × |C 5 | = 15, then, the nodes in HDN(B, 1, S) can be represented by (C 1 , U 1

Fig. 2 .Fig. 3 .

 23 Fig. 2. Four stages of total exchange on HDN(B, 1, S) (only shows the case of node 0)

Table 2 .

 2 Comparison of topological properties

	Network	# of nodes	Degree	Diameter
	3D Torus	x * y * z	6

Table 3 .

 3 CRw with w1 = w2 = 50% for some HDN(B, k, S)

	Network	n	d	D	CR
	10-cube	1,024	10	10	1.00
	3D-Tori(10)	1,000	6	15	1.05
	HDN(B, 1, (1))	1,800	7	10	0.79
	HDN(B, 1, (2))	900	7	9	0.82
	HDN(B, 1, (3))	600	7	9	0.87
	19-cube	524,288	19	19	1.00
	3D-Tori(80)	512,000	6	120	3.32
	HDN(B, 2, (2, 2))	810,000	8	19	0.69
	HDN(B, 2, (2, 5))	324,000	8	18	0.71
	HDN(B, 2, (5, 2))	129,600	8	17	0.74

Table 4 .

 4 Re-indexing examples

	Index in HDN(B, 1, S)	Index in HDN(B, 2, S)

  Algorithm 2: TOTAL EXCHANGE (HDN(B, k, S), M my id ) begin if k = 0 result ← TOTAL EXCHANGE (B, M my id ); else Divide M my id into two parts, M 1 my id and M 2 my id , where M 1 my id contains all messages to be sent to the nodes in the clusters of type my type, and M 2 my id contains the rest messages; partner ← the neighbor via cross-edge of level k; send message M 1 my id to partner; Receive message M 1partner from partner; M my id ← M 2 my id ∪ M 1partner;

(0000,0000) (0000,0001) (0000,0010) (0000,0011) (0000,0100) (0000,0101) (0000,0110) (0000,0111) (0000,1000) (0000,1001) (0000,1010) (0000,1011) (0000,1100) (0000,1101) (0000,1110) (0000,1111) (0000,0010) (0000,0011) (0000,0000) (0000,0001) (0000,1100) (0000,1101) (0000,1110) (0000,1111) (0000,0100) (0000,0101) (0000,0110) (0000,0111) (0000,1000) (0000,1001) (0000,1010) (0000,1011) (0000,0010) (0000,0011) (0000,0000) (0000,0001) (0000,1100) (0000,1101) (0000,1110) (0000,1111) (0000,0100) (0000,0101) (0000,0110) (0000,0111) (0000,1000) (0000,1001) (0000,1010) (0000,1011) (0000,0010) (0000,0011) (0000,0000) (0000,0001) (0000,1100) (0000,1101) (0000,1110) (0000,1111) (0000,0100) (0000,0101) (0000,0110) (0000,0111) (0000,1000) (0000,1001) (0000,1010) (0000,1011)

Algorithm 1: HDN ROUTING (HDN(B, k, S), u, v) input: HDN(B, k, S); input: node u = (u k 0 , u k 1 , u k 2 , u k 3 ) (the node representation of level k);