* Date: Tue, 21 Jul 92 15:49:01 -0700
* From: Cheryl Carey <carey@sccm.Stanford.EDU>
C*****************************************************************
C*****************************************************************
C           AN ITERATIVE BLOCK LANCZOS METHOD FOR THE
C              SOLUTION OF LARGE SPARSE SYMMETRIC
C                       EIGENPROBLEMS
C                            by
C                   Richard Ray Underwood
C
C
C     if further details concerning the content of the code are
C     required, you are advised to refer to R. R. Underwood's
C     PhD thesis - "An Iterative Block Lanczos Method for the
C     Solution of Large Sparse Symmetric Eigenproblems", 1975.
C*****************************************************************
C*****************************************************************
      IMPLICIT REAL*8 (A-H,O-Z)
C=================================================================
C     sample main program. The values for N, Q, PINIT, R, MMAX and
C     EPS are all problem dependent and should be altered
C     accordingly, as too should the dimensions of the arrays D
C     and X.
C=================================================================
      DIMENSION D(25),X(20000)
      EXTERNAL AX
      INTEGER Q,PINIT,R
      N = 300
      Q=12
      PINIT =4
      R=4
      MMAX=5000
      EPS=1D-03
      M=0
C
C------
C     MINVAL is used to compute the 4 least eigenvalues
C     of the matrix A = diag(-1,-1/2,...,-1/300) to an
C     approximate precision of 10**(-3). 12 vectors are
C     allowed for the block lanczos method and an initial
C     block size of 4 is chosen.
C------
      CALL MINVAL(N,Q,PINIT,R,MMAX,EPS,AX,M,D,X,IECODE)
C
      PRINT 1001,M,IECODE,(D(I),I=1,M)
 1001 FORMAT(/' =>M,IECODE=',I4,1X,I4/(' =>E ',5D23.15))
      STOP
      END
C
C
C
C
      SUBROUTINE MINVAL(N,Q,PINIT,R,MMAX,EPS,OP,M,D,X,IECODE)
C=================================================================
C     this subroutine is the main subroutine implementing the
C     iterative block lanczos method for computing the
C     eigenvalues and eigenvectors of symmetric matrices.
C=================================================================
      IMPLICIT REAL*8 (A-H,O-Z)
      INTEGER N,Q,PINIT,R,MMAX
      REAL*8 EPS
      EXTERNAL OP
      DIMENSION D(Q),X(N,Q)
      INTEGER IECODE
      DIMENSION E(25),C(25,25)
      DIMENSION U(1296),V(1296)
      INTEGER P,S,PS
C
C------
C     description of parameters:
C     N:     integer variable. The order of the symmetric
C            matrix A whose eigenvalues and eigenvectors are
C            being computed. The value of N should be less than
C            or equal to 1296 and greater than or equal to  2.
C
C     Q:     integer variable. The number of vectors of length
C            N contained in the array X. The value of Q should
C            be less than or equal to 25, at least one greater
C            than the value of R and less than or equal to N.
C
C     PINIT: integer variable. The initial block size to be used
C            in the block lanczos method. If PINIT is negative,
C            then -PINIT is used for the block size and columns
C            M+L,...,M+(-PINIT) of the array X are assumed to be
C            initialized to the initial matrix used to start the
C            block lanczos method. If the subroutine terminates
C            with a value of M less than R, then PINIT is assigned
C            a value -P where P is the final block size chosen.
C            In this circumstance, columns M+1,...M+P will contain
C            the most recent set of eigenvector approximations
C            which can be used to restart the subroutine if desired.
C
C     R:     integer variable. The number of eigenvalues and
C            eigenvectors being computed. That is, MINVAL attempts
C            to compute accurate approximations to the R least
C            eigenvalues and eigenvectors of the matrix A. The value
C            of R should be greater than zero and less than Q.
C
C     MMAX:  integer variable. The maximum number of matrix-vector
C            products A*X (where X is a vector) that are allowed
C            during one call of this subroutine to complete its task
C            of computing R eigenvalues and eigenvectors. Unless the
C            problem indicates otherwise, MMAX should be given a very
C            large value.
C
C     EPS:   REAL*8 variable. Initially, EPS should contain a value
C            indicating the relative precision to which MINVAL will
C            attempt to compute the eigenvalues and eigenvectors of A.
C            For eigenvalues less in modulus than 1, EPS will be an
C            absolute tolerance. Because of the way this method works,
C            it may happen that the later eigenvalues cannot be
C            computed to the same relative precision as those less in
C            value.
C
C     OP:    subroutine name. The actual argument corresponding to OP
C            should be the name of a subroutine used to define the
C            matrix A. This subroutine should have three arguments
C            N, U, and V, say, where N is an integer variable giving
C            the order of A, and U and V are two one-dimensional
C            arrays of length N. If W denotes the vector of order N
C            stored in U, then the statement
C                    CALL OP(N,U,V)
C            should result in the vector A*W being computed and stored
C            in V. The contents of U can be modified by this call.
C
C     M:     integer variable. M gives the number of eigenvalues and
C            eigenvectors already computed. Thus, initially, M should
C            be zero. If M is greater than zero, then columns one
C            through M of the array X are assumed to contain the
C            computed approximations to the M least eigenvalues and
C            eigenvectors of A. On exit, M contains a value equal to
C            the total number of eigenvalues and eigenvectors
C            computed including any already computed when MINVAL was
C            entered. Thus, on exit, the first M elements of D and the
C            first M columns of X will contain approximations to the
C            M least eigenvalues of A.
C
C     D:     REAL*8 array. D contains the computed eigenvalues. D
C            should be a one-dimensional array with at least Q
C            elements.
C
C     X:     REAL*8 array. X contains the computed eigenvectors. X
C            should be an array containing at least N*Q elements. X
C            is used not only to store the eigenvectors computed by
C            MINVAL, but also as working storage for the block lanczos
C            method. On exit, the first N*M elements of X contain the
C            eigenvector approximations - the first vector in the
C            first N elements, the second in the second N elements,
C            etc...
C
C
C     IECODE:integer variable. The value of IECODE indicates whether
C            MINVAL terminated successfully, and if not, the reason
C            why.
C
C               IECODE=0 : successful termination.
C               IECODE=1 : the value of N is less than 2.
C               IECODE=2 : the value of N exceeds 1296.
C               IECODE=3 : the value of R is less than 1.
C               IECODE=4 : the value of Q is less than or equal to R.
C               IECODE=5 : the value of Q is greater than 25.
C               IECODE=6 : the value of Q exceeds N.
C               IECODE=7 : the value of MMAX was exceeded before R
C                          eigenvalues and eigenvectors were
C                          computed.
C
C      Note that the subroutine has been designed to allow initial
C      approximations to the eigenvectors corresponding to the least
C      eigenvalues to be utilised if they are known (by storing them
C      in X and assigning PINIT minus the value of their number).
C      Furthermore, it has also been designed to allow restarting if
C      it stops with IECODE=7. Thus, the user of this program can
C      restart it after examining any partial results without loss of
C      previous work.
C------
C
C------
C     check that the initial values of the subroutine parameters
C     are in range.
C------
      IF (N.LT.2) GO TO 901
      IF (N.GT.1296) GO TO 902
      IF (R.LT.1) GO TO 903
      IF (Q.LE.R) GO TO 904
      IF (Q.GT.25) GO TO 905
      IF (Q.GT.N) GO TO 906
C
C------
C     choose initial values for the block size P, the number of
C     steps that the block lanczos method is carried out, and
C     choose an initial N-by-P orthonormal matrix X1 used to start
C     the block lanczos method.
C------
      P = PINIT
      IF (P.LT.0) P = -P
      S = (Q-M)/P
      IF (S.GT.2) GO TO 100
      S = 2
      P = Q/2
  100 IF (PINIT.LT.0) GO TO 150
      DO 120 K = 1,P
             CALL RANDOM(N,Q,K,X)
  120 CONTINUE
  150 IF (M.GT.0) GO TO 200
      CALL ORTHG(N,Q,M,P,C,X)
C
C------
C     rotate the initial N-by-P matrix X1 so that
C        X1'*A*X1=diag(D1,D2,...,DP)
C     where DI is stored in D(I), I=1,...,P.
C------
      CALL SECTN(N,Q,M,P,OP,X,C,D,U,V)
      ERRC = 0.D0
  200 ITER = 0
      IMM=0
C
C------
C     the main body of the subroutine starts here. IMM
C     counts the number of matrix-vector products computed,
C     which is the number of times the subroutine named by
C     OP is called. ERRC measures the accumulated error in
C     the eigenvalues and eigenvectors.
C------
C
  300 IF (M.GE.R) GO TO 900
      IF (IMM.GT.MMAX) GO TO 907
      ITER = ITER + 1
      PS = P*S
C
C------
C     BKLANC carries out the block lanczos method and stores
C     the resulting block tridiagonal matrix MS in C and the
C     N-by-PS orthonormal matrix XS in X. The initial N-by-P
C     orthonormal matrix is assumed to be stored in columns
C     M+1 through M+PS of X. The residuals for these vectors
C     and the eigenvalue approximations in D are computed and
C     stored in E.
C------
      CALL BKLANC(N,Q,M,P,S,OP,D,C,X,E,U,V)
C
C------
C     EIGEN solves the eigenproblem for MS, storing the eigenvalues
C     in elements M+1 through M+PS of D and the eigenvectors in the
C     first P*S rows and columns of C (overwriting MS, possibly).
C------
      CALL EIGEN(Q,M,P,PS,C,D)
C
C------
C     CNVTST determines how many of the eigenvalues and eigenvectors
C     have converged using the error estimates stored in E. The number
C     that have converged is stored in NCONV. If NCONV=0, then none
C     have converged.
C------
      CALL CNVTST(N,Q,M,P,ERRC,EPS,D,E,NCONV)
      write(6,*) nconv
C
C------
C     PCH chooses new values for P and S, the block size and the
C     number of steps for the block lanczos subprogram, respectively.
C------
      CALL PCH(N,Q,M,R,NCONV,P,S)
C
C------
C     ROTATE computes the eigenvectors of the restricted matrix
C     using XS stored in X and the eigenvectors of MS stored in C.
C     These vectors serve both as eigenvector approximations and
C     to form the matrix used to start the block lanczos method in
C     the next iteration.
C------
      CALL ROTATE(N,Q,M,PS,NCONV+P,C,X)
C
      M = M + NCONV
      IMM=IMM+P*S
      PRINT 1001,ITER,IMM,P,PS
 1001 FORMAT(' =>ITER,IMM,P,PS =',4I5)
      GO TO 300
C
C------
C     this is the end of the main body of the subroutine. Now set
C     the value of IECODE and EXIT.
C------
  900 IECODE = 0
      RETURN
  901 IECODE = 1
      RETURN
  902 IECODE = 2
      RETURN
  903 IECODE = 3
      RETURN
  904 IECODE = 4
      RETURN
  905 IECODE = 5
      RETURN
  906 IECODE = 6
      RETURN
  907 IECODE = 7
      PINIT = -P
C
      RETURN
      END
C
C
C
C
      SUBROUTINE BKLANC(N,Q,M,P,S,OP,D,C,X,E,U,V)
C====================================================================
C     this subroutine implements the block lanczos method
C     with reorthogonalization. BKLANC computes a block
C     tridiagonal matrix MS which it stores in rows and
C     columns M+1 through M+P*S of the array C, and an
C     orthonormal matrix XS which it stores in columns M+1
C     through M+P*S of the N-by-Q array X. MS is a symmetric
C     matrix with P-by-P symmetric matrices M(1),...,M(S) on
C     its diagonal and P-by-P upper triangular matrices
C     R(2),...,R(S) along its lower diagonal. Since MS is
C     symmetric and banded, only its lower triangle (P+1
C     diagonals) is stored in C. XS is composed of S N-by-P
C     orthonormal matrices X(1),...,X(S) where X(1) is given
C     and should be stored in columns M+1 through M+P of X.
C     Furthermore, X(1) is assumed to satisfy X(1)*A*X(1) =
C     diag(D(M+1),D(M+2),...,D(M+P)), and if M>0, then X(1) is
C     assumed to be orthogonal to the vectors stored in columns
C     1 through M of X. OP is the name of an external subroutine
C     used to define the matrix A. During the first step, the
C     subroutine ERR is called and the quantities EJ are computed
C     where EJ=||A*X1J-D(M+J)*X1J||, X1J is the J-th column of x(1),
C     and ||*|| denotes the Euclidean norm. EJ is stored in E(M+J),
C     J=1,2,...,P. U and V are auxilliary vectors used by OP.
C====================================================================
      IMPLICIT REAL*8 (A-H,O-Z)
      INTEGER N,Q,M,P,S
      DIMENSION D(Q),C(Q,Q),X(N,Q)
      DIMENSION E(Q),U(N),V(N)
C
      MP1 = M+1
      MPPS = M+P*S
      DO 90 L = 1,S
            LL = M+(L-1)*P+1
            LU = M + L*P
            DO 70 K = LL,LU
                  DO 10 I = 1,N
   10                   U(I) = X(I,K)
                  CALL OP(N,U,V)
                  IF (L.GT.1) GO TO 19
                  DO 12 I=K,LU
   12                   C(I,K) = 0.D0
                  C(K,K) = D(K)
                  DO 14 I = 1,N
   14                   V(I) = V(I) - D(K)*X(I,K)
                  GO TO 61
   19             DO 30 I = K,LU
                        T = 0
                        DO 20 J = 1,N
   20                         T = T+V(J)*X(J,I)
   30             C(I,K) = T
                  IT = K-P
                  DO 60 I = 1,N
                        T = 0
                        DO 40 J = IT,K
   40                         T = T + X(I,J)*C(K,J)
                  IF (K.EQ.LU) GO TO 60
                  KP1 = K+1
                  DO 50 J = KP1,LU
   50             T = T+X(I,J)*C(J,K)
   60             V(I) = V(I)-T
   61             IF (L.EQ.S) GO TO 70
                  DO 63 I = 1,N
   63                   X(I,K+P) = V(I)
   70    CONTINUE
         IF (L.EQ.1) CALL ERR(N,Q,M,P,X,E)
         IF (L.EQ.S) GO TO 90
         CALL ORTHG(N,Q,LU,P,C,X)
         IL = LU+1
         IT = LU
         DO 80 J = 1,P
               IT = IT+1
               DO 80 I = IL,IT
   80                C(I,IT-P) = C(I,IT)
   90 CONTINUE
C
      RETURN
      END
C
C
C
C
      SUBROUTINE PCH(N,Q,M,R,NCONV,P,S)
C====================================================================
C     based on the values of N, Q, M, R and NCONV, PCH chooses new
C     values for P and S, the block size and number of steps for the
C     block lanczos method. The strategy used here is to choose P to
C     be the smaller of the two following values:
C            1) the previous block size
C     and,   2) the number of values left to be computed. S is chosen
C     as large as possible subject to the constraints imposed by the
C     limits of storage. In any event, S is greater than or equal to
C     2. N is the order of the problem and Q is the number of vectors
C     available for storing eigenvectors and applying the block
C     lanczos method. M is the number of eigenvalues and eigenvectors
C     that have already been computed and R is the required number.
C     Finally, NCONV is the number of eigenvalues and eigenvectors
C     that have converged in the current iteration.
C====================================================================
      IMPLICIT REAL*8 (A-H,O-Z)
      INTEGER N,Q,M,R,NCONV,P,S
      INTEGER PT,ST
C
      MT = M+NCONV
      PT = R-MT
      IF (PT.GT.P) PT = P
      IF (PT.GT.0) GO TO 101
      P = 0
      RETURN
  101 CONTINUE
      ST = (Q-MT)/PT
      IF (ST.GT.2) GO TO 110
      ST = 2
      PT = (Q-MT)/2
  110 P = PT
      S = ST
C
      RETURN
      END
C
C
C
C
      SUBROUTINE ERR(N,Q,M,P,X,E)
C================================================
C     ERR computes the Euclidean lengths of the
C     vectors stored in the columns M+P+1 through
C     M+P+P of the N-by-Q array X and stores them
C     in elements M+1 through M+P of E.
C================================================
      IMPLICIT REAL*8 (A-H,O-Z)
      INTEGER N,Q,M,P
      DIMENSION X(N,Q),E(Q)
C
      MP1 = M+P+1
      MPP = M+P+P
      DO 200 K = MP1,MPP
             T = 0.D0
             DO 100 I = 1,N
  100               T = T+X(I,K)**2
  200        E(K-P) = DSQRT(T)
C
      RETURN
      END
C
C
C
C
      SUBROUTINE CNVTST(N,Q,M,P,ERRC,EPS,D,E,NCONV)
C===========================================================
C     CNVTST determines which of the P eigenvalues stored
C     in elements M+1 through M+P of D have converged. ERRC
C     is a measure of the accumulated error in the M
C     previously computed eigenvalues and eigenvectors. ERRC
C     is updated if more approximations have converged. The
C     norms of the residual vectors are stored in elements
C     M+1 through M+P of E. EPS is the precision to which we
C     are computing the approximations. Finally, NCONV is the
C     number that have converged. If NCONV=0, then none have
C     converged.
C============================================================
      IMPLICIT REAL*8 (A-H,O-Z)
      INTEGER Q,M,P
      REAL*8 EPS
      DIMENSION D(Q),E(Q)
      DATA CHEPS/2.22D-16/
C
      K = 0
      DO 100 I = 1,P
         T = DABS(D(M+I))
         IF (T.LT.1.D0) T = 1.D0
      IF (E(M+I).GT.T*(EPS+10.D0*N*CHEPS)+ERRC) GO TO 200
  100 K = I
  200 NCONV = K
      IF (K.EQ.0) RETURN
      T = 0.D0
      DO 300 I = 1,K
             T = T+E(M+I)**2
  300 CONTINUE
      ERRC = DSQRT(ERRC**2+T)
C
      RETURN
      END
C
C
C
C
      SUBROUTINE EIGEN(Q,M,P,PS,C,D)
C=======================================================
C     EIGEN solves the eigenproblem for the symmetric
C     matrix MS of order PS stored in rows and columns
C     M+1 through M+PS of C. The eigenvalues of MS are
C     stored in elements M+1 through M+PS of D and the
C     eigenvactors are stored in rows and columns 1
C     through PS of C possibly overwriting MS. EIGEN
C     simply re-stores MS in a manner acceptable to the
C     subroutines TRED2 and TQL2. These two routines are
C     available through eispack.
C=======================================================
      IMPLICIT REAL*8 (A-H,O-Z)
      INTEGER M,P,Q,PS,PP1
      DIMENSION C(Q,Q),D(Q)
      DIMENSION DD(25),V(25)
C
      DO 150 I=1,PS
             LIM = I-P
             IF (I.LE.P) LIM = 1
             IF (LIM.LE.1) GO TO 130
             LM1 = LIM-1
         DO 120 J = 1,LM1
  120           C(I,J) = 0.D0
  130    DO 140 J = LIM,I
  140           C(I,J) = C(I+M,J+M)
  150 CONTINUE
C
C------
C     call eispack routines here
C------
      CALL TRED2(Q,PS,C,DD,V,C)
      CALL TQL2(Q,PS,DD,V,C,IERR)
C
      PRINT 601,PS,(DD(I),I=1,PS)
  601 FORMAT(' => ORDER =',I4,/,('    EIGENVALUES =',10d10.4))
      PP1 = P+1
      DO 155 J = 1,PP1
  155 PRINT 602,J,(C(I,J),I=1,PS)
  602 FORMAT(' => J =',I4,/,('    EIGENVECTORS =',10d10.4))
      DO 160 I = 1,PS
  160        D(M+I) = DD(I)
C
      RETURN
      END
C
C
C
C
      SUBROUTINE SECTN(N,Q,M,P,OP,X,C,D,U,V)
C==============================================================
C     SECTN transforms the N-by-P orthonormal matrix X1,
C     say, stored in columns M+1 through M+P of the N-by-Q
C     array X so that X1'*A*X1 = diag(D1,D2,...,DP), where
C     ' denotes transpose and A is a symmetric matrix of
C     order N defined by the subroutine OP. The values D1,...
C     ,DP are stored in elements M+1 through M+P of D. SECTN
C     forms the matrix X1'*A*X1 = CP, storing CP in the array
C     C. The values D1,D2,...,DP and the eigenvectors QP of CP
C     are computed by EIGEN and stored in D and C respectively.
C     ROTATE then carries out the transformation X1<=X1*QP.
C==============================================================
      IMPLICIT REAL*8 (A-H,O-Z)
      INTEGER N,Q,M,P
      EXTERNAL OP
      DIMENSION X(N,Q),C(Q,Q),D(Q),U(N),V(N)
C
      ICOL1 = M
      DO 300 J = 1,P
             ICOL1 = ICOL1 + 1
         DO 100 I = 1,N
  100           U(I) = X(I,ICOL1)
                CALL OP(N,U,V)
                ICOL2 = M
             DO 300 I = 1,J
                    ICOL2 = ICOL2+1
                    T = 0.D0
                DO 200 K = 1,N
  200                  T = T + V(K)*X(K,ICOL2)
  300 C(ICOL1,ICOL2) = T
      CALL EIGEN(Q,M,P,P,C,D)
      CALL ROTATE(N,Q,M,P,P,C,X)
C
      RETURN
      END
C
C
C
C
      SUBROUTINE ROTATE(N,Q,M,PS,L,C,X)
C======================================================
C     ROTATE computes the first L columns of the matrix
C     XS*QS where XS is an N-by-PS orthonormal matrix
C     stored in columns M+1 through M+PS of the N-by-Q
C     array X and QS is a PS-by-PS orthonormal matrix
C     stored in rows and columns 1 through Ps of the
C     array C. The result is stored in columns M+1
C     through M+L of X overwriting part of XS.
C======================================================
      IMPLICIT REAL*8 (A-H,O-Z)
      INTEGER N,Q,M,PS,L
      DIMENSION C(Q,Q),X(N,Q)
      DIMENSION V(25)
C
      DO 300 I = 1,N
         DO 200 K = 1,L
                T = 0.D0
            DO 100 J = 1,PS
  100              T = T+X(I,M+J)*C(J,K)
  200       V(K) = T
         DO 300 K = 1,L
  300 X(I,M+K) = V(K)
C
      RETURN
      END
C
C
C
C
      SUBROUTINE ORTHG(N,Q,F,P,B,X)
C=========================================================
C     ORTHG reorthogonalizes the N-by-P matrix Z stored
C     in columns F+1 through F+P of the N-by-Q array X
C     with respect to the vectors stored in the first F
C     columns of X and then decomposes the resulting
C     matrix into the product of an N-by-P orthonormal
C     matrix XORTH, say, and a P-by-P upper triangular
C     matrix R. XORTH is stored over Z and the upper
C     triangle of R is stored in rows and columns F+1
C     through F+P of the Q-by-Q array B. A stable variant
C     of the Gram-Schmidt orthogonalization method is
C     utilised.
C=========================================================
      IMPLICIT REAL*8 (A-H,O-Z)
      INTEGER N,Q,F,P
      DIMENSION B(Q,Q),X(N,Q)
      INTEGER FP1,FPP
      LOGICAL ORIG
C
      IF (P.EQ.0) RETURN
      FP1 = F+1
      FPP = F+P
      DO 50 K = FP1,FPP
            ORIG = .TRUE.
            KM1 = K-1
   10       T = 0.D0
            IF (KM1.LT.1) GO TO 25
         DO 20 I = 1,KM1
               S = 0.D0
            DO 15 J = 1,N
   15          S = S+X(J,I)*X(J,K)
            IF (ORIG.AND.I.GT.F) B(I,K) = S
            T = T + S*S
         DO 20 J = 1,N
   20          X(J,K) = X(J,K) - S*X(J,I)
   25    S = 0.D0
         DO 30 J = 1,N
   30          S = S + X(J,K)*X(J,K)
         T = T+S
         IF (S.GT.T/100.D0) GO TO 40
         ORIG = .FALSE.
         GO TO 10
   40    S = DSQRT(S)
         B(K,K) = S
         IF (S.NE.0) S = 1.D0/S
         DO 50 J = 1,N
   50 X(J,K) = S*X(J,K)
C
      RETURN
      END
C
C
C
C
      SUBROUTINE RANDOM(N,Q,L,X)
C======================================================
C     RANDOM computes and stores a sequence of N
C     pseudo-random numbers in the L-th column of the
C     N-by-Q array X. RANDOM generates two sequences of
C     pseudo-random numbers, filling an array with one
C     sequence and using the second to access the array
C     in a random fashion.
C======================================================
      IMPLICIT REAL*8 (A-H,O-Z)
      INTEGER N,Q,L,FT,X1,F1,F2
      INTEGER A,C,X0
      DIMENSION X(N,Q)
      DIMENSION T(100)
      DATA F1/71416/,F2/27183/
      DATA A/6821/,C/5327/,X0/5328/
C
      DO 100 I = 1,100
            X1 = A *X0+C
            IF (X1.GE.10000) X1 = X1 - 10000
            T(I) = X1/9999.D0 - 0.5D0
  100 X0 = X1
      DO 200 I = 1,N
             FT = F1+F2
             IF (FT.GE.1000000) FT = FT-1000000
             F1 = F2
             F2 = FT
             K = FT/1.D6*100 + 1
             X(I,L) = T(K)
             X1 = A*X0 + C
             IF (X1.GE.10000) X1 = X1 - 10000
             T(K) = X1/9999D0 - 0.5D0
  200    X0 = X1
C
      RETURN
      END
C
C
C
C
      SUBROUTINE AX(N,U,V)
C======================================================
C     user dependent subroutine. It is defined to be
C     the actual argument corresponding to OP and is
C     used to define the matrix A. This subroutine
C     should have three arguments N, U and V, say,
C     where N is an integer variable giving the order
C     of A, and U and V are two one-dimensional arrays
C     of length N.
C======================================================
      IMPLICIT REAL*8(A-H,O-Z)
      INTEGER N
      DIMENSION U(N),V(N)
C
C------
C     in this example, the matrix A is defined to be
C         A = diag(-1,-1/2,...,-1/300).
C------
      DO 100 I=1,N
             V(I)=(-1.D0/I)*U(I)
  100 CONTINUE
C
      RETURN
      END
C
C
C
C