Mathematical modeling and the related computer experiment are now one of the main means of studying objects, processes and phenomena of various nature: in science, engineering, economics, society, etc. A significant improvement in the quality of mathematical modeling in many areas of science and engineering is possible only with the use of fundamentally new three-dimensional models, the transition from computer simulation of individual components and assemblies to the calculation and optimization of the product as a whole. It is obvious that the consideration of problems in such a formulation leads to discrete mathematical models of super-large sizes. Existing supercomputers of different parallel architectures make it possible to efficiently solve such problems. However, the time for solving problems on parallel computers consists of the time of the actual solution and the time of performing additional operations, that are necessary for the exchange of information between computing devices, that is overhead costs. This is especially true for problems of linear algebra with different structures of sparse matrices of large volumes, that arise in the mathematical modeling of processes. Sparse matrix compaction schemes, decomposition of data arrays between processors are one of the main factors for the effective solution of these problems on parallel computers. The paper considers efficient methods for processing sparse matrices of arbitrary structure for the purpose of effective mathematical modeling of structural strength problems on parallel computers. Various methods of regularization and decomposition of sparse matrices of arbitrary structure, efficient data storage schemes, technology for studying the conditionality of a matrix with approximate data on a computer are proposed. This way of using sparse matrices in mathematical modeling ensures more efficient use of computing resources and reliability of computer results. Problems of mathematical modeling are presented, where the considered methods of processing sparse matrices were effectively applied. Mathematical modeling, structural strength problems, sparse matrices, linear algebra algorithms, methods of structuring sparse matrices
Nowadays, one of the main directions of scientific and technical progress is the mathematical modeling of complex objects, processes, and phenomena, which is the basis of modern applied developments in various fields of science and technology. With the advent and constant development of supercomputers of various architectures, the problems of their effective use, the expansion of the range of tasks that need to be solved, and the assurance of the reliability of computer results become very relevant. As a practice and literary sources show, mathematical modeling of 80-85% of all scientific and technical problems is reduced to solving problems of linear algebra - systems of linear algebraic equations (SLAE) and algebraic eigenvalue problem (AEVP) for sparse matrices of various structures. The efficiency of solving the entire problem of mathematical modeling largely depends on the efficiency of solving these problems.
Despite the fact that there is a huge amount of software for solving linear algebra problems on various computer architectures on the international market, new applied problems are appearing that require new efficient algorithms. This is primarily due to the difference in the types of matrices (dense and sparse), structural properties (band, arbitrary structure, triangular and others), mathematical properties (singular and non-singular, positive definite, sign-definite etc.), and ways of processing them on different computer architectures.
When designing objects for various purposes, it often becomes necessary to carry out calculations of the strength of structures of the whole objects and their individual elements. Such problems arise, in particular, in various branches of mechanical engineering (ship building, aircraft building, rocket building, engine building, etc.), in industrial and civil engineering (calculation of individual structural elements or structures as a whole), etc.
Increasing requirements for the quality of design solutions, the use of new constructive materials lead to the development of refined three-dimensional mathematical models with huge amounts of processed information. Despite the availability of powerful supercomputers of different architectures, there is a need to use efficient methods of computer data processing, the volume of which is constantly growing. 2.1.
Mathematical problems of calculating the strength of structures using the principle of possible
displacements are reduced to the following variational problems [
It is necessary to find a vector function that integral identity for a static problem for dynamic problem
a ( u , v )+b ¿ u (t0 )=u(0) , u ( {t} rsub {0} )= {u} ^ {(1)} for the problem of eigen oscillations a ( u , v )=l ( f , v ) ; a(u , v )= λb(u , v ) ; (1) (2) (3) (4) where U 0 – is the infinite-dimensional functional space of possible displacements; symmetrical bilinear forms a ( u , v ), b ¿, c (u ' , v ) symmetrical bilinear forms are proportional in accordance with the potential, kinetic energy of deformation and braking energy, and the linear form l ( f , v ) is proportional to the work of the applied (external) forces under load. The first time derivative of the vector function u (t , χ ) is denoted u ', and the second u .
Only linear problems are considered here, since it is assumed that the solution of a nonlinear problem can be reduced to solving a sequence of linear problems.
The theoretical basis of most software tools for calculating the strength of structures, for
example, the SP LIRA software package designed to solve problems of calculating the strength of
building structures [
To obtain a discrete problem with the help of the FEM, the area occupied by the structure is distributed into finite elements, nodes and their degrees of freedom (displacement and rotation angles of nodes) are assigned. The degrees of freedom correspond to the basic (coordinate, approximating) vector-functions ϕi, of nonzeros only on elements containing a node corresponding to a given degree of freedom. In addition, the following relations hold for the degrees of freedom and basis functions (5)
L j( ϕi )=δij, where
δij – is the Kronecker symbol, and the result of the operation L j( ϕi ) is the value of the component component of the vector-function ϕi for the degrees of freedom L j.
Approximate solutions of the corresponding problems are sought in the finite-dimensional subspace U 0h of space U 0. Vector-functions from subspace Uh0 are piecewise polynomial and can be represented as a linear combination of basis vector-functions, which satisfies the main (kinetic) conditions for a linear combination of basis vector-functions
n (6) uh ( χ )=∑ xi ϕi( χ ), i=1 h where ϕi (i=1,2 ,... , n) – the piecewise polynomial basis mentioned above U 0.
Then the discrete problems, that should be solved have the form [
Ax=b; problems with initial conditions for dynamic problem (2), (3)
Bx ' ' ( t )+Cx ' (t )+ Ax (t )=b (t ), x (t 0)=x(0), x ' (t 0)=x(1); algebraic eigenvalue problem for the problem of eigen oscillations (4)
Ax= λh Bx . (7) (8) (9)
To investigate the sustainability of composite materials, the static method of the three-dimensional
linearized theory of stability is used within the framework of the second version of the theory of small
subcritical deformations. The subcritical state is determined from the solution of the plane problem of
the linear theory of elasticity of piecewise homogeneous bodies for various boundary conditions,
which correspond to a three-layer element of periodicity in the case of a material and a three-layer
element of a structure with lateral sides free from stresses in the case of a composite sample [
When using the static Euler method, the stability problem is reduced to a generalized
eigenvalue problem of the form (9), in which the minimum (in general, modulo) eigenvalue
determines the critical load, and the corresponding function determines uk (uk ,1 , uk ,2) – the form of
stability loss. The following are accepted as the main mathematical models: the equation of the linear
theory of elasticity of piecewise homogeneous environment to determine the subcritical stress-strain
state; equation of the three-dimensional linearized theory of stability at small deformations to
determine the critical parameters. The model of a piecewise homogeneous medium is used as a
mechanical model. Composite components are modeled by linearly elastic isotropic bodies [
Thus, for a plane deformation, we obtain an eigenvalue problem for a system of partial differential equations with variable coefficients depending on two variables. In this case (the investigation of the stability of a composite material under a compressive load), the conditions of selfadjointness of the operators of the problem of determining the initial stress-strain state and the generalized eigenvalue problem are satisfied, which corresponds to the problem of three-dimensional stability theory. In addition, the differential operators of the problems under consideration are positive-definite under the condition that there is no rigid displacement.
For the numerical solution of the problem of the three-dimensional linearized theory of stability
and the problem of the linear theory of elasticity, a discrete model is used, which is constructed by the
variational-difference method (finite differences) using the concept of basic schemes. In the
computational domain, with the help of straight lines, a difference grid nonuniform in each direction
is introduced so that the material is homogeneous within one grid cell. Thus, the discrete problem of
determining the subcritical stress-strain state (the problem of the linear theory of elasticity) is written
in the form of an SLAE of the form (7) with a sparse symmetric positive definite matrix. The discrete
problem for the stability problem is written in the form of a generalized AEVP of the form (9) of
sparse symmetric positive definite matrices [
In order to provide an effective computer solution of problems (7)–(9), it is necessary to take
into account the following features [
Therefore, for the effective use of computing resources of modern supercomputers in the mathematical modeling of structural strength problems with approximate data and ensuring the reliability of the obtained computer results, it is necessary: identify the original structure of the sparse matrix; determine the need for structural regularization and choose the most efficient (for a given computer architecture) regular structure of matrix; carry out structural regularization of sparse matrices - reducing an arbitrary initial structure of matrix to a regular block-sparse one (block-tape, block-profile, block-skyscraper, block-diagonal with bordering, etc.); choose or develop the most efficient (for a given computer architecture, taking into account the regularized structure of matrices) parallel algorithm for solving the problem; choose, in accordance with the algorithm, a decomposition of matrices scheme (the distribution between processor devices of a parallel computer); choose, in accordance with the algorithm, the scheme for storage of matrix elements in a computer; investigate of the correspondence the mathematical properties of the matrix, entered into the computer, with approximate data and the functional purpose of the algorithm for solving the problem.
Such a scheme for researching and solving SLAE and AEVP for sparse matrices of various
structures is implemented in the created parallel algorithms for modern computers with parallel
organization of computing, namely for MIMD-architecture computers with multi-core processors and
hybrid architecture [
The structure of non-zero elements (the structure of a sparse matrix) is determined by the numbering of unknown problems and can be either regular (for example, band, profile, block-diagonal with framing, etc.), or irregular (arbitrary). The structures of factorization matrices can be obtained by performing the so-called symbol factorization of the initial matrix. In most software tools, before forming the matrices of calculation problems, the unknowns are renumbered (i.e., the sparse structure is reorganized) in order to reduce the total number of floating-point arithmetic operations when solving calculation problems. Figure 1 shows the structures of matrices (actually factorized matrices) of calculation problems that arise when modeling the stress-strain state of structures.
Problems of solving linear algebra problems with sparse matrices and ways to solve them are
considered in [
Among these methods, the following should be mentioned:
• the Cuthill-McKee method and the factor tree method – the use of these methods ensures the
concentration of non-zero elements as close as possible to the main diagonal, which, as a rule,
makes it possible to represent the matrix in a band or in profile form;
• method of parallel sections – the use of this method makes it possible to represent the matrix
as a block-diagonal (with large diagonal blocks) with a border (the order of the diagonal frame
block should be much less than the orders of the main diagonal of blocks) [
Discrete problems with matrices of this type arise in the problems of calculating the stability of structures. For effective implementation in algorithms, for example, matrix-matrix operations of sparse matrices, it is advisable to first reorder (or immediately form) them so that the vast majority of nonzero elements are located in the most dense matrix blocks: i.e, determine zero blocks (all elements are identically equal to 0 ) and non-zero (have at least one non-zero element). It is also desirable that non-zero blocks are filled as much as possible. The block sizes are determined based on the architecture of the computer.
If in the future it is supposed to use the decomposition of the sparse matrix into the product of triangular matrices (e.g., LU , L LT , etc.), then it is necessary to determine the block-sparse structure of the decomposition matrices using the symbolic decomposition of the original matrix. After that, it is necessary to optimize the block-sparse structure of the original matrix or decomposition matrices, using one of the above algorithms and replacing the matrix elements in the algorithms with blocks. In many cases, it is expedient to combine adjacent (in a row or column) non-zero off-diagonal blocks of the optimized matrix structure into one rectangular block (tile), which will be further processed as a whole. In some cases, it is advisable to optimize the structure of diagonal blocks to reduce the number of arithmetic operations with individual blocks of the resulting block-sparse matrix.
Consider another type of sparse matrix – a bordered block-diagonal matrix (figure 3).
Many methods for solving linear algebra problems, that is used in modeling structural strength
problems, are based on algorithms for decomposing a matrix into a product of matrices of standard
types, for example, lower and upper triangular matrices. Such algorithms are characterized by a
gradual decrease from step to step in the size of the processed part of the matrix. Therefore, it is
important to ensure approximately the same amount of calculations, exchanges and synchronizations
performed by each process or thread on a parallel computing model, that is, to exclude the influence
of the Haydn effect [
As research has shown, parallel versions of these algorithms, in which the so-called cyclic
distribution schemes and processing matrices are created (see, for example, [
Quite often, one-dimensional block-cyclic schemes are used: rows or columns of non-zero blocks of the matrix are distributed cyclically between processes or threads of the upper level of parallelism, for example, if the q-th process is distributed elements of matrix rows with numbers sr+ 1 , ... , sr + r, then ( q+ 1) -th process will have rows with numbers i-th process will have rows with numbers s (r +1)+1 , ... , s (r +1)+ s, where s is the number of rows in the block (we can talk about the r-th and (r + 1)-th rows of square matrix blocks of order s).
In the case of a regular sparse matrix structure, parallel algorithms that use one-dimensional
block-cyclic schemes for allocating matrix elements can provide an approximately equal amount of
computations and exchanges performed by each parallel process or thread at each instant of time (see,
for example, [
Such a regular structure is primarily the band structure of the matrix under the obvious condition that the half-width of the matrix band m exceeds the product s× p , where p is the number of processes or threads. In the case of a regular profile structure of the system matrix, by varying the value of s and p , one can practically balance the load of parallel processes at each moment of time if ce / n> sp , where ce is the total number of subdiagonal elements in the matrix profile.
For block-diagonal bordered matrices, cyclic distribution schemes do not allow balancing the load of processes, and it is advisable to use block distribution schemes between computational elements. That means that each diagonal block, together with the corresponding off-diagonal border block, is distributed to a separate process (or thread). A diagonal border block can also be allocated in a separate process (or thread).
There are various storage schemes for sparse matrices, that differ in the way zero elements usage. In some cases, it is allowed to store a portion of zeros in exchange for simplifying the storage scheme; in others, all zeros of the matrix are used; thirdly, zeros are not used at all. The choice of storage format, of course, affects memory requests, and therefore significantly affects the efficiency of software implementations of algorithms for processing sparse matrices (for example, decomposition of matrices into a product of triangular matrices).
In the algorithms created for solving SLAE and AEVP, different formats for storing elements of
sparse matrices are used: coordinate format, sparse row or column format, ELLPACK format, hybrid
format, etc. [
Mathematical models describing applied problems always contain errors in the initial data [
An error in the solution of a mathematical problem, which is caused by an error in the initial
data set, is called a hereditary error. If the hereditary error of the solution of a mathematical problem
is large, then the obtained mathematical solution may not have a physical meaning, that is, such a
solution will not contain the solution of a physical problem. Therefore, it is necessary to determine the
scope of the mathematical model, because the hereditary error cannot be corrected by mathematical
methods of solving the problem. To reduce the hereditary error, you can either increase the accuracy
of the input data or reformulate the problem with respect to other parameters [
The investigation of the correctness of the problem statement in a computer is reduced to
checking two relations:
1.0 + γ ≠ 1.0,
(‖ΔA‖‖A‖) h ( A )<1;
(10)
where h ( A ) – is the condition number of the system matrix, γ =h−1 ( A ). The first condition, fulfilled in
floating-point arithmetic, means that the matrix is non-singular within the limits of machine accuracy (machine
non-singular), and the second one means that it is non-singular within the limits of the accuracy of the input data
[
Thus, when conditions (10) are met, the solution of the computer problem with approximate data exists, is unified and stable. Such a computer task should be considered as correctly posed within the accuracy of the initial data input. Otherwise, the matrix of the system may turn out to be a matrix of incomplete rank, and the computer problem should be considered as incorrectly posed. The fulfillment of condition (10) at an increased computer bit rate in accordance with the value of the condition number of the matrix indicates that the original problem is poorly conditioned relative to the previous bit rate and it is possible to obtain a single solution at an increased bit rate.
So, the computer problem, which should to be solved in the end, always has approximate data relative to the original problem (due to the hereditary error in the initial data, due to the discretization errors, due to the errors in obtaining numerical data in the computer).
In the created parallel algorithms for solving SLAE and AEVP for sparse matrices of various structures on modern high-performance computing systems with parallel organization of calculations, a research is provided in a computer environment of the mathematical properties of matrices using various computer bit rate (the technology is described in the work [16]).
Let’s consider some practical problems of mathematical modeling of structural strength, which are reduced to solving linear algebra problems with sparse matrices.
To analyze the strength of building structures on single-processor computers, Software package
Lira-Sapr is widely used – a multi-threaded version of PC Lira [17]. The developed parallel
algorithms for solving SLAE with sparse matrices of various structures were used as a processor
component of the parallel version of PC Lira – the software package Lira-cluster [
Consider the problem of analyzing the strength of a 27-storey building (figure 4). The sparse symmetric matrix obtained as a result of discretization by the finite element method (FEM-mesh) has an order of 661590 and a band half-width of 34242, density (the matrix tape is filled with non-zero elements by 5%). After optimizing the matrix structure by the minimum degree method, a skyscraper structure matrix was obtained (order – 661590, ribbon width – 541257, density <1%). The time for solving the problem using Lira-Sapr on a personal computer is 3254 sec, using Lira-cluster on the hybrid architecture of the SKIT-4 supercomputer is 57 sec. Thus, an acceleration of 23 times was achieved.
Supercomputer SKIT-4 of the V.M. Glushkov Institute of Cybernetics National Academy of Sciences of Ukraine has the following specifications: CPU – Intel(R) Xeon(R) CPU E5-26700, clock speed 2.60 GHz, speed 8 GT/s, cache memory – 20 MB, per node: 2 CPUs with 8 cores + Hyperthrefding =32 cores, Max Memory Size – 384 GB; GPU – Nvidia Tesla M2050, 3 GB memory, peak performance – 515 Gflops [19]. 4.1.
The created parallel algorithms for solving SLAE and AEVP for sparse matrices of various
structures were used in mathematical modeling of the problem of structural stability - the problem of
stability of a layered two-component composite material with a regular structure (Fig. 5, left) under a
uniform single-axis load on computers of various architectures. Since the stability problem is solved
for the computational domain of finite dimensions (Fig. 5, right) and the influence of the boundary
conditions on the lateral sides of the computational domain on the critical stability parameters is
studied, we distinguish between the loss of stability in a composite material and a composite sample
corresponding to this material [
Comparative results of solving the obtained partial generalized eigenvalue problem Ax =λ Bx in the problem of mathematical modeling of the stability of new composite materials are given for band symmetric positive definite matrices (the order of matrices A and B is 12282, the half-width band of matrix A is 6212, the half-width tape of the matrix B is 71).
Compared to the sequential version of the algorithm, the case of solving the AEVP of band symmetric matrices, a 7-fold acceleration was obtained by the parallel algorithm on the MIMD architecture; on a hybrid architecture - 18 and 33 times, using one and two GPUs respectively compared. In the case of the AEVP solution for sparse matrices reduced to the block-diagonal form with border, an acceleration of 43 times was obtained on the hybrid architecture using two GPUs compared to the sequential version of the algorithm [21].
Methods for regularization, decomposition and storage of sparse matrices are proposed, which allow to significantly increase the efficiency of solving problems of mathematical modeling of structural strength. New methods and computer block-cyclic and block algorithms for highperformance computing created on their basis for solving linear algebra problems (systems of linear algebraic equations, algebraic eigenvalue problem) with approximate data can be used for mathematical modeling of different processes on modern supercomputers.
[15] O. V. Popov (2013) About parallel algorithms for factorization of sparsed ed matrices. Computer
Mathematics: Collection of scientific papers. N. 2. P. 115–124. (in Ukrainian) [16] О. М. Khimich, T. V. Chistjakova, V. A. Sidoruk, P. S. Ershov (2021) Adaptive computer technologies for solving problems of computational and applied mathematics. Cybernetics and Systems Analysis. Vol. 57. N 6, P. 990–997. https://doi.org/10.1007/s10559-021-00424-z. [17] LIRA-SAPR [Online] Available from: http://www.liraland.ua. [Accessed: 1 June 2022] [18] Programs for calculation and design of structures [Online] Available from: http://www.lira.com.ua. [Accessed: 1 June 2022] [19] Supercomputer complex SKIT. [Online] Available from: http://icybcluster.org.ua/ [Accessed: 1
June 2022] [20] A. N. Khimich, A. V. Popov, A. V. Chistyakov (2017) Hybrid Algorithms for Solving the Algebraic Eigenvalue Problem with Sparse Matrices. Cybernetics and Systems Analysis. N 6.
P. 132–146. https://doi.org/10.1007/s10559-017-9996-5. [21] A. N. Khimich, A. V. Popov, O. V. Chistyakov, V. A. Sidoruk (2020) A Parallel Algorithm for Solving a Partial Eigenvalue Problem for Block-Diagonal Bordered Matrices. Cybernetics and Systems Analisys. V. 53. N. 6. P. 913–923. https://doi.org/10.1007/s10559-020-00311-z.