Scroll to:

# Implementation of Basic Operations for Sparse Matrices when Solving a Generalized Eigenvalue Problem in the ACELAN-COMPOS Complex

https://doi.org/10.23947/2687-1653-2023-23-2-121-129

### Abstract

**Introduction. **The widespread use of piezoelectric materials in various industries stimulates the study of their physical characteristics and determines the urgency of such research. In this case, modal analysis makes it possible to determine the operating frequency and the coefficient of electromechanical coupling of piezoelectric elements of various devices. These indicators are of serious theoretical and applied interest. The study was aimed at the development of numerical methods for solving the problem of determining resonance frequencies in a system of elastic bodies. To achieve this goal, we needed new approaches to the discretization of the problem based on the finite element method and the execution of the software implementation of the selected method in C# on the .net platform. Current solutions were created in the context of the ACELAN-COMPOS class library. The known methods of solving the generalized eigenvalue problem based on matrix inversion are not applicable to large-dimensional matrices. To overcome this limitation, the presented scientific work implemented the logic of constructing mass matrices and created software interfaces for exchanging data on eigenvalue problems with pre- and postprocessing modules.

**Materials and Methods. **A platform was used to implement numerical methods .net and the C# programming language. Validation of the research results was carried out through comparing the values found with solutions obtained in well-known SAE packages (computer-aided engineering). The created routines were evaluated in terms of performance and applicability for large-scale tasks. Numerical experiments were carried out to validate new algorithms in small-dimensional problems that were solved by known methods in MATLAB. Next, the approach was tested on tasks with a large number of unknowns and taking into account the parallelization of individual operations. To avoid finding the inverse matrix, a modified Lanczos method was programmatically implemented. We examined the formats for storing matrices in RAM: triplets, CSR, СSC, Skyline. To solve a system of linear algebraic equations (SLAE), an iterative symmetric LQ method adapted to these storage formats was used.

**Results. **New calculation modules integrated into the class library of the ACELAN-COMPOS complex were developed. Calculations were carried out to determine the applicability of various formats for storing sparse matrices in RAM and various methods for implementing operations with sparse matrices. The structure of stiffness matrices constructed for the same task, but with different renumbering of nodes of a finite element grid, was graphically visualized. In relation to the problem of the theory of electroelasticity, data on the time required to perform basic operations with stiffness matrices in various storage formats were summarized and presented in the form of a table. It has been established that the renumbering of grid nodes gives a significant increase in performance even without changing the internal structure of the matrix in memory. Taking into account the objectives of the study, the advantages and weaknesses of known matrix storage formats were named. Thus, CSR was optimal when multiplying a matrix by a vector, SKS was optimal when inverting a matrix. In problems with the number of unknowns of the order of 10^{3}, iterative methods for solving a generalized eigenvalue problem won in speed. The performance of the software implementation of the Lanczos method was evaluated. The contribution of all operations to the total solution time was measured. It has been found that the operation of solving SLAE takes up to 95% of the total time of the algorithm. When solving the SLAE by symmetric LQ method, the greatest computational costs were needed to multiply the matrix by a vector. To increase the performance of the algorithm, parallelization with shared memory was resorted to. When using eight threads, the performance gain increased by 40–50%.

**Discussion and Conclusion. **The software modules obtained as part of the scientific work were implemented in the ACELAN-COMPOS package. Their performance for model problems with quasi-regular finite element grids was estimated. Taking into account the features of the structures of the stiffness and mass matrices obtained through solving the generalized eigenvalue problem for an electroelastic body, the preferred methods for their processing were determined.

### Keywords

#### For citations:

Oganesyan P.A.,
Shtein O.O.
Implementation of Basic Operations for Sparse Matrices when Solving a Generalized Eigenvalue Problem in the ACELAN-COMPOS Complex. *Advanced Engineering Research (Rostov-on-Don)*. 2023;23(2):121-129.
https://doi.org/10.23947/2687-1653-2023-23-2-121-129

**Introduction. **Devices made of piezoelectric materials have been widely used, actively studied and improved for a long time. Medical ultrasound devices (diagnostic equipment, ultrasonic scalpels) [1–4] and mobile energy generators [5] should be noted separately. Paper [6] described the combination of photo- and piezoelectric effects to create efficient compact energy sources. New materials designed for the application under specific conditions are being studied in science and industry. In [7], the creation of a lead-free piezo-active composition suitable for operation at various temperatures was considered.

In the study on piezoelectric elements, a key role is played by the modal analysis stage, which enables to establish the resonance and antiresonance frequencies of the device. These data:

- are needed to find out the operating frequency of the device;
- provide determining the electromechanical coupling factor — an important performance indicator of the device;
- are input information in numerical experiments for problems on forced oscillations.

The study was aimed at the creation of numerical methods for solving the problem of determining resonance frequencies for a system of elastic bodies. Achieving the stated goal required solving two tasks. The first was to develop methods of discretization of the problem based on the finite element method (FEM). The second was to carry out a software implementation of the selected method in C# on the .net platform. All known programs take into account the context of the ACELAN-COMPOS class library [8]. When solving a generalized eigenvalue problem, methods based on matrix inversion are widely used. However, they are not applicable to large-dimensional matrices. In the presented research, this limitation was overcome as follows:

- the logic of constructing mass matrices was additionally implemented;
- software interfaces were created for exchanging data on eigenvalue tasks with pre- and postprocessing modules.

**Materials and Methods.** Principally, the proposed approach was designed to solve static problems of electroelasticity when implementing the averaging method [9], which was used to calculate the effective properties of piezo composites. In this regard, only stiffness matrices were presented at the stage of constructing global FEM matrices. In this study, we additionally implemented the logic of constructing mass matrices and developed software interfaces (application programming interface, API) for exchanging data on eigenvalue tasks with pre- and postprocessing modules. The developed routines were evaluated in terms of performance and applicability for large-scale tasks. Numerical experiments were carried out to validate the algorithms created for such small-dimensional problems that provide obtaining a solution by general methods in the MATLAB computing package. Next, testing was performed on tasks with a large number of unknowns and taking into account the parallelization of individual operations.

The mathematical model of the problem being solved consists of the defining relations [9]:

,,(1)

, , (2)

, (3)

Here, σ— stress tensor; ρ* _{j}* — body density; ε — strain tensor;

*u*— displacement vector;

*D*— electric displacement vector; E — electric-field vector;

— body force vector; — electric potential; , , — damping coefficients;

, , — tensors of elastic constants, piezoelectric modules and dielectric permittivity; index

*j*— body number in the model.

The discretization is performed by replacing:

Here, — shape function matrix for the displacement field; *N*_{ϕ} — shape function vector for electric potential;

— global vectors of the corresponding nodal degrees of freedom.

In this case, the original problem (1–3) takes the form:

. (4)

Here, matrices *M* and *K* are global matrices of mass and stiffness, respectively, and the vector is a general vector of unknowns:

a = [*U*, Ф].

In the problem of the theory of electroelasticity:

(5)

Matrices , and — symmetric. In the case of harmonic oscillations at natural frequency , it is possible to write:

, denoting the corresponding eigenvector by .

Consider free oscillations if . In this case, task (4) is represented as:

. (6)

Thus, the original problem is reduced to a generalized eigenvalue problem (6). For nonzero , inequality (6) is solved by finding the matrix inverse to . However, at the same time, the sparse matrix becomes full, i.e., the method is unsuitable for large matrices. Therefore, it is needed to use other methods that do not require finding the inverse matrix. To solve this problem, a modified Lanczos method was programmatically implemented in this paper [10]. The author of this modification is T. S. Martynova. The description of the development is not given in this article. Of the operations used in the method, the most expensive from the point of view of computational resources was the solution to a system of linear algebraic equations (SLAE), needed for performing a spectral transformation.

Matrices *M* and *K* — sparse, with a small number of nonzero elements. Several formats are used to store such matrices in RAM:

- triplets, or coordinate format;
- CSR (compressed sparse row);
- СSC (compressed sparse column);
- Skyline storage format (SKS method).

The coordinate format involves storing triples (triplets) of values (*i*,* j*, *k*), representing coordinates (*i*, *j*) and values (*k*) of nonzero elements. CSR is sometimes referred to as CRS or Yale format. It involves storing a sparse matrix in the form of three arrays. Consider matrix of *N* size with *NZ* nonzero elements. We describe the possible organization of its storage. All nonzero elements must be placed in one array of *NZ* size. The positions of these elements in the columns should be placed in another array of *NZ* size, and the third array of *N* size should be used to store the indices of the first elements of the rows. Similarly, the storage in CSV format is implemented.

The SKS format assumes the storage of a variable-width matrix band that includes all nonzero elements. In this case, zeros are allowed. The efficiency of this format depends on the renumbering of the matrix rows. Methods for reducing the size of the tape are described in [11]; however, their applicability to the stiffness matrix obtained when solving a three-dimensional problem using FEM requires a separate study.

To solve the SLAE (system of linear algebraic equations), an iterative symmetric LQ method (SYMLQ [12]) adapted to the storage formats listed above was used.

**Research Results.** At the beginning of the study, we chose the optimal storage format for sparse matrices. The coordinate format enabled to quickly add and change an element of the matrix. These operations were needed at the stage of assembling the global matrix and taking into account the boundary conditions. In addition, for ill-conditioned matrices, to which refers, a preliminary transformation is often used for normalization. It is also convenient to perform it in the coordinate format. However, this format is ineffective when it comes to algebraic operations.

CSR is ill-adapted for changing the structure of the matrix: by adding a nonzero element, you need to insert into two arrays. In this case, the matrix is multiplied by a vector quite easily and efficiently.

SKS has similar problems with the addition of nonzero elements and is highly dependent on the renumbering of unknowns in the problem. We focus on the example of a quasi-regular grid, which is used in the ACELAN-COMPOS package to work with representative volumes of composites. The width of the band containing all nonzero elements can be predetermined and depends on the number of nodes and the type of final element. In the general case of an arbitrary finite element grid, it is difficult to estimate the size of the band in advance.

Four methods of numbering unknowns were used in numerical experiments. Figure 1 shows the structure of stiffness matrices constructed for the same task, but with different renumbering of nodes of a finite element grid.

Fig. 1. Stiffness matrix structure with various methods of node numbering:*a *— unknowns are ordered by nodes; *b* — first, displacement nodes are ordered, then potentials; *c* — nodes are sorted by layers of the FE grid, and unknowns — by nodes;*b, d *— nodes are sorted by layers of the FE grid, and unknowns — as in example

Thus, the grid was a cube with regular partitioning by eight-node finite elements. A model matrix of 500 lines was used for the illustrations. The matrices shown in Figures 1 *a* and 1 *b*, were not subjected to additional renumbering of nodes and differed only in the numbering of degrees of freedom. In 1 *a*:

.

In 1 *b,* the unknowns responsible for the potential distribution were collected at the end of the vector:

.

The unknowns in matrices 1 *c* and 1 *d* were numbered similarly, but the nodes of the finite element grid were pre-numbered according to their coordinates through alternately sorting all nodes by each of the coordinates. This technique is widely used to build more efficient SLAE solution modules, as it enables to work with the matrix in a suitable band format, convenient for parallelization. Similar external modules were implemented for the ACELAN-COMPOS complex [13–15]; however, in this work, only formats for storing sparse matrices of a general form were used.

Table 1 summarizes the data on the time required to perform basic operations with matrices in various formats.

Table 1

Time to perform basic operations with the stiffness matrix in the problem of the theory of electroelasticity.

19,652 rows

Storage format |
Operation |
Elapsed time, ms |
|||

1 |
1 |
1 |
1 |
||

CSR |
Conversion from coordinate format |
123 |
132 |
97 |
117 |

CSR |
Multiplication by vector, 100 operations |
260 |
260 |
260 |
260 |

SKS |
Conversion from coordinate format |
690 |
703 |
124 |
268 |

SKS |
Multiplication by vector, 100 operations |
60,558 |
61,450 |
7,616 |
22,113 |

The experimental results showed that the conversion operation from a coordinate storage format to a compact one took little time. At the same time, the renumbering of grid nodes to form a block-tape matrix made it possible to get a noticeable increase in performance even without changing the internal structure of the matrix in memory. CSR format turned out to be optimal in terms of the efficiency of the matrix-vector multiplication operation. When the matrix was inverted, SKS format was more efficient, but for problems with the number of unknowns of the order of 103, iterative methods for solving a generalized eigenvalue problem worked noticeably faster.

Further, the performance of the software implementation of the Lanczos method was experimentally evaluated. The contribution of all operations to the total solution time was measured. As a result, it was found that the operation of solving SLAE took up to 95% of the total time of the algorithm. In the course of the algorithm, a Krylov subspace was constructed, and depending on its dimension, the number of SLAE that needed to be solved changed. Note that the dimension of the Krylov subspace was chosen based on heuristics with respect to the number of desired eigenvalues. Here, the SLAE differed only in the right part, so that the requirements for allocated memory remained low. Among the basic operations used in the course of solving SLAE by the SYMLQ method, the greatest computational costs were needed for multiplying the matrix by the vector.

To increase the performance of the algorithm, the simplest parallelization with shared memory was implemented. Blocks of rows were allocated for the CRS format. They were transmitted to separate threads that calculated the corresponding components of the resulting vector. The performance gain was 40–50% when using 8 threads. At the same time, for matrices of the order of 10^{3} elements, the increase was about 40%, and for matrices of the order of 10^{4} — about 50%.

**Discussion and Conclusions.** Within the framework of this study, a method for solving a generalized eigenvalue problem for matrices obtained by modeling electroelastic bodies was implemented. Software modules were created in C# for constructing mass matrices by the finite element method and performing auxiliary operations within the framework of the Lanczos method (working with Krylov subspace vectors, reorthogonalization, finding eigenvectors). The computational complexity was mainly due to the operations of multiplying sparse matrices by a vector. In this regard, numerical experiments were carried out to determine the optimal formats for storing matrices, the optimal structure of the matrix obtained as a result of renumbering the nodes of the FE grid and degrees of freedom in the nodes. A version of the SYMMLQ iterative algorithm using parallel computing was developed. The final scheme of work included three points. First, global matrices were constructed in coordinate format with a renumbering algorithm (Fig. 1 *c*). Secondly, the data was converted to CRS format. Thirdly, the data was processed by the Lanczos method, which included the SYMLQ method for solving SLAE. The results of the work were included in the ACELAN-COMPOS software package.

## References

1. Lisong Deng, Mingxiang Ling. Design and Integrated Stroke Sensing of a High-Response Piezoelectric Direct Drive Valve Enhanced by Push-Pull Compliant Mechanisms. Review of Scientific Instruments. 2022;93(3):035008. https://doi.org/10.1063/5.0067483

2. Urtnasan Erdenebayar, Jong-Uk Park, Pilsoo Jeong, et al. Obstructive Sleep Apnea Screening Using a Piezo-Electric Sensor. Journal of Korean Medical Science. 2017;32(6):893–899. https://doi.org/10.3346/jkms.2017.32.6.893 3. Skaliukh AS, Gerasimenko TE, Oganesyan PA, Solovieva AA. Effect of Geometric and Physical Parameters on Resonant Frequencies of Ultrasonic Vibrations of Elastic and Piezoelectric Element System. Vestnik of Don State Technical University. 2017;17(4):5–13. https://doi.org/10.23947/1992-5980-2017-17-4-5-13

3. Bulletti A, Capineri L, Floridia D. Automatic System to Measure the Impedance of Piezoelectric Actuators Used in Ultrasonic Scalpels. In book: Sensors and Microsystems. Cham: Springer; 2014. Vol. 268. P. 71–74. https://doi.org/10.1007/978-3-319-00684-0_14

4. Keli Li, Qisheng He, Jiachou Wang, et al. Wearable Energy Harvesters Generating Electricity from Low-Frequency Human Limb Movement. Microsystems & Nanoengineering. 2018;4:24. https://doi.org/10.1038/s41378-018-0024-3

5. Wenbo Peng, Chenhong Wang, Fangpei Li, et al. Piezo- and Photo-Voltage Field-Effect Transistor. Nano Energy. 2022;105:108025. https://doi.org/10.1016/j.nanoen.2022.108025

6. Tangyuan Li, Chang Liu, Peng Shi, et al. High-Performance Strain of Lead-Free Relaxor-Ferroelectric Piezoceramics by the Morphotropic Phase Boundary Modification. Advanced Functional Materials. 2022;32(32):2202307. https://doi.org/10.1002/adfm.202270184

7. Kurbatova NV, Nadolin DK, Nasedkin AV, et al. Finite Element Approach for Composite Magneto-Piezoelectric Materials Modeling in ACELAN-COMPOS Package. In book: Analysis and Modelling of Advanced Structures and Smart Systems. Singapore: Springer; 2018. Vol. 81. P. 69–88. https://doi.org/10.1007/978-981-10-6895-9_5

8. Belokon AV, Nasedkin AV, Soloviev AN. Novye skhemy konechnoehlementnogo dinamicheskogo analiza p'ezoehlektricheskikh ustroistv. Journal of Applied Mathematics and Mechanics. 2002;66(3):491–501. (In Russ.)

9. Zhongming Teng, Lei-Hong Zhang. A Block Lanczos Method for the Linear Response Eigenvalue Problem. Electronic Transactions on Numerical Analysis. 2017;46:505–523. https://doi.org/10.13140/RG.2.2.16369.68962

10. Chagas G, Oliveira SLGD. Metaheuristic-Based Heuristics for Symmetric-Matrix Bandwidth Reduction: A Systematic Review. Procedia Computer Science. 2015;51:211–220. https://doi.org/10.1016/j.procs.2015.05.229

11. Paige CC, Saunders MA. Solution of Sparse Indefinite Systems of Linear Equations. SIAM Journal on Numerical Analysis. 1975;12(4):617–629. https://doi.org/10.1137/0712047

12. Fassbender H, Ikramov K. SYMMLQ-like Procedure of Ax = b where A is a Special Normal Matrix. Calcolo. 2006;43(1):17–37. https://doi.org/10.1007/s10092-006-0112-x

13. Vasilenko A, Veselovskiy V, Metelitsa E, et al. Precompiler for the ACELAN-COMPOS Package Solvers. In: Proc. 16th Int. Conf.: Parallel Computing Technologies. Cham: Springer; 2021. Vol. 12942. P. 103–116. https://doi.org/10.1007/978-3-030-86359-3_8 15. Steinberg BYa, Vasilenko AA, Veselovskiy VV, et al. Solvers for Systems of Linear Algebraic Equations with Block-Band Matrices. Bulletin of the South Ural State University. Series “Mathematical Modeling, Programming & Computer Software”. 2021;14(3):106–112. https://doi.org/10.14529/mmp210309

### About the Authors

**P. A. Oganesyan**Russian Federation

Pavel A. Oganesyan, Cand.Sci. (Phys.-Math.), associate professor of the Mathematical Modeling Department, Vorovich Institute for Mathematics, Mechanics, and Computer Science

105/42, Bolshaya Sadovaya St., Rostov-on-Don, 344006, RF

**O. O. Shtein**Russian Federation

Olga O. Shtein, teaching assistant of the Applied Mathematics and Programming Department, Vorovich Institute for Mathematics, Mechanics, and Computer Science

105/42, Bolshaya Sadovaya St., Rostov-on-Don, 344006, RF

### Review

**Supporting agencies:**the authors would like to thank AN Solovyov and TS Martynova for their assistance in the development of numerical methods. Thanks are extended to the Russian Scientific Foundation for financial support of the research with grant no. 22–21–00318, https://rscf.ru/project/22-21-00318

#### For citations:

Oganesyan P.A.,
Shtein O.O.
Implementation of Basic Operations for Sparse Matrices when Solving a Generalized Eigenvalue Problem in the ACELAN-COMPOS Complex. *Advanced Engineering Research (Rostov-on-Don)*. 2023;23(2):121-129.
https://doi.org/10.23947/2687-1653-2023-23-2-121-129