Gangfeng Liu

Yunlan Wang
email: wangyl@nwpu.edu.cn

Tianhai Zhao

Jianhua Gu

Dongyang Li

mHLogGP: a Parallel Computation Model for CPU/GPU Heterogeneous Computing Cluster

Keywords: CPU/GPU Heterogeneous computing, Parallel Computation Model, Parallel Computing

CPU/GPU heterogeneous computing has become a tendency in scientific and engineering computing. The conventional computation models cannot be used to estimate the application running time under the CPU/GPU heterogeneous computing environment. In this paper, a new model named mHLogGP is presented on the basis of mPlogP, LogGP and LogP. In mHLogGP, the communication and memory access is abstracted based on the characteristic of CPU/GPU hybrid computing cluster. This model can be used to study the behavior of application, estimate the execution time and guide the optimization of parallel programs. The results show that the predicted running time approaches to the actual execution of program.

Introduction

The compute-capable graphics processing units (GPUs) are attractive as accelerators for high-performance parallel computing. According to the top500 supercomputer sites, three out of the top five fastest supercomputers in November 2011, employ GPUs. CPU and GPU have their own advantages respectively. CPU is suitable for complex arithmetical, logical, and input/output operations. GPU consists of hundreds of smaller cores which gives high compute performance. Therefore, heterogeneous CPU/GPU computing provides new and powerful computing architecture. The conventional parallel computation models are not suitable for the new emerging heterogeneous computing environment. The PRAM [START_REF] Steven | Parallelism in Random Access Machines[END_REF] is unrealistic because it assumes that all processors work synchronously. The bulk-synchronous parallel model (BSP) [START_REF] Leslie | A bridging model for parallel computation[END_REF] models asynchronously processors, communication latency and limited bandwidth. The LogP model [START_REF] David | Logp: Towards a realistic model of parallel computation[END_REF] is capable of predicting the communication time for short messages but cannot accurately predict the communication time for long messages. The LogGP [START_REF] Albert | LogGP: incorporating long messages into the logp modelá one step closer towards a realistic model for parallel computation[END_REF] extends it and adds an additional parameter for long messages, which captures the bandwidth constraints for long messages. The PlogP [START_REF] Thilo | Fast measurement of logp parameters for message passing platforms[END_REF] considers the varying communication overheads for different message sizes. The memory LogP [START_REF] Qasim | Modeling advanced collective communication algorithms on cell-based systems[END_REF], log 3 P and log n P [START_REF] Kirk | Xian-HeSun: log n P and log 3 P: accurate analytical models of pointto-point communication. In: distributed systems[END_REF] models draw characteristics of data movement in the tradi-tional parallel computer. However, all the models mentioned above ignore the new characteristics of communication and memory access in the heterogeneous computing system. MMGP [START_REF] Filip | Modeling multi-grain parallelism on heterogeneous multicore processors: A case study of the cell be[END_REF] is a multi-dimensional parallel computation model for steering the parallelization process on heterogeneous multi-core processors. However, it ignores the influence of data movement across computational cores. mPlogP [START_REF] Liang | mPlogP: a Parallel Computation Model for Heterogeneous Multi-core Computer[END_REF] is also a model for heterogeneous multi-core computer. It extends PLogP and adds another level on top of the conventional level to model behaviors of computational cores. However, the characteristics of CPU/GPU hybrid computing environment is not considered in this model.

In this paper, a new model named mHLogGP is proposed, which is on the basis of mPlogP,LogGP and LogP. In order to analyze the characteristics of CPU/GPU heterogeneous computing cluster, some related parameters are defined to estimate the time performance of the parallel program. The hardware abstraction of computer node which includes GPU and CPU is showed in Fig. 2(a). The host memory is shared by GPU and CPU. GPU exchanges information with CPU by accessing host memory.

The hardware abstraction of the widely used CPU/GPU computer cluster is shown in Fig. 2(b). A computer cluster consists of dozens to thousands computers that are connected by high performance network such as Infiniband, 10-Gigabit Ethernet, etc.

2.2

The parameters of the model o(Overhead): the time needed by processor to send or receive a message, o ∈{o s ,o r }. o s ∈{o sg ,o sc }, o sg is the time that GPU sends messages to host memory. o sc is the time that CPU sending messages. o sc ∈{o scd ,o scn }, o scd defined as the time that CPU sends messages to global memory. o scn is the time that CPU sending messages to another computer node. o r ∈{ o rgb , o rn , o rd }, o r is the time that a processor receive a message. o rgb is the time of initializing the GPU. o rn (o rd ) is defined as the time that CPU (GPU) engaged reception of each message.

l : latency, or delay, incurred in communicating a message from its sender to receiver. l ∈ { l g , l c , l dm , l nn }, l g denotes the latency between the parallel data cache and the global memory. l c denotes the latency from the core of CPU to the host memory of another node. l dm denotes the latency between the global memory and the host memory. l nn is defined as the delay from one node to another. G: the Gap per byte for long messages, which is defined as the time per byte for a long message. The reciprocal of G characterizes the available communication bandwidth per processor for long messages.

g: the gap, defined as the time interval between consecutive message transmission or reception. The reciprocal of g corresponds to the available communication bandwidth per-processor and the speed of accessing the memory. g ∈{ g g , g c ,}, g g means the time interval between consecutive message transmission or reception in the global memory of GPU. g c denotes the time interval for sending a message from the host memory to network cache.

P: the number of computer nodes. P is mostly used to model the collective communication.

3

Usage of the Model

Point-to-Point Communication

In this section, the performance of point-to-point communications is predicted. In the CPU/GPU hybrid computing environment, the cost of sending short message between two nodes can be predicted as o s1 + l +o r1 ( LogP), where l is the latency of communicating a message from its sender to receiver, o s1 , o r1 is the time of sending or receiving the first message, where G is the Gap per byte. For long message, the Point-to-Point Communication time can be predicted as follows:

1 1 1 ( 1) sr T o k G l o      (LogGP) (1)
As shown in Fig. 3, if the non-blocking communication is used to send 'm' consecutive messages, the transmission time is as follows:

1 1 1 ( 1) ( 1) 1 ( 1) max((( 1)
), ) o

m s i c i s i nni nnm rm i T o k G k G g o l l                (2)
Where 1 i k  is the size of the (i+1)-th message. o si , o ri is the time of sending or receiving the i-th message, It is obvious that network bandwidth and the speed of accessing the memory have a great influence on inter-node communication time.

Broadcast Communication

We take the MPICH2 as example. For the short message or the number of nodes less than 8, the broadcast communication is taken by the short message algorithm, namely, binomial tree. If the number of nodes is P and the size of the message is k. it needs lgP steps in broadcast communication. Each byte travels through the network for l nn . Under the CPU/GPU heterogeneous computing environment, the broadcast communication cost can be predicted as follows: T broP = lgP *o+ lgP *k*(G+l nn ).

(3) For medium-size messages and power-of-two number of nodes, the broadcast communication uses a recursive doubling algorithm. This costs approximately T broP = 2*lgP *o+ 2*k*(P-1)/P*(G+l nn ). For long messages or medium-size messages and non-power-of-two node situations the broadcast communication use a ring algorithm, which takes P-1 steps, so the total cost is T broP = (lgP +P-1)*o+2*k*(P-1)/P*(G+l nn ).

Validating the Model

Our experimental platform is the Dawning TC3600 Blade system with Intel X5650 processors and NVIDIA Tesla C2050 GPUs at the National Super Computer Center in Shenzhen, China. The computing network is 40Gb InfiniBand.

Table 1. Sender and receiver overheads between Global memory and host memory

It can be shown in Table 1 that o scd remains unchanged when the message size is less than 256 KB. It increases linearly when the message size is larger than 256 KB. So does o rg , o scn and o rn when the size is larger than 64KB. Since the GPU must be initialized, no matter the size of message short or large, the value of o rgb has no large variation. Because the speed of intra-node communication is very high, the value of l g , l c and l dm can be ignored. l nn denotes the intra-node delay. The experiments show that l nn is random from 100 s  to 700 s  when the size of a message is less than 1MB. It increases linearly when the message size is larger than the 1MB.

As shown in Table 2, G is not stable when the size of a message is less than 1MB. When the message size larger than 1MB, it approaches a constant. g g and g c has changed little when the size of the data varies.

As shown in Table 3, for the messages size larger than 64KB, the total cost of broadcast communication is accord with ring algorithm. There may be errors of the

Experiments Validation

The matrix multiply [START_REF] Karunadasa | Accelerating High Performance Applications with CUDA and MPI[END_REF] is adopted to validate mHLogGP. The implementation of algorithm mainly adopts the point-to-point communication and broadcast communication. Let A and B be two matrices, n represents the number of lines in matrices A. Matrices A is divided into p parts. The master process sends every part and broadcast matrices B to p slave processes. The slave processes receive the data and do the calculation by GPU and send the results to the master process. The master process receives the results and save the information. The overhead of program can be estimated based on the formulas in Section 3. The cost of sending messages (for the master process) is:

T ptop1 = 1 1 1 1 ( 1) ( 1) 1
( 1) max((( 1) ), ) o

p scn i c i scn i nni nnp rnp rgbp rdp i o k G k G g o l l o o                
The cost of receiving messages (for the master process) is: The experiments are done on a cluster with 8 nodes including GPU. The comparison of the predicted and actual running time is showed in Fig. 5. It is shown that our model could get competitive accuracy. With the increase of computation scale, the errors of the predicted running time increases. But, at worst, when matrix size becomes 8192*8192 on 8 nodes, the error of our model is 15-25%.

T ptop2 = 1 ( ( 1) o

Conclusion And Future Work

Nowadays, CPU/GPU heterogeneous computing is a tendency. The conventional computation models do not fit in with the emerging CPU/GPU heterogeneous computing environment. In this paper, a new model named mHLogGP is presented. In mHLogGP, communication and memory access is abstracted by considering the characteristic of the CPU/GPU heterogeneous computing environment. This parallel computation model can be used to estimate the communication time and running time of parallel program, to find the bottleneck and guide the optimization of parallel programs in heterogeneous computing environment. The evaluation experiments show that this computation model is valid. In the future, the computing time of parallel program will be concerned. We will also study the execution time of concurrent parallel programs under the heterogeneous computing environment with dynamic load.

Fig. 1 (Fig. 2 .

 12 Fig.1(a) shows GPU abstract architecture[START_REF] Sunpyo | An Analytical Model for a GPU Architecture with Memorylevel and Thread-level Parallelism Awareness[END_REF]. Each thread processors has private memory and parallel data cache. The thread processors visit each other through the parallel data cache and global memory. Fig.1(b) shows CPU abstract architecture, the CPU cores exchange message via the bus and the host memory.

Fig. 3 .Fig. 4 .

 34 Fig. 3. Sending and receiving messages

Fig. 5 .

 5 Fig. 5. The comparison of the predicted and actual execution time

Table 2 .

 2 Overheads per byte for a long message

	Message Size	1B	4KB	16KB	64KB	256KB	1MB	4MB	16MB	256MB
	oscd,ord( s  )	3	45	52	64	173	381	1243	4754	71799
	Org( s  )	20	22	26	50	174	824	1329	4803	70209
	Oscn( s  )	4	4	8	17	32	65	128	904	20426
	Orgb( s  )	102055	103537	109984	109455	114427	111911	101656	101410	103542
	orn( s  )	3	5	8	15	32	67	129	1002	21527
	lnn( s  )	141	149	225	268	339	755	2589	10461	161469

Table 3 .

 3 Overheads of broadcast communication

Table 4 .

 4 The comparison between computing cost and total cost

	Matrix size	64*64	128*128	256*256	512*512	1024*1024	2048*2048	4192*4192	8192*8192
	Total cost( s  )	110268	109819	115474	115553	119466	219825	964722	6999979
	Computing	164	166	241	300	423	1042	2854	11284
	cost( s  )

  ) The value of parameters is based on the experiments data in section 4. For single GPU, the computing cost and total cost is shown in Table4. The computing time in matrix multiply is very short compared to communication time, so it is ignored.

	The cost of the broadcast Communication (for the master process) is:
				T broP = (lg	P P  	1)* 2* *( 1) / *( l ) nn o k P P G   
	Message Size	1B	256B	1KB		16 kB		64kB		256kB	1MB	16MB	64MB
	G(us/B)	4.5	x 10 6	141	0.578125	0.145508	0.014710		0.004577 3 x 10 6	0.001537	0.000845	0.000737	0.000749
	gg( s  )	4		4 mHLogGP Acutal	4		4			4		4		7 mHLogGP Acutal	68	104	249
	gc( s  ) size 3.5 time(us) 1 1.5 2 2.5 3 microsecond (μs)		0 1B	0	256B	0 16KB		42 64KB	2.5 256kB 39 1 1.5 microsecond 2 (μs)	40 1MB	4MB	41	16MB	55 64MB	35 256MB
	nodes	0.5										0.5
	2 4	0		83667 132472 64*64 128*128 256*256 512*512 1024*1024 2048*2048 4192*4192 8192*8192 85284 97732 106578 142002 140794 188769 The matrix size	106578 213175 0	95334 197967 64*64 128*128 256*256 512*512 1024*1024 2048*2048 4192*4192 8192*8192 101449 140239 211425 234129 284768 392290 The matrix size	417945 805469
	8			207340		203330	291777	294295	299879	300807	472381	599026	822744	1804826
				p
				i  	sgi o o scni 	i k     nni G l	rni

Acknowledgment

This research was partially supported by the National High-Tech Research and Development Plan of China (Grant No. 2009AA01Z142).


Li-Jia Chiu

Norman Huang
email: nthuang@ntnu.edu.tw

The Difference Analysis Between Demographic Variables and Personal Attributes -The Case of Internal Auditors in Taiwan

Keywords: Internal auditors, personal attributes, demographic variables

Internal control plays a pivotal role in enhancing security operations and quality control in an organization. From the perspective of management, internal auditors work for ensuring effective operations and improving performance for all the management systems within the organization. They bear great responsibility and play a key role. This study aimed to explore the internal audit staff demographic variables of personality traits differences. This research selected samples of the internal auditors from domestic companies, i.e. members of the Institute of Internal Auditors, in Taipei,

Introduction

Background

The major business scandals occurring in the recent years, such as the Eron case in 2001, which filed for bankruptcy, and the WorldCom case in 2002, whose CEO, Bernard Ebbers, was forced to resign due to his loans to the company in personal name, have drawn worldwide attention to corporate governance, and revealed the importance of the roles and functions of internal auditors. The Five Factors Model (FFM) of the personality traits is the most influential model in personality theory. The perspective of characteristics provides the general principle of the thinking on human beings' behavior. After rigorous scientific testing, the model has been approved by the circle of personality psychology. FFM, also known as Big Five, is one kind of characteristic perspective. It is NEO-PI-R personality inventory developed by Costa & McCrae (1992a). According to past literature, job satisfaction is relatively influenced by the employees' personality traits. Thus, this study treats internal auditors of firms as the subjects, and analyzes the personality traits of internal auditors' different statistical variables. Based on the research results, the suggestions are proposed for the reference of the business circle. [START_REF] Guilford | Personality[END_REF] suggested that personality is the combination of personal characteristics, attributes and properties. Thus, individuals have persistent and unique characteristics different from others. According to [START_REF] Smith | Choosing the right pond: the impact of group membership on self-esteem and group-oriented behaviors[END_REF], the society includes the groups defined upon population and members are decided by personal bodies or social attributes. Demographic statistics reveal quality and characteristics which include personality traits and common psychological capabilities. It is the individual difference of psychological process and behavior [START_REF] Brand | Personality and general intelligence[END_REF].

Personality is the unique composition of personal characteristics, and it determines the interaction between human beings and environment. Personality also interacts with the situations. [START_REF] Allport | Personality: A Psychological Interpretation[END_REF][START_REF] Allport | Pattern and Growth in Personality[END_REF] suggested that personality is the dynamic organization in the individuals' psychological system and the unique form which determines a person's "adaption to external environment" and "thinking & behavior". Therefore, the individuals' behavior reflects the unique personality characteristics. When the characteristics continuously appear in different situations, they are called personality traits. However, in the studies on personality traits, FFM is more commonly used and accepted by the scholars. In recent years, the psychologists have generalized a more proper theoretical framework to describe the employee difference in personnel matter arrangement. [START_REF] Tupes | Recurrent Personality Factors Based on Trait Ratings[END_REF], [START_REF] Norman | Toward and Adequate Taxonomy of Personality Attributes: Replicated Factor Structure[END_REF] and other studies suggested that FFM is legitimate to be the classification criterion of personality construction. According to the literatures on personality, FFM can be the principle framework for the researchers. [START_REF] Mount | The Big Five Personality Dimensions: Implications for Research and Practice in Human Resources Management[END_REF] suggested that FFM can be applied to personality measurement, particularly the selection of employees. [START_REF] Goldberg | An Alternative Description of Personality: The Big Five Factor Structure[END_REF] demonstrated the FFM framework of Norman. FFM of [START_REF] Costa | Four Ways Five Factors Are Basic[END_REF] is identified the most: agreeableness, conscientiousness, extraversion, neuroticism and experience openness.

Research Purposes

Based on research background and motivations above, this study intends to find the relationship between internal auditors' demographic variables and personality traits.

The findings can serve as the references for the recruitment of internal auditors. The research purpose is to find if internal auditors' demographic variables have significant differences on personality traits. According to above research purpose and literature review, this study assumes that internal auditors' personality traits will be significantly different according to different demographic variables.

Methods

Sample

This study conducted a questionnaire survey on the internal auditors in the branches of four regions in the Institute of Internal Auditors, R.O.C., and distributed the questionnaires on site. A total of 567 questionnaires were distributed to the internal auditors who participated in the project lecture on the day of survey, and 376 copies were retrieved, with a return rate of 66%. After eliminating 104 invalid samples, there were 272 valid samples, with a valid return rate of 48%.

Instrument

This study measured personality traits by NEO Five-Factor Inventory (NEO-FFI) developed by Costa & McCrace (1992b), and adopted their Form S. After the revision, there were 35 items in this inventory. Cronbach's alpha of the dimensions is extraversion 0.7181, agreeableness 0.6268, neuroticism 0.7926 and conscientiousness 0.7343.

Results

Dimensions of demographic variables

The sample structure is based on gender, age, education level, subordination level, job title and working years. Among the valid samples, 33.1% are males and 66.5 % are females; most of them are 26~35 years old (45.6%), followed by 36~45 years old (39.0%), 46~55 years old (12.1%), above 55 years old (2.6%) and below 25 years old (0.4%). There are few subjects aged below 25 years old and above 55 years old, and thus, they are respectively allocated in "26~35" and "46~55" and are revised as "below 35" and "above 46". Regarding education level, most of the subjects are graduated from universities (65.8%), followed by college (17.3%), above master (14.7%) and senior high school (vocational school) (1.8%). Generally speaking, 80.5% subjects are at least graduated from universities. Since few samples are graduated below senior high schools (vocational schools), they are allocated in "college" and it is revised as below (including) college. As to subordination level of internal auditors, according to sample distribution, general manager (office) is the most (50.4%), followed by president (office) (27.6%) and the board of directors (12.9%). Some do not belong to above (8.1%) (such as vice general manager, financial unit, etc.). In terms of job title, basic auditors (34.2%) are the most, followed by experienced auditors (31.3%) and general auditors (19.5%). Since job titles for internal auditors in different firms are varied, "others" are 14.7%; in addition, as to total working years as internal auditors, "below 3 years" is the most (38.2%), followed by "3~5 years" (29.0%), "6~10 years" (23.9%), "11~15 years (4.8%) and above 16 years (2.2%). Since there are few subjects above 16 years, they are allocated in "11~15 years" and revised as "above 11 years".

Difference analysis between gender and personality traits

Analysis is conducted to find if gender has a significant difference on four dimensions of personality traits. Independent sample t test is conducted on males and females to find if male and female internal auditors have a significant difference on the dimensions of personality traits. According to the findings, gender does not have a significant difference on dimensions of personality traits, indicating that internal auditors are not influenced by gender as a variable. Most of the valid samples are females, and this suggests that the firms prefer hiring the females. The analytical results do not indicate the difference, as shown in Table 1.

Difference analysis between age, education level, subordination level, job title, working years and personality traits

Analysis is conducted to find if age, education level, subordination level, job title and working years have significant differences on four dimensions of personality traits. One-way ANOVA is conducted on variables with at least two categories to find if there is significant difference among the groups of variables. Scheffe post hoc comparison is conducted on the groups with significant differences to compare the difference of pair groups. According to Table 2, internal auditors are not influenced by education level in different variables. Among the subjects, most of them are graduated from universities (65.8%), indicating that education level is a necessary condition for internal auditing; however, it is not the key factor of the variables. Internal auditors of different subordination levels do not show significant differences in the dimensions of personality traits, suggesting that internal auditors are not influenced by different subordination levels in the variables.

Job title has a significant difference on "conscientiousness" in personality traits. Scheffe post hoc comparison finds that "others" in conscientiousness is higher than entry-level auditors. However, "others" cannot be classified; thus, the analytical result cannot be indicated. In addition, according to one-way ANOVA, F value of job title on "extraversion" in personality traits is significant. However, Scheffe post hoc comparison does not show that the mean of any pair groups is significantly different. Although extraversion is different on the cognition of subordination level, it is uncertain if the supervisors or non-supervisors show a higher degree. Working years has a significant difference on "conscientiousness" in personality traits. Scheffe post hoc comparison finds that the degree of total working years for 6~10 years is higher than total working years below 3 years. In other words, when internal auditors' working years are longer, their conscientiousness is stronger. In addition, one-way ANOVA finds that F value of working years on "extraversion" in personality traits is significant; however, Scheffe post hoc comparison does not show a significant difference of mean in any two groups.

Mean and standard deviation of dimensions of personality traits are shown in Table 2. Regarding personality traits, average score of each item is 3.4369 and it is medium. Among the dimensions, the score of conscientiousness is the highest (mean =3.9164). It shows that internal auditors' characteristics of competency, responsibility, selfrequirement, carefulness, persistence and goal orientation are more significant. Score of neuroticism is the lowest (mean =2.5409). It means that internal auditors' neuroticism is insignificant, indicating that internal auditors should have the characteristics to be calm under the pressure.

Discussion

The results showed that job title has a significant difference on extraversion and conscientiousness in personality traits; total working years has a significant difference on extraversion and conscientiousness in personality traits; job title and total working years have significant differences on conscientiousness in internal auditors' personality traits. Internal auditors with total working years above 6-10 years have the strongest conscientiousness. In other words, internal auditors' conscientiousness varies according to working years. The auditors with longer working years have stronger conscientiousness.

In business operation and management, the role of internal auditing will become more important and influential. In the future, the firms will value the function and benefits of internal auditing, thus making internal auditors' personality traits a critical factor in the recruitment process, so that they could fulfill the important and effective roles. This study suggests that when selecting internal auditors, the firms can use various "personality traits evaluation scales". They should select internal auditors with per-sonality traits such as "extraversion", "agreeableness" and "conscientiousness". These auditors are more competent, they enjoy internal auditing and will have higher job performance.

The findings show that internal auditors having total working years less than 3 years account for the majority (38.2%), indicating that their turnover rate is high. [START_REF] Holland | Making Vocational Choices: A Theory of Vocational Personalities and Work Environments[END_REF] found that when personality traits match the jobs, the employees will have the highest job satisfaction and work efficiency, as well as the lowest turnover rate.

Therefore, it is suggested that senior managers should value internal auditors' job satisfaction, identify the problems, and improve them by proper measures in order to reduce the turnover of internal auditors and ensure that they can function properly in the organizations, enjoy the work, and contribute to the firms.

Table 1 .

 1 Differences between the various dimensions of gender in personality

	Variable	Gen-der	Sample Size	M	SD	t	p
	Extraversion	M F	90 181	3.6611 3.6961	.58337 .54918	.127	.722
	Agreeableness	M F	90 181	3.8519 3.8131	.46234 .43485	1.901	.169
	Neuroticism	M F	90 181	2.5653 2.5311	.56533 .59130	.317	.574
	Conscientiousness	M F	90 181	3.8386 3.9517	.46375 .46668	.805	.370

Table 2 .

 2 The demographic variables in the personality differences between the various dimensions

	Varia-			N	M	SD	F	Scheffe
	ble							Post
								hoc compari-
								son
	Extra-	Age	(1) under the age of 35	125	3.6560	.59709	.343
	version		(2) 36-45 years old	106	3.7005	.51457
			(3) More than 46 years old	40	3.7313	.56440
		Education	(1) College below	52	3.7548	.40972	.621
		level	(2) University	179	3.6592	.58040
			(3) Master above	40	3.7063	.63521
		Subordina-	(1) Board of directors	35	3.7357	.51429	2.367
		tion level	(2) Chairman of the board	75	3.5633	.57992
			(3) General Manager	137	3.7591	.56017
			(4) Other	22	3.5795	.50819
		Job Title	(1) Chief Auditor	53	3.7972	.54609	3.123*	n. s.
			(2) Senior auditor	85	3.7676	.51291
			(3) Primary auditor	93	3.5538	.57678
			(4) Other	40	3.6625	.58984
		Working	(1) less than three years	104	3.6250	.58397	3.906**	n. s.
		years	(2) 3-5 years	79	3.5886	.57601
			(3) 6-10 years	65	3.8192	.49511
			(4) More than 11 years	19	3.9474	.46082
	Holistic				3.6847	.55891

Table 3 .

 3 The demographic variables in the personality differences between the various dimensions(Continue)

								Scheffe
	Variable			N	M	SD	F	Post hoc
								comparison
	Agreea-	Age	(1) under the age of 35	125	3.7760	.50184	1.491
	bleness		(2) 36-45 years old	106	3.8648	.38566
			(3) More than 46 years old	40	3.8792	.38302
		Education	(1) College below	52	3.8750	.36436	.396
		level	(2) University	179	3.8156	.45916
			(3) Master above	40	3.8083	.47223
		Subordi-	(1) Board of directors	35	3.7333	.35974	1.182
		nation	(2) Chairman of the board	75	3.8000	.37066
		level	(3) General Manager	137	3.8735	.51032
			(4) Other	22	3.7879	.28257
		Job Title	(1) Chief Auditor	53	3.8648	.44585	1.473
			(2) Senior auditor	85	3.8882	.35485
			(3) Primary auditor	93	3.7581	.52963
			(4) Other	40	3.8000	.37780
		Working	(1) less than three years	104	3.7981	.43117	1.367
		years	(2) 3-5 years	79A	3.7785	.49694
			(3) 6-10 years	65	3.9179	.40731
			(4) More than 11 years	19	3.8333	.39675


Do-Hoon Kim

Jungwook Wee
email: jwwee@keti.re.kr

Kyu-Sung Hwang
email: kshwang@kiu.ac.kr

A Novel Time Reversal-Least Sidelobe Scheme to Minimize ISI and MUI

Keywords: UWB communication, Time-reversal, Sidelobe suppression, Additional prefilter

In this paper, the scheme that combines a prefilter for sidelobe minimization with the time-reversal prefilter in the ultra wide band wireless communication system is proposed. The proposed scheme can be used in the variable symbol interval situation. When the system has the symbol interval which is larger than one-tap, the proposed scheme exploits the symbol interval to minimize the ISI. In the multi user environment, the proposed scheme is combined with space-time multiplexing scheme to alleviate multi user interference (MUI). Through numerical simulations, we examine the bit error rate (BER) performance of the proposed scheme in single user and two-user environment, and show that the proposed scheme can offer some BER performance gains comparing to the conventional TR scheme.

Introduction

Recently, in indoor wireless communication environment, ultra wideband (UWB) multimedia wireless communication systems that can transmit high capacity and high speed multimedia data has caught major attentions, because it is able to support high speed transmission of the huge data. However, the UWB operates in a frequency selective fading environment, such that intersymbol interference (ISI) arises and degrades overall performances. In order to overcome the ISI is the time-reversal (TR) scheme in [START_REF] Fink | Time-reversed Acoustic[END_REF] which applies time-reversed complex conjugated channel impulse response (CIR) as a prefilter at transmitter. The TR scheme focuses the signal energy which is scattered by the frequency selective fading to maximize SNR and reduce ISI. In addition, the TR scheme focuses spatial signal energy to the intended receiver; it exploits the channel state information (CSI) to the intended specific position as a prefilter. Therefore, the TR scheme also reduces the multi user interference (MUI). Due to the spatial signal energy focusing property, this scheme has been applied to rich scattering environment, such as a radar system and underwater communication systems [START_REF] Jin | Time Reversal Detection in Clutter : Additional Experimental Results[END_REF][START_REF] Edelmann | Underwater Acoustic Communications using Time Reversal[END_REF].

However, the TR scheme cannot perfectly eliminate ISI due to the residual interference that is remained in the form of sidelobe in the TR equivalent channel impulse response (CIR). In addition, the receiver at the unintended position receives MUI in the form of the cross correlation function between the transmitter to the intended receiver channel and the transmitter to the unintended receiver channel. To mitigate the residual ISI, several schemes have been proposed. In joint TR-ZF scheme [START_REF] Kyritsi | Time Reversal and Zero-Forcing Equalization for Fixed Wireless Access Channels[END_REF], zeroforcing (ZF) prefilter was designed in conjunction with the TR prefilter to eliminate the residual ISI and maximize SNR. The schemes in [START_REF] Nguyen | Time Reversal Technique for Multi-user Wireless Communication with Single Tap Receiver[END_REF][START_REF] Sigit | Shifted Time Reversal Technique for Two-user Wireless Communication Using Variable Rate Back-off[END_REF] exploit the symmetric property of the equivalent CIR's power to reduce the multi-user interference (MUI). In these schemes, each transmitter sends the signal, which is shifted to each other to avoid the peak-power interference from another transmitter.

In this paper, we propose a novel structure of an additional prefilter at the transmitter to suppress the residual ISI and reduce the complexity of the receiver. We consider a variable rate system that the symbol interval can be adjustable. The proposed scheme exploits the symbol interval to suppress ISI in conjunction with semidefinite relaxation (SDR) [START_REF] Luo | Semidefinite Relaxation of Quadratic Optimization Problems[END_REF] technique to solve optimization problem of minimizing ISI with some constraints. The following notations are used for description throughput this paper. The operators * , (•) H , (•) * , and ||•|| indicate convolution, Hermitian transpose, complex conjugation, and norm operators respectively.

System Models

Indoor wireless channel model

In this paper, we consider an exponential decay model in [START_REF] Vaughan | Channels Propagation and Antennas for Mobile Communications[END_REF] for an indoor wireless channel model in order to simplify the representation and the performance analysis which is useful to express a typical frequency selective fading channel. The power delay profile (PDP) of an exponential decay model indicates the mean power at the l th tap and can be expressed as

              r l ) 1 ( exp 2 1 (1)
where Δτ is an interval between adjacent taps, and   represents the mean RMS delay spread. The coefficient at the l th tap of the channel impulse response is an identically and independently distributed circularly complex Gaussian random variable with zero mean and 2 l  variance.

Single user TR-LS system

Firstly, we consider a single user system that employs a TR prefilter at the transmitter side. The main advantage of this scheme is that it maximizes the transmit signal to noise ratio (SNR) since the CSI of the receiver is available at the transmitter while also greatly reducing the receiver's complexity by only employing a single-tap equalizer. However, the drawback of the TR scheme is that it has a residual ISI. The TR scheme makes an equivalent channel impulse response which is the original CIR's autocorrelation function that has maximum peak point and symmetrical sidelobes. The TR scheme maximizes the SNR but does not remove ISI completely. Thus, the performance of the system will be degraded by the residual ISI. To suppress the residual ISI, this paper proposes a scheme that employs an additional prefilter at the transmitter. The transmitter has two filters: one is a conventional TR prefilter while the other is a least sidelobe (LS) filter that minimizes the residual ISI from the TR filter. From the assumption of TR-LS system, the received signal can be written as

) ( ) ( * ) ( * ) ( ) ( t n t s t g t f t y LS TR   (2)
where s(t) is a message signal at time t, g LS (t) is the LS filter with length 2L-1 and f TR (t) is the TR equivalent channel that can be expressed as

( ) ( )* ( ) TR TR f t h t g t  (3)
The TR prefilter g TR can be represented as * ( )

( ) || ( ) || TR h t g t h t   (4)
where h * (-t) is time-reversed and conjugated form of the h(t). In ( 4), the nominator means a conjugation of the time-reversed channel impulse response. Variable n(t) is an additive zero-mean complex Gaussian noise. We can rewrite (3) in its discrete representation

] [ * ] [ ] [ k g k h k f TR TR  (5)
where k is a discrete sampling instant. The entire transmit prefilter is a combination of the TR filter and the LS filter which can be represented as

] [ * ] [ ] [ k g k g k g TR LS LS TR   (6)
where g LS [k] is a prefilter proposed in this paper. The prefilter g LS [k] performs peak power reduction and sidelobe expansion to minimize residual ISI. We consider the LS filter with length of 2L-1, similar to the conventional TR equivalent channel impulse response's length. The equivalent entire system channel impulse response

f TR-LS [k] can be represented as [ ] [ ]* [ ]. TR LS TR LS f k h k g k    (7)

Two user TR-LS system

For the two user TR-LS system, we consider a two transmitter-receiver pairs in the indoor wireless communication environment which can cause interference to each other. We assume the CSI of the intended receiver are available at each transmitter. However, the CSI of unintended receiver is unavailable. Since the transmitters have only the CSI of the intended user, each transmitter generate the prefilters using only the intended receiver's channels, respectively. To avoid the MUI, we introduce a scheme which combines the shifted transmission technique in [START_REF] Nguyen | Time Reversal Technique for Multi-user Wireless Communication with Single Tap Receiver[END_REF] and the adjustment of the mainlobe's width.

Based on the proposed two user TR-LS system, the received signal at receiver1 can be written as are the least sidelobe prefilter at the transmitter 1 and 2, respectively, while s 1 and s 2 are the transmitted symbol at the transmitter 1 and 2, respectively. Similarly, the received signal at receiver 2 can be described as

] [ ] [ * ] [ * ] [ ] [ * ] [ * ] [ ] [ 1 2 2 21 1 1 11 1 k n k s k g k f k s k g k f k y LS TR LS TR    (
] [ ] [ * ] [ * ] [ ] [ * ] [ * ] [ ] [ 2 1 1 21 1 2 22 2 k n k s k g k f k s k g k f k y LS TR LS TR    (11) where ] [ * ] [ ] [ 2 22 22 k g k h k f TR TR  (9) and ] [ * ] [ ] [ 1 12 12 k g k h k f TR TR  (10)
In ( 8) and (11), the first term is the intended signal, the second term is MUI, and the third term is an additive zero-mean complex gaussian noise.

TR-LS scheme for ISI suppression

This section describes out proposed scheme that exploits an additional prefilter at the transmitter to minimize the residual ISI in the TR wireless communication system. To achieve this, following subsections explain how to design the additional prefilter that minimize the sidelobe of TR equivalent impulse response.

Design of the least sidelobe prefilter

Optimization Problem Firstly, we examine how to design the least sidelobe prefilter g LS (t) that is designed by an optimization theory. The prefilter g LS (t) can be described as the matrix form, LS g . The objective of the prefilter is the minimization of residual ISI. Therefore, the optimization problem can be represented as where TR F  is a matrix that represents ISI components to be minimize, and c is a peak power constraint that has a relative value to the peak power of original TR scheme. TR f is a TR equivalent channel impulse response and Rev f is a reversed vector of TR f . The matrix TR F  is a modification of the original matrix TR F . The matrix TR F represents a TR equivalent channel.

[1] 0 0

[2] [1] 0 [2 1] [2 2] [1] 0 [ 2 1 ] [ 2 ] 0 0 [ 2 1 ] TR TR TR TR TR TR TR TR TR TR f f f f L f L f f L f f L                              F          , ( 15
)
and LS g is a solution of this optimization problem. This solution is a vector form of the least sidelobe prefilter.

Focusing the residual ISI within symbol interval In (14), the solution minimizes ISI term. The ISI term is represented by TR F  . This term has to be minimized by LS g .

TR F  is a modified matrix of the equivalent TR channel matrix TR F . In ( 14), central rows are nulled with zeros. Rows of TR F  are multiplied with LS g . Then, we can get a vector that represents ISI. The optimization solution LS g minimizes ISI power. Thus, the multiplication of TR F  and LS g has to be minimized by the solution. This is expressed by ( 14). When we compose TR F  , we have to consider the focusing of ISI within the symbol interval. Thus, the solution minimizes only the symbol interval tap and outside of the symbol interval. To represent the outside of the symbol interval, we need nulling the central rows of TR F  . Then central taps of TR LS F g  are omitted in the minimization target of the solution. Thus, in an equivalent channel impulse response, the taps inside of the symbol interval are not minimized. But they are not our concerns, because the taps do not interfere sampling time instants. The ignoring taps inside of symbol interval leads increase of the degree of freedom. Thus, the solution can exploit the increase of the degree of freedom to the suppressing ISI power. If the sytem gets more taps inside of symbol interval, more ISI power can be suppressed.

Conversion to solvable form

This subsection deals with the process to convert the optimization problem to solvable form. The optimization problem which discussed in previous section can be converted to a quadratic form. ( 14) can be converted as minimize ( )

( ) ; 0 1 ( ) 1 LS LS H LS LS H LS LS H LS g g c c subject to g     g g A g B g C (16) where Re Re ( ) , ( ) ,and ( )
.

TR H TR v H v TR H TR    A F F B f f C G G  
Now, the optimization problem is a quadratically constrained quadratic programming (QCQP). However, despite the conversion, the optimization problem is nonconvex, because the second constraint in ( 16) is a quadratic form equality constraint that is a non-convex. We can use the semi-definite programming relaxation (SDR) [START_REF] Luo | Semidefinite Relaxation of Quadratic Optimization Problems[END_REF] to convexity of the optimization problem. First, the QCQP has to convert to semidefinite programming (SDP). Thus, (16) can be reformulated as minimize

This optimization problem can be solved by the convex SDP that can derive the solution as matrix form, G LS . However we want a solution of the vector form, G LS . The vector form solution can be extracted from G LS by a rank-1 approximation. To extract

g LS , let 1 , r LS H i i i i    G qq ( 19
)
where r is a rank of G LS , and the eigenvalues of G LS are λ 1 ≥ λ 2 ≥ … ≥ λ r > 0, and q 1 , …, q r are the eigenvectors of G LS . We choose maximum eigenvalue and its eigenvector, 1 1

LS   g q . This is a solution of the optimization problem.

Single user system

In single user system which doesn't have MUI, reduction of the sidelobe can be achieved by handling the mainlobe width. To minimize ISI, we exploit the symbol interval. In previous section, we examined the wider mainlobe leads to lower sidelobe. However, if the expanded mainlobe interferes the adjacent sampling instance, ISI is occurred severely. Thus, the mainlobe width will set to a little bit shorter than adjacent sampling instance to avoid ISI from the mainlobe, therefore the mainlobe doesn't interfere to the adjacent sampling instance. The mainlobe width can be set W m = 2T s -1 where W m is a mainlobe width and T s is a symbol interval.

Two user system

In Two user system, we have to pay attention to MUI which cannot handle by the TR-LS prefilter due to the prefilter does not have the CSI to unintended receiver. Thus, in this paper, the proposed scheme is combined with the space time multiplexing scheme that explained in the previous chapter. To overcome MUI, the optimization problem of the TR-LS prefilter has to set the high c constraint, because the transmitter cannot mitigate MUI. Thus, the intended channel has to secure a measure of the peak power to overcome MUI. Also, the shifted transmission scheme can alleviate MUI [START_REF] Nguyen | Time Reversal Technique for Multi-user Wireless Communication with Single Tap Receiver[END_REF].

Simulation Results

In this section, we provide the BER performance of the proposed TR-LS with QPSK and compare that with the conventional schemes. The simulation parameters Fig. 2 BER performance comparison between the conventional TR scheme and the proposed TR-LS scheme with W m = 2. used in this paper is given by the channel impulse response L = 38, Δτ = 0.811ns, and 100 symbols of block length. In two user environment, the space-time multiplexing scheme [START_REF] Nguyen | Time Reversal Technique for Multi-user Wireless Communication with Single Tap Receiver[END_REF] is applied and the minimum symbol interval of 2 is used.

Figure 1 shows the comparison of BER performance in a single user environment between the proposed scheme, the conventional TR, and MMSE-TR in [START_REF] Strohmer | Applica-tion of Time-Reversal with MMSE Equalizer to UWB Communications[END_REF] with the fixed symbol interval, and variable mainlobe width W m . As the mainlobe is broadening, the BER performance is improved. As shown in Figure 1, the proposed scheme has better BER performance than conventional schemes. Figure 2 gives the BER performance of the two user system. In this simulation, we set c'= 0.9. The result shows that TR-LS scheme has better BER performance than conventional TR scheme and shifted TR scheme both. Since TR-LS scheme is combined with shifted transmission scheme, BER performance is slightly improved. As shown in Figure 2, he proposed scheme reduces bit error by 26 percent as compared with conventional TR scheme.

Conclusion

In this paper, we proposed the TR-LS scheme to minimize ISI and MUI in UWB communication system. In the frequency selective fading channel environment, the proposed scheme exploits the symbol interval to minimize the residual ISI in the TR communication system. Since the system broadens the symbol interval, the TR-LS scheme can reduce more ISI. In addition, the proposed scheme expands the mainlobe width to reduce ISI. The expanded mainlobe width makes lower sidelobe power. Thus, the proposed scheme shows an improved signal to-ISI ratio and BER performance. From some selected simulation results, the proposed scheme shows outperform the conventional schemes.

  The equivalent channel impulse response (10) can be called an interference chanel.

  is a non-convex SDP, because the constraint rank(G LS ) =1 is not a convex. By the relaxation of this constraint, we can change the non-convex SDP to the convex SDP. The changed optimization problem can be described as minimize

Fig. 1 .

 1 Fig. 1. BER performance comparison between the variable mainlobe width W m = 1, 7, and conventional schemes.

K.-S. Hwang is a corresponding author.

Acknowledgements

This work was supported by the IT R&D program of MKE/KEIT.


Jintao M Eng

Jianrui Yuan
email: jr.yuan@siat.ac.cn

Jiefeng Cheng
email: jf.cheng@siat.ac.cn

Yanjie Wei
email: yj.wei@siat.ac.cn

Shengzhong Feng
email: sz.feng@siat.ac.cn

Small World Asynchronous Parallel Model for Genome Assembly

Keywords: parallel computing, De Bruijn graph, genome assembly

Large de bruijn graph based algorithm is widely used in genome assembly and metagenetic assembly. The scale of this kind of graphs -in some cases billions of vertices and edges -poses challenges to genome assembly problem. In this paper, a one-step bi-directed graph is used to abstract the problem of genome assembly. After that small world asynchronous parallel model (SWAP) is proposed to handle the edge merging operation predefined in the graph. SWAP aims at making use of the locality of computing and communication to explore parallelism for graph algorithm. Based on the above graph abstraction and SWAP model, an assembler is developed, and experiment results shows that a factor of 20 times speedup is achieved when the number of processors scales from 10 to 640 when testing on processing C.elegans data.

Introduction

Current sequencing technology (Illu mina Solexa [START_REF] Bennet | Solexa ltd[END_REF], Applied Biosystems SoLiD [START_REF] Pandey | Applied Biosystems SOLiDTM System: Ligation-Based Sequencing[END_REF], and Helicos Biosciences Heliscope [START_REF] Wire | Helicos biosciences enters molecular diagnostics collaboration with renowned research center to sequence cancer-associated genes[END_REF]) allo ws one to read millions of short 35 to 100 nucleotide sequences per hour. Due to experimental errors, gaps, and genomic repeats, a much higher coverage depth of 50-fold to 300-fo ld is needed for accurate assembly. These factors contribute to a 300-fold to 1000-fold increase in the number of reads, which means there are billions of reads need to be processed, and this significantly complicate the genome assembly problem.

De Bru ijn Assembler based on de-bruijn graph strategy [START_REF] Idury | A New Algorithm for DNA Sequence Assembly[END_REF][START_REF] Pevzner | An Eulerian path approach to DNA fragment assembly[END_REF] is well suitable for the current generation high throughput short reads assembly. In De Bruijn graph each vertex represents a length-k substring (k-mer) in a length-L read or its reverse complement. A directed edge connects two vertex u and v, if the k-1 length suffix of u is the same as the k-1 length prefix of v. Each input read is a path in the graph. By connecting such vertex pairs through edges, this approach will output the longest path without any branches as contigs. We denote the assemblers using De Bruijn graph strategy as De Bruijn assembler.

The first De Bruijn assembler, EULER assembler [START_REF] Pevzner | An Eulerian path approach to DNA fragment assembly[END_REF] was proposed by Pevzner, who had transformed the frag ment assembly problem into a variation of the classical Eu lerian path problem by divid ing reads into k-mers and then constructing k-mers into a path graph. This opens new possibilities for repeat resolution and generating error-free solutions of the large-scale frag ment assembly problem. Programs such as Velvet [START_REF] Zerbino | Velvet: algorithms for de novo short read assembly using De Bruijn graphs[END_REF], SOAPdenovo [START_REF] Li | De novo assembly of human genomes with massively parallel short read sequencing[END_REF], and IDBA [START_REF] Peng | IDBA-A Practical Iterative De Bruijn Graph De Novo Assembler[END_REF] imp licit ly use this framewo rk but are slightly different in local details. Velvet man ipulates these De Bru ijn graphs efficiently to both eliminate errors and resolve repeats by error correction algorith m. SOAPdenovo imp lement pre-assembler error correct ion on human genome assembly, after this operation the proportion of error free reads was improved fro m 64% to 70%, and nearly 60% percent of k-mers was filtered fro m the graph. IDBA also adopt preassembler error filtering technique, wh ich can save nearly 40-80% of memo ry compared with velvet. The second feature of IDBA is that it iterates fro m s mall k-mer to large k-mer to get longer contigs. So the quality of contigs is better than other tools.

The above assemble tools can only run on single machine, the human genome assembly with current sequencing technology needs about 1TB memo ry and takes weeks or even months on single server. The situation will be even worse for larger genome assembly or meta-genome assembly.

Parallel algorith m for sequencing assembly is an alternative to solve the problem. Parallel assemblers included ABySS [START_REF] Simpson | ABySS: a parallel assembler for short read sequence data[END_REF] and YA GA [START_REF] Jackson | Parallel Construction of Bidirected String Graphs for Genome Assembly[END_REF][START_REF] Jackson | Parallel short sequence assembly of transcriptomes[END_REF][START_REF] Jackson | Parallel de novo assembly of large genomes from high-throughput short reads[END_REF], are both based on De Bru ijn g raph strategy. ABySS distributes k-mers to mult i-servers to build a d istributed De Bruijn graph, and error removal and vertex merg ing were imp lemented over MPI communicat ion messages. YA GA constructs the distributed bi-directed De Bruijn graph by maintaining edge tuples in a community of servers. Unanimous chain compaction problem in YA GA was transformed to undirected list ranking, and then the authors designed a modified sparse ruling set algorith m for undirected lists. The computational comp lexity of YA GA is given by O(n/p) compute time, O(n/p) communicat ion volume, and O(log 2 n) communication rounds, where n is the number of nucleotides in all reads, p denotes the number of processors .

Efficient and scalable framewo rks or libraries for d istributed graphs are essential to parallel assembly based on De Bruijn graph. Existing works, such as BSPlib [START_REF] Richard | A Library for Bulk-Synchronous Parallel Programming[END_REF][START_REF] Goudreau | Portable and Effcient Parallel Computing Using the BSP M odel[END_REF][START_REF] Bonorden | The Paderborn University BSP (PUB) Library[END_REF], CGMgraph [START_REF] Chan | CGMGRAPH/CGM LIB: Implementing and Testing CGM Graph Algorithms on PC Clusters and Shared M emory M achines[END_REF], PBGL [START_REF] Gregor | The Parallel BGL: A Generic Library for Distributed Graph Computations[END_REF][START_REF] Gregor | Lifting Sequential Graph Algorithms for Distributed-M emory Parallel Computation[END_REF], Prejel [START_REF] Austern | Pregel: a system for large-scale graph processing[END_REF], are based on BSP [START_REF] Valiant | A bridging model for parallel computation[END_REF] model. The BSP model has advantage on simple co mputation-co mmunication programming model, whereas the barrel principle exists in the co mputation -communication phase and the synchronous phase over large clusters limits the scalability of the model. To our knowledge, the scalability these imp lementations under BSP model has not been evaluated beyond several hundreds of co mputers [START_REF] Austern | Pregel: a system for large-scale graph processing[END_REF]. No genome assembly tools have adapted the BSP library, although YA GA has used the BSP idea in its design on parallel list ranking algorith m implicitly. Another parallel programming model, MapReduce [START_REF] Dean | M apReduce: simplified data processing on large clusters[END_REF], has strength in loosely coupled work such as frequency statistics, sorting, indexing, and mach ine learning etc, and these works can be easily d istributed to clusters. However graph algorith m is a tight coupled work and dividing one graph into several meaningful sub-graphs is still a challenging problem.

This paper first demonstrates a one-step bi-directed graph for the problem of genome assembly. Genome can be recovered by merg ing semi-extended edges to fullextended edges. Then small wo rld asynchronous parallel (SWAP) model is proposed to realize edge merging over a distributed one-step bi-directed graph. Specially, we implement an assembler using the SWAP model. Given the nu mber of p rocesses p, the complexity of this problem is reduced to O(n/p) parallel co mpute time, O(n/p) communicat ion round, and O(nlog(n)/p) communication volu me, here n is total length of input sequences. Simu lation shows that Assembler has a factor of 20 t imes speedup when the number of processors scales from 10 to 640.

The rest of the paper is organized as fo llo ws: Section 2 abstracts the De Bru ijn graph based genome assembly problem; Sect ion 3 describes the SWAP model for large scale graphs with s mall world property, then an assembler, as SWAP's first application, is illustrated. Experimental results will be present in section 4. Finally section 5 concludes this paper.

Abstraction of De Bruijn Assembly

Let l sM  be a string of length L, where   , , , a t c g  . Any substring [ ] [ 1]... [ 1], s j s j s j k      01 j L k     is a k-mer o f s.
The set of all k- mers of a given string s is written as ( , ) sk , here k must be odd. The reverse complement of a k-mer  , denotes by   , is obtained by reversing  and complementing each base (

[ ] [ 1] i k i       ) by the following bijection o f  ,   : , , , a t t a c g g c      . Note that [ ] [ ] ii    and    . A k-mo lecule  is a pair of co mplementary k-mers { , }
 . Let  be the partial ordering relation among the string of equal length such that   indicates that the string  is lexicographically larger than  . We designate the lexicographically larger of the two co mplementary k-mers as the positive k-mer, denoted as   , and the lexicographically s maller one as the negative k-mer, denoted as   , here    . We choose the positive k-mer   as the representative k- mer  of the k-mo lecule. The set of all k-molecules of a given string s is called the k-spectrum of s and is written as E is defined as follows:

( , ) sk . Noted that ( , ) ( , ) s k s k   . The notation su f (a, l) (pre(a, l),
1 1 1 { ( , , , , ) | , ( , ), ( , 1)
( , 1) ( [ 1] ( ( , 1) ( , 1)))} s E e d d c s k suf k pre k k s k s k                           (2)
Equations ( 2) declares that any two overlapped k-molecu les can be connected with a 1-step bi-directed edge, if they are continuous in sequence s or its complementary.

Here  , there will be four possible connections, and for each type of connection exactly two equivalent 1-step bi-directed edge exist, 1.

d  is the direction of k-mer  , if    , '' d  , otherwise '' d  . Set 1 c  is initialized with one element [ 1] k   ,
1 1 1 1 ( , , , , ), ( , , , , ) e c e c
                              , 2. 1 1 1 1 ( , , , , ), ( , , , , ) e c e c
                              , 3. 1 1 1 1 ( , , , , ), ( , , , , ) e c e c                               , 4.
                             .
In each type of connection, the first bi-directed edge and the second one correspond to the same bi-d irected edge, but in different form. Within a distributed edge representation situation, the first bi-d irected edge in each type will be attached with kmo lecule  , and the second one will be with  .

(S) { , } { , } ii k S S s s i h i h G V E V E    



. The key property of this bi-d irected De Bru ijn graph

1 (S) k G
is that each read can be recovered by traversing the corresponding path in either direction, concatenating (k-1)-molecu le prefix o f the first node and the edge labels on the path. As all input reads of assembler are fro m chro mosomes, each chro mosome can now be seen as a long path in this graph. However because of read errors, and repeats in the sequence, we cannot expect to see continuity in sampling, our goal is to recover the genome as a large set of contigs by merging semi-extended edges into full-extended edges.

Assembler over SWAP

Vertices in large scale real world graph (such as social network, web lin k graph, et) always have limited number of neighbors, little computing work, and constant number of edges randomly connected to other vertices, this phenomenon is denoted as small world property. For a g iven vertex, its small world includes all its edges, neighbors and itself. Then any computing and communication work of a vertex can be done in its small world. As long as the work of a vertex on a graph does not interrupt others, we can run computational work of those vertices in parallel. Here we will present our work on pursuing parallelism in the computation of bi-directed graph for genome assembly.

Inspirited from CSMA/CA in wireless networks [START_REF] Tanenbaum | Computer Networks[END_REF], Small World Asynchronous Parallel model (SWAP) aims to improve parallelism on processing large scale graph problem with small world property. After having distributed graph over a network of processors, the main schedule of SWAP can be defined as a combination of following three steps:

1. Lock operation is applied to each vertex's small world, which includes itself and its neighbors. 2. Computati on and modification will be performed in each vertex's small world.. 3. Unlock operation will be triggered after each computation step. The basic schedule of SWAP is Lock-Co mputation-Unlock. Because of the locality of computing and communication in the small world, SWAP model utilizes local synchronization and global asynchronization mechanism to maximize underline parallelism for the graph algorithm.

An assembler over SWAP is the first application using SWAP model. In the following paragraphs we will describe its data structure on distributed de brujin graph, strategy on error removal, and the edge merging algorithm, respectively. Input sequences can be broken into overlapping k-mo lecules by sliding a window of length k along the input sequence. Each processor maintains a hash table to store k-mo lecules, and each k-molecule is represented as a base-4 number of its positive k-mer. Numerical values {0,1,2,3} are assigned to bases {A,C,G,T}. The location of a given k-molecule can be co mputed by take the mod of a large prime number and then take the mod of the number of processors. The large prime nu mber is used to evenly distribute k-molecules to all processors.

Parallel constuction of distributed one-step bi-directed De Bruijn graph

A single k-molecule can have up to eight edges, and each of them corresponds to a one-base extension, {A, C, G, T} in either direction. The adjacent k-mo lecule can be easily generated by adding the base extension in the edge set to the source k-molecule.

The construction of one-step bidirected De Bruijn graph can be achieved in O(n/p) parallel co mpute time, O( 1) round of all-to-all commun ication, and O(n/p) parallel communication volume.

Error removal

Sequencing errors make the assembly problem more co mplex. To identify erro rs, we assume that the errors are rando m, and they are unlikely to occur twice in the same base. As each base in the genome is sampled on average as many times as the overage number, the erroneous k-molecule will have lower frequency compared to the correct ones. According to this principle, we identify all k-molecules with low frequency as erroneous k-molecu les, and delete all of them fro m our vertex set of the graph. The complexity of this step is O(n/p) parallel compute time.

Edge merging

One step bidirected graph generated in the previous section will likely have many long chains, and each corresponds to a sequence that can be unambiguously assembled into a single contig. We will merge these chains into full-extend edges using Algorithm 1.

A subset of vertices and their associated edges are stored in processor i, and the set of semi-extended vertices stored in this processor is denoted as i V . We aim to delete all semi-extended vertices and merge their associated semi-extended edges to form full-extended edges.

As the bi-directed graph

Experimental Results

The assembler is written in C++ and MPI. The hardware and software architecture supporting this assembler is demonstrated in Fig 2 . We use Dawning 5000 as high performance cluster, which has 40 16-core servers with 32GB memo ry. The distributed file system is lustre. All the components are interconnected with infiniband 20Gbit Router.

Perl scripts [START_REF] Zhang | A practical comparison of de novo genome assembly software tools for next-generation sequencing technologies[END_REF] are used to generate the following two theoretical datasets: 50x coverage of Yeast chromosomes containing 17 million reads, and 50x coverage of C.elegans chromosomes containing 141 million reads. The error rate is set to be 1%, and the length of reads ranges fro m 36bp to 50bp. The p rimary goal of this experiment is to demonstrate the scalability of SWAP model on handling large-scale graphs using parallel system with distributed memory.

We first test the performance of this assembler on Yeast dataset. The runtime of assembler is displayed in figure 3, and the time is div ided into three phases, Parallel File I/ O, graph build ing, and edge merging. The first phase is the time spent on reading dataset from a distributed file system, the second phase is the time used to construct the one-step bi-directed graph over the cluster, and the last phase is the time cost on edge merg ing operations. The run time is dominated by the third phase, where hundreds of processors are sending messages to lock their neighbors, merging edges, deleting semi-extended nodes and edges, and unlocking their neighbors. Figure 3 shows that this phase has good scalability. The speedup is about 50 when the number of processor scales from 10 to 640 and the overall runtime of assembler on Yeast dataset is reduced by a factor of 30.

The C.elegans dataset is nearly ten times larger than Yeast, and its corresponding data on time usage is demonstrated in figure 4. In this figure, the time used in parallel file I/O is very short co mpared to the other two phases. The graph construction phase has a decreasing trend on its running time. This phase have a speedup of 35x. This speedup is slightly smaller than that of the yeast dataset. The total speedup of all three phases on C.elegans dataset is about 20.

Definition 1 .Definition 2 .

 12 respectively) is used to denote the length l suffix (prefix, respectively) of string a. Let the symbol denotes the concatenation operation between two strings, and the number of edges attached to k-molecule  is denoted as degree(  ). The number of edges pointing out fro m k-mo lecule  is denoted as deg ( ) ree  . The vertex set Vs is defined as k-spectrum of s, The 1-step bi-directed edge set 1 s

Figure ( 1 )

 1 illustrates four possible connections and examples of a 1-step bi-directed edge graph.

Fig. 1 . 11 ( 2 ( and 1 e

 11121 Fig. 1. The illustration of four possible connections. Definition 3. 1-step bi-directed de bruijn graph of order k for a string s can be: 11 (s) { , } k s s G V E 

  m sequences sampled fro m a genome of length g, the total length of all these sequence is n. we aim to construct a one-step bi-directed De Bruijn graph 1 (S) k G with O(n) vertexes and edges distributed among p processors such that each processor store O(n/p) vertices.

  over p processors, each processor will store a subset of semi-extended k-mo lecules V i and the average number of kmo lecules in V i is O(n/p). Then the expected computational complexity of each processor on edge merg ing is given by O(n/p) parallel co mpute time, O(n/p) communication round, and O(nlog(n)/p) communication volume.

Fig. 2

 2 Fig.2 Time usage analysis in three phase on Yeast dataset

Figure 1 .

 1 Figure 1. Hardware and software architecture of the assembler on SWAP.

Fig. 3

 3 Fig.3 Time usage analysis in three phase on C.elegans dataset

  and

			suf	1 c k , )    . (
	Property 1. Given t wo k-mo lecules	,	( , ) sk

Acknowledgements

This work is supported by NSFC (Grant No. 61103049) and Shenzhen Research Fund (Grant No.JC201005270342A). The author also thanks Bingqiang Wang from BGI, and Prof. Francis Y.L. Chin from HKU for their suggestions on this work.


Daniel Orozco

Elkin Garcia

Robert Pavel

Orlando Ayala

Lian-Ping Wang

Guang Gao
email: ggao@capsl.udel.edu

Demystifying Performance Predictions of Distributed FFT3D Implementations

This paper presents a comprehensive story of the development of simpler performance models for distributed implementations of the Fast Fourier Transform in 3 Dimensions (FFT3D). We start by providing an overview of several implementations and their performance models. Then, we present arguments to support the use of a simple power function instead of the full performance models proposed by other publications. We argue that our model can be obtained for a particular problem size with minimal experimentation while other models require significant tuning to determine their constants.

Our advocacy for simpler performance models is inspired by the difficulties found when estimating the performance of FFT3D programs. Correctly estimating how well large-scale programs (such as FFT3D) will work is one of the most challenging problems faced by scientists. The significant effort devoted to this problem has resulted in the appearance of numerous works on performance modeling.

The results produced by an exhaustive performance modeling study may predict the performance of a program with a reasonably good accuracy. However, those studies may become unusable because their aim for accuracy can make them so difficult and cumbersome to use that direct experimentation with the program may be preferable, defeating their original purpose. We propose an alternative approach in this paper that does not require a full, accurate, performance model. Our approach mitigates the problem of existing performance models, each one of the parameters and constants in the model has to be carefully measured and tuned, a process that is intrinsically harder than direct experimentation with the program at hand. Instead, we were able to simplify our approach by (1) building performance models that target particular applications in their normal operating conditions and (2) using simpler models that still produce good approximations for the particular case of a program's normal operating environment.

We have conducted experiments using the Blue Fire Supercomputer at the National Center for Atmospheric Research (NCAR), showing that our simplified model can predict the performance of a particular implementation with a high degree of accuracy and very little effort when the program is used in its intended operating range.

Finally, although our performance model does not cover extreme cases, we show that its simple approximation under the normal operating conditions of FFT3D is able to provide solid, useful approximations.

Introduction

The Fast Fourier Transform (FFT) [START_REF] Cooley | An algorithm for the machine calculation of complex fourier series[END_REF] is an important tool in physics and engineering, and it is relevant to numerous High Performance Computing (HPC) applications.

This paper addresses the issue of how to efficiently select an FFT algorithm for the particular case of large-scale distributed scientific simulations. Examples of largescale scientific applications that use FFT include Cloud Physics [START_REF] Ayala | hybrid approach for simulating turbulent collisions of hydrodynamically-interacting particles[END_REF], Molecular Dynamics [START_REF] Laasonen | \Car-parrinello molecular dynamics with van-derbilt ultrasoft pseudopotentials[END_REF], calculation of Coulomb Energies [START_REF] Bylaska | \Parallel implementation of the projector augmented plane wave method for charged systems[END_REF], Seismic Imaging [START_REF] Calandra | massively parallel implementation of the common azimuth pre-stack depth migration[END_REF] and Computational Geosciences [START_REF] Stellmach | efficient spectral method for the simulation of dynamos in Cartesian geometry and its implementation on massively parallel computers[END_REF].

Understanding the performance of FFT in large scale systems is important because FFT plays a significant role in many scientific programs. For example, the turbulent droplet collision simulation developed by Wang et al [START_REF] Wang | \Towards an integrated multiscale simulation of turbulent clouds on petascale computers[END_REF] uses a pseudo-spectral method for the flow simulation that requires multiple 3-Dimensional FFTs during each timestep of the computation. In fact, profiling of the flow simulation alone reveals that the 3D FFT alone can take up to 90\% of the total execution time in the program.

Because the FFT is not the objective of the computation in many cases, but rather a tool used to help it, it is of paramount importance to understand the effect of design choices on the overall program. For example, there may be a certain overhead in the passing of arguments to the FFT. The partitioning of the data may force the use of certain FFT strategies, or even the physical location of the compute nodes may greatly affect the performance of particular FFT strategies.

The necessity to understand the performance of FFT implementations and its relationship to the program using it has prompted researchers to build performance models that aim to predict the performance of particular implementations. These previous approaches [START_REF] Dmitruk | Scalable parallel fft for spectral simulations on a beowulf cluster[END_REF][START_REF] Li | 2decomp fft a highly scalable 2d decomposition library and fft interface[END_REF][START_REF] Pekurovsky | Ultrascalable fourier transfroms in three dimensions[END_REF][START_REF] Ayala | Parallel implementation and scalability analysis of 3d fast fourier transform using 2d domain decomposition[END_REF] have resulted in very detailed performance models that include many factors such as computation time, communication time, the presence of overlapping or pipelining, the location of the memory and so on.

The biggest drawback to building a performance model is that correctly developing each performance model requires significant effort: the algorithm must be considered, constants in the model have to be measured, and the model may only work for certain architectures or interconnects. The difficulties of building performance models are so great that when one is correctly found it is published in a recognized venue (e.g. Dimitruk published his model in the Parallel Computing Journal [START_REF] Dmitruk | Scalable parallel fft for spectral simulations on a beowulf cluster[END_REF] and Ayala's model [START_REF] Ayala | Parallel implementation and scalability analysis of 3d fast fourier transform using 2d domain decomposition[END_REF]is currently under review in the same journal).

Instead, we have looked at alternatives that would allow scientists to make decisions about 3-Dimensional FFT implementations without the burden of developing full performance models. We have approached the problem by studying the nature of the 3D-FFT computation itself and existing performance models for different 3D-FFT strategies to see if there is a simpler way to efficiently select a particular 3D-FFT implementation over a set of alternatives.

We have sought to understand the performance of different 3D-FFT implementations on distributed memory machines through preliminary experiments (Figure 1) conducted on the Bluefire Supercomputer from the National Center for Atmospheric Research (NCAR).

Figure 1 shows the results of using two very different implementations of FFT3D on data of the same size. The differences in performance between the two techniques shown are due to differences in the implementation such as communication strategy and algorithm used. Our experiments, which focus on large-scale systems, suggest an alternative approach to analyze FFT3D implementations: The running time of a particular FFT3D implementation can be approximated by a line in a logarithmic graph. This approximation corresponds to a power function. The power-function approach has the advantage of being significantly easier to use and build while still allowing the estimation (albeit with reduced accuracy when communication latency is high) of the performance of particular 3D-FFT implementations.

To test our hypothesis regarding the validity of our power-law approximation, we have assembled several 3D-FFT implementations. A brief overview for each one of them is given in Section 2, along with their performance models, when they are available.

We develop our prediction model in Section 3 and we conduct experiments to test the validity of our model in Section 4. Finally, we provide some conclusions in Section 5.

Background on FFT3D Algorithms

Many techniques have been investigated and developed to perform a distributed 3dimensional FFT. For the purpose of this study, we focused on approaches that did not change the underlying FFT algorithm but instead changed the decomposition of the data and the communication patterns. The algorithms presented here leverage the principle that a 3-dimensional FFT can be performed by executing a 1-dimensional FFT along each one of the three dimensions of the data. We will use the notation in Table 1 to aid us in our explanation of the FFT3D algorithms and models. The first decomposition we studied is the 1D domain decomposition, as seen in Figure 2. This decomposition is often referred to as a ``slab'' decomposition where each processor is responsible for a single slab. In the case of Figure 2, each processor initially holds a slab of size N*N*(N/P).

Because a slab must have a width of at least one element, the parallelism of 1D decomposition implementations is limited to using at most N processors. Each of the P processors is responsible for a single slab. First, each processor performs N/P 2D FFTs along the xy plane. Then, a transpose is performed such that the data in the z direction is now saved in the y direction and vice versa. With this transposed data, the required 1D FFTs are performed. Finally, a second transpose may be performed to restore the data to its original orientation, depending upon the needs of the application.

However, the essence of Dmitruk's work is in his implementation of the transpose. There are three methods that work on the principle of breaking the slab into ``blocks'', moving the blocks to the appropriate processors, and then locally rearranging the data. What differs is the communication scheme used. The first method involves a cyclic strategy: Each processor sends data to the next processor in the list while receiving data from the previous processor. This is repeated until every processor has traversed the entire list of processors. The second strategy is similar in that each processor sends to the processor that is currently receiving from via a temporary pairing strategy. This is repeated in parallel until every processor has paired with every other processor. Finally, MPI collectives in the form of the ``all-to-all'' operator were used for the third method. Of these three methods, the second method's pairwise strategy gave the best performance under Dmitruk's tests.

Equation 1 shows the full performance model presented by Dmitruk's paper [START_REF] Dmitruk | Scalable parallel fft for spectral simulations on a beowulf cluster[END_REF].

(1)

2D Domain Decomposition: "Pencils"

A 2D domain decomposition strategy is a natural extension to the 1D idea. The 2D domain decomposition can be done along X axis, referred to as X-pencil, Y axis (Ypencil) or Z axis (Z-pencil).

Assuming that the P processors are arranged in a grid so that P=P y *P z , Figure 2 shows a case where each pencil is of size N * (N/P y ) * (N/P z ) and each processor is responsible for a single pencil. The three libraries studied below follow the same general pattern. Three sets of 1D FFTs are performed along each axis with a transpose used to reorder data. In the case of Figure 3, this means 1D FFT in the x-direction, transpose, y-direction, transpose, and then z-direction. Depending on the application, a final transpose will be performed to restore the data. These libraries differ in how the transpose is handled.

2DECOMPFFT

2DECOMPFFT [START_REF] Li | 2decomp fft a highly scalable 2d decomposition library and fft interface[END_REF] is designed for applications using 3-dimensional structured mesh and spatially implicit numerical algorithms. It implements a general-purpose 2D pencil decomposition for data distribution on distributed memory platforms. The transpositions are done using the MPI ``all-to-all'' operator, but with a complex communication pattern according to the orientation of pencils and their associated memory pattern.

p3dfft.

The p3dfft algorithm [START_REF] Pekurovsky | Ultrascalable fourier transfroms in three dimensions[END_REF] is an algorithm developed by Pekurovsky that uses a 2D pencil decomposition with the objective of maximizing parallelism. p3dfft was writ-ten in Fortran and MPI and uses FFTW [START_REF] Frigo | Fftw: An adaptive software architecture for the fft[END_REF] as the underlying library for FFT. Communication operations are done using MPI ``all-to-all'' operations.

Pekurovsky's FFT3D requires a total of 3 transpositions to compute a transformation. Data is redistributed during each step using call to MPI All-to-All that requires internal redistribution of data. In each of the 3 communication steps, each processor transfers a block of size N 3 /P with a bandwidth of B W . Equation 2 summarizes these steps into a performance model. As with Dmitruk's work, the optimization of the transpose plays a large part in improving performance. In this case, Ayala chose to follow Dmitruk's cyclic strategy. This provides better performance because during a transpose for a 2D Decomposition, only processors in the same plane as the transpose need to communicate. For example, a transpose along the x and y axes will only require communication between processors in the same xy plane to communicate.

Equation 3 presents a performance model taken from Ayala and Wang's paper [START_REF] Ayala | Parallel implementation and scalability analysis of 3d fast fourier transform using 2d domain decomposition[END_REF].

(3)

3D Domain Decomposition: "Blocks".

Finally, we studied the 3D Domain Decomposition in which the data is broken up into blocks. The P processors that participate in the computation are arranged in a 3-Dimensional grid so that P =P x * P y * P z . The amount of data held by each processor (Figure 4) is (N/P x ) * (N/P y ) * (N/P z ). 3Decomp. 3Decomp computes FFT3D by executing four data transpositions as shown in Figure 4. In each case the communication is done using nonblocking point-to-point MPI calls. Due to the differences in topology between each step, each processor needs to communicate with several other processors. The local 1-Dimensional FFT transforms required are computed using the FFTW [START_REF] Frigo | Fftw: An adaptive software architecture for the fft[END_REF] library. The performance model of 3Decomp (Equation 4) is obtained by analyzing its data movement and the amount of computation per-formed: Four communication steps are used where all data held by each processor (N 3 /P in size) is reorganized locally and then sent to another processor. Each processor, in parallel, computes 3N 2 /P onedimensional FFTs, each of complexity N \ log(N). ( 4)

3

A simplified model approach for FFT3D

Section 2 presented a selection of algorithms for FFT3D and their performance models.

Correctly developing and using performance models such as those presented in Equations 1, 2, 3 and 4 require a significant effort: It requires a reasonable understanding of the implementation, and it requires extensive testing on the machine to identify the constants of the model.

Our previous experiments (Figure 1) showed a case where plotting execution time as a function of the number of processors resulted in a straight line when plotted on a logarithmic scale, indicating a power-function behavior. These results have motivated us to approximate the performance models of FFT3D as a power function as presented in Equation 5.

(5) Equation 5 is significantly simpler than other performance models, such as those presented in Section 2, and it is appropriate to represent the behavior observed in Figure 1.

The power-function approximation is intended to provide estimates of parallelism for large scale systems in a range of situations where there is enough parallelism available. i.e. Equation 5 will provide a reasonable approximation for a particular number of processors if the FFT3D implementation is able to scale reasonably well up to that number of processors.

It is worth noting that our approximation doesn't account for the startup latency of a message. Thus, as P approaches N, the accuracy of our model may decrease. The following section provides experiments that test our hypothesis and the usability of our simplified performance model, even under these circumstances.

Experiments

We conducted experiments to confirm the validity of our claims. Specifically, we wanted to test whether or not the use of the power-law model of Equation 5was sufficient to predict the performance of particular implementations of FFT3D.

Our experiments were conducted on the Bluefire Supercomputer at the National Center for Atmospheric Research (NCAR). The Bluefire supercomputer possesses 4096 Power6 processors running at 4.7GHz and having 2GB of memory. Bluefire's processors are connected by eight 4X infiniband DDR links with a total bandwidth of 2.5GB/s. Further information about Bluefire's hardware can be obtained from NCAR's website [START_REF][END_REF].

The experiments seek to establish whether or not it is possible to predict the performance of FFT3D implementations using the model of Equation 5. As explained in Section 3, the power-law model attempts to predict the performance of a particular implementation in common situations.

To test our hypothesis, we have ran the implementations described in Section 2. In our experiments, we chose several input sizes (N 3 =256 3 , N 3 =512 3 and N 3 =1024 3 ) for FFT3D that are typical of large-scale programs such as in cloud microphysics studies by direct numerical simulations. For each implementation, we run each problem size while the number of processors was changed. Due to memory, execution time, or processor count constraints, not all combinations of problem sizes and processor counts were possible.

The results of our experiments are reported in Figures 5,6 and 7. To confirm the validity of our model, we have built simplified performance models for each implementation in the style of Equation 5. For each case, the constants A N and B N found in Equation 5 have been found using a least squares regression model. The resulting models are presented in Table 2. The table shows the models as well as an indication (R 2 parameter) of how well the models are able to represent the experimental data obtained (R 2 =1 means perfect match).

As can be seen from Figures 5, 6 and 7 and Table 2, the power function approximation works remark-ably well. In most situations, the R 2 parameter is greater than 0.99, indicating that the model was able to fit the data with a high degree of accuracy.

Summary and Conclusions

We demonstrated that a simplified performance model is enough to capture the behavior of distributed-memory implementations of FFT3D in the particular case of large-scale systems. We have targeted our efforts at predicting the performance of implementations under their typical operating conditions. That is, when enough parallelism is available and the number of processors used is reasonably large.

Preliminary results showed that the performance of FFT3D implementations followed a power-function trend. Our full set of experiments, conducted on NCAR's Bluefire supercomputer has shown that the use of a power function model is adequate to represent the performance of FFT3D implementations on traditional distributedmemory machines.

The R 2 indicator, which provides a quantitative way to measure the fidelity of an approximation, shows, in all cases, a high degree of accuracy for our method. In most cases, the R 2 indicator was greater than 0.99, further validating the effectiveness of our approach.

An important feature of the power function model is that obtaining it re-quires very little experimentation. The constants of the model can be obtained with as little as two measurements, while still producing useful predictions.

Our current model has been designed to produce predictions related to the execution of programs over data of a particular size. To study data of a different size, experiments must be run to determine the adjusted constants. Our future work will focus on extending our model to be a function of problem size in addition to the number of available processors.

Fig. 1 .

 1 Fig. 1. Performance of two different FFT3D techniques. The relationship between time to execute and the number of processors appears to follow a straight line in this logarithmic plot.

Fig. 2 .

 2 Fig. 2. 1D domain (slab) decomposition for FFT 3D

Fig. 3 .

 3 Fig. 3. 2D domain (pencil) decomposition for FFT 3D

( 2 )

 2 Ayala's 2D Domain Decomposition. Ayala and Wang[START_REF] Ayala | Parallel implementation and scalability analysis of 3d fast fourier transform using 2d domain decomposition[END_REF]built upon the work of Dmitruk et al. As with the other libraries, a 2D Pencil Decomposition is used. Due to the target application of Direct Numerical Simulation, the last transpose is not necessary in Ayala's method.

Fig. 4 .

 4 Fig. 4. 3D domain (block) decomposition for FFT 3D

Fig. 5 .

 5 Fig. 5. FFT3D performance for N 3 =256 3 . Solid lines show actual experimental data. Dashed lines show a prediction using a power function.

Fig. 6 .

 6 Fig. 6. FFT3D performance for N 3 =512 3 . Solid lines show actual experimental data. Dashed lines show a prediction using a power function.

Fig. 7 .

 7 Fig. 7. FFT3D performance for N 3 =1024 3 . Solid lines show actual experimental data. Dashed lines show a prediction using a power function.

Figures 5 ,

 5 6 and 7 also give an intuitive indication of the quality of our models: The predicted models, represented with dashed lines, agree with the data gathered (plotted with solid lines).

Table 1 .

 1 Parameter and constant definitions

2.1

1D Domain Decomposition: "Slabs"

Acknowledgments

This research was made possible by the generous support of the NSF through grants CCF-0833122, CCF-0925863, CCF-0937907, CNS-0720531, and OCI-0904534.

Table 2. Power-function models for several FFT3D implementations

The high level of similarity between the solid lines (measured) and dashed lines (model) and the high value of the R 2 indicator lead us to conclude that our model is an effective way to predict the performance of FFT3D implementations in their normal operating range. This is an important conclusion because the simple, power-function model can be quickly obtained with as few as only two measurements, one at low processor counts, and one at high processor counts. This is significantly simpler than the other models presented in Section 2 while still presenting good accuracy.

Of particular interest was the issue of startup latency and congestion when the maximum parallelism of a problem is approached. This is most evident for Dmitruk's 1D Domain Decomposition when N=P. Figure 5 does show this limitation. However, as the problem size is increased, the impact of the startup latency is greatly decreased, and our approximation's accuracy remains high.

The experimental information obtained here supports our hypothesis: A simplified performance model for FFT3D is enough to make sufficiently accurate predictions when those predictions pertain normal operating conditions. Our model has the advantage of being significantly simpler while still being useful.