NASA Contractor Report

ICASE Report No. 89-37

181860

ICASE

BUNCH-KAUFMAN FACTORIZATION FOR REAL

SYMMETRIC INDEFINITE BANDED MATRICES

Mark T. Jones

Merrell L. Patrick

(N&S&-CR-181860| BU WCH--KAUFmAN

FACTOHIZATION FOR REAL SIB_BTRIC IgDE/IgITE

OA_DED M&TRICES Final Repoct (IC&SE) 15 p

CSCL 12&

G3/6q

N89-283_1

Unclas

022q027

Contract No. NAS1-18605

May 1989

Institute for Computer Applications in Science and Engineering

NASA Langley Research Center

Hampton, Virginia 23665-5225

Operated by the Universities Space Research Association

IW A

National Aeronaulics and

Space Adminislration

Langley Research Center

Hampton, Virginia 23665-5225

Recently, ICASE has begun differentiating between reports with a mathematical or applied science theme and reports whose main emphasis is some aspect of

computer science by producing the computer science reports with a yellow cover.

The blue cover reports will now emphasize mathematical research. In all other

aspects the reports will remain the same; in particular, they will continue to be

submitted to the appropriate journals or conferences for formal publication.

Bunch-Kaufman Factorization for Real

Symmetric Indefinite Banded Matrices

Mark T. Jones*and Merrell L. Patrick *t

Abstract

The Bunch-Kaufman algorithm for factoring symmetric indefinite

matrices has been rejected for banded matrices because it destroys the

banded structure of the matrix. Herein, it is shown that for a subclass of real symmetric matrices which arise in solving the generalized

eigenvalue problem using Lanczos's method, the Bunch-Kaufman algorithm does not result in major destruction of the bandwidth. Space

time complexities of the algorithm are given and used to show that

the Bunch-Kaufman algorithm is a significant improvement over LU

factorization.

*Department of Computer Science, Duke University, Durham, NC 27706

tThis research was supported by the National Aeronautics and Space Administration

under NASA contract Nos. NASI-18107 and NAS1-18605 and the Air Force Office of

Scientific Research under AFOSR grant No. 88-0117 while the authors were in residence

at ICASE. Additional support was provided by NASA grant No. NAG-I-466.
1 Introduction

The Bunch-Kaufman algorithm is considered one of the best methods for

factoring full, symmetric, indefinite matrices [BK77], [BG76]. It has also

been modified and successfully used to factor sparse matrices [DRMN79].

However, to date it has been rejected for banded, symmetric indefinite matrices because it destroys the banded structure of the matrix [BK77]. Herein

it is shown that for a subclass of real symmetric indefinite matrices, which

arise in solving the generalized eigenvalue problem using Lanczos's method,

the Bunch-Kaufman algorithm does not result in major destruction of the

bandwidth. Furthermore, for our class of problems, the Bunch-Kaufman

factorization algorithm is a significant improvement over LU factorization,

the standard of comparison for such methods [BK77]. In addition to taking

advantage of symmetry, the Bunch-Kaufman algorithm yields the inertia of

the matrix essentially for free [BK77], which is important in eigenvalue calculations. LU factorization does not yield the inertia as a by-product and

destroys the symmetry of the matrix, thus increasing storage requirements

for its implementation.

In section 2 we give one of the several variations of the Bunch-Kaufman

algorithm and in section 3 describe a subclass of matrices to which we apply

it. An efficient implementation of the method is described in section 4 and

the space/time complexity of the implementation is disussed in section 5.

Conclusions are drawn in section 6.

2 The Bunch-Kaufman Algorithm

The Bunch-Kaufman algorithm factors A, an r_ × n real symmetric indefinite

matrix, into LDL T while doing symmetric permutations on A to maintain

stability, resulting in the following equation:

PAR T = LDL T. (1)

Although several variations of the algorithm exist, algorithm D of the

Bunch-Kaufman paper is the least destructive of the banded structure

[BK77]. The algorithm is shown in figure 1.
1) for i---- 1, n

begin

2) if the previous step was not a 2x2 pivot then

begin

3) ._ = maxi=,+l,. ] ai,i]

4) set r to the row number of this value

5) if_a < ] a_,_] then

begin

6) perform a lxl pivot

end

else

begin

7) cr = maxj=,+l,n [ a.,i ]

8) if a.k _ < a ] ai,_l then

begin

9) perform a lxl pivot

end

else

begin

10) exchange rows and columns r and i + 1

11) perform a 2x2 pivot

end

end

end

12) end

13) if inertia is desired, then scan the D matrix

Figure 1: The Bunch-Kaufman Factorization Algorithm
The parameter, a, is chosen such that stability is maximized and has

been shown by Bunch and Kaufman to be approximately 0.525 [BK77].

The exchange of rows and columns in step 10 maintains the symmetry of

the matrix, unlike LU factorization which destroys the symmetry of the

matrix by permuting only rows.

3 Applicable Set of Matrices

Bunch and Kaufman show that, in general, if ra is the semi-bandwidth of

a matrix being factored, then a 2x2 pivot can increase the semi-bandwidth

from rn to (2rn) - 2 and that this can happen at every step thus resulting

in the complete destruction of the band structure due to fill-in outside the

band [BK77]. However, it will be shown in section four that for a subclass

of matrices this fill-in can be controlled. The number of 2x2 pivots is

bounded above by the number of negative eigenvalues of A, because each

2x2 pivot represents a positive-negative eigenvalue pair [BK77]. Also, the

increase of the semi-bandwidth from rn to (2rn) - 2 is a worst case that

in practice is not likely to occur. Therefore, for matrices with a small

number of negative eigenvalues (in relation to the size of the matrix), it is

possible to use Bunch-Kaufman factorization with very little fill-in. Such

matrices arise in eigenvalue calculations where the smallest eigenvalues are

sought. Methods such as inverse iteration and Lanczos's method are often

used to find the smallest eigenvalues of a matrix, A. To do so, they often

require the factorization of a matrix, (A - a/), where a is normally very

near the left end of A's spectrum, but may not be to the left of the smallest

eigenvalue, thus the matrix is indefinite [NOPT83] but has only a small

number of negative eigenvalues. These matrices can be banded, as they

are in structural mechanics [BH87]. The difficulty is that the location and

amount of the fill-in outside the band is not possible to predict a priori.

In the following section, a detailed algorithm which dynamically allows for

fill-in during factorization will be presented.

4 Efficient Implementation of the Algorithm
• • • • • X _ row _ _ • • • • • x

0 • • • • x x _rowk+l_ 0 • • • • xx

0 0 • • • XXX 0 0 • • • XX X

O00o • xxxx O00o* xxxx

0000- xxxxx O000o xxxxx

00000 x x x x x x _rowr_ 00000 x x x x x x

O00000xxxxxx O0000[xxxxxx

O000000xxxxxx O0000ffxxxxxx

O0000000xxxxxx O0000fffxxxxxx

O00000000xxxxxx O0000ffffxxxxxx

O000000000xxxxxx O0000fffffxxxxxx

O0000000000xxxxxx O0000000000xxxxxx

Figure 2: Example of Fill-in (Note: this is an example of worst case fill-in)

As the following algorithm is executed the original matrix is copied,

piece-wise, from one place in memory to another. This allows for dynamic

allocation of fill-in as well as only requiring part of the matrix to be in main

memory at any particular time. Fill-in only takes place in a small triangle

when a 2x2 pivot occurs. If a pivot occurs at step k, this triangle is of

the form shown in figure 2, where *'s represent eliminated elements in L,

x's are uneliminated non-zeros, O's are zeros outside the band for which no

storage is needed, and fs are areas where fill occurs. The triangle of fill is

from row r + 1 to row r + rn, where m is the semi-bandwidth (this area may

already contain non-zeros depending on the value of r, so no extra memory

may be needed). The algorithm is as follows:

0) set upto to 0

1) fori= 1, n

begin

2) if the previous step was not a 2x2 pivot then

begin

3) read rows upto to min(n,i + rn) of the matrix A into L,
4)

5)

6)

7)

S)

9)

C

C

C

C

10)

c

11)

12)

13)

14)

15)

16)

17)

18)

19)

20)

no extra space for fill needs to be added for these rows

set upto to min(n,i + m)

)t --- max$=i+l,upto I aj,i I

set r to the row number of this value

if,_a <1 ai,i I then

begin

go to 11

end

tr = maxj=i+l,upto l a,,i I

(it may be necessary to access some elements that are not

read in at this point, but the number of elements is small,

so they may be read into L or simply discarded,

this is only a concern if i/o is taking place)

ifa)_2 <_ a [ al,i I then

begin

perform a lxl pivot

set pl = i

set dl,i = ai,i

set di,i+l = 0.0

set ai, i = 1.0

for j = i+ 1, upto

begin

Vj = aj,i

vlj = aj,i/ai,i

aj,i = vii

end

for j = i+ 1, upto

begin

fork=i+l,j

begin

6
21)

C

22)

23)

24)

25)

26)

27)

28)

29)

30)

31)

32)

33)

34)

35)

36)

37)

3s)

39)

40)

aj,k = ai,k -vljvk

end

end

end

else

begin

permute the matrix and then perform a 2x2 pivot

read rows upto to min(n,r + m) of the matrix A into L,

and allocate space for the fill triangle

set upto to min(n,r + rn)

exchange rows and columns r and i + 1

set p_ = i

set Pi+l = r

set di,i = ai,i

set di+l.i+l = ai+l,i+l

set di,i+l = ai+l,i

set di+,,i+_ = 0.0

set determinant = (((d,,,d,+l,,+l)/di,,+,) - d,,,+,)d,,,+l

for j= i+2, upto

begin

v I - aj,i

v2 i = aj,i+l

vii = aj,idi+l,i+,- ai,i+ldi,i+l

v121 = -a£idi,i+, + aj,i+ldi,i

ai, _ = vlj

ai,i+ * = vl21

end

for j= i+2, upto

begin

for k = i + 2, j
begin

41) a_,_ = _,_ - (vt_v_ + vl2_v2_)

end

end

end

end

end

P is a vector representing the permutation matrices. The only time fill

outside the band occurs is in step 24 of the algorithm when a 2x2 pivot

occurs and then storage for the fill is allocated dynamically.

5 Speed and Storage Analysis

In this section we compare the space/time requirements of our implementation of the Bunch-Kaufman algorithm with LU factorization. The storage

requirements for both algorithms will be analyzed for two different situations: 1) when simply factoring a matrix that falls in the subclass described

in section 2, and 2) when factoring a matrix pencil such as (K- crM) where

K and M are symmetric, K is positive definite and o is near the left end

of K's spectrum.

In the first situation, the storage required by the algorithm presented

in section 4 is significantly less than that required by LU factorization for

the set of matrices that was described in section 2. The storage required

by LU factorization is approximately 3ran [BK77]. The storage needed by

this implementation of Bunch-Kaufman is ran for the original storage from

which the matrix is copied, plus ran for the locations to which the matrix

is copied, plus an additional amount C which is the amount of storage

necessary for the fill-in triangles. C is much less than ran, because of the

small number of negative eigenvalues. In addition, two vectors of length n

are needed for storing the D matrix giving a total of 2n(ra+l)+C. So when

C is small, approximately (ra - 2)n storage locations are saved factoring

matrices using the Bunch-Kaufman algorithm instead of LU factorization.

In the second situation (which arises in an efficient implementation of

Lanczos's method for solving Kz = )tMz), the shift a may change during
Method adds. mults, divisions sq. roots comps.

Chol. 44433080 48140336 1824 1824 0

B-K 48277445 48686784 1831 0 446326

LU 137241687 137648943 1823 0 409079

fill

0

2083

2rnr_

Figure 3: Operation Counts for Factorization: n=1824, m---240, 5 negative

eigenvalues

Method adds. mults, divisions sq. roots comps, fill

Chol. 44433080 48140336 1824 1824 0 0

B-K 52023663 52445756 1837 0 485452 14837

LU 137412094 137819350 1823 0 409079 2mn

Figure 4: Operation Counts for Factorization: n=1824, m=240, 19 negative

eigenvalues

execution of the algorithm, so K and M must be saved throughout the computation. In this situation, the storage requirements for LU factorization

is increased to (4ran), but the storage needed by Bunch-Kaufman remains

the same, namely 2n(rn + 1) + C making it even more attractive in this

Ease.

The operation counts for factorization are the same in both cases. The

operation count for Bunch-Kaufman is significantly less than that of LU

factorization because symmetry is exploited and the fill-in is limited. For

simplicity, the operations added by the fill-in during Bunch-Kaufman are

ignored, since the amount that is added is trivial. The high order term

in the operation counts for Bunch-Kaufman is approximately nm 2 arithmetic operations plus approximately nrn comparisons while the high order

term for the operation counts for LU is approximately 4nm 2 arithmetic

operations plus approximately nrn comparisons.

The Bunch-Kaufman method also vectorizes well if the semi-bandwidth

is large enough. The gains from vectorization are much the same as those

9
Method adds. mults, divisions sq. roots comps.

Chol. 3321051 3434180 1980 1980 0

B-K 3322513 3435664 1985 0 142978

LU 10342067 10455196 1979 0 115108

fill

0

22

2rnr_

Figure 5: Operation Counts for Factorization: n=1980, m--59, 5 negative

eigenvalues

Method adds. mults, divisions sq. roots comps.

Chol. 3321051 3434180 1980 1980 0

B-K 3324670 3437856 1985 0 152370

LU 10618321 10731450 1979 0 115108

fill

0

57

2ran

Figure 6: Operation Counts for Factorization: n--1980, m=59, 15 negative

eigenvalues

N M

1824 240

1824 240

1980 59

i980 59

No. of neg. Eigenvalues No. of 2x2 pivots

5 4

19 7

15 3

5 3

Figure 7: The number of 2x2 pivots for each problem

10
for Choleski factorization.

The operations counts for both types of factorization, as well as Choleski

factorization, when using Lanczos's method for solving the generalized

eigenvalue problem are given in figures 3, 4, 5, and 6. The fill-in during factorization is also shown in these figures. The amount of fill-in when

using Bunch-Kaufman can be seen to increase when the number of negative eigenvalues increases. The implementation of LU factorization that is

used for the comparison is sgbfa from the Linpack package [DBMS78]. The

measurements for Choleski factorization are given only as a reference point,

the matrices that were solved were shifted to make them positive definite

for the Choleski factorization runs, otherwise Choleski factorization would

have failed due to the indefiniteness of the system. These matrices arise

from a problem in a structural engineering application [BH87]. In figure 7

the number of 2x2 pivots that occurred in each problem can be examined.

The solution phase that occurs after factorization takes slightly longer

for Bunch-Kaufman than for LU factorization due to the fact that three

matrices, L, D, and L t, arise from Bunch-Kaufman (see equation 1) rather

than just two matrices, L and U, that arise from LU factorization. This

solution phase however takes much less time than factorization, so this is

not significant.

6 Conclusions

The Bunch-Kaufman method has been shown to be a more efficient factorization method than LU factorization in terms of time and storage for

banded real symmetric indefinite matrices with a small number of eigenvalues. An algorithm has been presented that greatly limits the fill needed for

factorization as well as taking advantage of the symmetry of the matrix.

This method has been shown to be nearly as stable as LU factorization by

Bunch [BK77].

References

[BG76] Victor Barwell and Alan George. A comparison of algorithms

for solving symmetric indefinite systems of linear equations.

11
[BH87]

[BK771

[DBMS78]

[DRMN79]

[NOPT83]

ACM Transactions on Mathematical Software, 2(3):242-251,

September 1976.

Charles P. Blankenship and Robert J. Hayduk. Potential supercomputer needs for structural analysis. Presentation at the

Second International Conference on Supercomputing, May 3-8

1987. Santa Clara, CA.

James R. Bunch and Linda Kaufman. Some stable methods

for calculating inertia and solving symmetric linear systems.

Mathematics of Computation, 31(137):163-179, January 1977.

J. Dongarra, J. R. Bunch, C. B. Moler, and G. W. Stewart.

LINPA CK User's Guide, 1978.

I. S. Duff, J. K. Reid, N. Munksgaard, and H. B. Nielsen. Direct solution of sets of linear equations whose matrix is sparse,

symmetric and indefinite. J. Inst. Maths. Applies., 23:235-250,

1979.

Bahram Nour-Omid, Beresford N. Parlett, and Robert L. Taylor. Lanczos versus subspace iteration for solution of eigenvalue

problems. International Journal for Numerical Methods in Engineering, 19:859-871, 1983.

12

Report Documentation Page

1. Report No.

NASA CR- 181860

ICASE Report No. 89-37

2. Government Accession No.

4. Titteand Subtitle

Bunch-Kaufman Factorization for Real Symmetric

Indefinite Banded Matrices

7, Author(s)

Mark T. Jones

Merrell L. Patrick

9. Performing Organization Name and Address

Institute for Computer Applications in Science

and Engineering

Mail Stop 132C, NASA Langley Research Center

Hampton, VA 23665-5225

12. Sponsoring Agency Name and Address

National Aeronautics and Space Administration

Langley Research Center

Hampton, VA 23665-5225

3. Recipient's Catalog No.

5. Report Date

May 1989

6. Performing Organization Code

8. Performing Organization Report No.

89-37

10. Work Unit No.

505-90-21-01

11. Contract or Grant No.

NASI-18605

13. Type of Report and Period Covered

Contractor Report

14. Sponsoring/_,gency Code

15. Supplementaw Notes

Langley Technical Monitor:

Richard W. Barnwell

Final Report

Journal of Approximation

Theory

16. Ab=ract

The Bunch-Kaufman algorithm for factoring symmetric indefinite matrices has been

rejected for banded matrices because it destroys the banded structure of the matrix.

Herein, it is shown that for a subclass of real symmetric matrices which arise in

solving the generalized eigenvalue problem using Lanczos's method, the Bunch-Kaufman

algorithm does not result in major destruction of the bandwidth. Space time

complexities of the algorithm are given and used to show that the Bunch-Kaufman

algorithm is a significant improvement over LU factorization.

17. Key Words(SuggestedbyAuthor(s))

symmetric, indefinite, banded matrices,

Bunch-Kaufman algorithm

19. SecuriW Cla_if. (of this repot)

Unclassified

NASA FORM 1626 OCT 86

18. Distribution Statement

64 - Numerical Analysis

Unclassified - Unlimited

20. SecuriW Cla_if. (of this page)

Unclassified

21, No, of pages

14

22. Price

A03

NASA-Langley, 1989