REAL FUNCTION  L7SVX(P, L, X, Y)
C
C  ***  ESTIMATE LARGEST SING. VALUE OF PACKED LOWER TRIANG. MATRIX L
C
C  ***  PARAMETER DECLARATIONS  ***
C
      INTEGER P
      REAL L(1), X(P), Y(P)
C     DIMENSION L(P*(P+1)/2)
C
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C  ***  PURPOSE  ***
C
C     THIS FUNCTION RETURNS A GOOD UNDER-ESTIMATE OF THE LARGEST
C     SINGULAR VALUE OF THE PACKED LOWER TRIANGULAR MATRIX L.
C
C  ***  PARAMETER DESCRIPTION  ***
C
C  P (IN)  = THE ORDER OF L.  L IS A  P X P  LOWER TRIANGULAR MATRIX.
C  L (IN)  = ARRAY HOLDING THE ELEMENTS OF  L  IN ROW ORDER, I.E.
C             L(1,1), L(2,1), L(2,2), L(3,1), L(3,2), L(3,3), ETC.
C  X (OUT) IF  L7SVX RETURNS A POSITIVE VALUE, THEN X = (L**T)*Y IS AN
C             (UNNORMALIZED) APPROXIMATE RIGHT SINGULAR VECTOR
C             CORRESPONDING TO THE LARGEST SINGULAR VALUE.  THIS
C             APPROXIMATION MAY BE CRUDE.
C  Y (OUT) IF  L7SVX RETURNS A POSITIVE VALUE, THEN Y = L*X IS A
C             NORMALIZED APPROXIMATE LEFT SINGULAR VECTOR CORRESPOND-
C             ING TO THE LARGEST SINGULAR VALUE.  THIS APPROXIMATION
C             MAY BE VERY CRUDE.  THE CALLER MAY PASS THE SAME VECTOR
C             FOR X AND Y (NONSTANDARD FORTRAN USAGE), IN WHICH CASE X
C             OVER-WRITES Y.
C
C  ***  ALGORITHM NOTES  ***
C
C     THE ALGORITHM IS BASED ON ANALOGY WITH (1).  IT USES A
C     RANDOM NUMBER GENERATOR PROPOSED IN (4), WHICH PASSES THE
C     SPECTRAL TEST WITH FLYING COLORS -- SEE (2) AND (3).
C
C  ***  SUBROUTINES AND FUNCTIONS CALLED  ***
C
C         V2NRM - FUNCTION, RETURNS THE 2-NORM OF A VECTOR.
C
C  ***  REFERENCES  ***
C
C     (1) CLINE, A., MOLER, C., STEWART, G., AND WILKINSON, J.H.(1977),
C         AN ESTIMATE FOR THE CONDITION NUMBER OF A MATRIX, REPORT
C         TM-310, APPLIED MATH. DIV., ARGONNE NATIONAL LABORATORY.
C
C     (2) HOAGLIN, D.C. (1976), THEORETICAL PROPERTIES OF CONGRUENTIAL
C         RANDOM-NUMBER GENERATORS --  AN EMPIRICAL VIEW,
C         MEMORANDUM NS-340, DEPT. OF STATISTICS, HARVARD UNIV.
C
C     (3) KNUTH, D.E. (1969), THE ART OF COMPUTER PROGRAMMING, VOL. 2
C         (SEMINUMERICAL ALGORITHMS), ADDISON-WESLEY, READING, MASS.
C
C     (4) SMITH, C.S. (1971), MULTIPLICATIVE PSEUDO-RANDOM NUMBER
C         GENERATORS WITH PRIME MODULUS, J. ASSOC. COMPUT. MACH. 18,
C         PP. 586-593.
C
C  ***  HISTORY  ***
C
C     DESIGNED AND CODED BY DAVID M. GAY (WINTER 1977/SUMMER 1978).
C
C  ***  GENERAL  ***
C
C     THIS SUBROUTINE WAS WRITTEN IN CONNECTION WITH RESEARCH
C     SUPPORTED BY THE NATIONAL SCIENCE FOUNDATION UNDER GRANTS
C     MCS-7600324, DCR75-10143, 76-14311DSS, AND MCS76-11989.
C
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C  ***  LOCAL VARIABLES  ***
C
      INTEGER I, IX, J, JI, JJ, JJJ, JM1, J0, PM1, PPLUS1
      REAL B, BLJI, SMINUS, SPLUS, T, YI
C
C  ***  CONSTANTS  ***
C
      REAL HALF, ONE, R9973, ZERO
C
C  ***  EXTERNAL FUNCTIONS AND SUBROUTINES  ***
C
      REAL  D7TPR,  V2NRM
      EXTERNAL  D7TPR,  V2NRM, V2AXY
C
C/6
C     DATA HALF/0.5E+0/, ONE/1.E+0/, R9973/9973.E+0/, ZERO/0.E+0/
C/7
      PARAMETER (HALF=0.5E+0, ONE=1.E+0, R9973=9973.E+0, ZERO=0.E+0)
C/
C
C  ***  BODY  ***
C
      IX = 2
      PPLUS1 = P + 1
      PM1 = P - 1
C
C  ***  FIRST INITIALIZE X TO PARTIAL SUMS  ***
C
      J0 = P*PM1/2
      JJ = J0 + P
      IX = MOD(3432*IX, 9973)
      B = HALF*(ONE + FLOAT(IX)/R9973)
      X(P) = B * L(JJ)
      IF (P .LE. 1) GO TO 40
      DO 10 I = 1, PM1
         JI = J0 + I
         X(I) = B * L(JI)
 10      CONTINUE
C
C  ***  COMPUTE X = (L**T)*B, WHERE THE COMPONENTS OF B HAVE RANDOMLY
C  ***  CHOSEN MAGNITUDES IN (.5,1) WITH SIGNS CHOSEN TO MAKE X LARGE.
C
C     DO J = P-1 TO 1 BY -1...
      DO 30 JJJ = 1, PM1
         J = P - JJJ
C       ***  DETERMINE X(J) IN THIS ITERATION. NOTE FOR I = 1,2,...,J
C       ***  THAT X(I) HOLDS THE CURRENT PARTIAL SUM FOR ROW I.
         IX = MOD(3432*IX, 9973)
         B = HALF*(ONE + FLOAT(IX)/R9973)
         JM1 = J - 1
         J0 = J*JM1/2
         SPLUS = ZERO
         SMINUS = ZERO
         DO 20 I = 1, J
              JI = J0 + I
              BLJI = B * L(JI)
              SPLUS = SPLUS +  ABS(BLJI + X(I))
              SMINUS = SMINUS +  ABS(BLJI - X(I))
 20           CONTINUE
         IF (SMINUS .GT. SPLUS) B = -B
         X(J) = ZERO
C        ***  UPDATE PARTIAL SUMS  ***
         CALL V2AXY(J, X, B, L(J0+1), X)
 30      CONTINUE
C
C  ***  NORMALIZE X  ***
C
 40   T =  V2NRM(P, X)
      IF (T .LE. ZERO) GO TO 80
      T = ONE / T
      DO 50 I = 1, P
 50      X(I) = T*X(I)
C
C  ***  COMPUTE L*X = Y AND RETURN SVMAX = TWONORM(Y)  ***
C
      DO 60 JJJ = 1, P
         J = PPLUS1 - JJJ
         JI = J*(J-1)/2 + 1
         Y(J) =  D7TPR(J, L(JI), X)
 60      CONTINUE
C
C  ***  NORMALIZE Y AND SET X = (L**T)*Y  ***
C
      T = ONE /  V2NRM(P, Y)
      JI = 1
      DO 70 I = 1, P
         YI = T * Y(I)
         X(I) = ZERO
         CALL V2AXY(I, X, YI, L(JI), X)
         JI = JI + I
 70      CONTINUE
       L7SVX =  V2NRM(P, X)
      GO TO 999
C
 80    L7SVX = ZERO
C
 999  RETURN
C  ***  LAST CARD OF  L7SVX FOLLOWS  ***
      END

.