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