*     G13DSF Example Program Text
*     Mark 15 Revised. NAG Copyright 1991.
*     .. Parameters ..
      INTEGER          NIN, NOUT
      PARAMETER        (NIN=5,NOUT=6)
      INTEGER          KMAX, IPMAX, IQMAX, NMAX, NPARMX, LWORK, LIW,
     +                 MMAX, LDRCM, LDCM
      PARAMETER        (KMAX=3,IPMAX=3,IQMAX=3,NMAX=100,
     +                 NPARMX=(IPMAX+IQMAX)*KMAX*KMAX+KMAX,LWORK=2000,
     +                 LIW=100,MMAX=20,LDRCM=MMAX*KMAX*KMAX,LDCM=NPARMX)
*     .. Local Scalars ..
      DOUBLE PRECISION CGETOL, CHI, RLOGL, SIGLEV
      INTEGER          I, IDF, IFAIL, IP, IPRINT, IQ, ISHOW, J, K, M,
     +                 MAXCAL, N, NITER, NPAR, OUTCHN
      LOGICAL          EXACT, MEAN
*     .. Local Arrays ..
      DOUBLE PRECISION CM(LDCM,NPARMX), G(NPARMX), PAR(NPARMX),
     +                 QQ(KMAX,KMAX), R(KMAX,KMAX,MMAX), R0(KMAX,KMAX),
     +                 RCM(LDRCM,MMAX*KMAX*KMAX), V(KMAX,NMAX),
     +                 W(KMAX,NMAX), WORK(LWORK)
      INTEGER          IW(LIW)
      LOGICAL          PARHLD(NPARMX)
*     .. External Subroutines ..
      EXTERNAL         G13DDF, G13DSF, X04ABF
*     .. Executable Statements ..
      WRITE (NOUT,*) 'G13DSF Example Program Results'
      OUTCHN = NOUT
*     Skip heading in data file
      READ (NIN,*)
      READ (NIN,*) K, N
*
      CALL X04ABF(1,OUTCHN)
*
      IF (K.GT.0 .AND. K.LE.KMAX .AND. N.GE.3 .AND. N.LE.NMAX) THEN
         DO 20 I = 1, K
            READ (NIN,*) (W(I,J),J=1,N)
   20    CONTINUE
         READ (NIN,*) IP, IQ, MEAN, M
         IF (IP.GE.0 .AND. IP.LE.IPMAX .AND. IQ.GE.0 .AND. IQ.LE.IQMAX)
     +       THEN
            NPAR = (IP+IQ)*K*K
            IF (MEAN) NPAR = NPAR + K
            IF (NPAR.LE.NPARMX) THEN
               DO 40 I = 1, NPAR
                  PAR(I) = 0.0D0
                  PARHLD(I) = .FALSE.
   40          CONTINUE
               DO 80 J = 1, K
                  DO 60 I = J, K
                     QQ(I,J) = 0.0D0
   60             CONTINUE
   80          CONTINUE
               PARHLD(3) = .TRUE.
               EXACT = .TRUE.
*              ** Set IPRINT > 0 to obtain intermediate output **
               IPRINT = -1
               CGETOL = 0.0001D0
               MAXCAL = 40*NPAR*(NPAR+5)
               ISHOW = 2
               IFAIL = 1
*
               CALL G13DDF(K,N,IP,IQ,MEAN,PAR,NPAR,QQ,KMAX,W,PARHLD,
     +                     EXACT,IPRINT,CGETOL,MAXCAL,ISHOW,NITER,RLOGL,
     +                     V,G,CM,LDCM,IFAIL)
*
               WRITE (NOUT,*)
               IF (IFAIL.LT.0) THEN
                  WRITE (NOUT,*)
                  WRITE (NOUT,99999) ' ** G13DDF returned with'//
     +              ' IFAIL = ', IFAIL
               ELSE
                  IF (IFAIL.NE.0) THEN
                     WRITE (NOUT,99999) 'G13DDF fails. IFAIL =', IFAIL
                     WRITE (NOUT,*)
                  END IF
                  IF ((IFAIL.EQ.0 .OR. IFAIL.GE.4) .AND. M.LE.MMAX) THEN
                     WRITE (NOUT,*) 'Output from G13DSF'
                     WRITE (NOUT,*)
                     ISHOW = 1
                     IFAIL = -1
*
                     CALL G13DSF(K,N,V,KMAX,IP,IQ,M,PAR,PARHLD,QQ,ISHOW,
     +                           R0,R,RCM,LDRCM,CHI,IDF,SIGLEV,IW,LIW,
     +                           WORK,LWORK,IFAIL)
*
                     IF (IFAIL.NE.0) WRITE (NOUT,99999)
     +                   'G13DSF fails. IFAIL =', IFAIL
                  END IF
               END IF
            END IF
         END IF
      END IF
*
99999 FORMAT (1X,A,I5)
      END