BIT Numerical Mathematics 44:181-188.2004. 181 2004 Kluwer Academic Publishers.Printed in the Netherlands. STABILITY OF THE MATRIX FACTORIZATION FOR SOLVING BLOCK TRIDIAGONAL SYMMETRIC INDEFINITE LINEAR SYSTEMS* JINXI ZHAO',WEIGUO WANG2 and WEIQING REN3.** 1State Key Laboratory for Novel Software Technology,Nanjing University,Nanjing 210093, P.R.China.email:jrzhao@nju.edu.cn 2Department of Mathematics,Ocean University of China,Qingdao,266071.P.R.China 3Department of Mathematics,Nanjing University,Nanjing 210008,P.R.China Abstract. In this paper,we propose a new factorization method for block tridiagonal symmetric indefinite matrices.We also discuss the stability of the factorization method.As a mea- surement of stability,an effective condition number is derived by using backward error analysis and perturbation analysis.It shows that under some suitable assumptions the solution obtained by this factorization method is acceptable.Numerical results demonstrate that the factorization is stable if its condition number is not too large. AMS subject classification (2000):65F05,65F35. Key words:indefinite systems,backward stability,condition number,rounding-error analysis,quasidefinite matrix. 1 Introduction. We are concerned in this paper with the stability of a factorization for the equation (1.1) Bx=b. where B is a block tridiagonal symmetric indefinite matrix.The block tridiagonal matrix is of the form -4 0 (1.2) B -AT -C G 0 GT D/ where K is a symmetric positive definite matrix of order m and A and G have dimensions m x n and n x I with full column rank,C and D are symmetric positive semidefinite.The matrix(1.2)is symmetric indefinite and nonsingular. The special linear system (1.1)has many applications,e.g.,the remaining (linearized)Euler-Lagrange equations can be obtained in the matrix form(1.2) (see[11): Received January 2002.Revised October 2003.Communicated by Axel Ruhe. *The work is partly supported by China NSF grant (60073029)
BIT Numerical Mathematics 44: 181–188, 2004. © 2004 Kluwer Academic Publishers. Printed in the Netherlands. 181 STABILITY OF THE MATRIX FACTORIZATION FOR SOLVING BLOCK TRIDIAGONAL SYMMETRIC INDEFINITE LINEAR SYSTEMS JINXI ZHAO1, WEIGUO WANG2 and WEIQING REN3, 1State Key Laboratory for Novel Software Technology, Nanjing University, Nanjing 210093, P.R.China. email: jxzhao@nju.edu.cn 2Department of Mathematics, Ocean University of China, Qingdao, 266071, P.R.China 3Department of Mathematics, Nanjing University, Nanjing 210008, P.R.China Abstract. In this paper, we propose a new factorization method for block tridiagonal symmetric indefinite matrices. We also discuss the stability of the factorization method. As a measurement of stability, an effective condition number is derived by using backward error analysis and perturbation analysis. It shows that under some suitable assumptions, the solution obtained by this factorization method is acceptable. Numerical results demonstrate that the factorization is stable if its condition number is not too large. AMS subject classification (2000): 65F05, 65F35. Key words: indefinite systems, backward stability, condition number, rounding-error analysis, quasidefinite matrix. 1 Introduction. We are concerned in this paper with the stability of a factorization for the equation (1.1) Bx = b, where B is a block tridiagonal symmetric indefinite matrix. The block tridiagonal matrix is of the form B = K −A 0 −A −C G 0 G D (1.2) , where K is a symmetric positive definite matrix of order m and A and G have dimensions m × n and n × l with full column rank, C and D are symmetric positive semidefinite. The matrix (1.2) is symmetric indefinite and nonsingular. The special linear system (1.1) has many applications, e.g., the remaining (linearized) Euler–Lagrange equations can be obtained in the matrix form (1.2) (see [11]). Received January 2002. Revised October 2003. Communicated by Axel Ruhe. The work is partly supported by China NSF grant (60073029)
182 JINXI ZHAO.WEIGUO WANG AND WEIQING REN In [7],with a coupled DEM-FEM formulation,the necessity to use Lagrange multipliers will induce a three-field mixed system that,in the case of an imper- vious porous material with an incompressible pore fluid,has the form(1.2). Direct methods with partial pivoting strategies for symmetric indefinite ma- trices producing a block diagonal matrix consisting of blocks of order 1 or 2 are given in 1,2.The stability analysis of the Cholesky factorization for quasidefi- nite systems is developed in [4].Using the stability analysis for the factorization LDMT of nonsymmetric positive definite matrix,the conditions are derived under which the Cholesky factorization is stable for quasidefinite systems.Fors- gren,Gill and Shinnerl give a rounding-error analysis of the symmetric indefinite factorization when applied to t-diagonally dominant systems(see [3]). The factorization method was presented in [10,12].This method inherited the advantage of Cholesky factorization with small storage and low computation costs.For this we have the following theorem [10]. THEOREM 1.Given any symmetric indefinite matrir as (1.2),then we have (1.3) B=LJLT, D11 (1.4 L21 0 L32L33 where L11∈Rmxm,L22∈Rnxn,L33∈Rixt are lower triangular matrir,L21∈ Rnxm,L32E Rixn.Im:In and I are identity matrir. Lij can be easily calculated from the following matrix equations: (1.5) K LnLh, (1.6) -A=L21H, (1.7) C+21L团=L22L, (1.8) -GT=L32, (1.9) D+L22=L33 In this paper,we discuss the backward rounding error and the stability of the factorization method(1.3)considering the special structures B and the pertur- bation E sufficiently.We will examine conditions under which the factorization method may be used reliably.In Section 2 we give a backward rounding-error analysis of (1.1)by means of the factorization method.The stability condition of the factorization is described in Section 3.In Section 4 we describe the sensitivity of the solution of this class of systems when the matrix is changed by perturbation.The analysis suggests an effective condition number for these equations and indicates that under suitable assumptions,the solution can be computed accurately.Finally,some numerical results are given to demonstrate the prediction
182 JINXI ZHAO, WEIGUO WANG AND WEIQING REN In [7], with a coupled DEM-FEM formulation, the necessity to use Lagrange multipliers will induce a three-field mixed system that, in the case of an impervious porous material with an incompressible pore fluid, has the form (1.2). Direct methods with partial pivoting strategies for symmetric indefinite matrices producing a block diagonal matrix consisting of blocks of order 1 or 2 are given in [1, 2]. The stability analysis of the Cholesky factorization for quasidefi- nite systems is developed in [4]. Using the stability analysis for the factorization LDM T of nonsymmetric positive definite matrix, the conditions are derived under which the Cholesky factorization is stable for quasidefinite systems. Forsgren, Gill and Shinnerl give a rounding-error analysis of the symmetric indefinite factorization when applied to t-diagonally dominant systems (see [3]). The factorization method was presented in [10, 12]. This method inherited the advantage of Cholesky factorization with small storage and low computation costs. For this we have the following theorem [10]. Theorem 1. Given any symmetric indefinite matrix as (1.2), then we have B = LJLT (1.3) , L = L11 0 L21 L22 0 L32 L33 , J = Im −In Il (1.4) , where L11 ∈ Rm×m, L22 ∈ Rn×n, L33 ∈ Rl×l are lower triangular matrix, L21 ∈ Rn×m, L32 ∈ Rl×n. Im, In and Il are identity matrix. Lij can be easily calculated from the following matrix equations: K = L11LT 11 (1.5) , −AT = L21LT 11 (1.6) , C + L21LT 21 = L22LT 22 (1.7) , −GT = L32LT 22 (1.8) , D + L32LT 32 = L33LT 33 (1.9) . In this paper, we discuss the backward rounding error and the stability of the factorization method (1.3) considering the special structures B and the perturbation E sufficiently. We will examine conditions under which the factorization method may be used reliably. In Section 2 we give a backward rounding-error analysis of (1.1) by means of the factorization method. The stability condition of the factorization is described in Section 3. In Section 4 we describe the sensitivity of the solution of this class of systems when the matrix is changed by perturbation. The analysis suggests an effective condition number for these equations and indicates that under suitable assumptions, the solution can be computed accurately. Finally, some numerical results are given to demonstrate the prediction
STABILITY OF THE MATRIX FACTORIZATION 183 2 Backward rounding error analysis In this section we give a backward error analysis of the solution of (1.1)by means of the factorization method.Throughout,we use the "standard model" of floating-point arithmetic in which the evaluation of an expression in floating- point arithmetic is denoted by fl(),with f(aob)=(aob)(1+6),l6川≤u,o=+,-,*,/ and f(x2/2)=x2/2(1+e),le≤1.00001u (see,for example,Higham 6).Here u is the unit round-off associated with the particular machine being used. In the next lemma,the backward error analysis is presented. LEMMA 2.Let B be the same as that in (1.2)and assume max{m +3, n+3u<0.01.L is the computed factorization factor of B.Then we have B+△B=iJiT, where△B satisfies (2.1) Lm+71.01创i. 1△B1≤1-3.00002u PRooF.Let denote the subblock of L in analogy with that in (1.4),which are computed from Equations(1.5)~(1.9)respectively.The backward error of the factorization is accumulated from each of the five steps.First,let K+△K=iH then we have the following inequality as it in [8], (2.2) 1aK1≤2at2L2 Similar results can be established for the computable factors L21 and L22: -(AT+△AT)=21h, (2.3) 1.01(m+3严1i2话l 1△A1≤1-2.02u and C+△C+i21话=i22五 (2.4) 1AC≤maxn+7,n+o1010i21li1+1i2liD. 1-3.00002u
STABILITY OF THE MATRIX FACTORIZATION 183 2 Backward rounding error analysis. In this section we give a backward error analysis of the solution of (1.1) by means of the factorization method. Throughout, we use the “standard model” of floating-point arithmetic in which the evaluation of an expression in floatingpoint arithmetic is denoted by fl(·), with fl(a ◦ b)=(a ◦ b)(1 + δ), |δ| ≤ u, ◦ = +, −, ∗, / and fl(x1/2) = x1/2(1 + ε), |ε| ≤ 1.00001u (see, for example, Higham [6]). Here u is the unit round-off associated with the particular machine being used. In the next lemma, the backward error analysis is presented. Lemma 2. Let B be the same as that in (1.2) and assume max{m + 3, n + 3}u ≤ 0.01. Lˆ is the computed factorization factor of B. Then we have B + ∆B = LJˆ LˆT, where ∆B satisfies |∆B| ≤ (m + 7)1.01u 1 − 3.00002u |Lˆ||LˆT (2.1) |. Proof. Let Lˆij denote the subblock of Lˆ in analogy with that in (1.4), which are computed from Equations (1.5)∼(1.9) respectively. The backward error of the factorization is accumulated from each of the five steps. First, let K + ∆K = Lˆ11LˆT 11 then we have the following inequality as it in [8], |∆K| ≤ 1.01(m + 2)u 1 − 3.00002u |Lˆ11||LˆT 11 (2.2) |. Similar results can be established for the computable factors Lˆ21 and Lˆ22: −(AT + ∆AT) = Lˆ21LˆT 11, |∆AT| ≤ 1.01(m + 3)u 1 − 2.02u |Lˆ21||LˆT 11| (2.3) and C + ∆C + Lˆ21LˆT 21 = Lˆ22LˆT 22, |∆C| ≤ max{m + 7, n + 6}1.01u 1 − 3.00002u (|Lˆ21||LˆT 21| + |Lˆ22||LˆT 22|). (2.4)
184 JINXI ZHAO,WEIGUO WANG AND WEIQING REN Finally,we establish for the computed factors L32 and L33.Let L22=vj],L32= [ijl,L33=[hijl,from (1.8)and (1.9),the uij and hij are computed from 1一1 (2.5) 4=-9i- yi,i=1,2,,l,j=1,2,.,n D- and i=j, (2.6) D+∑=1a-h品P, (Dg+∑=hk4k-∑h/h,i>i, i=1,2,,l,j=1,2,,i. We obtain (2.7) 1△c91≤101m+3业ie场 1-2.02u and (2.8) 1ADl≤maxn+7l+6101Ii2li语1+1ial2原D. 1-3.00002u Then, B+△B=iJiT From(2.2)~(2.8),we have (2.9) 1△B≤maxn+7+7,U+61,01 QEETD. 1-3.000024 For A and G are full column rank,i.e.,m >n >1.(2.1)is derived immediately from(2.9). 0 We consider the backward error resulting from solution of the triangular systems, Lx=f and Uy=g, where L E Rnxn is lower triangular and U E Rnxn upper triangular.It follows from [9]that the intermediate vectors i and i satisfy (L+AL)i=f and (0+△U)i=g,where△Land△U have the same element-.wise bound (2.10) a≤自d1a≤n In the next theorem we show that the computed solution of(1.1)is the exact solution of(B+E)=b,where E is an error matrix.The theorem is established by accumulating the backward error. THEOREM 3.Let B be as in (1.2),and i is the computed solution of Br=b by the factorization.Then i is the exact solution of (B+E)i=b with 3m+n+)ui21 1E到≤1-(m+n+0m
184 JINXI ZHAO, WEIGUO WANG AND WEIQING REN Finally, we establish for the computed factors Lˆ32 and Lˆ33. Let L22 = [νij ], L32 = [µij ], L33 = [hij ], from (1.8) and (1.9), the µij and hij are computed from µij = − gji − j−1 p=1 µipνjp (2.5) /νjj , i = 1, 2, . . .,l, j = 1, 2, . . . , n. and hij = (Dii + n k=1 µ2 ik − j−1 p=1 h2 ip)1/2, i = j, (Dij + n k=1 µikµjk − j−i p=1 hiphjp)/hjj , i > j, (2.6) i = 1, 2, . . .,l, j = 1, 2, . . ., i. We obtain |∆GT| ≤ 1.01(m + 3)u 1 − 2.02u |Lˆ32||LˆT 22 (2.7) | and |∆D| ≤ max{n + 7, l + 6}1.01u 1 − 3.00002u (|Lˆ32||LˆT 32| + |Lˆ33||LˆT 33 (2.8) |). Then, B + ∆B = LJˆ LˆT. From (2.2)∼(2.8), we have |∆B| ≤ max{n + 7, l + 7, l + 6}1.01u 1 − 3.00002u (|Lˆ||LˆT (2.9) |). For A and G are full column rank, i.e., m ≥ n ≥ l. (2.1) is derived immediately from (2.9). ✷ We consider the backward error resulting from solution of the triangular systems, Lx = f and Uy = g, where L ∈ Rn×n is lower triangular and U ∈ Rn×n upper triangular. It follows from [9] that the intermediate vectors ˆx and ˆy satisfy (Lˆ + ∆L)ˆx = f and (Uˆ + ∆U)ˆy = g, where ∆L and ∆U have the same element-wise bound |∆L| ≤ nu 1 − nu|Lˆ| and |∆U| ≤ nu 1 − nu|Uˆ (2.10) |. In the next theorem we show that the computed solution ˆx of (1.1) is the exact solution of (B +E)ˆx = b, where E is an error matrix. The theorem is established by accumulating the backward error. Theorem 3. Let B be as in (1.2), and xˆ is the computed solution of Bx = b by the factorization. Then xˆ is the exact solution of (B + E)ˆx = b with |E| ≤ 3(m + n + l)u 1 − (m + n + l)u |Lˆ||LˆT|
STABILITY OF THE MATRIX FACTORIZATION 185 PRooF.Let LJLT be the computed factorization of B.In the following we denoteILT as U,an upper triangular matrix.To solve Br=b,we must solve two triangular systems.Let the computed solution Y and i satisfy (L+AL)=b and(d+△U):=y,with bounds for△Land△given in(2.l0).Hence主 satisfies (i+△L)(0+△U)定=b, i.e., (B+Ei=6, where E=△B+△LU+L△U+△L△U Ignoring elements of order u2,from(2.1)and(2.10),we have IE≤I△B+I△LI01+I△川+|△L△U川 ≤ m+701.0ui0+1-(m+n+0u 1-3.00002u 2(m+n+0u1创0 3m+n+0ui21 ≤1-(m+n+)u In the last inequality,we assume (m+n+)≥1.01(m+7). 口 In the preceding Theorem 3,the error E is a symmetric error matrix with the special structure as B. 3 Stability of the factorization. In the backward error analysis of solving Bx=b,it is shown that the computed solution i is the exact solution of the perturbed system(B+E)i=b,where the size of E is bounded by an expression involving the size of the computed factor L. Algorithms that produce L of sufficiently bounded size are therefore considered stable.Based on the factorization LDMT for nonsymmetric positive definite matrix H,the stability analysis of the Cholesky factorization for symmetric quasidefinite matrix is given (see [4]).Here the nonsymmetric positive definite matrix is and the symmetric quasidefinite matrix is H=()=(e)(-） In the following,we consider directly the stability of the factorization as(1.2) without using the factorization LDMT
STABILITY OF THE MATRIX FACTORIZATION 185 Proof. Let LJˆ LˆT be the computed factorization of B. In the following we denote JLˆT as Uˆ, an upper triangular matrix. To solve Bx = b, we must solve two triangular systems. Let the computed solution Yˆ and ˆx satisfy (Lˆ + ∆L)ˆy = b and (Uˆ + ∆U)ˆx = y, with bounds for ∆L and ∆U given in (2.10). Hence ˆx satisfies (Lˆ + ∆L)(Uˆ + ∆U)ˆx = b, i.e., (B + E)ˆx = b, where E = ∆B + ∆LUˆ + Lˆ∆U + ∆L∆U. Ignoring elements of order u2, from (2.1) and (2.10), we have |E|≤|∆B| + |∆L||Uˆ| + |Lˆ||∆U| + |∆L||∆U| ≤ (m + 7)1.01u 1 − 3.00002u |Lˆ||Uˆ| + 2(m + n + l)u 1 − (m + n + l)u |Lˆ||Uˆ| ≤ 3(m + n + l)u 1 − (m + n + l)u |Lˆ||LˆT|. In the last inequality, we assume (m + n + l) ≥ 1.01(m + 7). ✷ In the preceding Theorem 3, the error E is a symmetric error matrix with the special structure as B. 3 Stability of the factorization. In the backward error analysis of solving Bx = b, it is shown that the computed solution ˆx is the exact solution of the perturbed system (B +E)ˆx = b, where the size of E is bounded by an expression involving the size of the computed factor Lˆ. Algorithms that produce Lˆ of sufficiently bounded size are therefore considered stable. Based on the factorization LDM T for nonsymmetric positive definite matrix H¯ , the stability analysis of the Cholesky factorization for symmetric quasidefinite matrix is given (see [4]). Here the nonsymmetric positive definite matrix is H¯ = K A −A C and the symmetric quasidefinite matrix is H = K −A −A −C = K A −A C In −Im . In the following, we consider directly the stability of the factorization as (1.2) without using the factorization LDM T
186 JINXI ZHAO.WEIGUO WANG AND WEIQING REN ASSUMPTION 4.For some scalar y of moderate size, IIZIIETIE≤lILILTIE. From Equation (1.5)~(1.9),we have =tr(K), L2=tr(ATK-1A), lL22ll=tr(ATK-1A+C). L32l=tr(GT(ATK-1A+C)-1G), L33=tr(GT(ATK-1A+C)-1G+D). Hence IIZIILTIE≤LIFLTIIE tr(K)+tr(C)+tr(D)+2tr(ATK-1A) +2tr(GT(ATK-1A+C)-G) =(1+w(B))(tr(K)+tr(C)+tr(D)), where (3.1) w(B)= 2tr(ATK-1A)+2tr(GT(ATK-1A+C)-G) tr(K)+tr(C)+tr(D) Combining Theorem 3 and Assumption 4,we have a2≤Br≤2Cr+国Xu(W+O+D THEOREM 5.If B is as in (1.2),the factorization B=LJLT is stable if w(B)is not too large. 4 The condition number of the augmented system. Let the computed solution i of Br=b satisfy the perturbed system (B+E)x=6. Then the usual sensitivity bound takes the form where (B). B\2 From Assumption 4 and(3.2),we have (4.1) 3(m+n+)utr()+r(C)+(D（1+w(B)K2(B). a≤1-(m+n+0)u B2 A simple calculation shows that (4.2) tr(K)+tr(C)+tr(D)≤(m+n+I)川Bll2. From (4.1)and (4.2),we obtain the following result showing that the relative error is bounded by quantity involving (1+w(B))K2(B)
186 JINXI ZHAO, WEIGUO WANG AND WEIQING REN Assumption 4. For some scalar γ of moderate size, |Lˆ||LˆT|F ≤ γ|L||LT|F . From Equation (1.5)∼(1.9), we have L112 F = tr(K), L212 F = tr(ATK−1A), L222 F = tr(ATK−1A + C), L322 F = tr(GT(ATK−1A + C) −1G), L332 F = tr(GT(ATK−1A + C) −1G + D). Hence |L||LT|F ≤ LF LTF = tr(K) + tr(C) + tr(D) + 2tr(ATK−1A) + 2tr(GT(ATK−1A + C) −1G) = (1 + ω(B))(tr(K) + tr(C) + tr(D)), where ω(B) = 2tr(ATK−1A) + 2tr(GT(ATK−1A + C)−1G) tr(K) + tr(C) + tr(D) (3.1) . Combining Theorem 3 and Assumption 4, we have E2 ≤ EF ≤ 3γ(m + n + l)u 1 − (m + n + l) (3.2) (1 + ω(B))(tr(K) + tr(C) + tr(D)). Theorem 5. If B is as in (1.2), the factorization B = LJLT is stable if ω(B) is not too large. 4 The condition number of the augmented system. Let the computed solution ˆx of Bx = b satisfy the perturbed system (B + E)x = b. Then the usual sensitivity bound takes the form x − xˆ x2 ≤ α 1 − α, where α = E2 B2 κ2(B). From Assumption 4 and (3.2), we have α ≤ 3γ(m + n + l)u 1 − (m + n + l)u tr(K) + tr(C) + tr(D) B2 (4.1) (1 + ω(B))κ2(B). A simple calculation shows that (4.2) tr(K) + tr(C) + tr(D) ≤ (m + n + l)B2. From (4.1) and (4.2), we obtain the following result showing that the relative error is bounded by quantity involving (1 + ω(B))κ2(B)
STABILITY OF THE MATRIX FACTORIZATION 187 THEOREM 6.Let B be symmetric indefinite as (1.2),and i is the computed solution of Bx =b,then s。 with 3y(m+n+1)2u a≤-m+n+(B), where (B)=(1+w(B))K2(B),w(B)is defined in (3.1). The theorem tells us that under Assumption 4,system (1.1)can be solved accurately as long as (B)is not too large.We therefore interpret (B)to be the effective condition number of generalized Cholesky factorization for solving the system (1.1). 5 Numerical experiments. We have used MATLAB in PC to implement the factorization of B=b, where B has the form(1.2).The results confirm our prediction. ExAMPLE 1.Let m=n=10,1=5;K diag(ki),k11=E,k22,...,kmm be positive and randomly generated by MATLAB,A and G have full column rank, C and D be zero matrices.We list some results corresponding to the example in Table 5.1. EXAMPLE 2.Let K,A,G be as in Example 1;let C and D be symmetric positive semidefinite random matrices.The results listed in Table 5.2. The results show that for various e,though the spectral condition number of B is almost invariant,the precision of the computed solution is different.This confirms our predictions:(i)for small w(B),the factorization is stable;(ii)for small (B),the computed solution i is reliable. Note that small w(B)and (B)are both sufficient,but not necessary,for ensuring an accurate factorization factor and the reliable solution of Br =b. For large w(B)and (B),it is difficult to draw any conclusion. Table 5.1:Numerical results for the Example 1. 2(B) w(B) (B) B-LJLTIF Iz-i2 2 102 7.1860e+02 6.4539e+00 5.3564e+03 1.6834e-14 9.0382e-14 100 4.2444e+01 3.3010e+02 1.4435e+03 4.7234e-15 5.0320e-15 10-2 4.2236e+01 2.3321e+02 9.8921e+03 3.3934e-14 3.8136e-14 10-4 4.2241e+01 1.9735e+04 8.3365e+05 3.0106e-12 1.9367e-12 10-6 4.2241e+01 1.9699e+06 8.3209e+07 2.8257e-10 1.8954e-10 108 4.2241e+01 1.9698e+08 8.3207e+09 2.7447e-08 2.2862e-08
STABILITY OF THE MATRIX FACTORIZATION 187 Theorem 6. Let B be symmetric indefinite as (1.2), and xˆ is the computed solution of Bx = b, then x − xˆ x2 ≤ α 1 − α, with α ≤ 3γ(m + n + l)2u 1 − (m + n + l)u φ(B), where φ(B) = (1 + ω(B))κ2(B), ω(B) is defined in (3.1). The theorem tells us that under Assumption 4, system (1.1) can be solved accurately as long as φ(B) is not too large. We therefore interpret φ(B) to be the effective condition number of generalized Cholesky factorization for solving the system (1.1). 5 Numerical experiments. We have used MATLAB in PC to implement the factorization of Bx = b, where B has the form (1.2). The results confirm our prediction. Example 1. Let m = n = 10, l = 5; K = diag(kii), k11 = ε, k22,...,kmm be positive and randomly generated by MATLAB, A and G have full column rank, C and D be zero matrices. We list some results corresponding to the example in Table 5.1. Example 2. Let K, A, G be as in Example 1; let C and D be symmetric positive semidefinite random matrices. The results listed in Table 5.2. The results show that for various ε, though the spectral condition number of B is almost invariant, the precision of the computed solution is different. This confirms our predictions: (i) for small ω(B), the factorization is stable; (ii) for small φ(B), the computed solution ˆx is reliable. Note that small ω(B) and φ(B) are both sufficient, but not necessary, for ensuring an accurate factorization factor and the reliable solution of Bx = b. For large ω(B) and φ(B), it is difficult to draw any conclusion. Table 5.1: Numerical results for the Example 1. ε κ2(B) ω(B) φ(B) B − LJLTF x−xˆ2 x2 102 7.1860e+02 6.4539e+00 5.3564e+03 1.6834e−14 9.0382e−14 100 4.2444e+01 3.3010e+02 1.4435e+03 4.7234e−15 5.0320e−15 10−2 4.2236e+01 2.3321e+02 9.8921e+03 3.3934e−14 3.8136e−14 10−4 4.2241e+01 1.9735e+04 8.3365e+05 3.0106e−12 1.9367e−12 10−6 4.2241e+01 1.9699e+06 8.3209e+07 2.8257e−10 1.8954e−10 10−8 4.2241e+01 1.9698e+08 8.3207e+09 2.7447e−08 2.2862e−08
188 JINXI ZHAO.WEIGUO WANG AND WEIQING REN Table 5.2:Numerical results for the Example 2. K2(B) w(B) φ(B) B-LJLTIF 带 10 1.6485e+02 2.2222e+00 5.3116e+02 5.0286e-15 4.7348e-15 100 1.6504e+02 2.8013e+00 6.2736e+02 6.3152e-15 4.0532e-15 10-2 1.6515e+02 2.3192e+01 3.9954e+03 4.5329e-14 3.1294e-14 10-4 1.6516e+02 2.0582e+03 3.4009e+05 5.7003e-12 2.3596e-12 10-6 1.6516e+02 2.0556e+05 3.3939e+07 4.5493e-10 2.7224e-10 10-8 1.6516e+02 2.0556e+07 3.3949e+09 3.8521e-08 3.3042e-08 Acknowledgement. The authors are indebted to Prof.Ake Bjorck and Prof.Axel Ruhe for their comments which improved the formulation and interpretation of several re- sults.Also the authors wish to thank the anonymous referees for their valuable suggestions that improved the paper. REFERENCES 1.J.R.Bunch and B.N.Parlett,Direct methods for solving symmetric indefinite systems of linear equations,SIAM J.Numer.Anal.,8 (1971),pp.639-655. 2.J.R.Bunch,Partial pivoting strategies for symmetric matrices,SIAM J.Numer.Anal, 11(1974),pp.521-528. 3.A.Forsgren,P.E.Gill,and J.R.Shinnerl,Stability of symmetric ill-conditioned systems arising in interior methods for constrained optimization,SIAM J.Matrix Anal.Appl., 17(1996),Pp.187-211. 4.P.E.Gill.M.A.Saunders,and J.R.Shinnerl,On the stability of Cholesky factorization for symmetric quasidefinite systems,SIAM J.Matrix Anal.Appl.,17 (1996),pp.35-46. G.H.Golub and C.F.Van Loan,Matriz Computation,3rd ed.,The John Hopkins University Press,Baltimore,1996. 6. N.J.Higham,Accuracy and Stability of Numerical Algorithms,SIAM Publications, Philadelphia,PA,1996. 7.H.Modaressi.A diffuse element-finite element technique for transient coupled analysis. Int.J.Numer.Methods Eng.,39 (1996),pp.3809-3838. Wei-qing Ren,The solution of augmented systems,manuscript,1996 9.J.Stoer and R.Bulirsch,Introduction to Numerical Analysis,2nd ed.,Springer-Verlag New York,1993. 10. Weiguo Wang,A matrir factorization method for solving block tridiagonal symmetric indefinite linear systems,J.Nanjing Univ.Math.Biqu.,14 (1997). 11.S.L.Weissman,High-accuraty low-order three-dimensional brick elements.Int.J.Numer. Methods Eng.,39(1996),pp.2337-2361. 12.Jinxi Zhao,The generalized Cholesky factorization method for saddle point problem Appl.Math.Comput.,92 (1998),pp.49-58
188 JINXI ZHAO, WEIGUO WANG AND WEIQING REN Table 5.2: Numerical results for the Example 2. ε κ2(B) ω(B) φ(B) B − LJLTF x−xˆ2 x2 10 1.6485e+02 2.2222e+00 5.3116e+02 5.0286e−15 4.7348e−15 100 1.6504e+02 2.8013e+00 6.2736e+02 6.3152e−15 4.0532e−15 10−2 1.6515e+02 2.3192e+01 3.9954e+03 4.5329e−14 3.1294e−14 10−4 1.6516e+02 2.0582e+03 3.4009e+05 5.7003e−12 2.3596e−12 10−6 1.6516e+02 2.0556e+05 3.3939e+07 4.5493e−10 2.7224e−10 10−8 1.6516e+02 2.0556e+07 3.3949e+09 3.8521e−08 3.3042e−08 Acknowledgement. The authors are indebted to Prof. ˚Ake Bj¨orck and Prof. Axel Ruhe for their comments which improved the formulation and interpretation of several results. Also the authors wish to thank the anonymous referees for their valuable suggestions that improved the paper. REFERENCES 1. J. R. Bunch and B. N. Parlett, Direct methods for solving symmetric indefinite systems of linear equations, SIAM J. Numer. Anal., 8 (1971), pp. 639–655. 2. J. R. Bunch, Partial pivoting strategies for symmetric matrices, SIAM J. Numer. Anal., 11 (1974), pp. 521–528. 3. A. Forsgren, P. E. Gill, and J. R. Shinnerl, Stability of symmetric ill-conditioned systems arising in interior methods for constrained optimization, SIAM J. Matrix Anal. Appl., 17 (1996), pp. 187–211. 4. P. E. Gill, M. A. Saunders, and J. R. Shinnerl, On the stability of Cholesky factorization for symmetric quasidefinite systems, SIAM J. Matrix Anal. Appl., 17 (1996), pp. 35–46. 5. G. H. Golub and C. F. Van Loan, Matrix Computation, 3rd ed., The John Hopkins University Press, Baltimore, 1996. 6. N. J. Higham, Accuracy and Stability of Numerical Algorithms, SIAM Publications, Philadelphia, PA, 1996. 7. H. Modaressi, A diffuse element-finite element technique for transient coupled analysis, Int. J. Numer. Methods Eng., 39 (1996), pp. 3809–3838. 8. Wei-qing Ren, The solution of augmented systems, manuscript, 1996. 9. J. Stoer and R. Bulirsch, Introduction to Numerical Analysis, 2nd ed., Springer-Verlag, New York, 1993. 10. Weiguo Wang, A matrix factorization method for solving block tridiagonal symmetric indefinite linear systems, J. Nanjing Univ. Math. Biqu., 14 (1997). 11. S. L. Weissman, High-accuraty low-order three-dimensional brick elements, Int. J. Numer. Methods Eng., 39 (1996), pp. 2337–2361. 12. Jinxi Zhao, The generalized Cholesky factorization method for saddle point problem, Appl. Math. Comput., 92 (1998), pp. 49–58
