c
c
c*****************************************************************************
c*                                                                           *
c*  Copyright Compaq Computer Corporation 1999.         All rights reserved. *
c*                                                                           *
c*****************************************************************************
c
c
C
C
C
         PROGRAM EXAMPLE_SIG_PROC
C
C ***** TO DEMONSTRATE THE USE OF THE SIGNAL PROCESSING ROUTINES.
C 
C ***** THIS PROGRAM ILLUSTRATES THE FOLLOWING:
C
C       (1) USE OF THE SINGLE PRECISION FFT ROUTINE.
C       (2) USE OF THE DOUBLE PRECISION FFT ROUTINE.
C       (3) USE OF THE SINGLE PRECISION, COMPLEX 3-D FFT ROUTINE.
C       (4) USE OF THE DOUBLE PRECISION, COMPLEX 3-D FFT ROUTINE.
C       (5) USE OF THE SINGLE PRECISION CONVOLUTION ROUTINE.
C
c
c ***** sample output follows:
c
c  ILLUSTRATING THE USE OF SIGNAL PROCESSING ROUTINES
c  ILLUSTRATING THE USE OF THE SINGLE PRECISION FFT USING RANDOM INPUT DATA:
c    PERFORM A FORWARD TRANSFORM (REAL INPUT, COMPLEX OUTPUT),
c    FOLLOWED BY AN INVERSE TRANSFORM (COMPLEX INPUT, REAL OUTPUT)
c       MAXIMUM ERROR BETWEEN INPUT AND OUTPUT = 0.23841858E-06
c    PERFORM A FORWARD TRANSFORM (REAL INPUT, REAL OUTPUT),
c    FOLLOWED BY AN INVERSE TRANSFORM (REAL INPUT, REAL OUTPUT)
c       MAXIMUM ERROR BETWEEN INPUT AND OUTPUT = 0.23841858E-06
c  ILLUSTRATING THE USE OF THE DOUBLE PRECISION FFT USING RANDOM INPUT DATA:
c    PERFORM A FORWARD TRANSFORM (REAL INPUT, COMPLEX OUTPUT),
c    FOLLOWED BY AN INVERSE TRANSFORM (COMPLEX INPUT, REAL OUTPUT)
c       MAXIMUM ERROR BETWEEN INPUT AND OUTPUT = 0.30531133E-15
c    PERFORM A FORWARD TRANSFORM (REAL INPUT, COMPLEX OUTPUT),
c    FOLLOWED BY AN INVERSE TRANSFORM (COMPLEX INPUT, REAL OUTPUT)
c       MAXIMUM ERROR BETWEEN INPUT AND OUTPUT = 0.30531133E-15
c  ILLUSTRATING THE USE OF THE SINGLE PRECISION COMPLEX 3-D FFT USING
c  RANDOM INPUT DATA:
c    PERFORM A FORWARD TRANSFORM (COMPLEX INPUT, COMPLEX OUTPUT),
c    FOLLOWED BY AN INVERSE TRANSFORM (COMPLEX INPUT, COMPLEX OUTPUT)
c       MAXIMUM ERROR BETWEEN INPUT AND OUTPUT = 0.38165595E-06
c    PERFORM A FORWARD TRANSFORM (COMPLEX INPUT, COMPLEX OUTPUT),
c    FOLLOWED BY AN INVERSE TRANSFORM (COMPLEX INPUT, COMPLEX OUTPUT)
c       MAXIMUM ERROR BETWEEN INPUT AND OUTPUT = 0.29802322E-06
c  ILLUSTRATING THE USE OF THE DOUBLE PRECISION COMPLEX 3-D FFT USING
c  RANDOM INPUT DATA:
c    PERFORM A FORWARD TRANSFORM (COMPLEX INPUT, COMPLEX OUTPUT),
c    FOLLOWED BY AN INVERSE TRANSFORM (COMPLEX INPUT, COMPLEX OUTPUT)
c       MAXIMUM ERROR BETWEEN INPUT AND OUTPUT = 0.35544480E-15
c    PERFORM A FORWARD TRANSFORM (REAL INPUT, REAL OUTPUT),
c    FOLLOWED BY AN INVERSE TRANSFORM (REAL INPUT, REAL OUTPUT)
c       MAXIMUM ERROR BETWEEN INPUT AND OUTPUT = 0.33306691E-15
c  ILLUSTRATING THE USE OF THE SINGLE PRECISION EXTENDED PERIODIC AND
c  EXTENDED NON-PERIODIC CONVOLUTION USING RANDOM INPUT DATA:
c    PERFORM A NON-PERIODIC CONVOLUTION FOLLOWED BY A PERIODIC CONVOLUTION
c       MAXIMUM ERROR BETWEEN THE PERIODIC AND NON PERIODIC CONVOLUTION =
c        0.57220459E-05
c
c
         INCLUDE 'CXMLDEF.FOR'

C
         REAL*4 FLOAT_TOLERANCE
         PARAMETER (FLOAT_TOLERANCE = 1.0E-4)
C
         REAL*8 DFLOAT_TOLERANCE
         PARAMETER (DFLOAT_TOLERANCE = 1.0E-11)
C
         INTEGER*4 MAX_SIZE
         PARAMETER (MAX_SIZE = 1000)
C
         REAL*8 A2(MAX_SIZE), B2(MAX_SIZE), C2(MAX_SIZE), D2(MAX_SIZE)
C
         REAL*4 A(MAX_SIZE), B(MAX_SIZE), C(MAX_SIZE), D(MAX_SIZE)
C
         COMPLEX*8 AC(MAX_SIZE), CC(MAX_SIZE)
	   COMPLEX*16 AC2(MAX_SIZE), CC2(MAX_SIZE)
C
         INTEGER*4 STATUS, I, TEMP_ARRAY(MAX_SIZE), IOUNIT
C
C
C ***** SETUP OUTPUT FILE
C
         IOUNIT = 7
         OPEN (UNIT=IOUNIT, FILE='OUTPUT.DATA', STATUS='UNKNOWN')
         REWIND IOUNIT
C
         WRITE (IOUNIT, 101)
C
         STATUS = 0
C
C ***** ILLUSTRATE USE OF SINGLE PRECISION FFT
C
         IF (STATUS .EQ. 0) THEN
             CALL CHECK_SFFT (A, B, MAX_SIZE, FLOAT_TOLERANCE, STATUS,
     $                        IOUNIT)
         END IF
C
C ***** ILLUSTRATE USE OF DOUBLE PRECISION FFT
C
         IF (STATUS .EQ. 0) THEN
             CALL CHECK_DFFT (A2, B2, MAX_SIZE/2, DFLOAT_TOLERANCE,
     $                        STATUS, IOUNIT)
         END IF
C
C ***** ILLUSTRATE USE OF 3-D FFT (COMPLEX SINGLE PRECISION)
C
         IF (STATUS .EQ. 0) THEN
             CALL CHECK_CFFT_3D (A, B, C, D, AC, CC, MAX_SIZE, 
     $                           FLOAT_TOLERANCE, STATUS, IOUNIT)
         END IF
C
C ***** ILLUSTRATE USE OF 3-D FFT (COMPLEX DOUBLE PRECISION)
C
         IF (STATUS .EQ. 0) THEN
             CALL CHECK_ZFFT_3D (A2, B2, C2, D2, AC2, CC2, MAX_SIZE/2,
     $                           DFLOAT_TOLERANCE, STATUS, IOUNIT)
         END IF
C
C ***** ILLUSTRATE USE OF CONVOLUTION
C
         IF (STATUS .EQ. 0) THEN
             CALL CHECK_SCONV (A, B, C, D, MAX_SIZE, FLOAT_TOLERANCE,
     $                         STATUS, IOUNIT)
         END IF
C
  101    FORMAT (/,2X, 'ILLUSTRATING THE USE OF SIGNAL PROCESSING',
     $                 ' ROUTINES',/)
C
         END
C
C ----------------------------------------------------------------------
C
         SUBROUTINE CHECK_SFFT (A, B, MAX_SIZE, TOLERANCE, STATUS, 
     $                          IOUNIT)
C
C ***** THIS ROUTINE ILLUSTRATES THE USE OF SINGLE PRECISION FFT.
C
         INCLUDE 'CXMLDEF.FOR'

C
         INTEGER*4 FFT_SIZE
         PARAMETER (FFT_SIZE = 64)
C
         INTEGER*4 MAX_SIZE
         REAL*4 A(0:MAX_SIZE-1), B(0:MAX_SIZE-1), TOLERANCE,
     $          MAX_ERROR
C
         INTEGER*4 STATUS, I, ISEED, ASTATUS, BSTATUS, IOUNIT
C
C
         WRITE (IOUNIT, 101)
C
         IF (FFT_SIZE+2 .GT. MAX_SIZE) THEN
            STATUS = 2
            RETURN
         ENDIF
C
         ISEED = 1000000001
         STATUS = 0
C
C ***** INITIALISE DATA
C
         DO I = 0, FFT_SIZE-1, 1
            A(I) = RAN(ISEED)
         ENDDO
C
C ***** PERFORM A FORWARD FFT USING REAL INPUT FORMAT AND COMPLEX OUTPUT
C       FORMAT, FOLLOWED BY AN INVERSE TRANSFORM USING COMPLEX INPUT AND 
C       REAL OUTPUT. CHECK THE ERROR BETWEEN THE ORIGINAL INPUT AND THE 
C       FINAL OUTPUT AFTER THE FOWARD AND INVERSE TRANSFORMS.
C
         IF (STATUS .EQ. 0) THEN
C
            WRITE (IOUNIT, 102)
C
            ASTATUS = SFFT('R', 'C', 'F', A, B, FFT_SIZE, 1)
            BSTATUS = SFFT('C', 'R', 'B', B, B, FFT_SIZE, 1)
C
            MAX_ERROR = FLOAT(0)
            DO I = 0, FFT_SIZE-1, 1
               MAX_ERROR = MAX ( MAX_ERROR, ABS(A(I)-B(I)) )
            ENDDO
C
            WRITE (IOUNIT, 104) MAX_ERROR
C
            IF ( (MAX_ERROR .GT. TOLERANCE) .OR. 
     $           (ASTATUS .NE. DXML_SUCCESS())  .OR. 
     $           (BSTATUS .NE. DXML_SUCCESS()) )  THEN
               STATUS = 1
               WRITE (IOUNIT, 105)
            END IF
C
         ENDIF
C
C ***** PERFORM A FORWARD FFT USING REAL INPUT FORMAT AND REAL OUTPUT
C       FORMAT, FOLLOWED BY AN INVERSE TRANSFORM USING REAL INPUT AND 
C       REAL OUTPUT. CHECK THE ERROR BETWEEN THE ORIGINAL INPUT AND THE 
C       FINAL OUTPUT AFTER THE FOWARD AND INVERSE TRANSFORMS.
C
         IF (STATUS .EQ. 0) THEN
C
            WRITE (IOUNIT, 103)
C
            ASTATUS = SFFT('R', 'R', 'F', A, B, FFT_SIZE, 1)
            BSTATUS = SFFT('R', 'R', 'B', B, B, FFT_SIZE, 1)
C
            MAX_ERROR = FLOAT(0)
            DO I = 0, FFT_SIZE-1, 1
               MAX_ERROR = MAX ( MAX_ERROR, ABS(A(I)-B(I)) )
            ENDDO
C
            WRITE (IOUNIT, 104) MAX_ERROR
C
            IF ( (MAX_ERROR .GT. TOLERANCE) .OR. 
     $           (ASTATUS .NE. DXML_SUCCESS()) .OR.
     $           (BSTATUS .NE. DXML_SUCCESS()) ) THEN
               STATUS = 1
               WRITE (IOUNIT, 105)
            END IF
C
         ENDIF
C
  101    FORMAT (///,2X, 'ILLUSTRATING THE USE OF THE SINGLE PRECISION',
     $                 ' FFT USING RANDOM INPUT DATA:',/)
  102    FORMAT (/,4X, 'PERFORM A FORWARD TRANSFORM',
     $                 ' (REAL INPUT, COMPLEX OUTPUT),'
     $           /,4X, 'FOLLOWED BY AN INVERSE TRANSFORM',
     $                 ' (COMPLEX INPUT, REAL OUTPUT)')
  103    FORMAT (/,4X, 'PERFORM A FORWARD TRANSFORM',
     $                 ' (REAL INPUT, REAL OUTPUT),'
     $           /,4X, 'FOLLOWED BY AN INVERSE TRANSFORM',
     $                 ' (REAL INPUT, REAL OUTPUT)')
  104    FORMAT (/,7X, 'MAXIMUM ERROR BETWEEN INPUT AND OUTPUT =',
     $                  E15.8)
  105    FORMAT (/,7X, 'ERROR IN SIGNAL PROCESSING ROUTINE')
C
         RETURN
         END
C
C ----------------------------------------------------------------------
C
         SUBROUTINE CHECK_DFFT (A, B, MAX_SIZE, TOLERANCE, STATUS,
     $                          IOUNIT)
C
C ***** THIS ROUTINE ILLUSTRATES THE USE OF DOUBLE PRECISION FFT.
C
         INCLUDE 'CXMLDEF.FOR'

C
         INTEGER*4 FFT_SIZE
         PARAMETER (FFT_SIZE = 64)
C
         INTEGER*4 MAX_SIZE
         REAL*8 A(0:MAX_SIZE-1), B(0:MAX_SIZE-1), TOLERANCE, 
     $          MAX_ERROR
C
         INTEGER*4 STATUS, I, ISEED, ASTATUS, BSTATUS, IOUNIT
C
C
         WRITE (IOUNIT, 101)
C
         ISEED = 1000000001
         STATUS = 0
C
C ***** INITIALISE DATA
C
         DO I = 0, FFT_SIZE-1, 1
            A(I) = RAN(ISEED)
         ENDDO
C
C ***** PERFORM A FORWARD FFT USING REAL INPUT FORMAT AND COMPLEX OUTPUT
C       FORMAT, FOLLOWED BY AN INVERSE TRANSFORM USING COMPLEX INPUT AND 
C       REAL OUTPUT. CHECK THE ERROR BETWEEN THE ORIGINAL INPUT AND THE 
C       FINAL OUTPUT AFTER THE FOWARD AND INVERSE TRANSFORMS.
C
C
         IF (STATUS .EQ. 0) THEN
C
            WRITE (IOUNIT, 102)
C
            ASTATUS = DFFT('R', 'C', 'F', A, B, FFT_SIZE, 1)
            BSTATUS = DFFT('C', 'R', 'B', B, B, FFT_SIZE, 1)
C
            MAX_ERROR = DFLOAT(0)
            DO I = 0, FFT_SIZE-1, 1
               MAX_ERROR = MAX ( MAX_ERROR, ABS(A(I)-B(I)) )
            ENDDO
C
            WRITE (IOUNIT, 104) MAX_ERROR
C
            IF ( (MAX_ERROR .GT. TOLERANCE) .OR. 
     $           (ASTATUS .NE. DXML_SUCCESS())  .OR.
     $           (BSTATUS .NE. DXML_SUCCESS()) ) THEN
                STATUS = 1
                WRITE (IOUNIT, 105)
            END IF
C
         ENDIF
C
C ***** PERFORM A FORWARD FFT USING REAL INPUT FORMAT AND REAL OUTPUT
C       FORMAT, FOLLOWED BY AN INVERSE TRANSFORM USING REAL INPUT AND 
C       REAL OUTPUT. CHECK THE ERROR BETWEEN THE ORIGINAL INPUT AND THE 
C       FINAL OUTPUT AFTER THE FOWARD AND INVERSE TRANSFORMS.
C
         IF (STATUS .EQ. 0) THEN
C
            WRITE (IOUNIT, 102)
C
            ASTATUS = DFFT('R', 'R', 'F', A, B, FFT_SIZE, 1)
            BSTATUS = DFFT('R', 'R', 'B', B, B, FFT_SIZE, 1)
C
            MAX_ERROR = DFLOAT(0)
            DO I = 0, FFT_SIZE-1, 1
               MAX_ERROR = MAX ( MAX_ERROR, ABS(A(I)-B(I)) )
            ENDDO
C
            WRITE (IOUNIT, 104) MAX_ERROR
C
            IF ( (MAX_ERROR .GT. TOLERANCE) .OR. 
     $           (ASTATUS .NE. DXML_SUCCESS()) .OR.
     $           (BSTATUS .NE. DXML_SUCCESS()) ) THEN
                STATUS = 1
                WRITE (IOUNIT, 105)
            END IF
C
         ENDIF
C
  101    FORMAT (///,2X, 'ILLUSTRATING THE USE OF THE DOUBLE PRECISION',
     $                 ' FFT USING RANDOM INPUT DATA:',/)
  102    FORMAT (/,4X, 'PERFORM A FORWARD TRANSFORM',
     $                 ' (REAL INPUT, COMPLEX OUTPUT),'
     $           /,4X, 'FOLLOWED BY AN INVERSE TRANSFORM',
     $                 ' (COMPLEX INPUT, REAL OUTPUT)')
  103    FORMAT (/,4X, 'PERFORM A FORWARD TRANSFORM',
     $                 ' (REAL INPUT, REAL OUTPUT),'
     $           /,4X, 'FOLLOWED BY AN INVERSE TRANSFORM',
     $                 ' (REAL INPUT, REAL OUTPUT)')
  104    FORMAT (/,7X, 'MAXIMUM ERROR BETWEEN INPUT AND OUTPUT =',
     $                  E15.8)
  105    FORMAT (/,7X, 'ERROR IN SIGNAL PROCESSING ROUTINE')
C
         RETURN
         END
C
C ----------------------------------------------------------------------
C
         SUBROUTINE CHECK_CFFT_3D (A, B, C, D, AC, BC, MAX_SIZE, 
     $                             TOLERANCE, STATUS, IOUNIT)
C
C ***** THIS ROUTINE ILLUSTRATES THE USE OF COMPLEX, SINGLE PRECISION
C       3-D FFT.
C
         INCLUDE 'CXMLDEF.FOR'

C
         INTEGER*4 FFT_SIZE_X, FFT_SIZE_Y, FFT_SIZE_Z
         PARAMETER (FFT_SIZE_X = 6, FFT_SIZE_Y = 4, FFT_SIZE_Z = 8)
C
         COMPLEX*8 AC(0:FFT_SIZE_X-1, 0:FFT_SIZE_Y-1, 0:FFT_SIZE_Z-1),
     $             BC(0:FFT_SIZE_X-1, 0:FFT_SIZE_Y-1, 0:FFT_SIZE_Z-1)
C
         REAL*4 A(0:FFT_SIZE_X-1, 0:FFT_SIZE_Y-1, 0:FFT_SIZE_Z-1),
     $          B(0:FFT_SIZE_X-1, 0:FFT_SIZE_Y-1, 0:FFT_SIZE_Z-1),
     $          C(0:FFT_SIZE_X-1, 0:FFT_SIZE_Y-1, 0:FFT_SIZE_Z-1),
     $          D(0:FFT_SIZE_X-1, 0:FFT_SIZE_Y-1, 0:FFT_SIZE_Z-1)
C
         REAL*4 MAX_ERROR, TOLERANCE
C
         INTEGER*4 MAX_SIZE, STATUS,
     $             I, J, K, ISEED, ASTATUS, BSTATUS, IOUNIT
C
C
         WRITE (IOUNIT, 101)
C
         IF ((FFT_SIZE_X*FFT_SIZE_Y*FFT_SIZE_Z) .GT. MAX_SIZE) THEN
            STATUS = 2
            RETURN
         ENDIF
C
         ISEED = 1000000001
         STATUS = 0
C
C ***** PERFORM A FORWARD FFT USING COMPLEX INPUT FORMAT AND COMPLEX OUTPUT
C       FORMAT, FOLLOWED BY AN INVERSE TRANSFORM USING COMPLEX INPUT AND 
C       COMPLEX OUTPUT. CHECK THE ERROR BETWEEN THE ORIGINAL INPUT AND THE 
C       FINAL OUTPUT AFTER THE FOWARD AND INVERSE TRANSFORMS.
C
C
         IF (STATUS .EQ. 0) THEN
C
            WRITE (IOUNIT, 102)
C
C ***** INITIALISE DATA
C
            DO K = 0, FFT_SIZE_Z-1, 1
               DO J = 0, FFT_SIZE_Y-1, 1
                  DO I = 0, FFT_SIZE_X-1, 1
                     AC(I,J,K) = CMPLX ( RAN(ISEED), RAN(ISEED) )
                  ENDDO
               ENDDO
            ENDDO
C
            ASTATUS = CFFT_3D('C', 'C', 'F', AC, BC, 
     $                         FFT_SIZE_X, FFT_SIZE_Y, FFT_SIZE_Z,
     $                         FFT_SIZE_X, FFT_SIZE_Y, 1, 1, 1)
C
            BSTATUS = CFFT_3D('C', 'C', 'B', BC, BC,
     $                         FFT_SIZE_X, FFT_SIZE_Y, FFT_SIZE_Z,
     $                         FFT_SIZE_X, FFT_SIZE_Y, 1, 1, 1)
C
            MAX_ERROR = FLOAT(0)
            DO K = 0, FFT_SIZE_Z-1, 1
               DO J = 0, FFT_SIZE_Y-1, 1
                  DO I = 0, FFT_SIZE_X-1, 1
                     MAX_ERROR = MAX ( MAX_ERROR, 
     $                                 ABS(AC(I,J,K)-BC(I,J,K)) )
                  ENDDO
               ENDDO
            ENDDO
C
            WRITE (IOUNIT, 104) MAX_ERROR
C
            IF ( (MAX_ERROR .GT. TOLERANCE) .OR. 
     $           (ASTATUS .NE. DXML_SUCCESS()) .OR.
     $           (BSTATUS .NE. DXML_SUCCESS()) ) THEN
               STATUS = 1
               WRITE (IOUNIT, 105)
            END IF
C
         ENDIF
C
C ***** PERFORM A FORWARD FFT USING REAL INPUT FORMAT AND REAL OUTPUT
C       FORMAT, FOLLOWED BY AN INVERSE TRANSFORM USING REAL INPUT AND 
C       REAL OUTPUT. CHECK THE ERROR BETWEEN THE ORIGINAL INPUT AND THE 
C       FINAL OUTPUT AFTER THE FOWARD AND INVERSE TRANSFORMS.
C
C
         IF (STATUS .EQ. 0) THEN
C
            WRITE (IOUNIT, 102)
C
C ***** INITIALISE DATA
C
            DO K = 0, FFT_SIZE_Z-1, 1
               DO J = 0, FFT_SIZE_Y-1, 1
                  DO I = 0, FFT_SIZE_X-1, 1
                     A(I,J,K) = RAN(ISEED)
                     B(I,J,K) = RAN(ISEED)
                  ENDDO
               ENDDO
            ENDDO
C
            ASTATUS = CFFT_3D('R', 'R', 'F', A, B, C, D,
     $                        FFT_SIZE_X, FFT_SIZE_Y, FFT_SIZE_Z,
     $                        FFT_SIZE_X, FFT_SIZE_Y, 1, 1, 1)
C
            BSTATUS = CFFT_3D('R', 'R', 'B', C, D, C, D,
     $                        FFT_SIZE_X, FFT_SIZE_Y, FFT_SIZE_Z,
     $                        FFT_SIZE_X, FFT_SIZE_Y, 1, 1, 1)
C
            MAX_ERROR = FLOAT(0)
            DO K = 0, FFT_SIZE_Z-1, 1
               DO J = 0, FFT_SIZE_Y-1, 1
                  DO I = 0, FFT_SIZE_X-1, 1
                     MAX_ERROR = MAX ( MAX_ERROR,
     $                                 ABS(A(I,J,K)-C(I,J,K)),
     $                                 ABS(B(I,J,K)-D(I,J,K)) )
                  ENDDO
               ENDDO
            ENDDO
C
            WRITE (IOUNIT, 104) MAX_ERROR
C
            IF ( (MAX_ERROR .GT. TOLERANCE) .OR. 
     $           (ASTATUS .NE. DXML_SUCCESS())  .OR.
     $           (BSTATUS .NE. DXML_SUCCESS()) ) THEN
               STATUS = 1
               WRITE (IOUNIT, 105)
            END IF
C
         ENDIF
C
C
  101    FORMAT (///,2X, 'ILLUSTRATING THE USE OF THE SINGLE PRECISION',
     $                   ' COMPLEX 3-D FFT USING',
     $             /,2X, 'RANDOM INPUT DATA:',/)
  102    FORMAT (/,4X, 'PERFORM A FORWARD TRANSFORM',
     $                 ' (COMPLEX INPUT, COMPLEX OUTPUT),'
     $           /,4X, 'FOLLOWED BY AN INVERSE TRANSFORM',
     $                 ' (COMPLEX INPUT, COMPLEX OUTPUT)')
  103    FORMAT (/,4X, 'PERFORM A FORWARD TRANSFORM',
     $                 ' (REAL INPUT, REAL OUTPUT),'
     $           /,4X, 'FOLLOWED BY AN INVERSE TRANSFORM',
     $                 ' (REAL INPUT, REAL OUTPUT)')
  104    FORMAT (/,7X, 'MAXIMUM ERROR BETWEEN INPUT AND OUTPUT =',
     $                  E15.8)
  105    FORMAT (/,7X, 'ERROR IN SIGNAL PROCESSING ROUTINE')
C
         RETURN
         END
C
C ----------------------------------------------------------------------
C
         SUBROUTINE CHECK_ZFFT_3D (A, B, C, D, AC, BC, MAX_SIZE, 
     $                             TOLERANCE, STATUS, IOUNIT)
C
C ***** THIS ROUTINE ILLUSTRATES THE USE OF DOUBLE PRECISION, COMPLEX
C       3D FFT.
C
         INCLUDE 'CXMLDEF.FOR'

         INTEGER*4 FFT_SIZE_X, FFT_SIZE_Y, FFT_SIZE_Z
         PARAMETER (FFT_SIZE_X = 6, FFT_SIZE_Y = 4, FFT_SIZE_Z = 8)
C
         COMPLEX*16 AC(0:FFT_SIZE_X-1, 0:FFT_SIZE_Y-1, 0:FFT_SIZE_Z-1),
     $              BC(0:FFT_SIZE_X-1, 0:FFT_SIZE_Y-1, 0:FFT_SIZE_Z-1)
C
         REAL*8 A(0:FFT_SIZE_X-1, 0:FFT_SIZE_Y-1, 0:FFT_SIZE_Z-1),
     $          B(0:FFT_SIZE_X-1, 0:FFT_SIZE_Y-1, 0:FFT_SIZE_Z-1),
     $          C(0:FFT_SIZE_X-1, 0:FFT_SIZE_Y-1, 0:FFT_SIZE_Z-1),
     $          D(0:FFT_SIZE_X-1, 0:FFT_SIZE_Y-1, 0:FFT_SIZE_Z-1)
C  
         REAL*8 MAX_ERROR, TOLERANCE
C
         INTEGER*4 I, J, K, ISEED, ASTATUS, BSTATUS,
     $             MAX_SIZE, STATUS, IOUNIT
C
C
         WRITE (IOUNIT, 101)
C
         IF ((FFT_SIZE_X*FFT_SIZE_Y*FFT_SIZE_Z) .GT. MAX_SIZE) THEN
            STATUS = 2
            RETURN
         ENDIF
C
         ISEED = 1000000001
         STATUS = 0
C
C ***** PERFORM A FORWARD FFT USING COMPLEX INPUT FORMAT AND COMPLEX OUTPUT
C       FORMAT, FOLLOWED BY AN INVERSE TRANSFORM USING COMPLEX INPUT AND 
C       COMPLEX OUTPUT. CHECK THE ERROR BETWEEN THE ORIGINAL INPUT AND THE 
C       FINAL OUTPUT AFTER THE FOWARD AND INVERSE TRANSFORMS.
C
C
         IF (STATUS .EQ. 0) THEN
C
            WRITE (IOUNIT, 102)
C
C ***** INITIALISE DATA
C
            DO K = 0, FFT_SIZE_Z-1, 1
               DO J = 0, FFT_SIZE_Y-1, 1
                  DO I = 0, FFT_SIZE_X-1, 1
                     AC(I,J,K) = CMPLX ( RAN(ISEED), RAN(ISEED) )
                  ENDDO
               ENDDO
            ENDDO
C
            ASTATUS = ZFFT_3D ('C', 'C', 'F', AC, BC, 
     $                          FFT_SIZE_X, FFT_SIZE_Y, FFT_SIZE_Z,
     $                          FFT_SIZE_X, FFT_SIZE_Y, 1, 1, 1)
C
            BSTATUS = ZFFT_3D ('C', 'C', 'B', BC, BC,
     $                         FFT_SIZE_X, FFT_SIZE_Y, FFT_SIZE_Z,
     $                         FFT_SIZE_X, FFT_SIZE_Y, 1, 1, 1)
C
            MAX_ERROR = DFLOAT(0)
            DO K = 0, FFT_SIZE_Z-1, 1
               DO J = 0, FFT_SIZE_Y-1, 1
                  DO I = 0, FFT_SIZE_X-1, 1
                     MAX_ERROR = MAX ( MAX_ERROR, 
     $                                 ABS(AC(I,J,K)-BC(I,J,K)) )
                  ENDDO
               ENDDO
            ENDDO
C
            WRITE (IOUNIT, 104) MAX_ERROR
C
            IF ( (MAX_ERROR .GT. TOLERANCE) .OR. 
     $           (ASTATUS .NE. DXML_SUCCESS()) .OR.
     $           (BSTATUS .NE. DXML_SUCCESS()) ) THEN
                STATUS = 1
                WRITE (IOUNIT, 105)
            END IF
C
         ENDIF
C
C ***** PERFORM A FORWARD FFT USING REAL INPUT FORMAT AND REAL OUTPUT
C       FORMAT, FOLLOWED BY AN INVERSE TRANSFORM USING REAL INPUT AND 
C       REAL OUTPUT. CHECK THE ERROR BETWEEN THE ORIGINAL INPUT AND THE 
C       FINAL OUTPUT AFTER THE FOWARD AND INVERSE TRANSFORMS.
C
C
         IF (STATUS .EQ. 0) THEN
C
            WRITE (IOUNIT, 103)
C
C ***** INITIALISE DATA
C
            DO K = 0, FFT_SIZE_Z-1, 1
               DO J = 0, FFT_SIZE_Y-1, 1
                  DO I = 0, FFT_SIZE_X-1, 1
                     A(I,J,K) = RAN(ISEED)
                     B(I,J,K) = RAN(ISEED)
                  ENDDO
               ENDDO
            ENDDO
C
            ASTATUS = ZFFT_3D ('R', 'R', 'F', A, B, C, D,
     $                         FFT_SIZE_X, FFT_SIZE_Y, FFT_SIZE_Z,
     $                         FFT_SIZE_X, FFT_SIZE_Y, 1, 1, 1)
C
            BSTATUS = ZFFT_3D ('R', 'R', 'B', C, D, C, D, 
     $                         FFT_SIZE_X, FFT_SIZE_Y, FFT_SIZE_Z,
     $                         FFT_SIZE_X, FFT_SIZE_Y, 1, 1, 1)
C
            MAX_ERROR = DFLOAT(0)
            DO K = 0, FFT_SIZE_Z-1, 1
               DO J = 0, FFT_SIZE_Y-1, 1
                  DO I = 0, FFT_SIZE_X-1, 1
                     MAX_ERROR = MAX ( MAX_ERROR,
     $                                 ABS(A(I,J,K)-C(I,J,K)),
     $                                 ABS(B(I,J,K)-D(I,J,K)) )
                  ENDDO
               ENDDO
            ENDDO
C
            WRITE (IOUNIT, 104) MAX_ERROR
C
            IF ( (MAX_ERROR .GT. TOLERANCE) .OR. 
     $           (ASTATUS .NE. DXML_SUCCESS()) .OR.
     $           (BSTATUS .NE. DXML_SUCCESS()) ) THEN
               STATUS = 1
               WRITE (IOUNIT, 105)
            END IF
C
         ENDIF
C
  101    FORMAT (///,2X, 'ILLUSTRATING THE USE OF THE DOUBLE PRECISION',
     $                 ' COMPLEX 3-D FFT USING',
     $             /,2X, 'RANDOM INPUT DATA:',/)
  102    FORMAT (/,4X, 'PERFORM A FORWARD TRANSFORM',
     $                 ' (COMPLEX INPUT, COMPLEX OUTPUT),'
     $           /,4X, 'FOLLOWED BY AN INVERSE TRANSFORM',
     $                 ' (COMPLEX INPUT, COMPLEX OUTPUT)')
  103    FORMAT (/,4X, 'PERFORM A FORWARD TRANSFORM',
     $                 ' (REAL INPUT, REAL OUTPUT),'
     $           /,4X, 'FOLLOWED BY AN INVERSE TRANSFORM',
     $                 ' (REAL INPUT, REAL OUTPUT)')
  104    FORMAT (/,7X, 'MAXIMUM ERROR BETWEEN INPUT AND OUTPUT =',
     $                  E15.8)
  105    FORMAT (/,7X, 'ERROR IN SIGNAL PROCESSING ROUTINE')
C
         RETURN
         END
C
C ----------------------------------------------------------------------
C
         SUBROUTINE CHECK_SCONV (A, B, C, D, MAX_SIZE, 
     $                           TOLERANCE, STATUS, IOUNIT)
C
C ***** ILLUSTRATE THE USE OF EXTENDED NON-PERIODIC CONVOLUTION AND
C       EXTENDED PERIODIC CONVOLUTION IN SINGLE PRECISION.
C
         INCLUDE 'CXMLDEF.FOR'

C
         INTEGER*4 DATA_SIZE
         PARAMETER (DATA_SIZE = 64)
C
         INTEGER*4 MAX_SIZE
         REAL*4 A(0:MAX_SIZE-1), B(0:MAX_SIZE-1),
     $          C(0:MAX_SIZE-1), D(0:MAX_SIZE-1)
C
         REAL*4 MAX_ERROR, TOLERANCE
C
         INTEGER*4 STATUS, I, ISEED, ASTATUS, BSTATUS, IOUNIT
C
C
         WRITE (IOUNIT, 101)
C
         IF (2*DATA_SIZE .GT. MAX_SIZE) THEN
            STATUS = 2
            RETURN
         ENDIF
C
         ISEED = 1000000001
         STATUS = 0
C
C ***** INITIALISE DATA
C
         DO I = 0, DATA_SIZE-1, 1
            A(I) = RAN(ISEED)
            B(I) = RAN(ISEED)
         ENDDO
C
C ***** CALL THE NONPERIODIC CONVOLUTION ROUTINE FOLLOWED BY THE
C       PERIODIC CONVOLUTION ROUTINE.
C
         WRITE (IOUNIT, 102)
C
         ASTATUS = SCONV_NONPERIODIC_EXT (A, 1, B, 1, C, 1,
     $                                    DATA_SIZE, DATA_SIZE, 0,
     $                                    2*DATA_SIZE-2,
     $                                    .FALSE., .FALSE., 0., 0)
C
         BSTATUS = SCONV_PERIODIC_EXT (A, 1, B, 1, D, 1,
     $                                 DATA_SIZE, 0, DATA_SIZE-1,
     $                                 .FALSE., .FALSE., 0., 0)
C
         MAX_ERROR = FLOAT(0)
         DO I = 0, DATA_SIZE-2, 1
            MAX_ERROR = MAX (MAX_ERROR, 
     $                       ABS((C(I)+C(DATA_SIZE+I))-D(I)) )
         ENDDO
C
         WRITE (IOUNIT, 104) MAX_ERROR
C
         IF ( (MAX_ERROR .GT. TOLERANCE) .OR. 
     $        (ASTATUS .NE. DXML_SUCCESS()) .OR. 
     $        (BSTATUS .NE. DXML_SUCCESS()) ) THEN
             STATUS = 1
             WRITE (IOUNIT, 105)
         END IF
C
  101    FORMAT (///,2X, 'ILLUSTRATING THE USE OF THE SINGLE ',
     $                   'PRECISION EXTENDED PERIODIC AND',
     $             /,2X, 'EXTENDED NON-PERIODIC CONVOLUTION USING',
     $                   ' RANDOM INPUT DATA:',/)
  102    FORMAT (/,4X, 'PERFORM A NON-PERIODIC CONVOLUTION',
     $                 ' FOLLOWED BY A PERIODIC CONVOLUTION')
  104    FORMAT (/,7X, 'MAXIMUM ERROR BETWEEN THE PERIODIC AND',
     $                 ' NON PERIODIC CONVOLUTION =',/,8X,E15.8)
  105    FORMAT (/,7X, 'ERROR IN SIGNAL PROCESSING ROUTINE')
C
         RETURN
         END

