*     F11DRF 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        CKDRF, CKXNF, NORM, PRECON, TRANS, WEIGHT
      CHARACTER*8      METHOD
*     .. Local Arrays ..
      COMPLEX *16      A(LA), B(NMAX), RDIAG(NMAX), WORK(LWORK), X(NMAX)
      DOUBLE PRECISION WGT(NMAX)
      INTEGER          ICOL(LA), IROW(LA), IWORK(LIWORK)
*     .. External Subroutines ..
      EXTERNAL         F11BRF, F11BSF, F11BTF, F11DRF, F11XNF
*     .. Intrinsic Functions ..
      INTRINSIC        MAX
*     .. Executable Statements ..
      WRITE (NOUT,*) 'F11DRF 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(121+N*(3+M)+M*(M+5),120+7*N,120+(2*N+M)*(M+2)+2*N,
     +           120+10*N)
         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 rhs vector b and initial approximate solution x
*
         READ (NIN,*) (B(I),I=1,N)
         READ (NIN,*) (X(I),I=1,N)
*
*        Call F11BRF to initialize solver
*
         WEIGHT = 'N'
         MONIT = 0
         IFAIL = 0
         CALL F11BRF(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.0D0) THEN
                     RDIAG(IROW(I)) = 1.0D0/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 F11BSF to solve the linear system
*
         IREVCM = 0
         CKXNF = 'C'
         CKDRF = 'C'
*
  100    CONTINUE
*
         CALL F11BSF(IREVCM,X,B,WGT,WORK,LWORK,IFAIL)
*
         IF (IREVCM.EQ.1) THEN
*
*           Compute matrix-vector product
*
            TRANS = 'N'
            CALL F11XNF(TRANS,N,NNZ,A,IROW,ICOL,CKXNF,X,B,IFAIL)
            CKXNF = 'N'
            GO TO 100
*
         ELSE IF (IREVCM.EQ.-1) THEN
*
*           Compute conjugate transposed matrix-vector product
*
            TRANS = 'T'
            CALL F11XNF(TRANS,N,NNZ,A,IROW,ICOL,CKXNF,X,B,IFAIL)
            CKXNF = 'N'
            GO TO 100
*
         ELSE IF (IREVCM.EQ.2) THEN
*
*           SSOR preconditioning
*
            TRANS = 'N'
            CALL F11DRF(TRANS,N,NNZ,A,IROW,ICOL,RDIAG,OMEGA,CKDRF,X,B,
     +                  IWORK,IFAIL)
            CKDRF = 'N'
            GO TO 100
*
         ELSE IF (IREVCM.EQ.4) THEN
*
*           Termination
*
            CALL F11BTF(ITN,STPLHS,STPRHS,ANORM,SIGMAX,WORK,LWORK,IFAIL)
*
            WRITE (NOUT,'(1X,A,I10,A)') 'Converged in', ITN,
     +        ' iterations'
            WRITE (NOUT,'(1X,A,1P,D16.3)') 'Matrix norm =', ANORM
            WRITE (NOUT,'(1X,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,'','',1P,D16.4,'')'')') X(I)
  120       CONTINUE
*
         END IF
*
  140    CONTINUE
*
      END IF
      STOP
      END