Johann Cervenka was born in Schwarzach, Austria, in 1968. He studied electrical engineering at the Technische Universität Wien, where he received the degree of Diplomingenieur in 1999. He then joined the Institute for Microelectronics at the Technische Universität Wien and received his PhD in 2004. His scientific interests include three-dimensional mesh generation, as well as algorithms and data structures in computational geometry.
Matrix Classes in C++ for Python
During the simulation of physical problems, frequently huge equation systems have to be assembled and solved afterwards. Fortunately, the resulting equation systems are often sparsely occupied, and with the utilization of specialized storage methods, these systems achieve practicable memory demands, even with up-to-date hardware.
When it comes to selecting the desired programming languages, a trade-off between convenience and execution speed has to be accounted for. There exist several libraries, which handle the solving procedure of the equation systems by interfaces to well-established, even parallel, fast solvers. As the assembly mechanism strongly depends on the problem under examination, flexible high-level programming languages, such as Python, are preferred, which may result in slower execution speeds compared to other languages, such as FORTRAN or C. However, it is often neglected that the complexity of assembling the equation systems is in the same regime as solving these equations. The combination of programming languages enables a wide range of activity. We implemented the Python-C++ bindings via Cppyy and developed a sparse matrix class for C++, which offers an easy access model.
As mentioned previously, basic dense two-dimensional arrays are extremely inefficient for usage in sparse systems in terms of memory limitations and execution times, as a vast amount of unnecessary zero entries have to be processed. There exist several techniques to store sparse entries. During assembly, attention has to be paid to arbitrary access methods, which can be provided by, for example, special (sorted) arrays, trees or hash tables.
Upon examination of these methods, it was found that a trade-off has to be considered where the cached memory access of the machines is concerned, which is best handled by arrays, but these have shortcomings when it comes to the performance of frequent search and insert operations. On the other hand, trees show flexible search and insert mechanisms, but good caching cannot be guaranteed.
The performance of the methods strongly depends on the occupation ratio of the matrices. Unexpectedly, the caching performance using arrays dominates over the linear search performance up to high fill rates. During development, different methods have been evaluated and applied to a Wigner framework, written in Python and interfaced with C++ for assembly and for solving the equation systems, which also require parallel execution due to their large size and, intrinsic to this method, even high fill grades.
The Wigner formalism describes transport processes in phase space, which means for one spatial dimension, two dimensions in phase space (x,k) have to be considered. The shape of a typical matrix during a Wigner simulation is shown in Fig. 1. To first order, a large sparse matrix can be seen, but when you look closer, you can see a great variety of different block shapes. Further development is ongoing to improve performance.
Fig. 1: Plot of the nonzero entries in a matrix of size 190000 x 190000 resulting from a Wigner simulation. The entries are located along the diagonal and show a width of 300 elements.