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