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