C*************************************************************************
C                                                                        *
C                                                                        *
C     SPARSE SVD VIA EIGENSYSTEM OF A'A SYMMETRIC MATRIX                 *
C                                                                        *
C                                                                        *
C*************************************************************************
C                                                                        *
C                         (C) COPYRIGHT 1991                             *
C                          MICHAEL W. BERRY                              *
C                         ALL RIGHTS RESERVED                            *
C                                                                        *
C           PERMISSION TO COPY ALL OR PART OF ANY OF THIS                *
C           SOFTWARE IS ONLY GRANTED UPON APPROVAL FROM THE              *
C           AUTHOR: MICHAEL W. BERRY, DEPT. OF COMPUTER                  *
C           SCIENCE, UNIVERSITY OF TENNESSEE,  107 AYRES HALL,           *
C           KNOXVILLE, TN,  37996-1301 (BERRY@CS.UTK.EDU)                *
C                                                                        *
C                                                                        *
C*************************************************************************
C
      PROGRAM SIS2
C
C.... THIS SAMPLE PROGRAM USES RITZIT TO COMPUTE SINGULAR TRIPLETS OF A VIA
C.... THE EQUIVALENT SYMMETRIC EIGENVALUE PROBLEM
C....
C....  B = A'A , AND A IS M (NROW) BY N (NCOL)  (NROW>>NCOL),
C....
C.... SO THAT {U,SQRT(LAMBDA),V} IS A SINGULAR TRIPLET OF A.
C.... (A' = TRANSPOSE OF A)
C....
C.... USER SUPPLIED ROUTINES: OPA,OPB
C....
C.... OPA(     X,Y) TAKES AN N-VECTOR X AND SHOULD RETURN A*X IN Y.
C.... OPB(NCOL,X,Y) TAKES AN N-VECTOR X AND SHOULD RETURN B*X IN Y.
C....
C.... USER SHOULD REPLACE CALLS TO ETIME WITH AN APPROPRIATE CALL TO
C.... AN INSTRINSIC TIMING ROUTINE THAT RETURNS ELAPSED USER TIME
C....
C.... NSIG SPECIFIES MAXIMUM NUMBER OF SINGULAR TRIPLETS DESIRED
C.... VECTORS SPECIFIES WHETHER OR NOT VECTORS ARE DESIRED ON OUTPUT
C....
      PARAMETER(NSIG=100)
C....
      PARAMETER(NMAX=500,NZMAX=2000)
      REAL*8    VALUE(NZMAX)
      INTEGER   POINTR(NMAX), ROWIND(NZMAX)
      COMMON /MATRIXA/ NCOL,NROW,VALUE,POINTR,ROWIND
C
C------------------------------------------------------------
C
      REAL*8 X(NMAX,NSIG),U(NMAX),W(NMAX,NSIG),
     1       D(NSIG),F(NSIG),CX(NMAX),RQ(NSIG),B(NSIG,NSIG),
     2       S(NSIG),EPS,XNORM,DDOT,TMP1,TMP2
C
      REAL    TIME, T0, TMP(2)
C
C------------------------------------------------------------
C
C ----- HARWELL-BOEING COLLECTION DATA STRUCTURE
C
C------------------------------------------------------------
C
      INTEGER  NROW, NCOL
      INTEGER  NNZERO, NRHS, RESULTS, MATRIX, P, EM, EM2, INF
C
      CHARACTER TITLE*72, KEY*8, TYPE*3, PTRFMT*16, INDFMT*16,
     1 VALFMT*20, RHSFMT*20
C
      CHARACTER*30 NAME
      LOGICAL VECTORS
C
C---------------------
C
      EXTERNAL INTROS,OPB,OPA
C
      COMMON /OUTPUT/ INF
      COMMON /COUNT/ MXVCOUNT
      INTEGER MXVCOUNT
C
      INPUT=1
      MATRIX=7
      RESULTS=2
      INF=5
C
C     OPEN FILES FOR INPUT/OUTPUT
C
      OPEN (MATRIX,FILE='LAI7')
      OPEN (RESULTS,FILE='SIO2')
      OPEN (INPUT,FILE='SII1')
      OPEN (INF,FILE='SIO5')
      WRITE(INF,9994)
 
9994  FORMAT('-----------------------------------------------',/
     *       'INTERMEDIATE OUTPUT PARMS:',//
     *       'M := CURRENT DEGREE OF CHEBYSHEV POLYNOMIAL',/
     *       'S := NEXT ITERATION STEP                   ',/
     *       'G := NUMBER OF ACCEPTED EIGENVALUES  OF B  ',/
     *       'H := NUMBER OF ACCEPTED EIGENVECTORS OF B  ',/
     *       'F := VECTOR OF ERRORS FOR EIGENPAIRS OF B  ',/
     *       '-----------------------------------------------',/
     *       '       M       S       G       H        F',/
     *       '-----------------------------------------------'/)

C
      READ (MATRIX,10) TITLE, KEY,
     1   TYPE, NROW, NCOL, NNZERO, NRHS, PTRFMT, INDFMT, VALFMT, RHSFMT
 10   FORMAT (A72, A8 // A3, 11X, 4I14 / 2A16, 2A20)
C
C LEAVE IF MATRIX TOO BIG ----
C
	IF (NCOL .GE. NMAX .OR. NNZERO .GT. NZMAX) THEN
		WRITE(*,*) ' SORRY, YOUR MATRIX IS TOO BIG '
		STOP
		ENDIF
C
C READ DATA...
C
      READ (MATRIX,PTRFMT) (POINTR (I), I = 1, NCOL+1)
      READ (MATRIX,INDFMT) (ROWIND (I), I = 1, NNZERO)
      READ (MATRIX,VALFMT) (VALUE  (I), I = 1, NNZERO)
C
C DEFINE LAST ELEMENT OF POINTR IN CASE IT IS NOT.
C
	POINTR(NCOL+1) = NNZERO + 1
C
      LDA = NMAX
      N=NCOL
C
      READ (INPUT,*) NAME,EM,NUMEXTRA,KM,EPS,VECTORS
C
C     NAME     := DATASET NAME
C     EM       := NUMBER OF DESIRED TRIPLETS
C     NUMEXTRA := ADDITIONAL VECTORS TO CARRY
C     KM       := MAXIMUM NUMBER OF ITERATIONS
C     EPS      := TOLERANCE FOR ACCEPTING SINGULAR VECTORS
C     VECTORS  := OUTPUT SINGULAR VECTORS (BOOLEAN)
C     
C
      P = EM + NUMEXTRA
      EM2=EM
C
      T0 = ETIME(TMP(1))
      CALL RITZIT(LDA,N,P,KM,EPS,OPB,INTROS,EM,X,D,U,W,B,F,CX,RQ,S,
     *            IMEM)
      TIME = ETIME(TMP(1))-T0
C
      WRITE(RESULTS,2000) N,KM,EM2,EM,P,IMEM,VECTORS,EPS,
     *              TITLE,NAME,NROW,NCOL,N
C
9998  FORMAT(/1X,'...... SISVD EXECUTION TIME=',1PE10.2)
9999  FORMAT(1X,'...... '
     *   /1X,'...... ',16X,' MXV  =',I12
     *   /1X,'...... '
     *   /1X,'...... ',4X,'COMPUTED SINGULAR VALUES',2X,
     *    '(','RESIDUAL NORMS',')'
     *   /1X,'...... '
     *   /(1X,'...... ',I3,3X,1PE22.14,2X,'(',1PE11.2,')'))
2000  FORMAT(
     *    1X,'... '
     *   /1X,'... SOLVE THE [A^TA]   EIGENPROBLEM'
     *   /1X,'... NO. OF EQUATIONS          =',I10
     *   /1X,'... MAX. NO. OF ITERATIONS    =',I10
     *   /1X,'... NO. OF DESIRED EIGENPAIRS =',I10
     *   /1X,'... NO. OF APPROX. EIGENPAIRS =',I10
     *   /1X,'... INITIAL SUBSPACE DIM.     =',I10
     *   /1X,'... MEMORY REQUIRED (BYTES)   =',I10
     *   /1X,'... WANT S-VECTORS?   [T/F]   =',L10
     *   /1X,'... TOLERANCE                 =',1PE10.2
     *   //1X,A50
     *   / 1X,A50
     *   /1X,'... NO. OF DOCUMENTS (ROWS)   = ',I8
     *   /1X,'... NO. OF TERMS     (COLS)   = ',I8
     *   /1X,'... ORDER OF MATRIX B         = ',I8
     *   /1X,'... ')
C
      IF(VECTORS) THEN
         IVEC=3
         OPEN(IVEC,FILE='SIO3',FORM='UNFORMATTED')
      ENDIF
C
      T0 = ETIME(TMP(1))
C
      DO 35 J = 1, EM2
         CALL OPB(N,X(1,J),CX)
         TMP1=DDOT(N,X(1,J),1,CX,1)
         CALL DAXPY(N,-TMP1,X(1,J),1,CX,1)
         TMP1=DSQRT(TMP1)
         XNORM=DSQRT(DDOT(N,CX,1,CX,1))
C
C........MULTIPLY BY MATRIX A TO GET (SCALED) LEFT S-VECTOR
C
         CALL OPA(X(1,J),U)
         TMP2=1.0D0/TMP1
         CALL DSCAL(NROW,TMP2,U,1)
         XNORM=XNORM*TMP2
         F(J)=XNORM
         D(J)=TMP1
         IF(VECTORS) WRITE(IVEC) (U(I),I=1,NROW)
         IF(VECTORS) WRITE(IVEC) (X(I,J),I=1,N)
35    CONTINUE
C
      TIME = TIME + ETIME(TMP(1))-T0
C
      WRITE(RESULTS,9998) TIME
      WRITE(RESULTS,9999)MXVCOUNT,(I,D(I),F(I),I = 1,EM2)
C
      STOP
      END
C
      SUBROUTINE INTROS(KS,KG,KH,F,M)
C
C     INTROS IS A SUBROUTINE USED TO OBTAIN
C     INFORMATION OR EXERT CONTROL DURING EXECUTION.  INTROS IS CALLED
C     WITH PARAMETERS (KS,KG,KH,F,M), WHERE KS IS THE NUMBER
C     OF THE NEXT ITERATION STEP, KG IS THE NUMBER OF ALREADY
C     ACCEPTED EIGENVECTORS, KH IS THE NUMBER OF ALREADY ACCEPTED
C     EIGENVALUES, AND F IS THE ARRAY OF ERROR QUANTITIES FOR
C     VECTORS OF X.  AN ELEMENT OF F HAS THE VALUE 4.0
C     UNTIL THE CORRESPONDING EIGENVALUE HAS BEEN ACCEPTED.
C     M IS THE DEGREE OF THE CURRENT CHEBYSHEV POLYNOMIAL.
C
      INTEGER KG,KH,INF,M,KS,L
      REAL*8 F(1)
      COMMON /OUTPUT/ INF
      L = KH + 1
      WRITE(INF,10) M,KS,KG,KH,F(1)
      IF(L.GE.2) WRITE(INF,11) (F(I),I=2,L)
   10 FORMAT(/4I8,E15.6)
   11 FORMAT( 32X,E15.6)
      RETURN
      END
C
      SUBROUTINE RITZIT(NM,N,KP,KM,EPS,OPB,INF,KEM,X,D,U,W,B,F,CX,
     *                  RQ,S,IMEM)
C
      INTEGER I,J,K,L,M,N,IG,II,IK,IP,KG,KH,KM,KP,KS,KZ,L1,M1,NM,KEM,
     *        KZ1,KZ2,JP,IMEM
      REAL*8 X(NM,KP),D(KP),U(NM),W(NM,KP),B(KP,KP),F(KP),CX(KP),
     *        RQ(KP)
      REAL*8 E,S(KP),T,E1,E2,XK,EE2,EPS,XKM,XKS,XM1
      REAL*8 DABS,DEXP,DFLOAT,DLOG,DMAX1,DSQRT,DDOT
C
      INTEGER IABS,IDINT,MAX0,MOD
C
C     CALLS SUBROUTINES OPB,INF,TRED2,TQL2,DGEMM,DGEMV,DAXPY
C     CALLS FUNCTION DDOT,PYTHAG,LSAME
C
C     THIS SUBROUTINE IS A TRANSLATION OF THE ALGOL PROCEDURE RITZIT,
C     NUM. MATH. 16, 205-223(1970) BY RUTISHAUSER.
C     HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 284-302(1971).
C
C     THIS SUBROUTINE DETERMINES THE ABSOLUTELY LARGEST EIGENVALUES
C     AND CORRESPONDING EIGENVECTORS OF A REAL SYMMETRIC MATRIX BY
C     SIMULTANEOUS ITERATION.
C
C     ON INPUT:
C
C        NM MUST BE SET TO THE ROW DIMENSION OF THE TWO-DIMENSIONAL
C          ARRAY PARAMETER X AS DECLARED IN THE CALLING PROGRAM
C          DIMENSION STATEMENT;
C
C        N IS THE ORDER OF THE MATRIX A (MATRIX B FOR SVD PROBLEM);
C
C        KP IS THE NUMBER OF SIMULTANEOUS ITERATION VECTORS;
C
C        KM IS THE MAXIMUM NUMBER OF ITERATION STEPS TO BE PERFORMED.
C          IF STARTING VALUES FOR THE ITERATION VECTORS ARE AVAILABLE,
C          KM SHOULD BE PREFIXED WITH A MINUS SIGN;
C
C        EPS IS THE TOLERANCE FOR ACCEPTING EIGENVECTORS;
C
C        OPB IS THE NAME OF THE SUBROUTINE THAT DEFINES THE MATRIX B.
C          OPB IS CALLED WITH PARAMETERS (N,U,W) AND MUST COMPUTE
C          W=BU WITHOUT ALTERING U FOR SVD PROBLEM;
C
C        INF IS THE NAME OF THE SUBROUTINE THAT MAY BE USED TO OBTAIN
C          INFORMATION OR EXERT CONTROL DURING EXECUTION.  INF IS CALLED
C          WITH PARAMETERS (KS,KG,KH,F,KM,KEM), WHERE KS IS THE NUMBER
C          OF THE NEXT ITERATION STEP, KG IS THE NUMBER OF ALREADY
C          ACCEPTED EIGENVECTORS, KH IS THE NUMBER OF ALREADY ACCEPTED
C          EIGENVALUES, AND F IS THE ARRAY OF ERROR QUANTITIES FOR
C          THE VECTORS OF X.  AN ELEMENT OF F HAS THE VALUE 4.0
C          UNTIL THE CORRESPONDING EIGENVALUE HAS BEEN ACCEPTED;
C
C        KEM IS THE NUMBER OF EIGENVALUES AND CORRESPONDING
C          EIGENVECTORS DESIRED.  KEM MUST BE LESS THAN KP;
C
C        X CONTAINS, IF KM IS NEGATIVE, THE STARTING VALUES FOR
C          THE ITERATION VECTORS.
C
C     ON OUTPUT:
C
C        KM IS RESET TO THE MAGNITUDE OF ITS INPUT VALUE;
C
C        KEM IS RESET TO THE NUMBER OF EIGENVALUES AND EIGENVECTORS
C          ACTUALLY ACCEPTED WITHIN THE LIMIT OF KM STEPS;
C
C        IMEM IS SET TO THE APPROXIMATE NUMBER OF BYTES NEEDED FOR
C          THIS INVOCATION.
C
C        X CONTAINS IN ITS FIRST KEM COLUMNS ORTHONORMAL EIGENVECTORS
C          OF A CORRESPONDING TO THE EIGENVALUES IN D.  THE REMAINING
C          COLUMNS CONTAIN APPROXIMATIONS TO FURTHER EIGENVECTORS;
C
C        D CONTAINS IN ITS FIRST KEM POSITIONS THE ABSOLUTELY LARGEST
C          EIGENVALUES OF A (B).  THE REMAINING POSITIONS CONTAIN
C          APPROXIMATIONS TO SMALLER EIGENVALUES;
C
C        U, W, B, F, CX, S, AND RQ ARE TEMPORARY STORAGE ARRAYS.
C
C     QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO B. S. GARBOW,
C     APPLIED MATHEMATICS DIVISION, ARGONNE NATIONAL LABORATORY.
C     MODERNIZATIONS  BY M. W. BERRY, UNIV. OF TENNESSEE, 1991.
C
C     ------------------------------------------------------------------
C
C.....GET AUXILLARY MEMORY COUNT (8-BYTE WORDS)
C     (ASSUME NM=N FOR COUNTS)
C
      IMEM =        N*KP + KP + KP
      IMEM = IMEM + N + N*KP + KP*KP + KP + KP + KP
      IMEM = 8*IMEM
C
      EE2 = 1.0D0 + 1.0D-1 * EPS
      E = 0.0D0
      KG = 0
      KH = 0
      KZ = 1367
      KZ1 = 0
      KZ2 = 0
      KS = 0
      M = 1
C
      DO 30 L = 1, KP
         F(L) = 4.0D0
         CX(L) = 0.0D0
         RQ(L) = 0.0D0
   30 CONTINUE
C
      IF (KM .LT. 0) GO TO 60
C     :::::::::: GENERATE RANDOM INITIAL ITERATION VECTORS ::::::::::
      DO 50 L = 1, KP
C
         DO 40 J = 1, N
            KZ = MOD(3125*KZ,65536)
            X(J,L) = DFLOAT(KZ-32768)
   40    CONTINUE
C
   50 CONTINUE
C
   60 KM = IABS(KM)
      L1 = 1
      II = 1
      GO TO 905
   65 IG = 1
      IP = KP - 1
C     :::::::::: JACOBI STEP MODIFIED ::::::::::
   70 DO 90 K = IG, KP
         CALL OPB(N,X(1,K),W(1,1))
C
         DO 85 J = 1, N
   85    X(J,K) = W(J,1)
C
   90 CONTINUE
C
      L1 = IG
      II = 2
      GO TO 905
  100 IF (KS .GT. 0) GO TO 130
C     :::::::::: MEASURES AGAINST UNHAPPY CHOICE OF
C                INITIAL VECTORS ::::::::::
      DO 120 K = 1, KP
         IF (B(K,K) .NE. 0.0D0) GO TO 120
C
         DO 110 J = 1, N
            KZ = MOD(3125*KZ,65536)
            X(J,K) = DFLOAT(KZ-32768)
  110    CONTINUE
C
         KS = 1
  120 CONTINUE
C
      IF (KS .NE. 1) GO TO 130
      L1 = 1
      II = 3
      GO TO 905
C
  130 DO 150 K = IG, KP
C
         DO 145 L = K, KP
            T = 0.0D0
C
            DO 140 I = L, KP
  140       T = T + B(I,K) * B(I,L)
C     :::::::::: NEGATE MATRIX TO REVERSE EIGENVALUE ORDERING ::::::::::
            B(L,K) = -T
  145    CONTINUE
C
  150 CONTINUE
C     :::::::::: TRED2, TQL2 ARE MEMBERS OF EISPACK PACKAGE ::::::::::
      J = KP - KG
      CALL TRED2(KP,J,B(IG,IG),D(IG),U,B(IG,IG))
      CALL TQL2(KP,J,D(IG),U,B(IG,IG),II)
C
      DO 190 K = IG, KP
  190 D(K) = DSQRT(DMAX1(-D(K),0.0D0))
C
C     DO 200 K = IG, KP
C        DO 200 J = 1, N
C 200    W(J,K) = 0.0D0
C
C     DO 205 K = IG, KP
C        DO 205 L = IG, KP
C           DO 205 J = 1, N
C 205       W(J,K) = W(J,K) + X(J,L) * B(L,K)
C
      CALL DGEMM('N','N',N,KP-IG+1,KP-IG+1,1.0D0,X(1,IG),NM,
     *           B(IG,IG),KP,0.0D0,W(1,IG),NM)
C
C
      DO 210 K = IG, KP
         DO 210 J = 1, N
  210    X(J,K) = W(J,K)
C
      KS = KS + 1
      XKS = DFLOAT(KS)
      IF (D(KP) .GT. E) E = D(KP)
C     :::::::::: RANDOMIZATION ::::::::::
      IF (KZ1 .GE. 3) GO TO 225
C
      DO 220 J = 1, N
         KZ = MOD(3125*KZ,65536)
         X(J,KP) = DFLOAT(KZ-32768)
  220 CONTINUE
C
      L1 = KP
      II = 4
      GO TO 905
C     :::::::::: COMPUTE CONTROL QUANTITIES CX ::::::::::
  225 DO 270 K = IG, IP
         T = (D(K) - E) * (D(K) + E)
         IF (T .GT. 0.0D0) GO TO 240
         CX(K) = 0.0D0
         GO TO 270
  240    IF (E .NE. 0.0D0) GO TO 260
         CX(K) = 1.0D3 + DLOG(D(K))
         GO TO 270
  260    CX(K) = DLOG((D(K)+DSQRT(T))/E)
  270 CONTINUE
C
C     :::::::::: ACCEPTANCE TEST FOR EIGENVALUES INCLUDING ADJUSTMENT
C                OF KEM AND KH SUCH THAT D(KEM) .GT. E, D(KH) .GT. E,
C                AND D(KEM) DOES NOT OSCILLATE STRONGLY ::::::::::
      DO 300 K = IG, KEM
         IF (D(K) .LE. E .OR.
     X       (KZ1 .GT. 1 .AND. D(K) .LE. 9.99D-1 * RQ(K))) GO TO 320
  300 CONTINUE
C
      GO TO 340
  320 KEM = K - 1
  340 IF (KEM .EQ. 0) GO TO 900
  350 K = KH + 1
      IF (D(K) .EQ. 0.0D0 .OR. D(K) .GT. EE2 * RQ(K)) GO TO 360
      KH = K
      GO TO 350
  360 IF (D(K) .LE. E) KH = K - 1
      K = K - 1
      IF (K .GT. KEM) GO TO 360
C     :::::::::: ACCEPTANCE TEST FOR EIGENVECTORS ::::::::::
      L = KG
      E2 = 0.0D0
C
      DO 560 K = IG, IP
         IF (K .NE. L + 1) GO TO 430
C     :::::::::: CHECK FOR NESTED EIGENVALUES ::::::::::
         L = K
         L1 = K
         IF (K .EQ. IP) GO TO 410
         IK = K + 1
         S(1) = 5.0D-1 / XKS
         T = 1.0D0 / DFLOAT(KS*M)
C
         DO 400 J = IK, IP
            IF (CX(J) * (CX(J) + S(1)) + T .LE. CX(J-1) * CX(J-1))
     X         GO TO 410
            L = J
  400    CONTINUE
C
  410    IF (L .LE. KH) GO TO 430
         L = L1 - 1
         GO TO 570
  430    CALL OPB(N,X(1,K),W(1,1))
         S(1) = 0.0D0
C
         DO 480 J = 1, L
            IF (DABS(D(J)-D(K)) .GE. 1.0D-2 * D(K)) GO TO 480
            T = 0.0D0
C
            DO 460 I = 1, N
  460       T = T + W(I,1) * X(I,J)
C
            DO 470 I = 1, N
  470       W(I,1) = W(I,1) - T * X(I,J)
C
            S(1) = S(1) + T * T
  480    CONTINUE
C
         T = 0.0D0
C
         DO 490 I = 1, N
  490    T = T + W(I,1) * W(I,1)
C
         IF (S(1) .NE. 0.0D0) GO TO 510
         T = 1.0D0
         GO TO 520
  510    T = DSQRT(T/(S(1)+T))
  520    IF (T .GT. E2) E2 = T
         IF (K .NE. L) GO TO 560
C     :::::::::: TEST FOR ACCEPTANCE OF GROUP OF EIGENVECTORS ::::::::::
         IF (L .GE. KEM .AND.
     X       D(KEM) * F(KEM) .LT. EPS * (D(KEM) - E)) KG = KEM
         IF (E2 .GE. F(L)) GO TO 555
C
         DO 550 J = L1, L
  550    F(J) = E2
C
  555    IF (L .LE. KEM .AND. D(L) * F(L) .LT. EPS * (D(L) - E)) KG = L
         IG = KG + 1
  560 CONTINUE
C     :::::::::: ADJUST M ::::::::::
  570 IF (E .GT. 4.0D-2 * D(1)) GO TO 590
      M = 1
      K = 1
      GO TO 600
  590 E2 = 2.0D0 / E
      E1 = 5.1D-1 * E2
      K = 2 * MAX0(IDINT(4.0D0/CX(1)),1)
      IF (M .GT. K) M = K
C     :::::::::: REDUCE KEM IF CONVERGENCE WOULD BE TOO SLOW ::::::::::
  600 XKM = DFLOAT(KM)
      IF (F(KEM) .EQ. 0.0D0 .OR. XKS .GE. 9.0D-1 * XKM) GO TO 660
      XK = DFLOAT(K)
      S(1) = XK * CX(KEM)
      IF (S(1) .GE. 5.0D-2) GO TO 640
      T = 5.0D-1 * S(1) * CX(KEM)
      GO TO 650
  640 T = CX(KEM) + DLOG(5.0D-1*(1.0D0+DEXP(-2.0D0*S(1)))) / XK
  650 S(1) = DLOG(D(KEM)*F(KEM)/(EPS*(D(KEM)-E))) / T
      IF ((XKM - XKS) * XKM .LT. S(1) * XKS) KEM = KEM - 1
  660 CALL INF(KS,KG,KH,F,M)
C
      IF (KG .GE. KEM .OR. KS .GE. KM) GO TO 900
C
      DO 670 K = IG, IP
  670 RQ(K) = D(K)
C
  680 IF (KS + M .LE. KM) GO TO 700
      KZ2 = -1
      IF (M .GT. 1) M = 2 * ((KM - KS + 1) / 2)
  700 M1 = M
C     :::::::::: SHORTCUT LAST INTERMEDIATE BLOCK IF ALL ERROR
C                QUANTITIES F ARE SUFFICIENTLY SMALL ::::::::::
      IF (L .LT. KEM) GO TO 750
      S(1) = D(KEM) * F(KEM) / (EPS * (D(KEM) - E))
      T = S(1) * S(1) - 1.0D0
      IF (T .LE. 0.0D0) GO TO 70
      S(1) = DLOG(S(1)+DSQRT(T)) / (CX(KEM) - CX(KH+1))
      M1 = 2 * IDINT(5.0D-1*S(1)+1.01D0)
      IF (M1 .LE. M) GO TO 740
      M1 = M
      GO TO 750
  740 KZ2 = -1
  750 XM1 = DFLOAT(M1)
C     :::::::::: CHEBYSHEV ITERATION ::::::::::
      IF (M .NE. 1) GO TO 790
C
      DO 780 K = IG, KP
         CALL OPB(N,X(1,K),W(1,1))
C
         DO 775 I = 1, N
  775    X(I,K) = W(I,1)
C
  780 CONTINUE
C
      GO TO 860
C     :::::::::: DEGREE .NE. ONE ::::::::::
C
  790 DO 850 K = IG, KP
         CALL OPB(N,X(1,K),W(1,1))
C
         DO 810 I = 1, N
  810    U(I) = E1 * W(I,1)
C
         CALL OPB(N,U(1),W(1,1))
C
         DO 820 I = 1, N
  820    X(I,K) = E2 * W(I,1) - X(I,K)
C
         IF (M1 .LT. 4) GO TO 850
C
         DO 840 J = 4, M1, 2
            CALL OPB(N,X(1,K),W(1,1))
C
            DO 830 I = 1, N
  830       U(I) = E2 * W(I,1) - U(I)
C
            CALL OPB(N,U,W(1,1))
C
            DO 835 I = 1, N
  835       X(I,K) = E2 * W(I,1) - X(I,K)
C
  840    CONTINUE
C
  850 CONTINUE
C
  860 L1 = IG
      II = 5
      GO TO 905
C     :::::::::: DISCOUNTING THE ERROR QUANTITIES F ::::::::::
  870 IF (KG .EQ. KH) GO TO 895
      IF (M .NE. 1) GO TO 885
C
      DO 880 K = IG, KH
  880 F(K) = F(K) * (D(KH+1) / D(K))
C
      GO TO 895
  885 T = DEXP(-XM1*CX(KH+1))
C
      DO 890 K = IG, KH
         S(1) = DEXP(-XM1*(CX(K)-CX(KH+1)))
         F(K) = S(1) * F(K) * (1.0D0 + T * T) /
     X          (1.0D0 + (S(1) * T) ** 2)
  890 CONTINUE
C
  895 KS = KS + M1
      KZ2 = KZ2 - M1
C     :::::::::: POSSIBLE REPETITION OF INTERMEDIATE STEPS ::::::::::
      IF (KZ2 .GE. 0) GO TO 680
      KZ1 = KZ1 + 1
      KZ2 = 2 * KZ1
      M = 2 * M
      GO TO 70
  900 KEM = KG
C     :::::::::: SOLVE EIGENVALUE PROBLEM OF PROJECTION
C                OF MATRIX A (MATRIX B FOR SVD PROBLEM) ::::::::::
      L1 = 1
      II = 6
C     :::::::::: IN-LINE PROCEDURE FOR EXTENDING ORTHONORMALIZATION
C                TO ALL KP COLUMNS OF X ::::::::::
  905 JP = KP
      IF (II .EQ. 6) JP = KP - 1
      DO 915 I = 1, JP
	 DO 906 K = L1, JP
  906    B(K,I) = 0.0D0
C
         IK = L1
         IF (I .LT. L1) GO TO 909
C
	 IK = I + 1
C
C        DO 907 J = 1, N
C 907    B(I,I) = B(I,I) + X(J,I) * X(J,I)
C        B(I,I) = DSQRT(B(I,I))
C
         B(I,I)=DSQRT(DDOT(N,X(1,I),1,X(1,I),1))
C
         T = 0.0D0
         IF (B(I,I) .NE. 0.0D0) T = 1.0D0 / B(I,I)
C
C        DO 908 J = 1, N
C 908    X(J,I) = T * X(J,I)
C
         CALL DSCAL(N,T,X(1,I),1)
C
  909    CONTINUE
C
         ILEN=JP-IK+1
C   
C        DO 912 K = IK, JP
C           DO 910 J = 1, N
C 910       B(K,I) = B(K,I) + X(J,I) * X(J,K)
C           DO 911 J = 1, N
C 911       X(J,K) = X(J,K) - B(K,I) * X(J,I)
C 912    CONTINUE
C
            CALL DGEMV('T',N,ILEN,1.0D0,X(1,IK),NM,
     *              X(1,I),1,0.0D0,B(IK,I),1)
            DO 912 K=IK,JP
               CALL DAXPY(N,-B(K,I),X(1,I),1,X(1,K),1)
912      CONTINUE
C
  915 CONTINUE
C
      GO TO (65,100,70,225,870,920), II
C
  920 DO 921 K = 1, IP
         DO 921 I = K, IP
  921    B(I,K) = 0.0D0
C
      DO 930 K = 1, IP
         CALL OPB(N,X(1,K),W(1,1))
C
         DO 925 I = 1, K
            DO 922 L = 1, N
  922       B(K,I) = B(K,I) - X(L,I) * W(L,1)
C
  925    CONTINUE
C
  930 CONTINUE
C
      CALL TRED2(KP,IP,B,D,U,B)
      CALL TQL2(KP,IP,D,U,B,II)
C     :::::::::: REORDERING OF EIGENVALUES AND EIGENVECTORS ACCORDING
C                TO THE MAGNITUDES OF THE FORMER ::::::::::
      DO 934 I = 1, IP
         IF (I .EQ. IP) GO TO 933
         K = I
         T = D(I)
         II = I + 1
C
         DO 931 J = II, IP
            IF (DABS(D(J)) .LE. DABS(T)) GO TO 931
            K = J
            T = D(J)
  931    CONTINUE
C
         IF (K .EQ. I) GO TO 933
         D(K) = D(I)
         D(I) = T
C
         DO 932 J = 1, IP
            S(J) = B(J,I)
            B(J,I) = B(J,K)
            B(J,K) = S(J)
  932    CONTINUE
C
  933    D(I) = -D(I)
  934 CONTINUE
C
C     USE LIBRARY KERNEL FOR MATRIX MULT.
C
C     DO 940 I = 1, IP
C        DO 940 J = 1, N
C 940    W(J,I) = 0.0D0
C
C     DO 945 I = 1, IP
C        DO 945 K = 1, IP
C           DO 945 J = 1, N
C 945       W(J,I) = W(J,I) + X(J,K) * B(K,I)
C
      CALL DGEMM('N','N',N,IP,IP,1.0D0,X,NM,B,KP,0.0D0,W,NM)
C
      DO 950 I = 1, IP
         DO 950 J = 1, N
  950    X(J,I) = W(J,I)
C
      D(KP) = E
      RETURN
C     :::::::::: LAST CARD OF RITZIT ::::::::::
      END

      SUBROUTINE DGEMM ( TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB,
     $                   BETA, C, LDC )
*     .. SCALAR ARGUMENTS ..
      CHARACTER*1        TRANSA, TRANSB
      INTEGER            M, N, K, LDA, LDB, LDC
      REAL*8   ALPHA, BETA
*     .. ARRAY ARGUMENTS ..
      REAL*8   A( LDA, * ), B( LDB, * ), C( LDC, * )
*     ..
*
*  PURPOSE
*  =======
*
*  DGEMM  PERFORMS ONE OF THE MATRIX-MATRIX OPERATIONS
*
*     C := ALPHA*OP( A )*OP( B ) + BETA*C,
*
*  WHERE  OP( X ) IS ONE OF
*
*     OP( X ) = X   OR   OP( X ) = X',
*
*  ALPHA AND BETA ARE SCALARS, AND A, B AND C ARE MATRICES, WITH OP( A )
*  AN M BY K MATRIX,  OP( B )  A  K BY N MATRIX AND  C AN M BY N MATRIX.
*
*  PARAMETERS
*  ==========
*
*  TRANSA - CHARACTER*1.
*           ON ENTRY, TRANSA SPECIFIES THE FORM OF OP( A ) TO BE USED IN
*           THE MATRIX MULTIPLICATION AS FOLLOWS:
*
*              TRANSA = 'N' OR 'N',  OP( A ) = A.
*
*              TRANSA = 'T' OR 'T',  OP( A ) = A'.
*
*              TRANSA = 'C' OR 'C',  OP( A ) = A'.
*
*           UNCHANGED ON EXIT.
*
*  TRANSB - CHARACTER*1.
*           ON ENTRY, TRANSB SPECIFIES THE FORM OF OP( B ) TO BE USED IN
*           THE MATRIX MULTIPLICATION AS FOLLOWS:
*
*              TRANSB = 'N' OR 'N',  OP( B ) = B.
*
*              TRANSB = 'T' OR 'T',  OP( B ) = B'.
*
*              TRANSB = 'C' OR 'C',  OP( B ) = B'.
*
*           UNCHANGED ON EXIT.
*
*  M      - INTEGER.
*           ON ENTRY,  M  SPECIFIES  THE NUMBER  OF ROWS  OF THE  MATRIX
*           OP( A )  AND OF THE  MATRIX  C.  M  MUST  BE AT LEAST  ZERO.
*           UNCHANGED ON EXIT.
*
*  N      - INTEGER.
*           ON ENTRY,  N  SPECIFIES THE NUMBER  OF COLUMNS OF THE MATRIX
*           OP( B ) AND THE NUMBER OF COLUMNS OF THE MATRIX C. N MUST BE
*           AT LEAST ZERO.
*           UNCHANGED ON EXIT.
*
*  K      - INTEGER.
*           ON ENTRY,  K  SPECIFIES  THE NUMBER OF COLUMNS OF THE MATRIX
*           OP( A ) AND THE NUMBER OF ROWS OF THE MATRIX OP( B ). K MUST
*           BE AT LEAST  ZERO.
*           UNCHANGED ON EXIT.
*
*  ALPHA  - REAL*8.
*           ON ENTRY, ALPHA SPECIFIES THE SCALAR ALPHA.
*           UNCHANGED ON EXIT.
*
*  A      - REAL*8 ARRAY OF DIMENSION ( LDA, KA ), WHERE KA IS
*           K  WHEN  TRANSA = 'N' OR 'N',  AND IS  M  OTHERWISE.
*           BEFORE ENTRY WITH  TRANSA = 'N' OR 'N',  THE LEADING  M BY K
*           PART OF THE ARRAY  A  MUST CONTAIN THE MATRIX  A,  OTHERWISE
*           THE LEADING  K BY M  PART OF THE ARRAY  A  MUST CONTAIN  THE
*           MATRIX A.
*           UNCHANGED ON EXIT.
*
*  LDA    - INTEGER.
*           ON ENTRY, LDA SPECIFIES THE FIRST DIMENSION OF A AS DECLARED
*           IN THE CALLING (SUB) PROGRAM. WHEN  TRANSA = 'N' OR 'N' THEN
*           LDA MUST BE AT LEAST  MAX( 1, M ), OTHERWISE  LDA MUST BE AT
*           LEAST  MAX( 1, K ).
*           UNCHANGED ON EXIT.
*
*  B      - REAL*8 ARRAY OF DIMENSION ( LDB, KB ), WHERE KB IS
*           N  WHEN  TRANSB = 'N' OR 'N',  AND IS  K  OTHERWISE.
*           BEFORE ENTRY WITH  TRANSB = 'N' OR 'N',  THE LEADING  K BY N
*           PART OF THE ARRAY  B  MUST CONTAIN THE MATRIX  B,  OTHERWISE
*           THE LEADING  N BY K  PART OF THE ARRAY  B  MUST CONTAIN  THE
*           MATRIX B.
*           UNCHANGED ON EXIT.
*
*  LDB    - INTEGER.
*           ON ENTRY, LDB SPECIFIES THE FIRST DIMENSION OF B AS DECLARED
*           IN THE CALLING (SUB) PROGRAM. WHEN  TRANSB = 'N' OR 'N' THEN
*           LDB MUST BE AT LEAST  MAX( 1, K ), OTHERWISE  LDB MUST BE AT
*           LEAST  MAX( 1, N ).
*           UNCHANGED ON EXIT.
*
*  BETA   - REAL*8.
*           ON ENTRY,  BETA  SPECIFIES THE SCALAR  BETA.  WHEN  BETA  IS
*           SUPPLIED AS ZERO THEN C NEED NOT BE SET ON INPUT.
*           UNCHANGED ON EXIT.
*
*  C      - REAL*8 ARRAY OF DIMENSION ( LDC, N ).
*           BEFORE ENTRY, THE LEADING  M BY N  PART OF THE ARRAY  C MUST
*           CONTAIN THE MATRIX  C,  EXCEPT WHEN  BETA  IS ZERO, IN WHICH
*           CASE C NEED NOT BE SET ON ENTRY.
*           ON EXIT, THE ARRAY  C  IS OVERWRITTEN BY THE  M BY N  MATRIX
*           ( ALPHA*OP( A )*OP( B ) + BETA*C ).
*
*  LDC    - INTEGER.
*           ON ENTRY, LDC SPECIFIES THE FIRST DIMENSION OF C AS DECLARED
*           IN  THE  CALLING  (SUB)  PROGRAM.   LDC  MUST  BE  AT  LEAST
*           MAX( 1, M ).
*           UNCHANGED ON EXIT.
*
*
*  LEVEL 3 BLAS ROUTINE.
*
*  -- WRITTEN ON 8-FEBRUARY-1989.
*     JACK DONGARRA, ARGONNE NATIONAL LABORATORY.
*     IAIN DUFF, AERE HARWELL.
*     JEREMY DU CROZ, NUMERICAL ALGORITHMS GROUP LTD.
*     SVEN HAMMARLING, NUMERICAL ALGORITHMS GROUP LTD.
*
*
*     .. EXTERNAL FUNCTIONS ..
      LOGICAL            LSAME
      EXTERNAL           LSAME
*     .. EXTERNAL SUBROUTINES ..
      EXTERNAL           XERBLA
*     .. INTRINSIC FUNCTIONS ..
      INTRINSIC          MAX
*     .. LOCAL SCALARS ..
      LOGICAL            NOTA, NOTB
      INTEGER            I, INFO, J, L, NCOLA, NROWA, NROWB
      REAL*8   TEMP
*     .. PARAMETERS ..
      REAL*8   ONE         , ZERO
      PARAMETER        ( ONE = 1.0D+0, ZERO = 0.0D+0 )
*     ..
*     .. EXECUTABLE STATEMENTS ..
*
*     SET  NOTA  AND  NOTB  AS  TRUE IF  A  AND  B  RESPECTIVELY ARE NOT
*     TRANSPOSED AND SET  NROWA, NCOLA AND  NROWB  AS THE NUMBER OF ROWS
*     AND  COLUMNS OF  A  AND THE  NUMBER OF  ROWS  OF  B  RESPECTIVELY.
*
      NOTA  = LSAME( TRANSA, 'N' )
      NOTB  = LSAME( TRANSB, 'N' )
      IF( NOTA )THEN
         NROWA = M
         NCOLA = K
      ELSE
         NROWA = K
         NCOLA = M
      END IF
      IF( NOTB )THEN
         NROWB = K
      ELSE
         NROWB = N
      END IF
*
*     TEST THE INPUT PARAMETERS.
*
      INFO = 0
      IF(      ( .NOT.NOTA                 ).AND.
     $         ( .NOT.LSAME( TRANSA, 'C' ) ).AND.
     $         ( .NOT.LSAME( TRANSA, 'T' ) )      )THEN
         INFO = 1
      ELSE IF( ( .NOT.NOTB                 ).AND.
     $         ( .NOT.LSAME( TRANSB, 'C' ) ).AND.
     $         ( .NOT.LSAME( TRANSB, 'T' ) )      )THEN
         INFO = 2
      ELSE IF( M  .LT.0               )THEN
         INFO = 3
      ELSE IF( N  .LT.0               )THEN
         INFO = 4
      ELSE IF( K  .LT.0               )THEN
         INFO = 5
      ELSE IF( LDA.LT.MAX( 1, NROWA ) )THEN
         INFO = 8
      ELSE IF( LDB.LT.MAX( 1, NROWB ) )THEN
         INFO = 10
      ELSE IF( LDC.LT.MAX( 1, M     ) )THEN
         INFO = 13
      END IF
      IF( INFO.NE.0 )THEN
         CALL XERBLA( 'DGEMM ', INFO )
         RETURN
      END IF
*
*     QUICK RETURN IF POSSIBLE.
*
      IF( ( M.EQ.0 ).OR.( N.EQ.0 ).OR.
     $    ( ( ( ALPHA.EQ.ZERO ).OR.( K.EQ.0 ) ).AND.( BETA.EQ.ONE ) ) )
     $   RETURN
*
*     AND IF  ALPHA.EQ.ZERO.
*
      IF( ALPHA.EQ.ZERO )THEN
         IF( BETA.EQ.ZERO )THEN
            DO 20, J = 1, N
               DO 10, I = 1, M
                  C( I, J ) = ZERO
   10          CONTINUE
   20       CONTINUE
         ELSE
            DO 40, J = 1, N
               DO 30, I = 1, M
                  C( I, J ) = BETA*C( I, J )
   30          CONTINUE
   40       CONTINUE
         END IF
         RETURN
      END IF
*
*     START THE OPERATIONS.
*
      IF( NOTB )THEN
         IF( NOTA )THEN
*
*           FORM  C := ALPHA*A*B + BETA*C.
*
            DO 90, J = 1, N
               IF( BETA.EQ.ZERO )THEN
                  DO 50, I = 1, M
                     C( I, J ) = ZERO
   50             CONTINUE
               ELSE IF( BETA.NE.ONE )THEN
                  DO 60, I = 1, M
                     C( I, J ) = BETA*C( I, J )
   60             CONTINUE
               END IF
               DO 80, L = 1, K
                  IF( B( L, J ).NE.ZERO )THEN
                     TEMP = ALPHA*B( L, J )
                     DO 70, I = 1, M
                        C( I, J ) = C( I, J ) + TEMP*A( I, L )
   70                CONTINUE
                  END IF
   80          CONTINUE
   90       CONTINUE
         ELSE
*
*           FORM  C := ALPHA*A'*B + BETA*C
*
            DO 120, J = 1, N
               DO 110, I = 1, M
                  TEMP = ZERO
                  DO 100, L = 1, K
                     TEMP = TEMP + A( L, I )*B( L, J )
  100             CONTINUE
                  IF( BETA.EQ.ZERO )THEN
                     C( I, J ) = ALPHA*TEMP
                  ELSE
                     C( I, J ) = ALPHA*TEMP + BETA*C( I, J )
                  END IF
  110          CONTINUE
  120       CONTINUE
         END IF
      ELSE
         IF( NOTA )THEN
*
*           FORM  C := ALPHA*A*B' + BETA*C
*
            DO 170, J = 1, N
               IF( BETA.EQ.ZERO )THEN
                  DO 130, I = 1, M
                     C( I, J ) = ZERO
  130             CONTINUE
               ELSE IF( BETA.NE.ONE )THEN
                  DO 140, I = 1, M
                     C( I, J ) = BETA*C( I, J )
  140             CONTINUE
               END IF
               DO 160, L = 1, K
                  IF( B( J, L ).NE.ZERO )THEN
                     TEMP = ALPHA*B( J, L )
                     DO 150, I = 1, M
                        C( I, J ) = C( I, J ) + TEMP*A( I, L )
  150                CONTINUE
                  END IF
  160          CONTINUE
  170       CONTINUE
         ELSE
*
*           FORM  C := ALPHA*A'*B' + BETA*C
*
            DO 200, J = 1, N
               DO 190, I = 1, M
                  TEMP = ZERO
                  DO 180, L = 1, K
                     TEMP = TEMP + A( L, I )*B( J, L )
  180             CONTINUE
                  IF( BETA.EQ.ZERO )THEN
                     C( I, J ) = ALPHA*TEMP
                  ELSE
                     C( I, J ) = ALPHA*TEMP + BETA*C( I, J )
                  END IF
  190          CONTINUE
  200       CONTINUE
         END IF
      END IF
*
      RETURN
*
*     END OF DGEMM .
*
      END
      SUBROUTINE DGEMV ( TRANS, M, N, ALPHA, A, LDA, X, INCX,
     $                   BETA, Y, INCY )
*     .. SCALAR ARGUMENTS ..
      REAL*8   ALPHA, BETA
      INTEGER            INCX, INCY, LDA, M, N
      CHARACTER*1        TRANS
*     .. ARRAY ARGUMENTS ..
      REAL*8   A( LDA, * ), X( * ), Y( * )
*     ..
*
*  PURPOSE
*  =======
*
*  DGEMV  PERFORMS ONE OF THE MATRIX-VECTOR OPERATIONS
*
*     Y := ALPHA*A*X + BETA*Y,   OR   Y := ALPHA*A'*X + BETA*Y,
*
*  WHERE ALPHA AND BETA ARE SCALARS, X AND Y ARE VECTORS AND A IS AN
*  M BY N MATRIX.
*
*  PARAMETERS
*  ==========
*
*  TRANS  - CHARACTER*1.
*           ON ENTRY, TRANS SPECIFIES THE OPERATION TO BE PERFORMED AS
*           FOLLOWS:
*
*              TRANS = 'N' OR 'N'   Y := ALPHA*A*X + BETA*Y.
*
*              TRANS = 'T' OR 'T'   Y := ALPHA*A'*X + BETA*Y.
*
*              TRANS = 'C' OR 'C'   Y := ALPHA*A'*X + BETA*Y.
*
*           UNCHANGED ON EXIT.
*
*  M      - INTEGER.
*           ON ENTRY, M SPECIFIES THE NUMBER OF ROWS OF THE MATRIX A.
*           M MUST BE AT LEAST ZERO.
*           UNCHANGED ON EXIT.
*
*  N      - INTEGER.
*           ON ENTRY, N SPECIFIES THE NUMBER OF COLUMNS OF THE MATRIX A.
*           N MUST BE AT LEAST ZERO.
*           UNCHANGED ON EXIT.
*
*  ALPHA  - REAL*8.
*           ON ENTRY, ALPHA SPECIFIES THE SCALAR ALPHA.
*           UNCHANGED ON EXIT.
*
*  A      - REAL*8 ARRAY OF DIMENSION ( LDA, N ).
*           BEFORE ENTRY, THE LEADING M BY N PART OF THE ARRAY A MUST
*           CONTAIN THE MATRIX OF COEFFICIENTS.
*           UNCHANGED ON EXIT.
*
*  LDA    - INTEGER.
*           ON ENTRY, LDA SPECIFIES THE FIRST DIMENSION OF A AS DECLARED
*           IN THE CALLING (SUB) PROGRAM. LDA MUST BE AT LEAST
*           MAX( 1, M ).
*           UNCHANGED ON EXIT.
*
*  X      - REAL*8 ARRAY OF DIMENSION AT LEAST
*           ( 1 + ( N - 1 )*ABS( INCX ) ) WHEN TRANS = 'N' OR 'N'
*           AND AT LEAST
*           ( 1 + ( M - 1 )*ABS( INCX ) ) OTHERWISE.
*           BEFORE ENTRY, THE INCREMENTED ARRAY X MUST CONTAIN THE
*           VECTOR X.
*           UNCHANGED ON EXIT.
*
*  INCX   - INTEGER.
*           ON ENTRY, INCX SPECIFIES THE INCREMENT FOR THE ELEMENTS OF
*           X. INCX MUST NOT BE ZERO.
*           UNCHANGED ON EXIT.
*
*  BETA   - REAL*8.
*           ON ENTRY, BETA SPECIFIES THE SCALAR BETA. WHEN BETA IS
*           SUPPLIED AS ZERO THEN Y NEED NOT BE SET ON INPUT.
*           UNCHANGED ON EXIT.
*
*  Y      - REAL*8 ARRAY OF DIMENSION AT LEAST
*           ( 1 + ( M - 1 )*ABS( INCY ) ) WHEN TRANS = 'N' OR 'N'
*           AND AT LEAST
*           ( 1 + ( N - 1 )*ABS( INCY ) ) OTHERWISE.
*           BEFORE ENTRY WITH BETA NON-ZERO, THE INCREMENTED ARRAY Y
*           MUST CONTAIN THE VECTOR Y. ON EXIT, Y IS OVERWRITTEN BY THE
*           UPDATED VECTOR Y.
*
*  INCY   - INTEGER.
*           ON ENTRY, INCY SPECIFIES THE INCREMENT FOR THE ELEMENTS OF
*           Y. INCY MUST NOT BE ZERO.
*           UNCHANGED ON EXIT.
*
*
*  LEVEL 2 BLAS ROUTINE.
*
*  -- WRITTEN ON 22-OCTOBER-1986.
*     JACK DONGARRA, ARGONNE NATIONAL LAB.
*     JEREMY DU CROZ, NAG CENTRAL OFFICE.
*     SVEN HAMMARLING, NAG CENTRAL OFFICE.
*     RICHARD HANSON, SANDIA NATIONAL LABS.
*
*
*     .. PARAMETERS ..
      REAL*8   ONE         , ZERO
      PARAMETER        ( ONE = 1.0D+0, ZERO = 0.0D+0 )
*     .. LOCAL SCALARS ..
      REAL*8   TEMP
      INTEGER            I, INFO, IX, IY, J, JX, JY, KX, KY, LENX, LENY
*     .. EXTERNAL FUNCTIONS ..
      LOGICAL            LSAME
      EXTERNAL           LSAME
*     .. EXTERNAL SUBROUTINES ..
      EXTERNAL           XERBLA
*     .. INTRINSIC FUNCTIONS ..
      INTRINSIC          MAX
*     ..
*     .. EXECUTABLE STATEMENTS ..
*
*     TEST THE INPUT PARAMETERS.
*
      INFO = 0
      IF     ( .NOT.LSAME( TRANS, 'N' ).AND.
     $         .NOT.LSAME( TRANS, 'T' ).AND.
     $         .NOT.LSAME( TRANS, 'C' )      )THEN
         INFO = 1
      ELSE IF( M.LT.0 )THEN
         INFO = 2
      ELSE IF( N.LT.0 )THEN
         INFO = 3
      ELSE IF( LDA.LT.MAX( 1, M ) )THEN
         INFO = 6
      ELSE IF( INCX.EQ.0 )THEN
         INFO = 8
      ELSE IF( INCY.EQ.0 )THEN
         INFO = 11
      END IF
      IF( INFO.NE.0 )THEN
         CALL XERBLA( 'DGEMV ', INFO )
         RETURN
      END IF
*
*     QUICK RETURN IF POSSIBLE.
*
      IF( ( M.EQ.0 ).OR.( N.EQ.0 ).OR.
     $    ( ( ALPHA.EQ.ZERO ).AND.( BETA.EQ.ONE ) ) )
     $   RETURN
*
*     SET  LENX  AND  LENY, THE LENGTHS OF THE VECTORS X AND Y, AND SET
*     UP THE START POINTS IN  X  AND  Y.
*
      IF( LSAME( TRANS, 'N' ) )THEN
         LENX = N
         LENY = M
      ELSE
         LENX = M
         LENY = N
      END IF
      IF( INCX.GT.0 )THEN
         KX = 1
      ELSE
         KX = 1 - ( LENX - 1 )*INCX
      END IF
      IF( INCY.GT.0 )THEN
         KY = 1
      ELSE
         KY = 1 - ( LENY - 1 )*INCY
      END IF
*
*     START THE OPERATIONS. IN THIS VERSION THE ELEMENTS OF A ARE
*     ACCESSED SEQUENTIALLY WITH ONE PASS THROUGH A.
*
*     FIRST FORM  Y := BETA*Y.
*
      IF( BETA.NE.ONE )THEN
         IF( INCY.EQ.1 )THEN
            IF( BETA.EQ.ZERO )THEN
               DO 10, I = 1, LENY
                  Y( I ) = ZERO
   10          CONTINUE
            ELSE
               DO 20, I = 1, LENY
                  Y( I ) = BETA*Y( I )
   20          CONTINUE
            END IF
         ELSE
            IY = KY
            IF( BETA.EQ.ZERO )THEN
               DO 30, I = 1, LENY
                  Y( IY ) = ZERO
                  IY      = IY   + INCY
   30          CONTINUE
            ELSE
               DO 40, I = 1, LENY
                  Y( IY ) = BETA*Y( IY )
                  IY      = IY           + INCY
   40          CONTINUE
            END IF
         END IF
      END IF
      IF( ALPHA.EQ.ZERO )
     $   RETURN
      IF( LSAME( TRANS, 'N' ) )THEN
*
*        FORM  Y := ALPHA*A*X + Y.
*
         JX = KX
         IF( INCY.EQ.1 )THEN
            DO 60, J = 1, N
               IF( X( JX ).NE.ZERO )THEN
                  TEMP = ALPHA*X( JX )
                  DO 50, I = 1, M
                     Y( I ) = Y( I ) + TEMP*A( I, J )
   50             CONTINUE
               END IF
               JX = JX + INCX
   60       CONTINUE
         ELSE
            DO 80, J = 1, N
               IF( X( JX ).NE.ZERO )THEN
                  TEMP = ALPHA*X( JX )
                  IY   = KY
                  DO 70, I = 1, M
                     Y( IY ) = Y( IY ) + TEMP*A( I, J )
                     IY      = IY      + INCY
   70             CONTINUE
               END IF
               JX = JX + INCX
   80       CONTINUE
         END IF
      ELSE
*
*        FORM  Y := ALPHA*A'*X + Y.
*
         JY = KY
         IF( INCX.EQ.1 )THEN
            DO 100, J = 1, N
               TEMP = ZERO
               DO 90, I = 1, M
                  TEMP = TEMP + A( I, J )*X( I )
   90          CONTINUE
               Y( JY ) = Y( JY ) + ALPHA*TEMP
               JY      = JY      + INCY
  100       CONTINUE
         ELSE
            DO 120, J = 1, N
               TEMP = ZERO
               IX   = KX
               DO 110, I = 1, M
                  TEMP = TEMP + A( I, J )*X( IX )
                  IX   = IX   + INCX
  110          CONTINUE
               Y( JY ) = Y( JY ) + ALPHA*TEMP
               JY      = JY      + INCY
  120       CONTINUE
         END IF
      END IF
*
      RETURN
*
*     END OF DGEMV .
*
      END
      LOGICAL          FUNCTION LSAME( CA, CB )
*
*  -- LAPACK AUXILIARY ROUTINE --
*     ARGONNE NATIONAL LABORATORY
*     OCTOBER 11, 1988
*
*     .. SCALAR ARGUMENTS ..
      CHARACTER          CA, CB
*     ..
*
*  PURPOSE
*  =======
*
*     LSAME RETURNS .TRUE. IF CA IS THE SAME LETTER AS CB REGARDLESS
*     OF CASE.
*
*  N.B. THIS VERSION OF THE ROUTINE IS ONLY CORRECT FOR ASCII CODE.
*       INSTALLERS MUST MODIFY THE ROUTINE FOR OTHER CHARACTER-CODES.
*
*       FOR EBCDIC SYSTEMS THE CONSTANT IOFF MUST BE CHANGED TO -64.
*       FOR CDC SYSTEMS USING 6-12 BIT REPRESENTATIONS, THE SYSTEM-
*       SPECIFIC CODE IN COMMENTS MUST BE ACTIVATED.
*
*  ARGUMENTS
*  =========
*
*  CA     - CHARACTER*1
*  CB     - CHARACTER*1
*           ON ENTRY, CA AND CB SPECIFY CHARACTERS TO BE COMPARED.
*           UNCHANGED ON EXIT.
*
*  AUXILIARY ROUTINE FOR LEVEL 2 BLAS.
*
*  -- WRITTEN ON 20-JULY-1986
*     RICHARD HANSON, SANDIA NATIONAL LABS.
*     JEREMY DU CROZ, NAG CENTRAL OFFICE.
*
*
*     .. PARAMETERS ..
      INTEGER            IOFF
      PARAMETER        ( IOFF = 32 )
*     ..
*     .. INTRINSIC FUNCTIONS ..
      INTRINSIC          ICHAR
*     ..
*     .. EXECUTABLE STATEMENTS ..
*
*     TEST IF THE CHARACTERS ARE EQUAL
*
      LSAME = CA.EQ.CB
*
*     NOW TEST FOR EQUIVALENCE
*
      IF( .NOT.LSAME ) THEN
         LSAME = ICHAR( CA ) - IOFF.EQ.ICHAR( CB )
      END IF
      IF( .NOT.LSAME ) THEN
         LSAME = ICHAR( CA ).EQ.ICHAR( CB ) - IOFF
      END IF
*
      RETURN
*
*  THE FOLLOWING COMMENTS CONTAIN CODE FOR CDC SYSTEMS USING 6-12 BIT
*  REPRESENTATIONS.
*
*     .. PARAMETERS ..
*     INTEGER            ICIRFX
*     PARAMETER        ( ICIRFX=62 )
*     .. SCALAR ARGUMENTS ..
*     CHARACTER*1        CB
*     .. ARRAY ARGUMENTS ..
*     CHARACTER*1        CA(*)
*     .. LOCAL SCALARS ..
*     INTEGER            IVAL
*     .. INTRINSIC FUNCTIONS ..
*     INTRINSIC          ICHAR, CHAR
*     .. EXECUTABLE STATEMENTS ..
*
*     SEE IF THE FIRST CHARACTER IN STRING CA EQUALS STRING CB.
*
*     LSAME = CA(1) .EQ. CB .AND. CA(1) .NE. CHAR(ICIRFX)
*
*     IF (LSAME) RETURN
*
*     THE CHARACTERS ARE NOT IDENTICAL. NOW CHECK THEM FOR EQUIVALENCE.
*     LOOK FOR THE 'ESCAPE' CHARACTER, CIRCUMFLEX, FOLLOWED BY THE
*     LETTER.
*
*     IVAL = ICHAR(CA(2))
*     IF (IVAL.GE.ICHAR('A') .AND. IVAL.LE.ICHAR('Z')) THEN
*        LSAME = CA(1) .EQ. CHAR(ICIRFX) .AND. CA(2) .EQ. CB
*     END IF
*
*     RETURN
*
*     END OF LSAME.
*
      END
      SUBROUTINE XERBLA( SRNAME, INFO )
*
*  -- LAPACK AUXILIARY ROUTINE --
*     ARGONNE NATIONAL LABORATORY
*     NOVEMBER 16, 1988
*
*     .. SCALAR ARGUMENTS ..
      CHARACTER*6        SRNAME
      INTEGER            INFO
*     ..
*
*  PURPOSE
*  =======
*
*     XERBLA  IS AN ERROR HANDLER FOR THE LAPACK ROUTINES.
*     IT IS CALLED BY AN LAPACK ROUTINE IF AN INPUT PARAMETER HAS AN
*     INVALID VALUE.  A MESSAGE IS PRINTED AND EXECUTION STOPS.
*
*     INSTALLERS MAY CONSIDER MODIFYING THE STOP STATEMENT IN ORDER TO
*     CALL SYSTEM-SPECIFIC EXCEPTION-HANDLING FACILITIES.
*
*  PARAMETERS
*  ==========
*
*  SRNAME - CHARACTER*6.
*           ON ENTRY, SRNAME SPECIFIES THE NAME OF THE ROUTINE WHICH
*           CALLED XERBLA.
*
*  INFO   - INTEGER.
*           ON ENTRY, INFO SPECIFIES THE POSITION OF THE INVALID
*           PARAMETER IN THE PARAMETER-LIST OF THE CALLING ROUTINE.
*
*
      WRITE( *, FMT = 9999 )SRNAME, INFO
*
      STOP
*
 9999 FORMAT( ' ** ON ENTRY TO ', A6, ' PARAMETER NUMBER ', I2, ' HAD ',
     $      'AN ILLEGAL VALUE' )
*
*     END OF XERBLA
*
      END
C
      REAL*8 FUNCTION PYTHAG(A,B)
      REAL*8 A,B
C
C     FINDS DSQRT(A**2+B**2) WITHOUT OVERFLOW OR DESTRUCTIVE UNDERFLOW
C
      REAL*8 P,R,S,T,U
C
      P = DMAX1(DABS(A),DABS(B))
      IF (P .EQ. 0.0D0) GO TO 20
      R = (DMIN1(DABS(A),DABS(B))/P)**2
   10 CONTINUE
         T = 4.0D0 + R
         IF (T .EQ. 4.0D0) GO TO 20
         S = R/T
         U = 1.0D0 + 2.0D0*S
         P = U*P
         R = (S/U)**2 * R
      GO TO 10
   20 PYTHAG = P
      RETURN
      END
      SUBROUTINE TRED2(NM,N,A,D,E,Z)
C
      INTEGER I,J,K,L,N,II,NM,JP1
      DOUBLE PRECISION A(NM,N),D(N),E(N),Z(NM,N)
      DOUBLE PRECISION F,G,H,HH,SCALE
C
C     THIS SUBROUTINE IS A TRANSLATION OF THE ALGOL PROCEDURE TRED2,
C     NUM. MATH. 11, 181-195(1968) BY MARTIN, REINSCH, AND WILKINSON.
C     HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 212-226(1971).
C
C     THIS SUBROUTINE REDUCES A REAL SYMMETRIC MATRIX TO A
C     SYMMETRIC TRIDIAGONAL MATRIX USING AND ACCUMULATING
C     ORTHOGONAL SIMILARITY TRANSFORMATIONS.
C
C     ON INPUT
C
C        NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL
C          ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM
C          DIMENSION STATEMENT.
C
C        N IS THE ORDER OF THE MATRIX.
C
C        A CONTAINS THE REAL SYMMETRIC INPUT MATRIX.  ONLY THE
C          LOWER TRIANGLE OF THE MATRIX NEED BE SUPPLIED.
C
C     ON OUTPUT
C
C        D CONTAINS THE DIAGONAL ELEMENTS OF THE TRIDIAGONAL MATRIX.
C
C        E CONTAINS THE SUBDIAGONAL ELEMENTS OF THE TRIDIAGONAL
C          MATRIX IN ITS LAST N-1 POSITIONS.  E(1) IS SET TO ZERO.
C
C        Z CONTAINS THE ORTHOGONAL TRANSFORMATION MATRIX
C          PRODUCED IN THE REDUCTION.
C
C        A AND Z MAY COINCIDE.  IF DISTINCT, A IS UNALTERED.
C
C     QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW,
C     MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY
C
C     THIS VERSION DATED AUGUST 1983.
C
C     ------------------------------------------------------------------
C
      DO 100 I = 1, N
C
         DO 80 J = I, N
   80    Z(J,I) = A(J,I)
C
         D(I) = A(N,I)
  100 CONTINUE
C
      IF (N .EQ. 1) GO TO 510
C     .......... FOR I=N STEP -1 UNTIL 2 DO -- ..........
      DO 300 II = 2, N
         I = N + 2 - II
         L = I - 1
         H = 0.0D0
         SCALE = 0.0D0
         IF (L .LT. 2) GO TO 130
C     .......... SCALE ROW (ALGOL TOL THEN NOT NEEDED) ..........
         DO 120 K = 1, L
  120    SCALE = SCALE + DABS(D(K))
C
         IF (SCALE .NE. 0.0D0) GO TO 140
  130    E(I) = D(L)
C
         DO 135 J = 1, L
            D(J) = Z(L,J)
            Z(I,J) = 0.0D0
            Z(J,I) = 0.0D0
  135    CONTINUE
C
         GO TO 290
C
  140    DO 150 K = 1, L
            D(K) = D(K) / SCALE
            H = H + D(K) * D(K)
  150    CONTINUE
C
         F = D(L)
         G = -DSIGN(DSQRT(H),F)
         E(I) = SCALE * G
         H = H - F * G
         D(L) = F - G
C     .......... FORM A*U ..........
         DO 170 J = 1, L
  170    E(J) = 0.0D0
C
         DO 240 J = 1, L
            F = D(J)
            Z(J,I) = F
            G = E(J) + Z(J,J) * F
            JP1 = J + 1
            IF (L .LT. JP1) GO TO 220
C
            DO 200 K = JP1, L
               G = G + Z(K,J) * D(K)
               E(K) = E(K) + Z(K,J) * F
  200       CONTINUE
C
  220       E(J) = G
  240    CONTINUE
C     .......... FORM P ..........
         F = 0.0D0
C
         DO 245 J = 1, L
            E(J) = E(J) / H
            F = F + E(J) * D(J)
  245    CONTINUE
C
         HH = F / (H + H)
C     .......... FORM Q ..........
         DO 250 J = 1, L
  250    E(J) = E(J) - HH * D(J)
C     .......... FORM REDUCED A ..........
         DO 280 J = 1, L
            F = D(J)
            G = E(J)
C
            DO 260 K = J, L
  260       Z(K,J) = Z(K,J) - F * E(K) - G * D(K)
C
            D(J) = Z(L,J)
            Z(I,J) = 0.0D0
  280    CONTINUE
C
  290    D(I) = H
  300 CONTINUE
C     .......... ACCUMULATION OF TRANSFORMATION MATRICES ..........
      DO 500 I = 2, N
         L = I - 1
         Z(N,L) = Z(L,L)
         Z(L,L) = 1.0D0
         H = D(I)
         IF (H .EQ. 0.0D0) GO TO 380
C
         DO 330 K = 1, L
  330    D(K) = Z(K,I) / H
C
         DO 360 J = 1, L
            G = 0.0D0
C
            DO 340 K = 1, L
  340       G = G + Z(K,I) * Z(K,J)
C
            DO 360 K = 1, L
               Z(K,J) = Z(K,J) - G * D(K)
  360    CONTINUE
C
  380    DO 400 K = 1, L
  400    Z(K,I) = 0.0D0
C
  500 CONTINUE
C
  510 DO 520 I = 1, N
         D(I) = Z(N,I)
         Z(N,I) = 0.0D0
  520 CONTINUE
C
      Z(N,N) = 1.0D0
      E(1) = 0.0D0
      RETURN
      END

      SUBROUTINE TQL2(NM,N,D,E,Z,IERR)
C
      INTEGER I,J,K,L,M,N,II,L1,L2,NM,MML,IERR
      DOUBLE PRECISION D(N),E(N),Z(NM,N)
      DOUBLE PRECISION C,C2,C3,DL1,EL1,F,G,H,P,R,S,S2,TST1,TST2,PYTHAG
C
C     THIS SUBROUTINE IS A TRANSLATION OF THE ALGOL PROCEDURE TQL2,
C     NUM. MATH. 11, 293-306(1968) BY BOWDLER, MARTIN, REINSCH, AND
C     WILKINSON.
C     HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 227-240(1971).
C
C     THIS SUBROUTINE FINDS THE EIGENVALUES AND EIGENVECTORS
C     OF A SYMMETRIC TRIDIAGONAL MATRIX BY THE QL METHOD.
C     THE EIGENVECTORS OF A FULL SYMMETRIC MATRIX CAN ALSO
C     BE FOUND IF  TRED2  HAS BEEN USED TO REDUCE THIS
C     FULL MATRIX TO TRIDIAGONAL FORM.
C
C     ON INPUT
C
C        NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL
C          ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM
C          DIMENSION STATEMENT.
C
C        N IS THE ORDER OF THE MATRIX.
C
C        D CONTAINS THE DIAGONAL ELEMENTS OF THE INPUT MATRIX.
C
C        E CONTAINS THE SUBDIAGONAL ELEMENTS OF THE INPUT MATRIX
C          IN ITS LAST N-1 POSITIONS.  E(1) IS ARBITRARY.
C
C        Z CONTAINS THE TRANSFORMATION MATRIX PRODUCED IN THE
C          REDUCTION BY  TRED2, IF PERFORMED.  IF THE EIGENVECTORS
C          OF THE TRIDIAGONAL MATRIX ARE DESIRED, Z MUST CONTAIN
C          THE IDENTITY MATRIX.
C
C      ON OUTPUT
C
C        D CONTAINS THE EIGENVALUES IN ASCENDING ORDER.  IF AN
C          ERROR EXIT IS MADE, THE EIGENVALUES ARE CORRECT BUT
C          UNORDERED FOR INDICES 1,2,...,IERR-1.
C
C        E HAS BEEN DESTROYED.
C
C        Z CONTAINS ORTHONORMAL EIGENVECTORS OF THE SYMMETRIC
C          TRIDIAGONAL (OR FULL) MATRIX.  IF AN ERROR EXIT IS MADE,
C          Z CONTAINS THE EIGENVECTORS ASSOCIATED WITH THE STORED
C          EIGENVALUES.
C
C        IERR IS SET TO
C          ZERO       FOR NORMAL RETURN,
C          J          IF THE J-TH EIGENVALUE HAS NOT BEEN
C                     DETERMINED AFTER 30 ITERATIONS.
C
C     CALLS PYTHAG FOR  DSQRT(A*A + B*B) .
C
C     QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW,
C     MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY
C
C     THIS VERSION DATED AUGUST 1983.
C
C     ------------------------------------------------------------------
C
      IERR = 0
      IF (N .EQ. 1) GO TO 1001
C
      DO 100 I = 2, N
  100 E(I-1) = E(I)
C
      F = 0.0D0
      TST1 = 0.0D0
      E(N) = 0.0D0
C
      DO 240 L = 1, N
         J = 0
         H = DABS(D(L)) + DABS(E(L))
         IF (TST1 .LT. H) TST1 = H
C     .......... LOOK FOR SMALL SUB-DIAGONAL ELEMENT ..........
         DO 110 M = L, N
            TST2 = TST1 + DABS(E(M))
            IF (TST2 .EQ. TST1) GO TO 120
C     .......... E(N) IS ALWAYS ZERO, SO THERE IS NO EXIT
C                THROUGH THE BOTTOM OF THE LOOP ..........
  110    CONTINUE
C
  120    IF (M .EQ. L) GO TO 220
  130    IF (J .EQ. 30) GO TO 1000
         J = J + 1
C     .......... FORM SHIFT ..........
         L1 = L + 1
         L2 = L1 + 1
         G = D(L)
         P = (D(L1) - G) / (2.0D0 * E(L))
         R = PYTHAG(P,1.0D0)
         D(L) = E(L) / (P + DSIGN(R,P))
         D(L1) = E(L) * (P + DSIGN(R,P))
         DL1 = D(L1)
         H = G - D(L)
         IF (L2 .GT. N) GO TO 145
C
         DO 140 I = L2, N
  140    D(I) = D(I) - H
C
  145    F = F + H
C     .......... QL TRANSFORMATION ..........
         P = D(M)
         C = 1.0D0
         C2 = C
         EL1 = E(L1)
         S = 0.0D0
         MML = M - L
C     .......... FOR I=M-1 STEP -1 UNTIL L DO -- ..........
         DO 200 II = 1, MML
            C3 = C2
            C2 = C
            S2 = S
            I = M - II
            G = C * E(I)
            H = C * P
            R = PYTHAG(P,E(I))
            E(I+1) = S * R
            S = E(I) / R
            C = P / R
            P = C * D(I) - S * G
            D(I+1) = H + S * (C * G + S * D(I))
C     .......... FORM VECTOR ..........
            DO 180 K = 1, N
               H = Z(K,I+1)
               Z(K,I+1) = S * Z(K,I) + C * H
               Z(K,I) = C * Z(K,I) - S * H
  180       CONTINUE
C
  200    CONTINUE
C
         P = -S * S2 * C3 * EL1 * E(L) / DL1
         E(L) = S * P
         D(L) = C * P
         TST2 = TST1 + DABS(E(L))
         IF (TST2 .GT. TST1) GO TO 130
  220    D(L) = D(L) + F
  240 CONTINUE
C     .......... ORDER EIGENVALUES AND EIGENVECTORS ..........
      DO 300 II = 2, N
         I = II - 1
         K = I
         P = D(I)
C
         DO 260 J = II, N
            IF (D(J) .GE. P) GO TO 260
            K = J
            P = D(J)
  260    CONTINUE
C
         IF (K .EQ. I) GO TO 300
         D(K) = D(I)
         D(I) = P
C
         DO 280 J = 1, N
            P = Z(J,I)
            Z(J,I) = Z(J,K)
            Z(J,K) = P
  280    CONTINUE
C
  300 CONTINUE
C
      GO TO 1001
C     .......... SET ERROR -- NO CONVERGENCE TO AN
C                EIGENVALUE AFTER 30 ITERATIONS ..........
 1000 IERR = L
 1001 RETURN
      END
C
      SUBROUTINE OPB(N, X, Y)
C
C.....MULTIPLICATION OF MATRIX B BY A VECTOR X , WHERE
C
C     B =   A'A, WHERE A IS NROW BY NCOL (NROW >> NCOL)
C
C     THE PARAMETER NRMAX SHOULD BE GREATER THAN OR EQUAL TO NROW
C
C     HENCE, B  IS OF ORDER N:=NCOL (Y STORES PRODUCT VECTOR)
C
      PARAMETER(NMAX=500,NRMAX=500,NZMAX=2000)
      REAL*8    VALUE(NZMAX)
      INTEGER POINTR(NMAX), ROWIND(NZMAX)
      COMMON /MATRIXA/ NCOL,NROW,VALUE,POINTR,ROWIND
C
      COMMON /COUNT/ MXVCOUNT
      INTEGER N,MXVCOUNT
      REAL*8 X(1), Y(1), ZTEMP(NRMAX)
C
C---------------------
C
      MXVCOUNT=MXVCOUNT+2
C
      DO 55 I=1,NCOL
         Y(I)=0.0
55    CONTINUE
C
      DO 85 I=1,NROW
         ZTEMP(I)=0.0
85    CONTINUE
C
C.....MULTIPLY BY SPARSE C
C
      DO 15 I=1,NCOL
C
         DO 10 J=POINTR(I),POINTR(I+1)-1
            ZTEMP(ROWIND(J))=ZTEMP(ROWIND(J))+VALUE(J)*X(I)
10       CONTINUE
C
15    CONTINUE
C
C.....MULTIPLY BY SPARSE C'
C
      DO 25 I =1,NCOL
C
         DO 20 J=POINTR(I),POINTR(I+1)-1
            Y(I)=Y(I)+VALUE(J)*ZTEMP(ROWIND(J))
20       CONTINUE
25    CONTINUE
C
      RETURN
      END
C
      SUBROUTINE OPA (X, Y)
C
C.....MULTIPLICATION OF MATRIX A BY VECTOR X , WHERE
C
C     A IS NROW BY NCOL (NROW >> NCOL)
C
C     (Y STORES PRODUCT VECTOR)
C
      PARAMETER(NMAX=500,NZMAX=2000)
      REAL*8    VALUE(NZMAX)
      INTEGER POINTR(NMAX), ROWIND(NZMAX)
      COMMON /MATRIXA/ NCOL,NROW,VALUE,POINTR,ROWIND
C
      COMMON /COUNT/ MXVCOUNT
      INTEGER MXVCOUNT
      REAL*8 X(1), Y(1)
C
C---------------------
C
      MXVCOUNT=MXVCOUNT+1
C
      DO 55 I=1,NROW
         Y(I)=0.0
55    CONTINUE
C
      DO 15 I=1,NCOL
C
         DO 10 J=POINTR(I),POINTR(I+1)-1
            Y(ROWIND(J))=Y(ROWIND(J))+VALUE(J)*X(I)
10    CONTINUE
C
15    CONTINUE
C
      RETURN
      END

.