1      SUBROUTINE REDUC2(NM,N,A,B,DL,IERR)
2C***BEGIN PROLOGUE  REDUC2
3C***DATE WRITTEN   760101   (YYMMDD)
4C***REVISION DATE  830518   (YYMMDD)
5C***CATEGORY NO.  D4C1C
6C***KEYWORDS  EIGENVALUES,EIGENVECTORS,EISPACK
7C***AUTHOR  SMITH, B. T., ET AL.
8C***PURPOSE  Reduces certain generalized symmetric eigenproblems
9C            standard symmetric eigenproblem, using Cholesky
10C            factorization.
11C***DESCRIPTION
12C
13C     This subroutine is a translation of the ALGOL procedure REDUC2,
14C     NUM. MATH. 11, 99-110(1968) by Martin and Wilkinson.
15C     HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 303-314(1971).
16C
17C     This subroutine reduces the generalized SYMMETRIC eigenproblems
18C     ABx=(LAMBDA)x OR BAy=(LAMBDA)y, where B is POSITIVE DEFINITE,
19C     to the standard symmetric eigenproblem using the Cholesky
20C     factorization of B.
21C
22C     On Input
23C
24C        NM must be set to the row dimension of two-dimensional
25C          array parameters as declared in the calling program
26C          dimension statement.
27C
28C        N is the order of the matrices A and B.  If the Cholesky
29C          factor L of B is already available, N should be prefixed
30C          with a minus sign.
31C
32C        A AND B contain the real symmetric input matrices.  Only the
33C          full upper triangles of the matrices need be supplied.  If
34C          N is negative, the strict lower triangle of B contains,
35C          instead, the strict lower triangle of its Cholesky factor L.
36C
37C        DL contains, if N is negative, the diagonal elements of L.
38C
39C     On Output
40C
41C        A contains in its full lower triangle the full lower triangle
42C          of the symmetric matrix derived from the reduction to the
43C          standard form.  The strict upper triangle of A is unaltered.
44C
45C        B contains in its strict lower triangle the strict lower
46C          triangle of its Cholesky factor L.  The full upper
47C          triangle of B is unaltered.
48C
49C        DL contains the diagonal elements of L.
50C
51C        IERR is set to
52C          Zero       for normal return,
53C          7*N+1      if B is not positive definite.
54C
55C     Questions and comments should be directed to B. S. Garbow,
56C     APPLIED MATHEMATICS DIVISION, ARGONNE NATIONAL LABORATORY
57C     ------------------------------------------------------------------
58C***REFERENCES  B. T. SMITH, J. M. BOYLE, J. J. DONGARRA, B. S. GARBOW,
59C                 Y. IKEBE, V. C. KLEMA, C. B. MOLER, *MATRIX EIGEN-
60C                 SYSTEM ROUTINES - EISPACK GUIDE*, SPRINGER-VERLAG,
61C                 1976.
62C***ROUTINES CALLED  (NONE)
63C***END PROLOGUE  REDUC2
64C
65      INTEGER I,J,K,N,I1,J1,NM,NN,IERR
66      REAL A(NM,N),B(NM,N),DL(N)
67      REAL X,Y
68C
69C***FIRST EXECUTABLE STATEMENT  REDUC2
70      IERR = 0
71      NN = IABS(N)
72      IF (N .LT. 0) GO TO 100
73C     .......... FORM L IN THE ARRAYS B AND DL ..........
74      DO 80 I = 1, N
75         I1 = I - 1
76C
77         DO 80 J = I, N
78            X = B(I,J)
79            IF (I .EQ. 1) GO TO 40
80C
81            DO 20 K = 1, I1
82   20       X = X - B(I,K) * B(J,K)
83C
84   40       IF (J .NE. I) GO TO 60
85            IF (X .LE. 0.0E0) GO TO 1000
86            Y = SQRT(X)
87            DL(I) = Y
88            GO TO 80
89   60       B(J,I) = X / Y
90   80 CONTINUE
91C     .......... FORM THE LOWER TRIANGLE OF A*L
92C                IN THE LOWER TRIANGLE OF THE ARRAY A ..........
93  100 DO 200 I = 1, NN
94         I1 = I + 1
95C
96         DO 200 J = 1, I
97            X = A(J,I) * DL(J)
98            IF (J .EQ. I) GO TO 140
99            J1 = J + 1
100C
101            DO 120 K = J1, I
102  120       X = X + A(K,I) * B(K,J)
103C
104  140       IF (I .EQ. NN) GO TO 180
105C
106            DO 160 K = I1, NN
107  160       X = X + A(I,K) * B(K,J)
108C
109  180       A(I,J) = X
110  200 CONTINUE
111C     .......... PRE-MULTIPLY BY TRANSPOSE(L) AND OVERWRITE ..........
112      DO 300 I = 1, NN
113         I1 = I + 1
114         Y = DL(I)
115C
116         DO 300 J = 1, I
117            X = Y * A(I,J)
118            IF (I .EQ. NN) GO TO 280
119C
120            DO 260 K = I1, NN
121  260       X = X + A(K,J) * B(K,I)
122C
123  280       A(I,J) = X
124  300 CONTINUE
125C
126      GO TO 1001
127C     .......... SET ERROR -- B IS NOT POSITIVE DEFINITE ..........
128 1000 IERR = 7 * N + 1
129 1001 RETURN
130      END
131