Analyzing GPU Tensor Core Potential for Fast Reductions

The Nvidia GPU architecture has introduced new computing elements such as the \textit{tensor cores}, which are special processing units dedicated to perform fast matrix-multiply-accumulate (MMA) operations and accelerate \textit{Deep Learning} applications. In this work we present the idea of using tensor cores for a different purpose such as the parallel arithmetic reduction problem, and propose a new GPU tensor-core based algorithm as well as analyze its potential performance benefits in comparison to a traditional GPU-based one. The proposed method, encodes the reduction of $n$ numbers as a set of $m\times m$ MMA tensor-core operations (for Nvidia's Volta architecture $m=16$) and takes advantage from the fact that each MMA operation takes just one GPU cycle. When analyzing the cost under a simplified GPU computing model, the result is that the new algorithm manages to reduce a problem of $n$ numbers in $T(n) = 5\log_{m^2}(n)$ steps with a speedup of $S = \frac{4}{5}\log_2(m^2)$.

I. INTRODUCTION
In the present day, the amount of data generated by all sorts of media grows at an exponential rate, making each year more difficult to process these data sets entirely. On top of this, the rise of applications such as autonomous vehicle platforms, pattern recognition in images, disease diagnosis, weather forecasting, financial modeling, robotic and language translation, among others impose an ever greater challenge which is to perform in real-time. One important computational pattern in many applications is the arithmetic reduction which computes a single value from a set of n elements. Arithmetic reductions are employed in many of the applications recently mentioned as well as in Monte Carlo methods, physical simulations such as the gravitational n-body problem and ray tracing, among many others. Therefore, it is of high interest to count with a fast reduction method and give the possibility to reach real-time performance in these applications.

A. The GPU Programming Model
Today GPU Computing [12] is a useful computational tool for attempting the solution to large problems in the big data era, as it offers massively parallel computation for data-parallel problems. This means that elements in a large data-set can actually be processed by a parallel algorithm if the problem can be defined as data-parallel. GPU Computing has been used successfully in the recent years, giving important performance speedups in physics [3], [13], all-pairs problems [11], medicine [8], image processing [4], deep learning [18] and computer graphics [5], among other fields.
When working with GPU Computing, its programming model plays an important role in the design and development of GPU accelerated programs, and it is governed by a threelevel hierarchy, which defines the parallel space of the GPU. This hierarchy corresponds to the thread, block and grid (See Figure 1). A thread is an abstract compute element that is in charge of executing the kernel once. A block is a container of threads, and has the property that all of its threads can synchronize among themselves and can share information via cache. Synchronization among blocks is not possible in the current programming model, unless the kernel is terminated. The grid is the last construct, and it contains all the blocks, spatially identified, that will be employed in the execution of the kernel. Lastly, the kernel is the function that will execute in GPU by all threads. With the help of programming tools such as OpenCL or CUDA, one can implement a parallel GPU algorithm that will accelerate the data-parallel computations.

B. GPU Tensor Cores
With the latest rise of Machine Learning applications, and more specifically the fast adoption of Deep Learning in multiple fields of science and technology, CPU and GPU Companies have started including application specific integrated circuits (ASICs) to their processors to further accelerate the computational tasks involved in the phases of training and inference in Deep Learning applications. This change has led to the inclusion of Tensor Core (TC) units in recent Nvidia GPUs, which are special purpose processing units that sit next to the GPU CUDA cores in the streaming multi-processors (SM) of the chip, as shown in Figure 2. Fig. 2. A processing group of the Nvidia Tesla V100 which has a total of 640 tensor cores. Image inspired from the Nvidia CUDA programming guide [16].
The aspect that make tensor cores an attractive feature is the performance it can offer in comparison to the traditional GPU cores. Today, the Nvidia Volta GPU Tesla V100, Quadro V100 and Titan V all include around 640 tensor cores, and they can offer up to 120 TFLOPS in mixed FP16-FP32 precision. In comparison, the traditional CUDA cores, which are 5120 in total for the GPUs recently mentioned, offer up to ∼ 15 TFLOPS of performance in FP32 precision and around ∼ 7 TFLOPS in FP64 precision. Indeed, the fast performance of tensor cores is an attractive opportunity to explore possible applications that can take advantage of this new technology. However, the use of tensor cores does not come free of restrictions and potential issues. The first restriction is that the programming of tensor cores is much more restrictive than the programming of traditional GPU cores. In fact, the only operation allowed when using tensor cores is the matrixmultiply-accumulate (MMA) in matrices 1 of 4 × 4, i.e., which takes just one GPU cycle. Therefore, the only way to make an algorithm to take advantage of tensor core acceleration is to redesign it as a set of many 4×4×4 MMA operations that can run in parallel.
The second restriction of tensor cores is the numerical precision. Although the resulting D matrix can be stored in FP32 precision, currently the operations A × B + C are done in FP16 precision. This mixed mode of operation may have negative effects in some applications. The reason of why tensor cores operate at FP16 is because in general Deep Learning applications do not suffer from numerical precision when doing training and inference in FP16.
Given the importance of arithmetic reductions in the era of big data and the recent advances in GPU Tensor Cores for fast MMA computations, we formulate the following research question: Can tensor core performance be exploited to accelerate arithmetic reductions and speedup applications from science and technology? This work aims at answering this question by analyzing the potential performance improvement one could encounter by doing arithmetic reductions with tensor cores. For this, we propose a new parallel reduction strategy, based on a hierarchy of MMA operations. Our findings show that in theory a tensor core based reduction can be significantly faster than a traditional GPU one, as it only takes T (n) = 5 log m 2 (n) time steps, with m being the linear size of the MMA matrices involved.
The rest of the manuscript presents the related work (Section II), an overview of the parallel reduction (Section III), the formulation of the new tensor core algorithm and its analysis (Section IV), and a discussion of the results obtained in Section V.

II. RELATED WORK
The parallel reduction has been implemented using different frameworks. In the case of OpenMP there are high level abstract constructs that allow the programmer to express a parallel reduction via OpenMP pragma commands [1]. In the case of GPUs, the parallel reduction has been addressed by Nickolls, Buck and Garland [14]. The authors propose a parallel sum reduction of O(log 2 (n)) time, where each thread loads one element of the input array, and then adds pairs of values in parallel as where p is the number of threads. The loop in this kernel implicitly builds a summation tree over the input elements, and at the end of this loop, the first data slot of each thread-block holds the reduction result of one iteration, i.e., x B [0] = p i=0 x i . More kernels are executed until the input problem is one value.
Harris optimized for CUDA an algorithm for parallel reduction [7] tree-based as well. In his work, the author illustrates seven different optimizations that are relevant to the parallel reduction, achieving a final version that is up to 30x times faster than the initial GPU version presented. Harris mentions that although the time complexity of the parallel reduction is indeed O(log 2 (n)), the cost is not efficient if p = n/2 threads are used. The author shows, with the help of Brent's Theorem, that using p = n/ log 2 (n) threads leads to a parallel efficient cost algorithm.
Another alternative, simpler in design complexity but still competitive in performance, is to use the intrinsic Atomic Add [15] GPU instruction. The operation is atomic in the sense that it is guaranteed to be performed without interference from other threads, but can still run in parallel if threads are using the atomic operation on different data addresses. The atomic add reads a number at some address in global or shared memory, adds another number to it, and writes the result back to the same address, in one transaction.
In the case of distributed computing, there are two levels of parallelism that must occur for the reduction to be completed; (1) local reduction and (2) distributed reduction. For the local reduction, the process may be carried with multi-core CPU or GPU computation as recently described. For the case of distributed computation, the results of different compute nodes must be merged with message passing tools such as MPI [17]. The result is an hybrid OpenMP-MPI or GPU-MPI reduction for massive scale systems. Recent tools such as MapReduce [6] also offer a higher-level of abstraction to accomplish parallel reduction in cluster environments.
The most recent and relevant work about CUDA GPU tensor core programming is the one by Markidis et al. [10], in which they studied current approaches to program NVIDIA Tensor Cores, as well as the performance and the precision loss due to computation in mixed precision. The authors show how NVIDIA CUDA provided three ways of programming the matrix-multiply-accumulate (MMA): CUDA Warp MMA (WMMA) API, CUTLASS, and cuBLAS GEMM. The tensor core programming is analyzed in different aspects such as programmability, performance and precision. The authors report that the maximum performance obtained was with the cuBLAS GEMM implementation, where they achieved 83 T F LOP S in their test environment (approximately 74% of the theoretical performance), followed by CUTLASS with 62 T f lops/s. The WWMA implementation did not provide any performance improvement, however the authors realized they did not use shared memory in the process, which is an important aspect in order to accomplish efficient tensor core computation. They also observed that when the size of the input matrix increases, the error may increase significantly in some cases. For this, the authors include mechanisms to increase precision such as the Kahan summation and iterative precision refinement.
In more general terms, the Google Tensor Processing Unit (TPU), deployed for data-centers in 2015, is another processor comparable to the tensor cores found in Nvidia GPUs, that accelerates the inference phase of neural networks [9]. Norman et al. compared the TPU to a server-class Intel Haswell CPU and an Nvidia K80 GPU. They used workloads written with the TensorFlow framework. The results showed that the TPU was on average about 15× -30× faster than its contemporary GPU or CPU, and about 30× -80× higher in TOPS/Watt (Tensor Operations per Second per Watt).

III. OVERVIEW OF THE CLASSIC PARALLEL REDUCTION
Given an array of n elements, X = {x 1 , x 2 , ..., x n }, the arithmetic reduction R(X) is defined as On a serial processor, one would write a simple loop with a single accumulator variable to construct the sum of all elements. Such algorithm is the best sequential one (as in principle no element can be left out) and it costs Θ(n). A parallel reduction algorithm cannot solve the problem in one time step, as there are dependencies and race conditions to be fulfilled within the reduction process. Nevertheless, it is known that a parallel reduction can sum pairs of values in parallel at each time step, leading to a parallel cost of T p=n/2 (n) = O(log 2 (n)) using p = n/2 processors. At each time step k, we have a problem of size n 2 k−1 and n 2 k threads prepared to do work. Each parallel thread sums a pair of values, i.e, the i-th thread sums the elements x i and x i+n/2 k and stores the partial sum in X[i]. Such parallel step costs O(1) time and cuts the problem in half. Consecutive applications of this process lead to the recurrence subject to T (2) = O(1), which by the master theorem it can be shown that it is T p=n/2 (n) = O(log 2 (n)).
When considering the parallel cost, which is defined as C p = T p (n) · p, with p the number of processors employed, one can realize that using n/2 processors leads to a parallel cost of C p (n) = log 2 (n) which makes the algorithm cost inefficient as it is greater than the O(n) cost of the sequential algorithm. To improve on this cost, one can use Brent's Theorem [2] that states the following inequality where in the case of the arithmetic reduction we have that T 1 (n) = O(n) and T ∞ (n) = log 2 (n). With this, one can choose the number of processors as p = n log 2 (n) without sacrificing parallel time, i.e, With this change, the cost of the parallel algorithm is now efficient, i.e, It is important to consider that there is an asymptotic lower bound of Ω(log 2 (n)) for the parallel reduction. Nonetheless, one can still improve in the constants involved as in the base of the logarithm, which can make an important difference in experimental performance. Such is the case of reducing with tensor cores, where a large number of arithmetic operations can be encoded into MMA operations and executed in just one GPU cycle, which is equivalent to one time unit.

IV. A NEW REDUCTION ALGORITHM BASED ON MMA
In this section we present the new algorithm for parallel arithmetic reductions using tensor cores operations and analyze its parallel time in both asymptotic and non-asymptotic forms.

A. Formulation
The tensor core programming model exposes a single operation to the programmer, the matrix-multiply-accumulate (MMA). That is, given three matrices A, B, C, the MMA operation computes In one GPU cycle. The tensor core computing model allows many MMA operations to occur simultaneously in parallel.
It is interesting to note that in the programming model the tensor core MMA operation is exposed in terms of 16 × 16 matrices to the programmer, even when the actual operation at hardware level is carried in terms of 4×4 matrices. The process of splitting the 16 × 16 workload into smaller 4 × 4 works is done automatically by the GPU scheduler, but splitting a large problem of size n into several 16×16 matrices is not automatic and the partition must be designed manually. This last aspect is the one important for the research, as it is related to the research question of wether a reduction problem of size n can be encoded into multiple MMA operations. The presentation of the new reduction algorithm will proceed in terms of m × m MMA matrices, as it favors the analysis in the next subsection.
The main intuition behind a tensor core MMA based reduction is to produce many partial summations of groups of m 2 numbers, in parallel. To achieve this, we employ two MMA operations. For the first MMA operation, m 2 elements of the input array X[ ] are inserted in A m×m in parallel such that A m,m is actually the m 2 -th element of the group. Then we set B m×m as an all-ones matrix, also in parallel, and C is a zero-matrix, also set in parallel. When the MMA operation is executed on A, B and C, the result is a matrix, namely D m×m = A m×m × B m×m + C m×m , of the form where each column k,j holds the whole set of partial summations from the group of m 2 elements. The second MMA operation is in charge of reducing the partial summations found in the columns of D, now into a single value. It is relevant to notice that dealing with one column of D is sufficient has it holds all the partial summations. The other columns of D hold copies of the result. At the same time, given the rigidness of the MMA operation, it takes less time to process the whole matrix within the MMA operation instead of filtering or doing extra considerations to just process a single column. As long as the result of one column is not compromised, doing a full MMA operation is still the most convenient approach in this contexts as we preserve the one GPU cycle cost per MMA.
The second MMA operation proceeds by changing the order of the multiplying matrices, and re-uses the output matrix D in the position of B, while using B in the position of A. With these changes, the second MMA operation becomes the following expression The resulting matrix D ′ contains the reduction of the m 2 numbers, replicated in all of its elements. Once the second MMA operation finishes, the thread in charge proceeds to write the reduction result, e.g, D ′ 1,1 , into global memory as part of the new array of partial sums.
The global idea of the algorithm is to subdivide the domain of n numbers into n m 2 blocks of m 2 elements and execute the 2-step MMA reduction proposed for each group in parallel. This reduces the problem size by a factor of m 2 . In the next iteration, the idea is to take the smaller version of the problem, of size n ′ = n m 2 and reduce it again with the 2-step MMA operations for all the possible groups of m 2 elements. This process is carried iteratively until the reduction set fits in just one group of m 2 elements, for which a final tensor core reduction returns the result of the whole reduction problem of size n.
This new tensor core MMA based reduction, namely R tc (X), with X = {x 1 , . . . , x n }, can be described by the following recurrence where M (...) is the tensor core based MMA reduction and the initial condition is defined as In the following subsection we perform a fine analysis (considering constants) of the performance of the tensor core based reduction algorithm.

B. Analysis of Performance
In order to analyze the cost of the tensor core based reduction algorithm, we utilize a simplified GPU Computing model similar to the PRAM but with extra restrictions related to the GPU architecture. In this simplified model, the costs for the different types of parallel operations are: In the case of a GPU reduction, it is possible to produce coalesced memory acceses, by not using a stride in the access pattern, and instead make threads access the array X in blocks without gaps among them. Once threads read their corresponding data, which takes one unit of time, they can load this data into the tensor core matrices, which reside in cache memory near the tensor cores. This loading of information takes another unit of time. Then, the algorithm proceeds and executes the 2-step tensor core MMA reduction simultaneously for all the m 2 groups that can be made for the array X[1..n]. Lastly, once the reduction of a group is done, the thread in charge writes this result in its corresponding location in X, taking another unit of time. Combining these costs into one expression leads to the following recurrence that stops at T tc (m 2 ) = 5. Solving the recurrence leads to the final cost of T tc (n) = 5 log m 2 (n) For comparison, a typical GPU reduction, under the same cost model, would take unit of time for reading the first element of each thread, another unit of time for reading the second element, a unit of time for adding two numbers, and another unit of time for storing the value back to memory. The recurrence for such algorithm would be T (n) = 4 + T (n/2) which solves into T (n) = 4 log 2 (n). Based on this cost model, the potential speedup for the tensor core based reduction would become S = T (n) T tc (n) = 4 log 2 (n) 5 log m 2 (n) = 4 5 log 2 (m 2 ).
A value of m ≥ 2, which is the minimum value m could take in order for the reduction to work, is already sufficient to provide a result of S > 1.

V. DISCUSSION AND CONCLUSIONS
The main contribution of this work is the presentation and analysis of a new tensor core based reduction based on matrix-multiply-accumulate operations. The intuition behind the algorithm lies in defining each iteration of the parallel reduction as a sequence of two MMA operations that act together to provide the summation for groups of m 2 values. In the first MMA operation the elements are inserted in the A matrix, while B is all-one and C is a zero matrix. In the second MMA operation, A is an all-one matrix, B is the result of the previous MMA, and C is zero again.
The main results obtained from the analysis are the time cost of the algorithm, which is T tc (n) = 5 log m 2 (n), and the speedup with respect to a classic parallel GPU reduction, in which we obtained a factor of S = (4/5) log 2 (m 2 ). It is interesting to note that the minimum value that makes the reduction work, i.e. m = 2, already makes S > 1 and leads to a favorable result. Current GPUs perform 4 × 4 × 4 MMA operations in one cycle at hardware level, therefore there is an important chance that the potential speedup one could observe with current GPUs, such as Tesla V100 or Titan V, would be S ≈ 3.2 if we consider m = 4, and S ≈ 6.4 if we consider m = 16 which is the value exposed to the programmer. In either case, the performance speedup is significant and would contribute greatly to the problem of parallel reductions.
The results obtained in this work have shown that there is an important amount of potential performance that can be exploited from tensor cores, with applications in fields that are not necessarily related to machine learning. Indeed, it remains unknown what is the level of precision loss by performing reductions in FP16, and future work should include experimental results validating both the speedup values and the % of precision loss observed in different types of reductions.


Online Packaging Method of Filament Cathode for Intelligent Manufacturing

It is challenging to automatically package spiral shaped filament cathodes that are widely used in worldwide microwave ovens. The difficulties include a fast packaging speed of over one parts per second and the dexterity required to handle tangling and embedding between spiral filaments’ multiple circles. In order to convert conventional labor intensive packaging lines of filament cathodes into automatic ones, an online fast automatic packaging method is reported. The method converts the fast incoming sequential flow of filament spirals into a batch of twenty-parallel queues with two stopping levers. The proposed method using two levers to release only the front twenty filaments at a time, one from each queue to fill one complete row of the plastic packaging tray. In this way, it speeds up packaging speed twenty times over conventional sequential packaging method. The inner diameter of the queue tube matches the outer diameter of the filament. It is this dimeter matching to eliminate the embedding problem occurred between neighboring filaments. The proposed method uses vibrating brush to decouple tangled filaments. Large number of validating experiments indicate that the proposed method could automatically package 200 filaments spirals within 2 minutes, with positional accuracy of over 97%. Thus, the proposed packaging method could meet the manufacturing requirement of filament cathodes, releasing four labors out of every production line at every work-shift. Further, this method also simultaneously improves the surface quality of the filaments by eliminating scratch caused by their sharp fresh-cutting cross sections.

Introduction
Cathode magnetron is widely used in both military applications 1,2 as well as household microwaves 3,4 . Cathode magnetron or more specifically its cathode filament is typically made of thorium tungsten or barium-tungsten alloy. When excited with a high voltage, the filament of a magnetron cathode within a microwave oven would generate Tera Hertz frequency microwaves causing water molecules of food to move, vibrate, and bump into other food molecules. In this way, the filament inside a cathode converts electromagnetic energy into heat, which is quickly absorbed by the food. Therefore, the filament within a cathode magnetron is the core component of a microwave oven. The chemical composition of the cathode filament determines the escape power, emission stability, and lifespan of a microwave oven. However, when the chemical composition is optimized and fixed, the manufacturing quality of the cathode filament determines the level of noise, stability, efficiency, and lifetime of a microwave oven 3 . The filament cathode within a microwave oven is usually spring-spiral or spiral shaped. Over three quarters of worldwide microwave-oven used filament cathodes are made in China. This is because this country is rich in rare earth raw materials of tungsten and thorium. The factories find it is challenging to automatically package the spiral shaped filaments. Existing spring spiral packaging method does not satisfy the speed requirement of online packaging of filament cathodes 5 .
The manufacturing of filament spirals includes multiple processes. The first step is to pull the molten thorium tungsten material into a long-wire. Secondly, bend and wound the long-wire into a spiral tube. Thirdly, cut the spiral tube into multiple filament spirals. Fourthly, packaging the filament spirals intro trays. Then the fifth or the last step is quality inspection under a microscope. The typical length of a filament spiral is about 25mm, with outer-circle diameter of about 5mm. The spiral contains about 10 circles of winding filament, with pitch size between two adjacent circles is about 2mm. Since the spirals are cut out of a long-winded tube line, its both ends may either have a roundor sharp-shaped cutting cross-section, depending on the cutting tool. With the advancement of automatic manufacturing of various products, the first three steps of the filament spiral processing lines have been successfully converted from labor-intensive type into semi-automatic ones. With at most one worker per machine. For example, with the spiral cutting process, one worker stands beside each automatic cutting machine to continuously feed it with spiral-like wires. However, for the packaging process, filament spirals produced by each single cutting machine need four workers to do the packaging, since it is challenging to make this process automatic. One problem comes from the fact that they contain multiple circles, causing them enbedding and tangling with each other. Fig. 1 shows some spiral-shaped thorium tungsten filaments used in worldwide microwave ovens. Fig. 1 also illustrates the ways how the spiral-shaped filaments tangling and embedding into each other. This makes it is challenging to automatically separate them into different pits of a packaging tray. Separating them needs dexterity of human fingers with careful eye-hand coordinating. Another problem is the speed of filament coming from the manufacturing line at over one part with a second. Existing The most important goals pursued by the tungsten molybdenum deep processing enterprises is the delivery time, quality, and manufacturing cost to guarantee a reasonable profit. In order to ensure that such premises are fulfilled under the situation of sharply increasing labor cost and quickly dropping working forces, it is necessary that these companies have automatic packaging machines as wells intelligent quality inspection system to release manual operators for more attractive and personalized tasks. It is possible that well-designed automatic packaging machine may simultaneously reach higher production quantity, reduce the number of workers and increase their motivation and productivity, and improve the quality of the products. In this work, a novel automatic packaging method is designed to fulfill all these goals. In a related paper, an automatic quality inspection method is reported separately.

Design Requirements for Automatic Packaging
The challenge of the packaging process of filament spiral comes from multiple aspects. First, the speed of the cutting process is quick, which is over 60 spirals (typically 68 spirals) over every minute. This is far too fast for any human to pick it up and place it orderly in a tray. Hence, the spirals are piled up in a plastic basket (Fig. 1, right column). This causes a second challenge, i.e., the piling of the filament spiral cause in embedding and tangling of the spirals into each other, making it difficult and time-consuming to separate them. Conventionally, filament spirals produced by one cutting machine usually need about four workers to manually place them in a tidy manner into a plastic tray (Fig.2). Further, the piling of the filament spiral may cause the surface scratched by the sharp-cutting cross section, which have to be removed away from the packaging process. In order to automatic this packaging process with higher efficiency and better quality, the design requirements of the packaging devices include: ① Fast, about to package on average over 60 filament spirals per minute. Use the same existing disposable plastic tray that is conventionally used in the manual packing process, each tray shall contain 200 filament spirals, with 10 rows by 20 columns. The length of the tray is about 25cm, and its width is about 21.5cm. This means about 3 minutes per tray. Leaving sufficient time to move full tray away from and empty tray into the output port of the cutting machine, this means about 2 minutes to package a full tray.
② There may exist unexpected interrupts of the cutting process due to many possible reasons. The major challenge to automatic the packaging process is the required speed. Conventional method using motorized XY stage to move the tray relative to the cutting machine. However, such approach would cause a dilemma between a high moving speed of the tray and accurate placement of the filament spiral into a vacant spot. The reasoning is described as below.
We first theoretically analyze the feasibility of using a conventional XY motorized stage. To sequentially fill in the 200 vacant spots of a tray, the tray needs to travel a distance of 250mmⅹ20 columns +210mm within 2 minutes. The means a maximum speed of about 41.8mm/s if the tray is moving at a constant speed without stop nor acceleration/deacceleration. This speed is far too high to allow the filament spiral to fall into the vacant spot without bumping out of it. The bumping force exerts on the spiral would be 41.8mm/sⅹ1gram)/δt. The contacting time is very short, and hard to measure. However, it is reasonable to estimate that it is of the order of microseconds. Hence, the bumping force would cause the spiral bump and jump out of the tray randomly.
To reduce the possibility of the spiral bumping, it is necessary to reduce the speed at the instant when the spiral touches the tray down to near zero. Typical processes to place a number of spirals into the vacancies of the tray involve the following continuous steps. First, move the tray horizontally so one of its vacant spot is just placed underneath the output port of the cutting machine at the time before the spiral is ready to drop into it. Then deaccelerate to let the spiral drop into the spot. Then accelerate and move the tray again so its next vacant spot is placed underneath the output port of the cutting machine, wait for the next spiral to come. This is repeated ten times to feed the whole row. After one row finishes, move the tray vertically to place the second vacant rows into place. Then continuously feed the second row. However, this approach would not work considering the limited acceleration and de-acceleration of any two-dimensional motorized stage, and the vibration caused by periodic acceleration and de-acceleration of the mechanical system. The feasible packaging speed using the XY motorized stage is typically 22 minutes per tray, which is about 10 times slower than the requirement of 2 minutes per tray.
Hence, a sequential packaging approach is proposed to solve the problem of filling in the 200 vacant spots of the disposable plastic tray within 2 minutes.

Design of the Packaging Device
The idea is to do parallel packaging by decoupling the tray's XY two dimensional movements into two one-dimensional movements of different objects. The major componets of the designed parallel packaging device include a flow redistribution disk, parallel filling device, and a tray grabber ( Figure   3). The designed flow redistribution disk is shown in Figure 44a. Its major componets include one fixed entry cup, one rotational arm, and 20 parallel output paths. The entry cup is beneath the output port of the cutting machine. Fast sequential flow of filament spirals coming at a speed of about 70 spirals per minutes from the cutting machine enters the disk via the entry cup. The entry cup shall not be flat, instead, it should be like a deep bowl shape to allow bumpped filament spiral fall back into the cup. Its material should be soft plastics to extend the contacting time between the falling spiral and the cup, in order to reduce the bumping force used to reduce the momentum of the filament down to zero. The rotational arm then bridge the filament to one of the 20 parallel paths. Special attention is needed to guarantee roughlly equal amount of filament move to each of the 20 parallel-paths. This is done by an infrared sensor to detect whether or not an actual filament is moved to a path. The 20 parallel output paths of the redistribution disk are connected to a total of 20 queuing tubes of a parallel-filling device (Figure 4b).
The parallel filling device consists of twenty queueing tubes, a releasing mechanism, and a motorized stage on which the tray is placed (Figure 4b)). The queuing tubes save as a reservoir of the filament spirals. Each queue may hold about 20 filament spirals, all together, the queues may hold upto two full trays of spirals. It is the parallel queuing tubes that decouple the manufacturing speed and the packaging speed. This is important because either procedure may be interrupt at a random time and it is impossible to synchronize them at the the time span of a single spiral duration. The diameter of the parallel queuing tubes is slighly bigger than the outter diameter of the filament spiral such that at any section of the tube, only one filament spiral is there. In such way the queuing tube avoid circles from one spiral enbed into its adjacent neighbors's. However, it is still possible that the front end of the newly arriving spiral may tangle the terminal end of the last spiral waiting in the queue. The releasing mechanism is rather complex. It contains two motorized stopping levers A and B, placed at a distance of about 25mm away. This distance is adjustable according to the length of the filament to be packaged. The releasing mechanism is placed between the exit end of the parallel queuing tubes and the plastic tray, with the lever B closer to the exit. When the tray is properly placed under the exit of the parallel queuing tube, then a message is sent by the control unit of the packaging system to signal the releasing mechanism to lift up its motorized stopping lever B and keep its lever A down. So when the stopping lever B is lifted, only the first filament spiral of each queuing tubes could be released, blocking the second filament spiral in each queue by lever A. Hence all spirals entering each queue after its second spiral will also be blocked by lever A. Each time when lever B is lifted up, one filament will slide from each of the twenty ques into one vancant spot of the plastic tray, all together twenty queus will fill in the full twenty-collumns in each row of the plastic tray. In such a way, the parallel packaging is achieved.
Once one row of the plastic tray is filled, the motorized stage is controled to move one step forward, then stays there staticaly. The step size of the stage is the distance between two raws of the tray. While the tray is carried by the stage to move forward to place its next raw in position to wait for the next batch of 20 filaments, the releasing mechanism is controlled to put down its lever B in order to close the door of all queues and to simultaneously lift up its lever A. This results in the full stream of filaments within a queue to move one step forward to lever B. The step size of the stream is the length of the filament. Once the forward stepping of the stream of filament is completed, then lever A is controled to put down to block the second filament and every one behind it, while lever B is lifted up so a batch of 20 filaments slide parallely into the raw of the tray, and so on till 10 raws are filled.
One trick worth to mention here is that when lever B is up, it is possible some of the queues will not let its front filament go freely. This may due to the possibility that two adjacent filaments tangle together by there spiral ends. A servior tangling case is shown in the right column of Figure 1, while the degree of tangling phenomenon within the reservoir is much lighter. Howerver, potential tangleing between adjacent filament will reduce the robusty of the packaging system. To solve this, one back-and-forth vibrating brush is designed and placed just between levers A and B. Each time after the lever B is lifted up, the brush will vibrating clockwise and anticlockwise to ensure all front filament be released. Once the 10 raws of a plastic tray are all filled ( Figure 5), the motorized stage will first carry the tray to slightly vibrate back and forth to cause filaments that are still on its plateu to fall into the nearby valley. Then carries it back away from the queuing tubes and enter the automatic tray grabbing system ( Figure 6). Then vacumed tray grabber will succk it up, move it, and place it on the top of a pile of trays that are full of spirals. The tray grabber will then suck an empty try and carry it to the now empty motoried stage. After receive the empty tray, the motorized stage will then move forward and place the first row of the empty tray underneath the exit of the queuing tubes. Then start the procedure of packaging the next tray.
It takes roughly from half to one minute for the automatic vacummed grabber to replace the full-packaged tray by an empty one. During this process, the cutting machine may keep on working without being interruptted, even though no tray is receiving its output filament spirals. All parts it produces enter the queuing tubes, waiting there till the tray is ready for them to flow in a batched manner into the tray.

Evaluating Results
The proposed parallel packaging system was prototyped and tested in a filament spiral-manufacturing factory. Twenty trays were packaged with the prototyped device. The experimental results are shown in table 1 below. The time to package one tray is 2 minutes on average, satisfy the design requirement, and matches the manufacturing speed of the cutting machine. The total number of filament spirals put into each tray is 200, thus the accuracy of packaging number is 100% correct. Though there were still occasionally jumped spirals that were not in the expected pits of the tray, the achieved positioning accuracies were all above 97%.

Summary and Conclusion
An automatic packaging method was proposed and prototyped to simultaneously deal with the problems of high-throughput packaging of spiral shaped filament spirals at a speed of over one part per second, accurate placement in terms of positioning as well as orientation of the spiral spirals, reduce tangling and embedding among filament spirals, and eliminate scratching by their sharp cutting ends. The experimental results indicate that the efficiency of the prototyped system is to package a 200-spiral tray with 2 minutes. Other performance aspects include 100% accuracy concerning the total number of filament spirals per tray, and over 97% position/orientation accuracy. The proposed automatic packaging facilitates the factory to convert its semi-automatic manufacturing lines into full automatic ones, reduce the manual working load of the most labor-intensive packaging process by one order of magnitude. Saving cost while improving its product quality.


extfieldsofstudy:
Physics
Computer Science

text:
Software Implementation of the discrete ordinates method in the problem of the Poiseuille flow in a plane channel with infinite parallel walls

The algorithm for calculating the gas macroparameters in a channel using the discrete ordinates method was considered. The system of finite-difference equations that determine the general solution of a homogeneous integro-differential equation after separation of variables was written in a matrix form. The coefficient matrix in the resulting system of equations was reduced to a symmetric one and then to a tridiagonal matrix using the Householder algorithm. The search for the eigenvalues of the tridiagonal matrix was performed using a QL algorithm with implicit shifts. The determination of the coefficients in the expansion of the solution of the required boundary value problem with allowance for the boundary conditions on the channel walls was reduced to a system of linear algebraic equations. The Gauss method with a permutation of the rows was used to solve it. The software implementation of the algorithm was carried out using the C++ programming language. The mass velocity profile of the gas and the flow rate through the cross section of the channel in the Poiseuille flow problem were calculated as an application. The comparison with similar results presented in the open press was performed.

Introduction
Technologies related to the design, manufacture and operation of micro and nano devices develop actively and require conducting systematic fundamental research of the processes taking place in such devices. Given that existing microelectromechanical systems (MEMS) have dimensions from a micrometer to several millimeters and nanoelectromechanical systems (NEMS) have submicron dimensions, conducting experimental studies is often difficult or even impossible. Therefore, the main direction of research is the use of mathematical modeling methods. The complexity of applying mathematical modeling methods to the description of transport processes in micro and nano devices is due to the fact that in practice most of them operate in a sufficiently wide range of Knudsen numbers. The latter circumstance leads to either using several models of gas flow for various flow regimes with «docking» them with each other, or using of universal approaches capable of describing all flow regimes [1]. Such universal approaches include the modeling of gas flows based on the Boltzmann kinetic equation or model kinetic equations using numerical methods and, in particular, the modified discrete ordinates method, which was implemented in [2] for the case of rarefied gas flows in channels with infinite parallel walls. Considering the high efficiency of the modified discrete ordinates method, it is of interest to generalize it for three-dimensional flows. The aim of the present work is to develop a software implementation of the discrete ordinates method and its application in the problem of the Poiseuille flow in a plane channel with infinite parallel walls.

Statement of the problem and mathematical model
Consider the flow of gas in a channel with walls located on the planes '' xa  in a rectangular Cartesian coordinate system. Assume that the temperature of the channel walls is constant everywhere and equals T , and the motion of the gas is due to a constant pressure gradient ' We direct the axis ' Oz of the Cartesian coordinate system along the pressure gradient. Then the BGK model of the Boltzmann kinetic equation in the chosen coordinate system is written in the form [1] v v ( ) Here ( ', ', ) f x y v -the gas molecules distribution function in coordinates and velocities, v x and v ythe projections of the gas molecules velocities, p и  -the pressure and the coefficient of dynamic viscosity of the gas, ( ', ', ) eq f x y v -the local equilibrium Maxwellian with parameters given on the channel walls. The relative pressure drop over the mean free path of the gas molecules g l is assumed to be small. This allows us to solve the problem in a linearized form. The distribution function may be represented in the form  r -the radius vector of the points of the channel walls,  -the accommodation coefficient of the tangential momentum of the gas molecules by the channel walls, n -the normal vector to the channel walls directed toward the gas [3]. Substituting (2) in (1) and (3) for finding ( ; ) x Z x C we arrive at the boundary-value problem ( The general solution of (4) may be written in the form (6) Substituting (6) into (4) to find ( , ) Yx we arrive at the homogeneous equation (7) Taking this into account, equation (7) can be rewritten in the form (9) Following [2], we pass from (8), (9) to the finite-difference equations in the velocity space Here k  -the weight coefficients of the quadrature formula used to find the value of the integrals in the right-hand sides of (8) and (9). In the present paper 6-point Newton-Cotes quadrature formulas were used to calculate the integrals.
. Similarly to the Case method [3], we seek the solution (10), (11) in the form (12) Here  -spectral parameter. The substitution of (12) into (10) and (11) leads to two matrix equations Here  (13) and (14), we obtain Transforming equation (17), we successively find   Thus, finding the values of the spectral parameter  , from expression (12) reduces to finding the eigenvalues of the system of equations (18). Since the matrix of the system of equations (18) is symmetric, its eigenvalues are real and different. In order to find them we have used the Householder algorithm, which made it possible to move from a symmetric matrix to a tridiagonal one. The search for the eigenvalues of the tridiagonal matrix was carried out using a QL algorithm with implicit shifts. Hence the problem of finding the spectral parameter values is solved. Now we turn to finding the eigenfunctions corresponding to the found spectral parameters. We use the normalization condition Substituting (12) into (10) and (11) with mentioned normalization condition, we find ( , ) Then the solutions of equations (8) and (9) can be written as a linear combination of the constructed solutions (20) A direct substitution shows that the constant A and the function ) ( (21) By virtue of condition (19), the values of the spectral parameter tend to zero for an infinite increase N and, as a consequence, the particular solutions corresponding to them in (20) and do not contribute With this in mind when writing (21), we excluded from (20) particular solutions that correspond to the minimum of the found values of the spectral parameter. Substituting (6), (21) into the boundary conditions (5), we arrive at a system of linear algebraic equations for determining the We have used the Gauss method with a permutation of the rows to solve the system of equations (22), (23). A solution of the system of equations (22), (23) completes the solution of the boundary value problem (4), (5). Proceeding from the statistical meaning of the distribution function [3], we finally find the mass velocity profile of the gas in the channel and then the flow rate through the cross section of the channel Q [1] Substituting the obtained results in (24) and (25), we finally find ) 3

Results and discussion
The software implementation for calculating the values of () qx and Q have been performed using the C ++ programming language. Since the integrand in (24) contains a factor proportional to

Conclusions
The algorithm for calculating the gas macroparameters in a channel using the discrete ordinates method is considered on the example of the Poiseuille flow problem. Software implementation using C ++ programming language is also presented. The BGK model of the Boltzmann kinetic equation is used as the basic equation describing the process kinetics, and Maxwell's diffusive reflection model is adapted as the boundary condition on the channel walls. The mass velocity profile of the gas and the flow rate through the cross section of the channel are calculated as an application.