We present an ongoing research in the area of transforming existing sequential Fortran programs into their parallel equivalents. Our approach is to use rewriting rules technique in order to automate the transformation process. Sequential source code is transformed into parallel code for shared-memory systems, such as multicore processors. Parallelizing and optimizing transformations are formally described as rewriting rules which facilitates their reuse. Using high-level algebraic models allows describing program transformations in a more concise manner. Performance measurements demonstrate high efficiency of obtained parallel programs.
Despite being one of the first programming languages, Fortran is still widely used, in
particular for solving scientific and engineering computation-intensive problems. Its
popularity is due to its relative simplicity and lack of complex facilities (e.g.
pointers), closeness to mathematical description of problem and efficiency of
generated binary code. Another reason for continued use of Fortran is that in more
than 50 years of its existence a vast repository of programs, libraries and routines for
solving different scientific problems has been created. Algorithms implemented in
such programs are still valuable, however there is a need to adapt this legacy code to
new parallel computational platforms. Furthermore, due to size and complexity of
existing code, manual adaptation is not a practical option: there is a need for
automated tools to facilitate conversion of legacy code to modern parallel platforms
[
In this paper we describe an ongoing research on parallelizing Fortran programs using rewriting rules technique. Sequential source code is transformed into parallel code for shared-memory parallel platform (such as multicore processors) using automated transformations. Parallelizing and optimizing transformations are formally described as rewriting rules which facilitates their reuse. Such approach is aimed at two main goals: to improve runtime efficiency of programs and to increase developer’s productivity. We illustrate our approach on two sample programs: a simple Gauss elimination algorithm and an applied problem of calculating electron density from the field of quantum chemistry.
There has been an extensive research in the area of parallelizing existing sequential
code, in particular for multicore architectures. Some approaches require manual code
modification and provide facilities that help a developer express parallelism. Such
approaches include parallel libraries [13], parallel extensions to existing languages
[14] and new parallel languages [16]. Another research direction is interactive
parallelization [11], when a developer manually selects the loops to be parallelized,
and the tool applies transformation automatically (our approach also belongs to this
category). Finally there are numerous approaches to automated parallelization, mostly
implemented as parallelizing compilers [
This paper continues our research on automation of process of designing and
development of efficient parallel programs, started in [
Below we describe our approach in more detail, provide examples of parallelizing transformations and illustrate them with parallelization and evaluation of two programs: small example program (Gauss elimination) and applied quantum chemistry problem (electron density). 2
As in our previous works [
After high-level program model is created, we use parallelizing transformations to implement a parallel version of the program on a given platform. Transformations are represented as rewriting rules and therefore can be applied in automated manner. (Selection of loops that could be transformed is performed manually.) The declarative nature of rewriting rules technique simplifies adding new transformations. Also transformations work with high-level model elements (on any level of abstraction), which means they are language-independent.
Usage of high-level algebraic models also allows proving correctness of the
developed transformations [
To automate program transformations we use the rewriting rules system Termware
[
Here source is a source term (a pattern for match), condition is a condition of rule application, destination is a transformed term, action is additional action that is performed when rule fires. Each of 4 components of a rule can contain variables (denoted as $var), so that rules are more generally applicable. Components condition and action are optional. They can execute any procedural code, in particular use the additional data on the program.
Termware supports a number of evaluation strategies, including TopDown (used in this paper), BottomUp and a possibility to implement additional strategies. Termware system itself doesn’t check that transformation process terminates, however the rules used in this paper are designed in such way that each model element is processed at most once, therefore the transformation process is guaranteed to terminate.
In addition to rewriting system, our tools include parsers and generators for target
languages that perform transformation between source code and low-level (syntax)
program model, which is represented as Termware term. We have previously
developed such tools for C# language [
In this section we describe the process of parallelizing sequential Fortran programs for parallel systems with shared memory, such as multicore processors. We parallelize source code of Fortran programs by replacing suitable loops with parallel loop constructs. To create multithreaded Fortran program we use OpenMP framework [14]. OpenMP PARALLEL DO directives are used to parallelize loops. For simple loops, just addition of such directive can produce quite efficient parallel code. In this case there is additional advantage of keeping transformed parallel code similar to existing sequential code. In more complex cases (when there is data dependency between iterations) there is a need of more significant transformations, such as using OpenMP library subroutines for advanced thread management. In such cases, the transformed source code contains significant changes. However, usage of high-level algebraic models allows describing these changes in concise and understandable form.
We will describe the details of our approach using as an example a Fortran program implementing Gauss elimination algorithm for solving systems of linear algebraic equations. The Fortran source code was transformed into a low-level syntax model using developed parser, then into a high-level algebraic model using Termware patterns. When working with legacy code, we found it useful to apply several levels of patterns. First we used generic linear algebra patterns, such as vector and matrix operations. The obtained algebraic model was language-independent, but still quite detailed. Then we applied patterns specific to the problem in question. In this way we obtained schematic representation of algorithm useful for its understanding and deciding where parallelizing transformations should be applied.
The high-level model of relevant fragment of program has the following form: DoCnt(K,1,N-1,
We will parallelize only two of the operators present in program, namely FindMaxElement and UpdateElements. Other operators have less computational complexity, therefore their parallelization is less effective.
Out of two operators, the simplest is UpdateElements, responsible for calculating new values for elements of submatrix:
UpdateElements = DoCnt(I,K+1,N, Assign(S,A(I,K)),
Here, DoCnt denotes common DO loop with counter. The iterations of the outer loop are independent, so this fragment is easily parallelized. We use the following rewriting rule:
The loop to be transformed is marked with _MARK_Parallel symbol to enable rule application. ParallelDoCnt operator is high-level model element responsible for parallel loop. In particular, for OpenMP platform it is transformed into OmpParallelDo operator that describes OpenMP directive represented in Fortran as a pair of special comments: !$OMP PARALLEL DO … !$OMP END PARALLEL DO.
Notice that for C language the same operator is represented as a single pragma statement: #pragma omp parallel for. Therefore, using multiple levels of patterns allows us to provide operators that are common for given platforms, use these generic operators in most rewriting rules and then specialize them only when transforming program model back into source code.
While UpdateElements operator can be parallelized by simple application of OpenMP directive, the other operator FindMaxElement is more complex. It also has the form of loop, but iterations of the loop update the same set of variables (value of the maximum element in submatrix and its indices). This is the case of reduction, when some local values are calculated on each iteration and then merged into one global value. OpenMP supports such cases with REDUCTION clause, however only a set of predefined reduction operators are supported: while finding just maximum value can be accomplished using OpenMP directives, finding maximum value and indices where it occurs is not directly supported.
Therefore we need to provide transformations that parallelize the loop in general case of reduction. We represent FindMaxElement as following combination: FindMaxElement=FindMaxElLoc1*…*FindMaxElLocTN*FindMaxEl
On each thread we execute local version of operator (FindMaxElLoc1,…, FindMaxElementLocTN), and then execute reduction operator FindMaxElReduct that combines local values into one global value.
Both already described parallelizing transformations are aimed at high-level
structure of algorithm. However, as we observed in [
In the Gaussian elimination program we have observed the same effect. We noticed that for certain sizes of input matrix (N=256*M) there was a sudden increase of execution time. We attribute this increase to the peculiarities of memory access: namely, caching adjacent matrix elements. For such matrix size, the adjacent matrix elements were put into the same cache items, therefore increasing the number of cache misses and greatly reducing overall performance. To overcome this peculiarity, we declare the matrix size as N+1 instead of N. The extra elements are not used in calculations, but they change location of elements and improve efficiency of memory access. The transformation is implemented with the following rules:
-> Declaration(A,Array(Double,[MN,MN+1]))
-> [Parameter(N,Integer,In):[Parameter(MN,Integer,In): $next]] 6. Call($name,[N:[A:$next]])
-> Call($name,[N:[MN:[A:$next]]])
The rule 1 adds new parameter, MN, denoting declared matrix size. The rule 2 specifies for which values of matrix size the transformation should be applied. The rule 3 modifies matrix declaration to use new size MN instead of N. Rules 4-6 propagate new parameter to all procedures, procedure parameters and procedure calls.
Notice that rules 4-6 are applied multiple times in a single program: for each procedure definition (rules 4-5) and for each procedure call (rule 6). One of the advantages of rewriting rules technique is that single rule can describe changes in multiple places, reducing effort to make the changes and preventing mistakes possible when applying such changes manually. Notice also that rules 1-6 work on lower level of abstraction compared with previously described rules. The ability to describe transformations on different model levels is another advantage of proposed approach and it allows describing different types of transformations with the same tools.
To evaluate effects of developed transformations, we have measured the performance of different versions of initial program of Gauss elimination. We have compared performance of 4 versions: initial sequential program (SEQ), parallel program with UpdateElements operator parallelized (PAR1), parallel program with both UpdateElements and FindMaxElement operators parallelized (PAR2), and program with both operators parallelized and memory optimization applied (MEM). The measurements were performed on 4-core parallel system, for matrix sizes from 256 to 2048. Obtained speedup (compared with SEQ program) is shown on fig. 1.
As can be seen from the diagram, all transformations result in some performance increase, although their effect differs. For small matrix sizes, both PAR1 and PAR2 show some noticeable speedup, while MEM is not very effective and is very close to PAR2. However, for larger matrix sizes (N>1024), the situation changes. PAR1 and PAR2 become less efficient, close to SEQ. However, MEM becomes much more efficient and demonstrates speedup of more than 10x. Therefore both high-level transformations of algorithms and taking care of low-level implementation details is necessary to obtain efficient parallel programs. Measurement results also demonstrate complex dependency of execution time on real parallel systems, as compared to ideal theoretical models that suggest simple O(N3) dependency.
After developing our tools on sample problem (Gauss elimination) we have tried them on real-world program in area of quantum chemistry. The program calculates electron and spin density in atoms of polycyclic aromatic hydrocarbons on N*N grid [12]. The size of the program is 1680 lines of Fortran code. Source code is not well structured – actual calculations are mixed with I/O operations, debug code and some hardcoded data. Also it contains mix of features from different versions of language – from Fortran 77 to Fortran 95. Therefore usage of high-level algebraic models helped us to understand this legacy code and apply parallelizing transformations in most efficient way.
We were able to reuse parallelizing transformations developed for Gauss elimination program also in electron density program. Only the first, most simple loop transformation was applied. However, the challenge was to select the most suitable loop for this transformation, as the program contained 54 loops and trying all of them was not a feasible option. We have used a profiler tool, Intel VTune Amplifier [10], to find hotspots in source code. Then we applied rewriting rules technique to detect all loops enclosing such code fragments. Thus the number of candidate loops was significantly reduced from 54 to 6. Out of these 6 loop, we applied transformation to second outermost loop (as the outermost loop contained too few iterations, and parallelizing inner loops was less efficient because of repeated cost of creating and synchronizing threads each time inner loop was executed).
We have compared execution time of initial sequential program (SEQ) and parallelized program (PAR) for grid dimensions N from 200 to 800 (see fig. 2).
Applying transformations has resulted in quite significant speedup – from 3.3X to 3.6X (depending on problem size) on 4-core system.
In this paper we have described our approach for parallelizing Fortran programs by applying formalized program transformations to existing sequential Fortran code. Using rewriting rules technique automates application of transformations and prevents mistakes that can appear when applying changes to source code manually. High-level algebraic models simplify understanding of legacy programs and their transformations, and enable transformation on different levels of abstraction. We have applied our approach both to simple program implementing Gauss elimination algorithm and real-world quantum chemistry problem (calculating electron density). Performance measurements demonstrate significant speedup for both programs.
Further research directions include development of the same approach for transforming legacy Fortran applications to target distributed-memory systems and GPUs. Our future plans also include extension to Grid and cloud-based platforms. Also we are planning to improve support for large and complex Fortran programs, in particular automate selection of most suitable place of application for transformations.