Reusing Preconditioners in Projection Based Model Order Reduction Algorithms

Dynamical systems are pervasive in almost all engineering and scientific applications. Simulating such systems is computationally very intensive. Hence, Model Order Reduction (MOR) is used to reduce them to a lower dimension. Most of the MOR algorithms require solving large sparse sequences of linear systems. Since using direct methods for solving such systems does not scale well in time with respect to the increase in the input dimension, efficient preconditioned iterative methods are commonly used. In one of our previous works, we have shown substantial improvements by reusing preconditioners for the parametric MOR (Singh et al. 2019). Here, we had proposed techniques for both, the non-parametric and the parametric cases, but had applied them only to the latter. We have three main contributions here. First, we demonstrate that preconditioners can be reused more effectively in the non-parametric case as compared to the parametric one. Second, we show that reusing preconditioners is an art via detailed algorithmic implementations in multiple MOR algorithms. Third and final, we demonstrate that reusing preconditioners for reducing a real-life industrial problem (of size 1.2 million), leads to relative savings of up to 64 % in the total computation time (in absolute terms a saving of 5 days).


I. INTRODUCTION
Dynamical systems arise in many engineering and scientific applications such as weather prediction, machine design, circuit simulation, biomedical engineering, etc.Generally, dynamical systems corresponding to real-world applications are extremely large in size.A set of equations describing a parametric nonlinear second-order dynamical system is represented as g( ẍ(t), p) = f ( ẋ(t), p) + h(x(t), p, u(t)), where t is the time variable, x(t) : R → R n is the state, p = (p 1 , p 2 , . . ., p k ) is the set of parameters (with p j ∈ R; for j= 1, . . ., k), u(t) : R → R m is the input, y(t) : R → R q is the output, C T ∈ R q×n is the output matrix, and g(•) : R n+k → R n , f (•) : R n+k → R n and h(•) : R n+k+m → R n are some nonlinear functions [1]- [6].If m and q both are equal to one, then we have a Single-Input Single-Output (SISO) system.Otherwise, it is called a Multi-Input Multi-Output (MIMO) (m and q > 1) system.The functions g(•), f (•), and h(•) are usually simplified as [2], [6] g( ẍ(t), p) = k ∑ j=1 g j (p)g( ẍ(t)), VOLUME 4, 2016 where g j (•), f j (•), h j (•) : R k → R are scalar-valued functions while g(•), f(•) : R n → R n , and h(•) : R n+m → R n are vector-valued.Next, we look at simplifications to (1) based upon the three predicates; the presence of parameters; the degree of non-linearity, and the order of the system.
• If g j (p), f j (p), and h j (p) are independent of the parameters, then (1) becomes a non-parametric dynamical system.
• Bilinear systems are one of the common types of nonlinear dynamical systems.Here, there is a product between the state variables and the input variables.Another important class of nonlinear dynamical systems is the quadratic systems.Here, there is product among the state variables.If g(•) and f(•) are linear functions of the state variables, and h(•) is a linear function of the state and the input variables, then (1) is called a linear dynamical system.• Finally, if the second derivative term in (1) is not present, then (1) becomes a first-order dynamical system.Simulation of large dynamical systems can be unmanageable due to high demands on computational resources.These large systems can be reduced into a smaller dimension by using Model Order Reduction (MOR) techniques [4], [7]- [11].The reduced system has approximately the same characteristics as the original system but it requires significantly less computational effort in simulation.MOR can be done in many ways such as balanced truncation, Hankel approximations, and Krylov projection [4], [7], [8], [11].Among these, the projection methods are quite popular, and hence, we focus on them.
Some of the commonly used projection-based MOR algorithms for different types of dynamical systems are summarized in Table 1.
In the above mentioned MOR algorithms, sequences of very large and sparse linear systems arise during the model reduction process.Solving such linear systems is the main computational bottleneck in efficient scaling of these MOR algorithms for reducing extremely large dynamical systems.Preconditioned iterative methods are commonly used for solving such linear systems [25], [26].In most of the above listed MOR algorithms, the change from one linear system to the next is usually very small, and hence, the applied preconditioner could be reused.
Next, we briefly summarize the past work that has been done in the field of reusing preconditioners.This technique was first applied in the QMC context, where it was referred to as recycling preconditioner [27], [28].In the optimization context, this approach was applied in [29], where it was termed as the preconditioner update.Such a technique was first applied in the MOR context in [12], and more recently in [30], where the focus was mostly on MOR of non-parametric linear first-order dynamical systems (part of the first category above).
The main goal of this paper is to demonstrate the reuse of preconditioners in the remainder of the algorithms for the first category above (MOR of non-parametric linear second-order dynamical systems) as well as the algorithms for the second category above (MOR of non-parametric bilinear/ bilinearquadratic dynamical systems).
In one of our recent works [31], we had proposed a general framework for reuse of preconditioners during MOR of both non-parametric and parametric dynamical systems.However, in [31] we had demonstrated application of this framework for the parametric case only.That is, the third category above (MOR of parametric linear dynamical systems).We are currently (and separately) working on the algorithms for the fourth category above as well (MOR of parametric bilinear/ bilinear-quadratic dynamical systems).
To summarize, in this paper we broadly demonstrate the application of our above mentioned framework for MOR of non-parametric dynamical systems.We have four contributions as below, which have not been catered in any of the above cited papers.
(i) We demonstrate that the reuse of preconditioners can be done more effectively in the non-parametric case as compared to the parametric case because of the lack of parameters in the former.(ii) We show that as the underlying MOR algorithms get more intricate, the reuse of preconditioners needs to be fine-tuned.(iii) We highlight that there are multiple pitfalls in the algorithmic implementation of reusing preconditioners, which if not done efficiently, could actually increase the computational complexity of the underlying MOR algorithms instead of reducing it.(iv) We experiment on a massively large and real-life industrial problem (BMW disc brake model), which is of size 1.2 million.We reduce the total computation time from 197 hours to about 72 hours (approximately), leading to a saving of 64%.The paper has four more sections.We discuss MOR techniques in Section II.The theory of reusing preconditioners is described in Section III.We support our theory with numerical experiments in Section IV.Finally, conclusions and future works are discussed in Section V.For the rest of this paper, • f denotes the Frobenius norm, • denotes the Euclidean norm for vectors and the induced spectral norm for matrices, ⊗ refers to the Kronecker product (i.e. an operation on two matrices of arbitrary size), vec(•) signifies the vectorization of a matrix, and I denotes the Identity matrix.
• Approximating the reduced state vector x(t) using V as where r(t) is the residual after projection.• Enforcing the residual r(t) to be orthogonal to V or V T r(t) = 0 leads to the reduced system given as follows: where M = V T MV, D = V T DV, K = V T KV, F = V T F, and ĈT = C T V .To compute this projection matrix V , AIRGA matches the moments of the original system transfer function and the reduced system transfer function.We briefly summarize AIRGA in Algorithm II.1, where parts relevant to solving linear systems are only listed.BIRKA [16] is a Petrov-Galerkin projection based algorithm for MOR of the bilinear first-order dynamical systems, which for the MIMO case are represented as where span two r-dimension subspaces (where, as earlier, r n ).In principle, the Petrov-Galerkin projection method involves the steps below.
• Approximating the reduced state vector x(t) using V as where r(t) is the residual after projection.
Output: M, D, K, F, Ĉ. 1: z = 1 2: while (no convergence) do 3: end for while (no convergence) do 9: for i = 1, . . ., do 10: end for 13: end while 15: "All the given set of expansion points (i.e.s 1 , s 2 , . . ., s ) are updated" 16: • Enforcing the residual r(t) to be orthogonal to W or W T r(t) = 0 leads to the reduced system given by , and (W T V ) −1 is assumed to be invertible.Here, V and W are computed by using interpolation, where the original system transfer function and its derivative are respectively matched with the reduced system transfer function and its derivative at a set of points.We briefly summarize BIRKA in Algorithm II.2, where again, only parts related to solving linear systems are listed.
• As before, approximating the reduced state vector x(t) using V as where r(t) is the residual after projection.• Enforcing the residual r(t) to be orthogonal to W or W T r(t) = 0 leads to the reduced system given by Here, V and W are computed by matching the moments of the original system transfer function and the reduced system transfer function.We briefly summarize QB-IHOMM in Algorithm II.3, where as earlier, only parts related to solving linear systems are listed.

III. PROPOSED WORK
Here, we discuss preconditioned iterative methods in Section III-A.In Section III-B, we revisit the theory of reusing preconditioners from [31].Finally, we discuss application of reusing preconditioners to the earlier discussed algorithms in Section III-C.

A. PRECONDITIONED ITERATIVE METHODS
Krylov subspace based methods are very popular class of iterative methods [32], [33].Let Ax = b be a linear system, with A ∈ R n×n , b ∈ R n , x 0 the initial solution and r 0 (where r 0 = b −Ax 0 ) the initial residual.We find the solution of a linear system in K k (A, r 0 ) = span{r 0 , Ar 0 , A2 r 0 , . . ., A k−1 r 0 }, where K k (•, •) represents the Krylov subspace.
Often iterative methods are slow or fail to converge, and hence, preconditioning is used to accelerate them.If P is a non-singular matrix that approximates the inverse of A (that is, P ≈ A −1 ), then the preconditioned system becomes AP x = b with x = P x 2 .We expect that the preconditioned iterative solves would find a solution in less amount of time as compared to the unpreconditioned ones.
For most of the input dynamical systems (as mentioned here), the Krylov subspace methods fail to converge (see Numerical Experiments section).Hence, we use a preconditioner.
The goal is to find a preconditioner that is cheap to compute as well as apply.There exist many preconditioning techniques [34]- [37], like incomplete factorizations, Sparse Approximate Inverse (SPAI) etc. SPAI preconditioners are known to work in the most general setting and can be easily parallelized.Hence, we use them.
For constructing a preconditioner P corresponding to a coefficient matrix A, we focus on methods for finding approximate inverse of A by minimizing the Frobenius norm of the residual matrix I − AP.This minimization problem can be rewritten as [36] min Here, the columns of residual matrix I − AP can be computed independently, which is an important property that can be exploited.Hence, the solution of ( 6) can be separated into n independent least square problems as where e i and p i are the i-th column of I and P, respectively.The above minimization problem can be implemented in parallel and one can efficiently obtain the explicit approximate inverse P of A.

B. THEORY OF REUSING PRECONDITIONERS
In general, the linear systems of equations generated by lines 4 and 10 of Algorithm II.1; lines 4 and 5 of Algorithm II.2; and lines 4 and 10 of Algorithm II.3 have the following form: Let P 1 be a good preconditioner for A 1 , that is, computed by min Now, we need to find a good preconditioner P 2 corresponding to A 2 .Using the standard SPAI theory, this means solving min If we are able to enforce A 1 P 1 = A 2 P 2 , then P 2 will be an equally good preconditioner for A 2 as much as P 1 is a good preconditioner for A 1 (since the Spectrum of A 2 P 2 would be same as that of A 1 P 1 , on which convergence of any Krylov subspace method depends).Since P 2 is unknown here, we have a degree of freedom in choosing how to form it.Without loss of generality, we assume that P 2 = Q 2 P 1 , where Q 2 is an unknown matrix.Here, we need to enforce Thus, instead of solving the minimization problem (8), we can solve min Note that P 2 here is never explicitly formed by multiplying two matrices Q 2 and P 1 .Rather, always a matrix-vector product is done to apply the preconditioner.Next, we apply a similar argument for finding a good preconditioner P i corresponding to A i .For this we refer to one of our recent works [31], which focused on MOR of parametric linear dynamical systems (category three from the Introduction).We can obtain P i by enforcing either A 1 P 1 = A i P i or A i−1 P i−1 = A i P i .For these two cases, P i would be as effective preconditioner for A i as P 1 is for A 1 or P i−1 is for A i−1 , respectively.These two approaches are summarized in Table 2.

First approach
Second approach In [31], we have conjectured (with evidence) the following two results: (a) In the parametric case, the first approach is more beneficial.This is because, in this case although the two approaches have a similarly hard minimization problem (attributed to slowly varying parameters, and in-turn, slowly changing matrices), the computation of P i from P 1 in the first approach leads to a preconditioner with less approximation errors, and hence, a one which is more accurate.(b) In the non-parametric case, the second approach is more suited.This is because in this case the minimization problem of the second approach is much easier to solve as compared to the first approach (attributed to rapidly changing expansion/ interpolation points, and in-turn, rapidly changing matrices).The computation of P i from P i−1 in this case (rather than P 1 as above) does have the drawback of accumulated approximation errors, however, solving the minimization problem efficiently is a bigger bottleneck for scaling to large problems.
As mentioned in the Introduction, in [31] we have extensively experimented for the parametric case (again, category three earlier) using the first approach.The focus here is to do a similar experimentation for the non-parametric case (first two categories earlier) using the second approach.

C. APPLICATION OF REUSING PRECONDITIONER
Here, we first discuss the application of the above presented theory of reusing preconditioners to the AIRGA algorithm.If we closely observe Algorithm II.1, as mentioned earlier, linear systems are solved at lines 4 and 10.Computation of preconditioners is done only at line 4 because at line 10, matrices do not change, only the right-hand sides do.Hence, we only focus on reusing preconditioners for line 4.
Delving further into the complexity of such linear systems, we observe that the matrices change with the index of outer while loop (line 2) as well as with the index of the for loop corresponding to the expansion points (line 3).Hence, we denote such matrices not only with a subscript as in previous subsection but also with a superscript.That is, A i D + K, where z = 1, . . ., z (until covergence) and i = 1, . . ., .As the matrix A (z) i changes with respect to two different indices, we can reuse preconditioners in many ways.However, here we use the second approach as discussed in the previous subsection.This approach is diagrammatically represented in Figure 1.
( Next, we show how the new preconditioners are computed for both, the horizontal direction and the vertical direction.While looking at the horizontal route, let, i−1 and s (z) i , respectively, with i = 2, . . ., .Using the above theory, we enforce A i in Figure 2. Thus, we eventually enforce A i−1 and solve the minimization problem min

This gives us the new preconditioner
This minimization is again performed for n independent least square problems as in (7).Similar steps are followed for reusing preconditioners along the rest of the horizontal directions, i.e. for all z = 1, . . ., z.Now, applying this technique for the vertical direction, we have for z = 2, . . ., z Following the steps as for the horizontal direction, here, we solve the minimization problem min This gives us the new preconditioner P (z) Again, this is solved as n independent least square problems as in (7).
AIRGA with an efficient implementation of the above discussed theory of reusing preconditioners is given in Algorithm III.1.If we closely look at line 4 of Algorithm II.1, the solution vector is denoted by X (0) (s i ), where the superscript "0" refers to the index of the inner while loop (line 8).We do not bother about this index because, as earlier, matrix does not change inside this inner loop.Rather, we need to capture the change because of the outer while loop indexed with z.Hence, we denote the solution vector as X (z) (s i ) in Algorithm III.1 (lines 8, 11, 19 & 22).It is important to emphasize again that preconditioners are never computed explicitly.Rather, they are obtained using matrix-vector product (please see line numbers 11, 19 & 22 of Algorithm III.1).
For sake of brevity, reusing preconditioners in BIRKA (Algorithm II.2) is discussed as part of Appendix A. Similarly, applying this theory to QB-IHOMM (Algorithm II.3) is discussed in Appendix B.

IV. NUMERICAL EXPERIMENTS
For supporting our proposed preconditioned iterative solver theory using the AIRGA algorithm [15], we perform experiments on two models.The first is a macroscopic equations of motion model (i.e.academic disk brake M 0 ) [38], and is discussed in Section IV-A.The second is also a similar model, however, this is a real-life industrial problem (i.e.industrial disk brake M 1 ) [38].The experiments on this model are discussed in Section IV-B.These models are described by the following set of equations [38]: where .
FIGURE 2: Expressing one linear system matrix in terms of the other.
(Algorithm II.1) for comparison.In Algorithms II.1 and III.1, at line 2 the overall iteration (while-loop) terminates when the change in the reduced model (computed as H 2error between the reduced models at two consecutive AIRGA iterations) is less than a certain tolerance.We take this tolerance as 10 −04 based upon the values in [15].There is one more stopping criteria in Algorithms II.1 at line 8 (also in Algorithm III.1 but not listed here).This checks the H 2error between two temporary reduced models.We take this tolerance as 10 −06 , again based upon the values in [15].Since this is an adaptive algorithm, the optimal size of the reduced model is determined by the algorithm itself, and is denoted by r.
The linear systems that arise here have non-symmetric matrices.There are many iterative methods available for solving such linear systems.We use the Generalized Minimal Residual (GMRES) method [32] because it is very popular [40].The stopping tolerance in GMRES is taken as 10 −06 , which is a common standard.As mentioned in Introduction, for both the given models, we observe that unpreconditioned GMRES fails to converge.Hence, we use the SPAI preconditioner as described above (without and with reuse).We use Modified Sparse Approximate Inverse (MSPAI 1.0) proposed in [39] as our preconditioner.This is because MSPAI uses a linear algebra library for solving sparse least square problems that arise here.We use standard initial settings of MSPAI i.e. tolerance (ep) of 10 −04 .We perform our numerical experiments on a machine with the following configuration: Intel Xeon (R) CPU E5-1620 V3 @ 3.50 GHz., frequency 1200 MHz., 8 CPU and 64 GB RAM.All the codes are written in MATLAB (2016b) (including AIRGA, GMRES) except SPAI and reusable SPAI.MATLAB is used because of ease of rapid prototyping.Computing SPAI and reusable SPAI in MATLAB is expensive, therefore, we use C++ version of these (SPAI is from MSPAI and reusable SPAI is written by us).MSPAI further uses BLAS, LAPACK and ATLAS libraries.Whenever a preconditioner has to be computed, we first compute SPAI and reusable SPAI separately (in-parallel) and save them.Then, we run MATLAB code along with the saved preconditioner matrices (i.e.SPAI and reusable SPAI).

A. ACADEMIC DISK BRAKE MODEL
This model is of size 4, 669.Based upon experience, the maximum reduced system size (r max ) is taken as 20.As mentioned earlier, however, due to the adaptive nature of the AIRGA algorithm, we obtain a reduced system of size r = 13.For this model, the AIRGA algorithm takes two outer iterations (line 2 of Algorithms II.1 and III.1) to converge (i.e.z = 2).
Reusing the SPAI preconditioner is beneficial when the values of I − A (z) i f / I f is large, and the values of f are small, which is true in this case (see Table 3).In this table, columns 1 and 2 list the AIRGA iterations and the four expansion points, respectively.The above three quantities are listed in columns 3, 4 and 5, respectively.For the first AIRGA iteration and the first expansion point, SPAI preconditioner cannot be reused because there is no earlier preconditioner (mentioned as NA in table).From the second expansion point (and the first AIRGA iteration), we perform horizontal reuse of preconditioner (see Figure 1).This is the same for the second AIRGA iteration as well.Vertical reuse of preconditioner is done only for the first expansion point (and the second AIRGA iteration; again see Figure 1).
In Table 4, we compare the SPAI and the reusable SPAI timings.As for Table 3, here columns 1 and 2 list the AIRGA iterations and the four expansion points, respectively.SPAI and reusable SPAI computation times are given in columns 3 and 4, respectively.At the first AIRGA iteration and the first expansion point, both SPAI and reusable SPAI take the same computation time.This is because, as above, reusing of SPAI preconditioner is not applicable here.From the second expansion point of the first AIRGA iteration, we see substantial savings because of the reuse of the SPAI preconditioner (approximately 68%).
Table 5 provides the iteration count and the computation time of GMRES.Here, we only provide GMRES execution Algorithm III.1 : AIRGA with reuse of SPAI preconditioner 1: z = 1 2: while no convergence do for i = 1, . . ., do 5: Compute initial P 1 by solving min i by solving min 1 ]X (1) end for for i = 1, . . ., do 16: 1 by solving min end if   details since the computation time of preconditioner has been discussed above.In this table, column 1 lists the AIRGA iterations.The number of linear solves and average GMRES iterations per linear solve are given in columns 2 and 3, respectively.Finally, columns 4 and 5 list the computation times of GMRES when using SPAI and reusable SPAI, respectively.We notice from this table that solving linear systems by GMRES with SPAI takes less computation time as compared to solving them by GMRES with reusable SPAI.This is because when we reuse the SPAI preconditioner in GMRES, additional matrix-vector products are performed, however, this extra cost is almost negligible when compared to the savings in the preconditioner computation time for the latter case (as evident in Table 3 above; also see total GMRES and preconditioner time below).Table 6 gives the computation time of GMRES plus SPAI (column 2) and GMRES plus reusable SPAI (column 3) at each AIRGA iteration (column 1).As evident from this table, reusing the SPAI preconditioner leads to about 60% savings in total time required for solving all the linear systems.
f are small, which is true in this case (see Table 7).The structure of this table is same as Table 3.As earlier, for the first AIRGA iteration and the first expansion point, SPAI preconditioner cannot be reused because there is no earlier preconditioner (mentioned as NA in table).From the second expansion point (and the first AIRGA iteration), we perform horizontal reuse of preconditioner (see Figure 1).This is the same for the second, the third and the fourth AIRGA iterations as well.Vertical reuse of preconditioner is done only for the first expansion point (and the second, the third, and the fourth AIRGA iterations; again see Figure 1).
In Table 8, we compare the SPAI and the reusable SPAI timings.The structure of this table is same as that of Table 4.As before, at the first AIRGA iteration and the first expansion point, both SPAI and reusable SPAI take the same computation time.This is because, as above, reusing of SPAI preconditioner is not applicable here.From the second expansion point of the first AIRGA iteration, we see substantial savings because of the reuse of the SPAI preconditioner (from 160 hours to 26 hrs 30 minutes; approximately 83%).
Table 9 provides the iteration count and the computation time of GMRES.Here, again we have only provided GMRES execution details since the computation time of the preconditioner has already been discussed above.The structure of this table is same as that of Table 5.As earlier, we notice from this table that solving linear systems by GMRES with SPAI takes less computation time as compared to solving them by GMRES with reusable SPAI.This is again because of additional matrix-vector products in the reusable SPAI case.Here also, this extra cost is almost negligible when compared to the savings in the preconditioner computation time (as evident in Table 8; also see the total GMRES and preconditioner time below).
Table 10 gives the computation time of GMRES plus SPAI (column 2) and GMRES plus reusable SPAI (column 3) at each AIRGA iteration (column 1).As before, it is evident from this table, reusing the SPAI preconditioner leads to about 64% savings in total time (from 197 hours 28 minutes to 72 hours 06 minutes).
To demonstrate the quality of the reduced system, we plot the relative H 2 error between the transfer function of the original system and the reduced system with respect to the different expansion points (in Figure 3).The reduced system considered here is obtained by using GMRES with reusable SPAI.These expansion points, denoted by S, are computed as 2π f , where the frequency variable f is linearly spaced between 1 and 500.As evident from this figure, the obtained reduced system is good (the error is very small).Further, we also observe from this figure that the reduced model is most accurate in 7-10 range of the expansion points.This is because the final expansion points, upon the convergence of the AIRGA algorithm, lie in this range.

V. CONCLUSIONS & FUTURE WORK
In this work, we have focused on MOR of non-parametric dynamical systems, specifically on the following three algorithms: AIRGA, BIRKA, and QB-IHOMM.Since solving  Specifically, we have demonstrated the following: ex-ploitation of the simplicity because of the lack of parameters in reusing preconditioners, multiple ways of reusing preconditioners within the algorithm, efficient implementation to ensure that the savings because of reusing preconditioners are not negated by bad coding, and experimentation on a massively large industrial problem.Numerical experiments show the effectiveness of our approach, where for a problem of size 1.2 million, we save upto 64% in the computation time.In absolute terms, this gives a saving of 5 days.
In future, we plan to explore two directions.The first is use of the randomized preconditioners in solving linear systems arising in MOR.This is giving promising results.The second is to use the spiking neural networks to optimize the parameters inside the preconditioners. .
case of proportionally damped system; as needed for AIRGA) with commonly used parameter values as Ω = 2π, α = 5 × 10 −02 , and β = 5 × 10 −06 .Further, F ∈ R n and C T ∈ R n are taken as [1 0 • • • 0] T , which is the most frequently used choice.We take four expansion points linearly spaced between 1 and 500 based upon experience.Although our purpose is to just reuse SPAI in AIRGA (Algorithm III.1), we also execute original SPAI in AIRGA

FIGURE 3 :
FIGURE 3:  Relative error between the original and reduced system for the industrial disk brake model.

TABLE 3 :
SPAI and reusable SPAI analysis for the academic disk brake model.

TABLE 4 :
SPAI and reusable SPAI computation time for the academic disk brake model.

TABLE 5 :
GMRES computation time for the academic disk brake model.

TABLE 6 :
GMRES with SPAI and reusable SPAI computation time for the academic disk brake model.This model is of size 1.2 million.Based upon experience, the maximum reduced system size (r max ) is taken as 100.As mentioned earlier, however, due to the adaptive nature of the AIRGA algorithm, we obtain a reduced system of size r = 52.For this model, the AIRGA algorithm takes four outer iterations (line 2 of Algorithms II.1 and III.1) to converge (i.e.z = 4).

TABLE 7 :
SPAI and reusable SPAI analysis for the industrial disk brake model.

TABLE 8 :
SPAI and reusable SPAI computation time for the industrial disk brake model.hrs 30 mins § All times given here differ in seconds (not evident because of the rounding to the nearest minute).

TABLE 9 :
GMRES computation time for the industrial disk brake model.

TABLE 10 :
GMRES with SPAI and reusable SPAI computation time for the industrial disk brake model.
large and sparse linear systems is a bottleneck in scaling these MOR algorithms for reduction of large sized dynamical systems, we have proposed reusing of the SPAI preconditioner.