*     F11DDF Example Program Text
*     Mark 20 Revised. NAG Copyright 2001.
*     .. Parameters ..
      INTEGER          NIN, NOUT
      PARAMETER        (NIN=5,NOUT=6)
      INTEGER          NMAX, LA, LIWORK, LWORK
      PARAMETER        (NMAX=1000,LA=10000,LIWORK=2*NMAX+1,LWORK=10000)
*     .. Local Scalars ..
      DOUBLE PRECISION ANORM, OMEGA, SIGMAX, STPLHS, STPRHS, TOL
      INTEGER          I, IFAIL, IREVCM, ITERM, ITN, LWNEED, LWREQ, M,
     +                 MAXITN, MONIT, N, NNZ
      CHARACTER        CKDDF, CKXAF, NORM, PRECON, TRANS, WEIGHT
      CHARACTER*8      METHOD
*     .. Local Arrays ..
      DOUBLE PRECISION A(LA), B(NMAX), RDIAG(NMAX), WGT(NMAX),
     +                 WORK(LWORK), X(NMAX)
      INTEGER          ICOL(LA), IROW(LA), IWORK(LIWORK)
*     .. External Subroutines ..
      EXTERNAL         F11BDF, F11BEF, F11BFF, F11DDF, F11XAF
*     .. Intrinsic Functions ..
      INTRINSIC        MAX
*     .. Executable Statements ..
      WRITE (NOUT,*) 'F11DDF Example Program Results'
      WRITE (NOUT,*)
*     Skip heading in data file
      READ (NIN,*)
*
*     Read algorithmic parameters
*
      READ (NIN,*) N
      IF (N.LE.NMAX) THEN
         READ (NIN,*) NNZ
         READ (NIN,*) METHOD
         READ (NIN,*) PRECON, NORM, ITERM
         READ (NIN,*) M, TOL, MAXITN
         READ (NIN,*) ANORM, SIGMAX
         READ (NIN,*) OMEGA
*
*        Check size of workspace
*
         LWREQ = MAX(N*(M+3)+M*(M+5)+101,7*N+100,(2*N+M)*(M+2)+N+100,
     +           10*N+100)
         IF (LWORK.LT.LWREQ) THEN
            WRITE (NOUT,*) 'LWORK must be at least', LWREQ
            STOP
         END IF
*
*        Read the matrix A
*
         DO 20 I = 1, NNZ
            READ (NIN,*) A(I), IROW(I), ICOL(I)
   20    CONTINUE
*
*        Read right-hand side vector b and initial approximate solution x
*
         READ (NIN,*) (B(I),I=1,N)
         READ (NIN,*) (X(I),I=1,N)
*
*        Call F11BDF to initialize solver
*
         WEIGHT = 'N'
         MONIT = 0
         IFAIL = 0
         CALL F11BDF(METHOD,PRECON,NORM,WEIGHT,ITERM,N,M,TOL,MAXITN,
     +               ANORM,SIGMAX,MONIT,LWNEED,WORK,LWORK,IFAIL)
*
*        Calculate reciprocal diagonal matrix elements if necessary
*
         IF (PRECON.EQ.'P' .OR. PRECON.EQ.'p') THEN
*
            DO 40 I = 1, N
               IWORK(I) = 0
   40       CONTINUE
*
            DO 60 I = 1, NNZ
               IF (IROW(I).EQ.ICOL(I)) THEN
                  IWORK(IROW(I)) = IWORK(IROW(I)) + 1
                  IF (A(I).NE.0.D0) THEN
                     RDIAG(IROW(I)) = 1.D0/A(I)
                  ELSE
                     WRITE (NOUT,*)
     +                 'Matrix has a zero diagonal element'
                     GO TO 140
                  END IF
               END IF
   60       CONTINUE
*
            DO 80 I = 1, N
               IF (IWORK(I).EQ.0) THEN
                  WRITE (NOUT,*)
     +              'Matrix has a missing diagonal element'
                  GO TO 140
               END IF
               IF (IWORK(I).GE.2) THEN
                  WRITE (NOUT,*)
     +              'Matrix has a multiple diagonal element'
                  GO TO 140
               END IF
   80       CONTINUE
*
         END IF
*
*        Call F11BEF to solve the linear system
*
         IREVCM = 0
         CKXAF = 'C'
         CKDDF = 'C'
*
  100    CONTINUE
*
         CALL F11BEF(IREVCM,X,B,WGT,WORK,LWORK,IFAIL)
*
         IF (IREVCM.EQ.1) THEN
*
*           Compute matrix-vector product
*
            TRANS = 'N'
            CALL F11XAF(TRANS,N,NNZ,A,IROW,ICOL,CKXAF,X,B,IFAIL)
            CKXAF = 'N'
            GO TO 100
*
         ELSE IF (IREVCM.EQ.-1) THEN
*
*           Compute transposed matrix-vector product
*
            TRANS = 'T'
            CALL F11XAF(TRANS,N,NNZ,A,IROW,ICOL,CKXAF,X,B,IFAIL)
            CKXAF = 'N'
            GO TO 100
*
         ELSE IF (IREVCM.EQ.2) THEN
*
*           SSOR preconditioning
*
            TRANS = 'N'
            CALL F11DDF(TRANS,N,NNZ,A,IROW,ICOL,RDIAG,OMEGA,CKDDF,X,B,
     +                  IWORK,IFAIL)
            CKDDF = 'N'
            GO TO 100
*
         ELSE IF (IREVCM.EQ.4) THEN
*
*           Termination
*
            CALL F11BFF(ITN,STPLHS,STPRHS,ANORM,SIGMAX,WORK,LWORK,IFAIL)
*
            WRITE (NOUT,'(A,I10,A)') ' Converged in', ITN,
     +        ' iterations'
            WRITE (NOUT,'(A,1P,D16.3)') ' Matrix norm =', ANORM
            WRITE (NOUT,'(A,1P,D16.3)') ' Final residual norm =',
     +        STPLHS
            WRITE (NOUT,*)
*
*           Output x
*
            WRITE (NOUT,*) '           X'
            DO 120 I = 1, N
               WRITE (NOUT,'(1X,1P,D16.4)') X(I)
  120       CONTINUE
*
         END IF
*
  140    CONTINUE
*
      END IF
      STOP
      END