path: root/src/hyper/pages/nagc.ht

\begin{page}{manpageXXc02}{NAG On-line Documentation: c02}
\beginscroll
\begin{verbatim}



     C02(3NAG)         Foundation Library (12/10/92)         C02(3NAG)



          C02 -- Zeros of Polynomials                   Introduction -- C02
                                    Chapter C02
                               Zeros of Polynomials

          1. Scope of the Chapter

          This chapter is concerned with computing the zeros of a
          polynomial with real or complex coefficients.

          2. Background to the Problems

          Let f(z) be a polynomial of degree n with complex coefficients
          a :
           i

                            n    n-1    n-2
                   f(z)==a z +a z   +a z   +...+a   z+a ,  a /=0.
                          0    1      2          n-1   n    0

          A complex number z  is called a zero of f(z) (or equivalently a
                            1
          root of the equation f(z)=0), if:

                                      f(z )=0.
                                         1

          If z  is a zero, then f(z) can be divided by a factor (z-z ):
              1                                                     1

                                f(z)=(z-z )f (z)                        (1)
                                         1  1

          where f (z) is a polynomial of degree n-1. By the Fundamental
                 1
          Theorem of Algebra, a polynomial f(z) always has a zero, and so
          the process of dividing out factors (z-z ) can be continued until
                                                  i
          we have a complete factorization of f(z)

                           f(z)==a (z-z )(z-z )...(z-z ).
                                  0    1     2        n

          Here the complex numbers z ,z ,...,z  are the zeros of f(z); they
                                    1  2      n
          may not all be distinct, so it is sometimes more convenient to
          write:

                                   m       m          m
                                    1       2          k
                     f(z)==a (z-z )  (z-z )  ...(z-z )  ,  k<=n,
                            0    1       2          k

          with distinct zeros z ,z ,...,z  and multiplicities m >=1. If
                               1  2      k                     i
          m =1, z  is called a single zero, if m >1, z  is called a
           i     i                              i     i
          multiple or repeated zero; a multiple zero is also a zero of the
          derivative of f(z).

          If the coefficients of f(z) are all real, then the zeros of f(z)
          are either real or else occur as pairs of conjugate complex
          numbers x+iy and x-iy. A pair of complex conjugate zeros are the
                                                 2
          zeros of a quadratic factor of f(z), (z +rz+s), with real
          coefficients r and s.

          Mathematicians are accustomed to thinking of polynomials as
          pleasantly simple functions to work with. However the problem of
          numerically computing the zeros of an arbitrary polynomial is far
          from simple. A great variety of algorithms have been proposed, of
          which a number have been widely used in practice; for a fairly
          comprehensive survey, see Householder [1]. All general algorithms
          are iterative. Most converge to one zero at a time; the
          corresponding factor can then be divided out as in equation (1)
          above - this process is called deflation or, loosely, dividing
          out the zero - and the algorithm can be applied again to the
          polynomial f (z). A pair of complex conjugate zeros can be
                      1
          divided out together - this corresponds to dividing f(z) by a
          quadratic factor.

          Whatever the theoretical basis of the algorithm, a number of
          practical problems arise: for a thorough discussion of some of
          them see Peters and Wilkinson [2] and Wilkinson [3]. The most
          elementary point is that, even if z  is mathematically an exact
                                             1
          zero of f(z), because of the fundamental limitations of computer
          arithmetic the computed value of f(z ) will not necessarily be
                                              1
          exactly 0.0. In practice there is usually a small region of
          values of z about the exact zero at which the computed value of
          f(z) becomes swamped by rounding errors. Moreover in many
          algorithms this inaccuracy in the computed value of f(z) results
          in a similar inaccuracy in the computed step from one iterate to
          the next. This limits the precision with which any zero can be
          computed. Deflation is another potential cause of trouble, since,
          in the notation of equation (1), the computed coefficients of
          f (z) will not be completely accurate, especially if z  is not an
           1                                                    1
          exact zero of f(z); so the zeros of the computed f (z) will
                                                            1
          deviate from the zeros of f(z).

          A zero is called ill-conditioned if it is sensitive to small
          changes in the coefficients of the polynomial. An ill-conditioned
          zero is likewise sensitive to the computational inaccuracies just
          mentioned. Conversely a zero is called well-conditioned if it is
          comparatively insensitive to such perturbations. Roughly speaking
          a zero which is well separated from other zeros is well-
          conditioned, while zeros which are close together are ill-
          conditioned, but in talking about 'closeness' the decisive factor
          is not the absolute distance between neighbouring zeros but their
          ratio: if the ratio is close to 1 the zeros are ill-conditioned.
          In particular, multiple zeros are ill-conditioned. A multiple
          zero is usually split into a cluster of zeros by perturbations in
          the polynomial or computational inaccuracies.

          2.1. References

          [1]   Householder A S (1970) The Numerical Treatment of a Single
                Nonlinear Equation. McGraw-Hill.

          [2]   Peters G and Wilkinson J H (1971) Practical Problems Arising
                in the Solution of Polynomial Equations. J. Inst. Maths
                Applics. 8 16--35.

          [3]   Wilkinson J H (1963) Rounding Errors in Algebraic Processes,
                Chapter 2. HMSO.

          3. Recommendations on Choice and Use of Routines

          3.1. Discussion

          Two routines are available: C02AFF for polynomials with complex
          coefficients and C02AGF for polynomials with real coefficients.

          C02AFF and C02AGF both use a variant of Laguerre's Method due to
          BT Smith to calculate each zero until the degree of the deflated
          polynomial is less than 3, whereupon the remaining zeros are
          obtained using the 'standard' closed formulae for a quadratic or
          linear equation.

          The accuracy of the roots will depend on how ill-conditioned they
          are. Peters and Wilkinson [2] describe techniques for estimating
          the errors in the zeros after they have been computed.

          3.2. Index

           Zeros of a complex polynomial                             C02AFF
           Zeros of a real polynomial                                C02AGF


          C02 -- Zeros of Polynomials                       Contents -- C02
          Chapter C02

          Zeros of Polynomials

          C02AFF  All zeros of complex polynomial, modified Laguerre method

          C02AGF  All zeros of real polynomial, modified Laguerre method

\end{verbatim}
\endscroll
\end{page}
\begin{page}{manpageXXc02aff}{NAG On-line Documentation: c02aff}
\beginscroll
\begin{verbatim}



     C02AFF(3NAG)      Foundation Library (12/10/92)      C02AFF(3NAG)



          C02 -- Zeros of Polynomials                                C02AFF
                  C02AFF -- NAG Foundation Library Routine Document

          Note: Before using this routine, please read the Users' Note for
          your implementation to check implementation-dependent details.
          The symbol (*) after a NAG routine name denotes a routine that is
          not included in the Foundation Library.

          1. Purpose

          C02AFF finds all the roots of a complex polynomial equation,
          using a variant of Laguerre's Method.

          2. Specification

                 SUBROUTINE C02AFF (A, N, SCALE, Z, W, IFAIL)
                 INTEGER          N, IFAIL
                 DOUBLE PRECISION A(2,N+1), Z(2,N), W(4*(N+1))
                 LOGICAL          SCALE

          3. Description

          The routine attempts to find all the roots of the nth degree
          complex polynomial equation

                               n    n-1    n-2
                       P(z)=a z +a z   +a z   +...+a   z+a =0.
                             0    1      2          n-1   n

          The roots are located using a modified form of Laguerre's Method,
          originally proposed by Smith [2].

          The method of Laguerre [3] can be described by the iterative
          scheme

                                             -n*P(z )
                                                   k
                          L(z )=z   -z = ----------------,
                             k   k+1  k             
                                         P'(z )+-  /H(z )
                                             k   \/    k

                                           2
          where H(z )=(n-1)*[(n-1)*(P'(z )) -n*P(z )P''(z )], and z  is
                   k                    k         k      k         0
          specified.

          The sign in the denominator is chosen so that the modulus of the
          Laguerre step at z , viz. |L(z )|, is as small as possible. The
                            k           k
          method can be shown to be cubically convergent for isolated roots
          (real or complex) and linearly convergent for multiple roots.
          The routine generates a sequence of iterates z , z , z ,..., such
                                                        1   2   3
          that |P(z   )|<|P(z )| and ensures that z   +L(z   ) 'roughly'
                   k+1       k                     k+1    k+1
          lies inside a circular region of radius |F| about z  known to
                                                             k
          contain a zero of P(z); that is, |L(z   )|<=|F|, where F denotes
                                               k+1
          the Fejer bound (see Marden [1]) at the point z . Following Smith
                                                         k
          [2], F is taken to be min(B,1.445*n*R), where B is an upper bound
          for the magnitude of the smallest zero given by

                                                         1/n
                      B=1.0001*min(\/n*L(z ),|r |,|a /a |   ),
                                          k    1    n  0

          r  is the zero X of smaller magnitude of the quadratic equation
           1

                                           2
                    2(P''(z )/(2*n*(n-1)))X +2(P'(z )/n)X+P(z )=0
                           k                       k         k

          and the Cauchy lower bound R for the smallest zero is computed
          (using Newton's Method) as the positive root of the polynomial
          equation

                         n      n-1      n-2
                    |a |z +|a |z   +|a |z   +...+|a   |z-|a |=0.
                      0      1        2            n-1     n

          Starting from the origin, successive iterates are generated
          according to the rule z   =z +L(z ) for k = 1,2,3,... and L(z )
                                 k+1  k    k                           k
          is 'adjusted' so that |P(z   )|<|P(z )| and |L(z   )|<=|F|. The
                                    k+1       k           k+1
          iterative procedure terminates if P(z   ) is smaller in absolute
                                               k+1
          value than the bound on the rounding error in P(z   ) and the
                                                           k+1
          current iterate z =z    is taken to be a zero of P(z). The
                           p  k+1
                              ~
          deflated polynomial P(z)=P(z)/(z-z ) of degree n-1 is then
                                            p
          formed, and the above procedure is repeated on the deflated
          polynomial until n<3, whereupon the remaining roots are obtained
          via the 'standard' closed formulae for a linear (n = 1) or
          quadratic (n = 2) equation.

          To obtain the roots of a quadratic polynomial, C02AHF(*) can be
          used.

          4. References

          [1]   Marden M (1966) Geometry of Polynomials. Mathematical
                Surveys. 3 Am. Math. Soc., Providence, RI.

          [2]   Smith B T (1967) ZERPOL: A Zero Finding Algorithm for
                Polynomials Using Laguerre's Method. Technical Report.
                Department of Computer Science, University of Toronto,
                Canada.

          [3]   Wilkinson J H (1965) The Algebraic Eigenvalue Problem.
                Clarendon Press.

          5. Parameters

           1:  A(2,N+1) -- DOUBLE PRECISION array                     Input
               On entry: if A is declared with bounds (2,0:N), then A(1,i)
               and A(2,i) must contain the real and imaginary parts of a
                                                                        i
                                          n-i
               (i.e., the coefficient of z   ), for i=0,1,...,n.
               Constraint: A(1,0) /= 0.0 or A(2,0) /= 0.0.

           2:  N -- INTEGER                                           Input
               On entry: the degree of the polynomial, n. Constraint: N >=
               1.

           3:  SCALE -- LOGICAL                                       Input
               On entry: indicates whether or not the polynomial is to be
               scaled. See Section 8 for advice on when it may be
               preferable to set SCALE = .FALSE. and for a description of
               the scaling strategy. Suggested value: SCALE = .TRUE..

           4:  Z(2,N) -- DOUBLE PRECISION array                      Output
               On exit: the real and imaginary parts of the roots are
               stored in Z(1,i) and Z(2,i) respectively, for i=1,2,...,n.

           5:  W(4*(N+1)) -- DOUBLE PRECISION array               Workspace

           6:  IFAIL -- INTEGER                                Input/Output
               On entry: IFAIL must be set to 0, -1 or 1. For users not
               familiar with this parameter (described in the Essential
               Introduction) the recommended value is 0.

               On exit: IFAIL = 0 unless the routine detects an error (see
               Section 6).

          6. Error Indicators and Warnings

          Errors detected by the routine:

          If on entry IFAIL = 0 or -1, explanatory error messages are
          output on the current error message unit (as defined by X04AAF).

          IFAIL= 1
               On entry A(1,0) = 0.0 and A(2,0) = 0.0,

               or       N < 1.

          IFAIL= 2
               The iterative procedure has failed to converge. This error
               is very unlikely to occur. If it does, please contact NAG
               immediately, as some basic assumption for the arithmetic has
               been violated. See also Section 8.

          IFAIL= 3
               Either overflow or underflow prevents the evaluation of P(z)
               near some of its zeros. This error is very unlikely to
               occur. If it does, please contact NAG immediately. See also
               Section 8.

          7. Accuracy

          All roots are evaluated as accurately as possible, but because of
          the inherent nature of the problem complete accuracy cannot be
          guaranteed.

          8. Further Comments

          If SCALE = .TRUE., then a scaling factor for the coefficients is
          chosen as a power of the base B of the machine so that the
                                                                EMAX-P
          largest coefficient in magnitude approaches THRESH = B      .
          Users should note that no scaling is performed if the largest
          coefficient in magnitude exceeds THRESH, even if SCALE = .TRUE..
          (For definition of B, EMAX and P see the Chapter Introduction
          X02.)

          However, with SCALE = .TRUE., overflow may be encountered when
          the input coefficients a ,a ,a ,...,a  vary widely in magnitude,
                                  0  1  2      n
                                                    (4*P)
          particularly on those machines for which B      overflows. In
          such cases, SCALE should be set to .FALSE. and the coefficients
          scaled so that the largest coefficient in magnitude does not
                  (EMAX-2*P)
          exceed B          .

          Even so, the scaling strategy used in C02AFF is sometimes
          insufficient to avoid overflow and/or underflow conditions. In
          such cases, the user is recommended to scale the independent
          variable (z) so that the disparity between the largest and
          smallest coefficient in magnitude is reduced. That is, use the
          routine to locate the zeros of the polynomial d*P(cz) for some
          suitable values of c and d. For example, if the original
                               -100   100 20                   -10
          polynomial was P(z)=2    i+2   z  , then choosing c=2    and
             100                                                     20
          d=2   , for instance, would yield the scaled polynomial i+z  ,
          which is well-behaved relative to overflow and underflow and has
                           10
          zeros which are 2   times those of P(z).

          If the routine fails with IFAIL = 2 or 3, then the real and
          imaginary parts of any roots obtained before the failure occurred
          are stored in Z in the reverse order in which they were found.
          Let n  denote the number of roots found before the failure
               R
          occurred. Then Z(1,n) and Z(2,n) contain the real and imaginary
          parts of the 1st root found, Z(1,n-1) and Z(2,n-1) contain the
          real and imaginary parts of the 2nd root found, ..., Z(1,n ) and
                                                                    R
          Z(2,n ) contain the real and imaginary parts of the n th root
               R                                               R
          found. After the failure has occurred, the remaining 2*(n-n )
                                                                     R
          elements of Z contain a large negative number (equal to
                         

          -1/(X02AMF().\/2)).

          9. Example

                                                 5    4    3    2
          To find the roots of the polynomial a z +a z +a z +a z +a z+a =0,
                                               0    1    2    3    4   5
          where a =(5.0+6.0i), a =(30.0+20.0i), a =-(0.2+6.0i),
                 0              1                2
          a =(50.0+100000.0i), a =-(2.0-40.0i) and a =(10.0+1.0i).
           3                    4                   5

          The example program is not reproduced here. The source code for
          all example programs is distributed with the NAG Foundation
          Library software and should be available on-line.

\end{verbatim}
\endscroll
\end{page}
\begin{page}{manpageXXc02agf}{NAG On-line Documentation: c02agf}
\beginscroll
\begin{verbatim}



     C02AGF(3NAG)      Foundation Library (12/10/92)      C02AGF(3NAG)



          C02 -- Zeros of Polynomials                                C02AGF
                  C02AGF -- NAG Foundation Library Routine Document

          Note: Before using this routine, please read the Users' Note for
          your implementation to check implementation-dependent details.
          The symbol (*) after a NAG routine name denotes a routine that is
          not included in the Foundation Library.

          1. Purpose

          C02AGF finds all the roots of a real polynomial equation, using a
          variant of Laguerre's Method.

          2. Specification

                 SUBROUTINE C02AGF (A, N, SCALE, Z, W, IFAIL)
                 INTEGER          N, IFAIL
                 DOUBLE PRECISION A(N+1), Z(2,N), W(2*(N+1))
                 LOGICAL          SCALE

          3. Description

          The routine attempts to find all the roots of the nth degree real
          polynomial equation

                               n    n-1    n-2
                       P(z)=a z +a z   +a z   +...+a   z+a =0.
                             0    1      2          n-1   n

          The roots are located using a modified form of Laguerre's Method,
          originally proposed by Smith [2].

          The method of Laguerre [3] can be described by the iterative
          scheme

                                             -n*P(z )
                                                   k
                          L(z )=z   -z = ----------------,
                             k   k+1  k             
                                         P'(z )+-  /H(z )
                                             k   \/    k

                                           2
          where H(z )=(n-1)*[(n-1)*(P'(z )) -n*P(z )P''(z )], and z  is
                   k                    k         k      k         0
          specified.

          The sign in the denominator is chosen so that the modulus of the
          Laguerre step at z , viz. |L(z )|, is as small as possible. The
                            k           k
          method can be shown to be cubically convergent for isolated roots
          (real or complex) and linearly convergent for multiple roots.
          The routine generates a sequence of iterates z , z , z ,..., such
                                                        1   2   3
          that |P(z +1)|<|P(z )| and ensures that z   +L(z   ) 'roughly'
                   k         k                     k+1    k+1
          lies inside a circular region of radius |F| about z  known to
                                                             k
          contain a zero of P(z); that is, |L(z   )|<=|F|, where F denotes
                                               k+1
          the Fejer bound (see Marden [1]) at the point z . Following Smith
                                                         k
          [2], F is taken to be min(B,1.445*n*R), where B is an upper bound
          for the magnitude of the smallest zero given by

                                                         1/n
                      B=1.0001*min(\/n*L(z ),|r |,|a /a |   ),
                                          k    1    n  0

          r  is the zero X of smaller magnitude of the quadratic equation
           1

                                           2
                    2(P''(z )/(2*n*(n-1)))X +2(P'(z )/n)X+P(z )=0
                           k                       k         k

          and the Cauchy lower bound R for the smallest zero is computed
          (using Newton's Method) as the positive root of the polynomial
          equation

                         n      n-1      n-2
                    |a |z +|a |z   +|a |z   +...+|a   |z-|a |=0.
                      0      1        2            n-1     n

          Starting from the origin, successive iterates are generated
          according to the rule z   =z +L(z ) for k=1,2,3,... and L(z ) is
                                 k+1  k    k                         k
                                 k+1       k           k+1
          iterative procedure terminates if P(z   ) is smaller in absolute
                                               k+1
          value than the bound on the rounding error in P(z   ) and the
                                                           k+1
          current iterate z =z    is taken to be a zero of P(z) (as is its
                           p  k-1
                    

          conjugate z  if z  is complex). The deflated polynomial
                     p     p
          ~
          P(z)=P(z)/(z-z ) of degree n-1 if z  is real
                        p                    p
           ~                   
          (P(z)=P(z)/((z-z )(z-z )) of degree n-2 if z  is complex) is then
                          p     p                     p
          formed, and the above procedure is repeated on the deflated
          polynomial until n<3, whereupon the remaining roots are obtained
          via the 'standard' closed formulae for a linear (n = 1) or
          quadratic (n = 2) equation.

          To obtain the roots of a quadratic polynomial, C02AJF(*) can be
          used.

          4. References

          [1]   Marden M (1966) Geometry of Polynomials. Mathematical
                Surveys. 3 Am. Math. Soc., Providence, RI.

          [2]   Smith B T (1967) ZERPOL: A Zero Finding Algorithm for
                Polynomials Using Laguerre's Method. Technical Report.
                Department of Computer Science, University of Toronto,
                Canada.

          [3]   Wilkinson J H (1965) The Algebraic Eigenvalue Problem.
                Clarendon Press.

          5. Parameters

           1:  A(N+1) -- DOUBLE PRECISION array                       Input
               On entry: if A is declared with bounds (0:N), then A(i)
                                                          n-i
               must contain a  (i.e., the coefficient of z   ), for
                             i
               i=0,1,...,n. Constraint: A(0) /= 0.0.

           2:  N -- INTEGER                                           Input
               On entry: the degree of the polynomial, n. Constraint: N >=
               1.

           3:  SCALE -- LOGICAL                                       Input
               On entry: indicates whether or not the polynomial is to be
               scaled. See Section 8 for advice on when it may be
               preferable to set SCALE = .FALSE. and for a description of
               the scaling strategy. Suggested value: SCALE = .TRUE..

           4:  Z(2,N) -- DOUBLE PRECISION array                      Output
               On exit: the real and imaginary parts of the roots are
               stored in Z(1,i) and Z(2,i) respectively, for i=1,2,...,n.
               Complex conjugate pairs of roots are stored in consecutive
               pairs of elements of Z; that is, Z(1,i+1) = Z(1,i) and
               Z(2,i+1)=-Z(2,i).

           5:  W(2*(N+1)) -- DOUBLE PRECISION array               Workspace

           6:  IFAIL -- INTEGER                                Input/Output
               On entry: IFAIL must be set to 0, -1 or 1. For users not
               familiar with this parameter (described in the Essential
               Introduction) the recommended value is 0.
               On exit: IFAIL = 0 unless the routine detects an error (see
               Section 6).

          6. Error Indicators and Warnings

          Errors detected by the routine:

          If on entry IFAIL = 0 or -1, explanatory error messages are
          output on the current error message unit (as defined by X04AAF).

          IFAIL= 1
               On entry A(0) = 0.0,

               or       N < 1.

          IFAIL= 2
               The iterative procedure has failed to converge. This error
               is very unlikely to occur. If it does, please contact NAG
               immediately, as some basic assumption for the arithmetic has
               been violated. See also Section 8.

          IFAIL= 3
               Either overflow or underflow prevents the evaluation of P(z)
               near some of its zeros. This error is very unlikely to
               occur. If it does, please contact NAG immediately. See also
               Section 8.

          7. Accuracy

          All roots are evaluated as accurately as possible, but because of
          the inherent nature of the problem complete accuracy cannot be
          guaranteed.

          8. Further Comments

          If SCALE = .TRUE., then a scaling factor for the coefficients is
          chosen as a power of the base B of the machine so that the
                                                                EMAX-P
          largest coefficient in magnitude approaches THRESH = B      .
          Users should note that no scaling is performed if the largest
          coefficient in magnitude exceeds THRESH, even if SCALE = .TRUE..
          (For definition of B, EMAX and P see the Chapter Introduction
          X02.)

          However, with SCALE = .TRUE., overflow may be encountered when
          the input coefficients a ,a ,a ,...,a  vary widely in magnitude,
                                  0  1  2      n
                                                    (4*P)
          particularly on those machines for which B      overflows. In
          such cases, SCALE should be set to .FALSE. and the coefficients
          scaled so that the largest coefficient in magnitude does not
                  (EMAX-2*P)
          exceed B          .

          Even so, the scaling strategy used in C02AGF is sometimes
          insufficient to avoid overflow and/or underflow conditions. In
          such cases, the user is recommended to scale the independent
          variable (z) so that the disparity between the largest and
          smallest coefficient in magnitude is reduced. That is, use the
          routine to locate the zeros of the polynomial d*P(cz) for some
          suitable values of c and d. For example, if the original
                               -100  100 20                   -10
          polynomial was P(z)=2    +2   z  , then choosing c=2    and
             100                                                     20
          d=2   , for instance, would yield the scaled polynomial 1+z  ,
          which is well-behaved relative to overflow and underflow and has
                           10
          zeros which are 2   times those of P(z).

          If the routine fails with IFAIL = 2 or 3, then the real and
          imaginary parts of any roots obtained before the failure occurred
          are stored in Z in the reverse order in which they were found.
          Let n  denote the number of roots found before the failure
               R
          occurred. Then Z(1,n) and Z(2,n) contain the real and imaginary
          parts of the 1st root found, Z(1,n-1) and Z(2,n-1) contain the
          real and imaginary parts of the 2nd root found, ..., Z(1,n ) and
                                                                    R
          Z(2,n ) contain the real and imaginary parts of the n th root
               R                                               R
          found. After the failure has occurred, the remaining 2*(n-n )
                                                                     R
          elements of Z contain a large negative number (equal to
                         

          -1/(X02AMF().\/2)).

          9. Example

          To find the roots of the 5th degree polynomial
           5   4   3   2
          z +2z +3z +4z +5z+6=0.

          The example program is not reproduced here. The source code for
          all example programs is distributed with the NAG Foundation
          Library software and should be available on-line.

\end{verbatim}
\endscroll
\end{page}
\begin{page}{manpageXXc05}{NAG On-line Documentation: c05}
\beginscroll
\begin{verbatim}



     C05(3NAG)         Foundation Library (12/10/92)         C05(3NAG)



          C05 -- Roots of One or More Transcendental Equations
                                                         Introduction -- C05
                                    Chapter C05
                   Roots of One or More Transcendental Equations

          1. Scope of the Chapter

          This chapter is concerned with the calculation of real zeros of
          continuous real functions of one or more variables. (Complex
          equations must be expressed in terms of the equivalent larger
          system of real equations.)

          2. Background to the Problems

          The chapter divides naturally into two parts.

          2.1. A Single Equation

          The first deals with the real zeros of a real function of a
          single variable f(x).

          At present, there is only one routine with a simple calling
          sequence. This routine assumes that the user can determine an
          initial interval [a,b] within which the desired zero lies, that
          is f(a)*f(b)<0, and outside which all other zeros lie. The
          routine then systematically subdivides the interval to produce a
          final interval containing the zero. This final interval has a
          length bounded by the user's specified error requirements; the
          end of the interval where the function has smallest magnitude is
          returned as the zero. This routine is guaranteed to converge to a
          simple zero of the function. (Here we define a simple zero as a
          zero corresponding to a sign-change of the function.) The
          algorithm used is due to Bus and Dekker.

          2.2. Systems of Equations

          The routines in the second part of this chapter are designed to
          solve a set of nonlinear equations in n unknowns


                                                           T
                   f (x)=0,  i=1,2,...,n,  x=(x ,x ,...,x )             (1)
                    i                          1  2      n

          where T stands for transpose.

          It is assumed that the functions are continuous and
          differentiable so that the matrix of first partial derivatives of
          the functions, the Jacobian matrix J  (x)=ddf /ddx  evaluated at
                                              ij       i    j
          the point x, exists, though it may not be possible to calculate
          it directly.

          The functions f  must be independent, otherwise there will be an
                         i
          infinity of solutions and the methods will fail. However, even
          when the functions are independent the solutions may not be
          unique. Since the methods are iterative, an initial guess at the
          solution has to be supplied, and the solution located will
          usually be the one closest to this initial guess.

          2.3. References

          [1]   Gill P E and Murray W (1976) Algorithms for the Solution of
                the Nonlinear Least-squares Problem. NAC 71 National
                Physical Laboratory.

          [2]   More J J, Garbow B S and Hillstrom K E (1974) User Guide for
                Minpack-1. ANL-80-74 Argonne National Laboratory.

          [3]   Ortega J M and Rheinboldt W C (1970) Iterative Solution of
                Nonlinear Equations in Several Variables. Academic Press.

          [4]   Rabinowitz P (1970) Numerical Methods for Nonlinear
                Algebraic Equations. Gordon and Breach.

          3. Recommendations on Choice and Use of Routines

          3.1. Zeros of Functions of One Variable

          There is only one routine (C05ADF) for solving a single nonlinear
          equation. This routine is designed for solving problems where the
          function f(x) whose zero is to be calculated, can be coded as a
          user-supplied routine.

          C05ADF may only be used when the user can supply an interval
          [a,b] containing the zero, that is f(a)*f(b)<0.

          3.2. Solution of Sets of Nonlinear Equations

          The solution of a set of nonlinear equations

                        f (x ,x ,...,x )=0,  i=1,2,...,n                (2)
                         i  1  2      n

          can be regarded as a special case of the problem of finding a
          minimum of a sum of squares

                           m
                           /                    2
                     s(x)= |  [f (x ,x ,...,x )]   (m>=n).              (3)
                           /    i  1  2      n
                           i=1

          So the routines in Chapter E04 of the Library are relevant as
          well as the special nonlinear equations routines.

          There are two routines (C05NBF and C05PBF) for solving a set of
          nonlinear equations. These routines require the f  (and possibly
                                                           i
          their derivatives) to be calculated in user-supplied routines.
          These should be set up carefully so the Library routines can work
          as efficiently as possible.

          The main decision which has to be made by the user is whether to
                                  ddf
                                     i
          supply the derivatives  ----. It is advisable to do so if
                                  ddx
                                     j
          possible, since the results obtained by algorithms which use
          derivatives are generally more reliable than those obtained by
          algorithms which do not use derivatives.

          C05PBF requires the user to provide the derivatives, whilst
          C05NBF does not. C05NBF and C05PBF are easy-to-use routines. A
          routine, C05ZAF, is provided for use in conjunction with C05PBF
          to check the user-provided derivatives for consistency with the
          functions themselves. The user is strongly advised to make use of
          this routine whenever C05PBF is used.

          Firstly, the calculation of the functions and their derivatives
          should be ordered so that cancellation errors are avoided. This
          is particularly important in a routine that uses these quantities
          to build up estimates of higher derivatives.

          Secondly, scaling of the variables has a considerable effect on
          the efficiency of a routine. The problem should be designed so
          that the elements of x are of similar magnitude. The same comment
          applies to the functions, all the f  should be of comparable
                                             i
          size.

          The accuracy is usually determined by the accuracy parameters of
          the routines, but the following points may be useful:

          (i)   Greater accuracy in the solution may be requested by
                choosing smaller input values for the accuracy parameters.
                However, if unreasonable accuracy is demanded, rounding
                errors may become important and cause a failure.

          (ii)  Some idea of the accuracies of the x  may be obtained by
                                                    i
                monitoring the progress of the routine to see how many
                figures remain unchanged during the last few iterations.

          (iii) An approximation to the error in the solution x, given by e
                where e is the solution to the set of linear equations

                J(x)e=-f(x)

                                                  T
                where f(x)=(f (x),f (x),...,f (x))  (see Chapter F04).
                             1     2         n

          (iv)  If the functions f (x) are changed by small amounts
                                  i
                (epsilon) , for i=1,2,...,n, then the corresponding change
                         i
                in the solution x is given approximately by (sigma), where
                (sigma) is the solution of the set of linear equations

                J(x)(sigma)=-(epsilon), (see Chapter F04).

                Thus one can estimate the sensitivity of x to any
                uncertainties in the specification of f (x), for
                                                       i
                i=1,2,...,n.

          3.3. Index

          Zeros of functions of one variable:
               Bus and Dekker algorithm                              C05ADF
          Zeros of functions of several variables:
               easy-to-use                                           C05NBF
               easy-to-use, derivatives required                     C05PBF
          Checking Routine:
               Checks user-supplied Jacobian                         C05ZAF


          C05 -- Roots of One or More Transcendental Equations
                                                             Contents -- C05
          Chapter C05

          Roots of One or More Transcendental Equations

          C05ADF  Zero of continuous function in given interval, Bus and
                  Dekker algorithm

          C05NBF  Solution of system of nonlinear equations using function
                  values only

          C05PBF  Solution of system of nonlinear equations using 1st
                  derivatives

          C05ZAF  Check user's routine for calculating 1st derivatives

\end{verbatim}
\endscroll
\end{page}
\begin{page}{manpageXXc05adf}{NAG On-line Documentation: c05adf}
\beginscroll
\begin{verbatim}



     C05ADF(3NAG)      Foundation Library (12/10/92)      C05ADF(3NAG)



          C05 -- Roots of One or More Transcendental Equations       C05ADF
                  C05ADF -- NAG Foundation Library Routine Document

          Note: Before using this routine, please read the Users' Note for
          your implementation to check implementation-dependent details.
          The symbol (*) after a NAG routine name denotes a routine that is
          not included in the Foundation Library.

          1. Purpose

          C05ADF locates a zero of a continuous function in a given
          interval by a combination of the methods of linear interpolation,
          extrapolation and bisection.

          2. Specification

                 SUBROUTINE C05ADF (A, B, EPS, ETA, F, X, IFAIL)
                 INTEGER          IFAIL
                 DOUBLE PRECISION A, B, EPS, ETA, F, X
                 EXTERNAL         F

          3. Description

          The routine attempts to obtain an approximation to a simple zero
          of the function f(x) given an initial interval [a,b] such that
          f(a)*f(b)<=0. The zero is found by calls to C05AZF(*) whose
          specification should be consulted for details of the method used.

          The approximation x to the zero (alpha) is determined so that one
          or both of the following criteria are satisfied:

               (i) |x-(alpha)|<EPS,

               (ii) |f(x)|<ETA.

          4. References

          None.

          5. Parameters

           1:  A -- DOUBLE PRECISION                                  Input
               On entry: the lower bound of the interval, a.

           2:  B -- DOUBLE PRECISION                                  Input
               On entry: the upper bound of the interval, b. Constraint: B
               /= A.

           3:  EPS -- DOUBLE PRECISION                                Input
               On entry: the absolute tolerance to which the zero is
               required (see Section 3). Constraint: EPS > 0.0.

           4:  ETA -- DOUBLE PRECISION                                Input
               On entry: a value such that if |f(x)|<ETA, x is accepted as
               the zero. ETA may be specified as 0.0 (see Section 7).

           5:  F -- DOUBLE PRECISION FUNCTION, supplied by the user.
                                                    External Procedure
               F must evaluate the function f whose zero is to be
               determined.

               Its specification is:

                      DOUBLE PRECISION FUNCTION F (XX)
                      DOUBLE PRECISION XX

                1:  XX -- DOUBLE PRECISION                            Input
                    On entry: the point at which the function must be
                    evaluated.
               F must be declared as EXTERNAL in the (sub)program from
               which C05ADF is called. Parameters denoted as Input
               must not be changed by this procedure.

           6:  X -- DOUBLE PRECISION                                 Output
               On exit: the approximation to the zero.

           7:  IFAIL -- INTEGER                                Input/Output
               Before entry, IFAIL must be assigned a value. For users not
               familiar with this parameter (described in the Essential
               Introduction) the recommended value is 0.

               Unless the routine detects an error (see Section 6), IFAIL
               contains 0 on exit.

          6. Error Indicators and Warnings

          Errors detected by the routine:

          IFAIL= 1
               On entry EPS <= 0.0,

               or       A = B,

               or       F(A)*F(B)>0.0.

          IFAIL= 2
               Too much accuracy has been requested in the computation,
               that is, EPS is too small for the computer being used. The
               final value of X is an accurate approximation to the zero.

          IFAIL= 3
               A change in sign of f(x) has been determined as occurring
               near the point defined by the final value of X. However,
               there is some evidence that this sign-change corresponds to
               a pole of f(x).

          IFAIL= 4
               Indicates that a serious error has occurred in C05AZF(*).
               Check all routine calls. Seek expert help.

          7. Accuracy

          This depends on the value of EPS and ETA. If full machine
          accuracy is required, they may be set very small, resulting in an
          error exit with IFAIL = 2, although this may involve more
          iterations than a lesser accuracy. The user is recommended to set
          ETA = 0.0 and to use EPS to control the accuracy, unless he has
          considerable knowledge of the size of f(x) for values of x near
          the zero.

          8. Further Comments

          The time taken by the routine depends primarily on the time spent
          evaluating F (see Section 5).

          If it is important to determine an interval of length less than
          EPS containing the zero, or if the function F is expensive to
          evaluate and the number of calls to F is to be restricted, then
          use of C05AZF(*) is recommended. Use of C05AZF(*) is also
          recommended when the structure of the problem to be solved does
          not permit a simple function F to be written: the reverse
          communication facilities of C05AZF(*) are more flexible than the
          direct communication of F required by C05ADF.

          9. Example

                                                            -x
          The example program below calculates the zero of e  -x within the
          interval [0,1] to approximately 5 decimal places.

          The example program is not reproduced here. The source code for
          all example programs is distributed with the NAG Foundation
          Library software and should be available on-line.

\end{verbatim}
\endscroll
\end{page}
\begin{page}{manpageXXc05nbf}{NAG On-line Documentation: c05nbf}
\beginscroll
\begin{verbatim}



     C05NBF(3NAG)      Foundation Library (12/10/92)      C05NBF(3NAG)



          C05 -- Roots of One or More Transcendental Equations       C05NBF
                  C05NBF -- NAG Foundation Library Routine Document

          Note: Before using this routine, please read the Users' Note for
          your implementation to check implementation-dependent details.
          The symbol (*) after a NAG routine name denotes a routine that is
          not included in the Foundation Library.

          1. Purpose

          C05NBF is an easy-to-use routine to find a solution of a system
          of nonlinear equations by a modification of the Powell hybrid
          method.

          2. Specification

                 SUBROUTINE C05NBF (FCN, N, X, FVEC, XTOL, WA, LWA, IFAIL)
                 INTEGER          N, LWA, IFAIL
                 DOUBLE PRECISION X(N), FVEC(N), XTOL, WA(LWA)
                 EXTERNAL         FCN

          3. Description

          The system of equations is defined as:

                        f (x ,x ,...,x )=0, for i=1,2,...,n.
                         i  1  2      n

          C05NBF is based upon the MINPACK routine HYBRD1 (More et al [1]).
          It chooses the correction at each step as a convex combination of
          the Newton and scaled gradient directions. Under reasonable
          conditions this guarantees global convergence for starting points
          far from the solution and a fast rate of convergence. The
          Jacobian is updated by the rank-1 method of Broyden. At the
          starting point the Jacobian is approximated by forward
          differences, but these are not used again until the rank-1 method
          fails to produce satisfactory progress. For more details see
          Powell [2].

          4. References

          [1]   More J J, Garbow B S and Hillstrom K E User Guide for
                MINPACK-1. Technical Report ANL-80-74. Argonne National
                Laboratory.

          [2]   Powell M J D (1970) A Hybrid Method for Nonlinear Algebraic
                Equations. Numerical Methods for Nonlinear Algebraic
                Equations. (ed P Rabinowitz) Gordon and Breach.

          5. Parameters

           1:  FCN -- SUBROUTINE, supplied by the user.
                                                    External Procedure
               FCN must return the values of the functions f  at a point x.
                                                            i

               Its specification is:

                      SUBROUTINE FCN (N, X, FVEC, IFLAG)
                      INTEGER          N, IFLAG
                      DOUBLE PRECISION X(N), FVEC(N)

                1:  N -- INTEGER                                      Input
                    On entry: the number of equations, n.

                2:  X(N) -- DOUBLE PRECISION array                    Input
                    On entry: the components of the point x at which the
                    functions must be evaluated.

                3:  FVEC(N) -- DOUBLE PRECISION array                Output
                    On exit: the function values f (x) (unless IFLAG is
                                                  i
                    set to a negative value by FCN).

                4:  IFLAG -- INTEGER                           Input/Output
                    On entry: IFLAG > 0. On exit: in general, IFLAG should
                    not be reset by FCN. If, however, the user wishes to
                    terminate execution (perhaps because some illegal point
                    X has been reached), then IFLAG should be set to a
                    negative integer. This value will be returned through
                    IFAIL.
               FCN must be declared as EXTERNAL in the (sub)program
               from which C05NBF is called. Parameters denoted as
               Input must not be changed by this procedure.

           2:  N -- INTEGER                                           Input
               On entry: the number of equations, n. Constraint: N > 0.

           3:  X(N) -- DOUBLE PRECISION array                  Input/Output
               On entry: an initial guess at the solution vector. On
               exit: the final estimate of the solution vector.

           4:  FVEC(N) -- DOUBLE PRECISION array                     Output
               On exit: the function values at the final point, X.

           5:  XTOL -- DOUBLE PRECISION                               Input
               On entry: the accuracy in X to which the solution is
               required. Suggested value: the square root of the machine
               precision. Constraint: XTOL >= 0.0.

           6:  WA(LWA) -- DOUBLE PRECISION array                  Workspace

           7:  LWA -- INTEGER                                         Input
               On entry: the dimension of the array WA. Constraint:
               LWA>=N*(3*N+13)/2.

           8:  IFAIL -- INTEGER                                Input/Output
               On entry: IFAIL must be set to 0, -1 or 1. For users not
               familiar with this parameter (described in the Essential
               Introduction) the recommended value is 0.

               On exit: IFAIL = 0 unless the routine detects an error (see
               Section 6).

          6. Error Indicators and Warnings

          Errors detected by the routine:

          If on entry IFAIL = 0 or -1, explanatory error messages are
          output on the current error message unit (as defined by X04AAF).

          IFAIL< 0
               The user has set IFLAG negative in FCN. The value of IFAIL
               will be the same as the user's setting of IFLAG.

          IFAIL= 1
               On entry N <= 0,

               or       XTOL < 0.0,

               or       LWA<N*(3*N+13)/2.

          IFAIL= 2
               There have been at least 200*(N+1) evaluations of FCN.
               Consider restarting the calculation from the final point
               held in X.

          IFAIL= 3
               No further improvement in the approximate solution X is
               possible; XTOL is too small.

          IFAIL= 4
               The iteration is not making good progress. This failure exit
               may indicate that the system does not have a zero, or that
               the solution is very close to the origin (see Section 7).
               Otherwise, rerunning C05NBF from a different starting point
               may avoid the region of difficulty.

          7. Accuracy

             ^
          If x is the true solution, C05NBF tries to ensure that

                                   ^             ^
                              ||x-x||<=XTOL*||x||.

                                                     -k
          If this condition is satisfied with XTOL=10  , then the larger
          components of x have k significant decimal digits. There is a
          danger that the smaller components of x may have large relative
          errors, but the fast rate of convergence of C05NBF usually avoids
          this possibility.

          If XTOL is less than machine precision, and the above test is
          satisfied with the machine precision in place of XTOL, then the
          routine exits with IFAIL = 3.

          Note: this convergence test is based purely on relative error,
          and may not indicate convergence if the solution is very close to
          the origin.

          The test assumes that the functions are reasonably well behaved.
          If this condition is not satisfied, then C05NBF may incorrectly
          indicate convergence. The validity of the answer can be checked,
          for example, by rerunning C05NBF with a tighter tolerance.

          8. Further Comments

          The time required by C05NBF to solve a given problem depends on n
          , the behaviour of the functions, the accuracy requested and the
          starting point. The number of arithmetic operations executed by
                                                            2
          C05NBF to process each call of FCN is about 11.5*n . Unless FCN
          can be evaluated quickly, the timing of C05NBF will be strongly
          influenced by the time spent in FCN.

          Ideally the problem should be scaled so that at the solution the
          function values are of comparable magnitude.

          9. Example

          To determine the values x ,..., x  which satisfy the tridiagonal
                                   1       9
          equations:

                                  (3-2x )x -2x =-1
                                       1  1   2

                       -x -1+(3-2x )x -2x   =-1,  i=2,3,...,8
                         i        i  i   i+1

                                  -x +(3-2x )x =-1.
                                    8      9  9

          The example program is not reproduced here. The source code for
          all example programs is distributed with the NAG Foundation
          Library software and should be available on-line.

\end{verbatim}
\endscroll
\end{page}
\begin{page}{manpageXXc05pbf}{NAG On-line Documentation: c05pbf}
\beginscroll
\begin{verbatim}



     C05PBF(3NAG)      Foundation Library (12/10/92)      C05PBF(3NAG)



          C05 -- Roots of One or More Transcendental Equations       C05PBF
                  C05PBF -- NAG Foundation Library Routine Document

          Note: Before using this routine, please read the Users' Note for
          your implementation to check implementation-dependent details.
          The symbol (*) after a NAG routine name denotes a routine that is
          not included in the Foundation Library.

          1. Purpose

          C05PBF is an easy-to-use routine to find a solution of a system
          of nonlinear equations by a modification of the Powell hybrid
          method. The user must provide the Jacobian.

          2. Specification

                 SUBROUTINE C05PBF (FCN, N, X, FVEC, FJAC, LDFJAC, XTOL,
                1                   WA, LWA, IFAIL)
                 INTEGER          N, LDFJAC, LWA, IFAIL
                 DOUBLE PRECISION X(N), FVEC(N), FJAC(LDFJAC,N), XTOL, WA
                1                 (LWA)
                 EXTERNAL         FCN

          3. Description

          The system of equations is defined as:

                          f (x ,x ,...,x )=0, i=1,2,...,n.
                           i  1  2      n

          C05PBF is based upon the MINPACK routine HYBRJ1 (More et al [1]).
          It chooses the correction at each step as a convex combination of
          the Newton and scaled gradient directions. Under reasonable
          conditions this guarantees global convergence for starting points
          far from the solution and a fast rate of convergence. The
          Jacobian is updated by the rank-1 method of Broyden. At the
          starting point the Jacobian is calculated, but it is not
          recalculated until the rank-1 method fails to produce
          satisfactory progress. For more details see Powell [2].

          4. References

          [1]   More J J, Garbow B S and Hillstrom K E User Guide for
                MINPACK-1. Technical Report ANL-80-74. Argonne National
                Laboratory.

          [2]   Powell M J D (1970) A Hybrid Method for Nonlinear Algebraic
                Equations. Numerical Methods for Nonlinear Algebraic
                Equations. (ed P Rabinowitz) Gordon and Breach.

          5. Parameters

           1:  FCN -- SUBROUTINE, supplied by the user.
                                                    External Procedure
               Depending upon the value of IFLAG, FCN must either return
               the values of the functions f  at a point x or return the
                                            i
               Jacobian at x.

               Its specification is:

                      SUBROUTINE FCN (N, X, FVEC, FJAC, LDFJAC, IFLAG)
                      INTEGER          N, LDFJAC, IFLAG
                      DOUBLE PRECISION X(N), FVEC(N), FJAC(LDFJAC,N)

                1:  N -- INTEGER                                      Input
                    On entry: the number of equations, n.

                2:  X(N) -- DOUBLE PRECISION array                    Input
                    On entry: the components of the point x at which the
                    functions or the Jacobian must be evaluated.

                3:  FVEC(N) -- DOUBLE PRECISION array                Output
                    On exit: if IFLAG = 1 on entry, FVEC must contain the
                    function values f (x) (unless IFLAG is set to a
                                     i
                    negative value by FCN). If IFLAG = 2 on entry, FVEC
                    must not be changed.

                4:  FJAC(LDFJAC,N) -- DOUBLE PRECISION array         Output
                    On exit: if IFLAG = 2 on entry, FJAC(i,j) must contain
                                  ddf
                                     i
                    the value of  ---- at the point x, for i,j=1,2,...,n
                                  ddx
                                     j
                    (unless IFLAG is set to a negative value by FCN).

                    If IFLAG = 1 on entry, FJAC must not be changed.

                5:  LDFJAC -- INTEGER                                 Input
                    On entry: the first dimension of FJAC.

                6:  IFLAG -- INTEGER                           Input/Output
                    On entry: IFLAG = 1 or 2:
                         if IFLAG = 1, FVEC is to be updated;

                         if IFLAG = 2, FJAC is to be updated.
                    On exit: in general, IFLAG should not be reset by FCN.
                    If, however, the user wishes to terminate execution
                    (perhaps because some illegal point x has been reached)
                    then IFLAG should be set to a negative integer. This
                    value will be returned through IFAIL.
               FCN must be declared as EXTERNAL in the (sub)program
               from which C05PBF is called. Parameters denoted as
               Input must not be changed by this procedure.

           2:  N -- INTEGER                                           Input
               On entry: the number of equations, n. Constraint: N > 0.

           3:  X(N) -- DOUBLE PRECISION array                  Input/Output
               On entry: an initial guess at the solution vector. On
               exit: the final estimate of the solution vector.

           4:  FVEC(N) -- DOUBLE PRECISION array                     Output
               On exit: the function values at the final point, X.

           5:  FJAC(LDFJAC,N) -- DOUBLE PRECISION array              Output
               On exit: the orthogonal matrix Q produced by the QR
               factorization of the final approximate Jacobian.

           6:  LDFJAC -- INTEGER                                      Input
               On entry:
               the first dimension of the array FJAC as declared in the
               (sub)program from which C05PBF is called.
               Constraint: LDFJAC >= N.

           7:  XTOL -- DOUBLE PRECISION                               Input
               On entry: the accuracy in X to which the solution is
               required. Suggested value: the square root of the machine
               precision. Constraint: XTOL >= 0.0.

           8:  WA(LWA) -- DOUBLE PRECISION array                  Workspace

           9:  LWA -- INTEGER                                         Input
               On entry: the dimension of the array WA. Constraint:
               LWA>=N*(N+13)/2.

          10:  IFAIL -- INTEGER                                Input/Output
               On entry: IFAIL must be set to 0, -1 or 1. For users not
               familiar with this parameter (described in the Essential
               Introduction) the recommended value is 0.

               On exit: IFAIL = 0 unless the routine detects an error (see
               Section 6).

          6. Error Indicators and Warnings

          Errors detected by the routine:

          If on entry IFAIL = 0 or -1, explanatory error messages are
          output on the current error message unit (as defined by X04AAF).

          IFAIL< 0
               A negative value of IFAIL indicates an exit from C05PBF
               because the user has set IFLAG negative in FCN. The value of
               IFAIL will be the same as the user's setting of IFLAG.

          IFAIL= 1
               On entry N <= 0,

               or       LDFJAC < N,

               or       XTOL < 0.0,

               or       LWA<N*(N+13)/2.

          IFAIL= 2
               There have been 100*(N+1) evaluations of the functions.
               Consider restarting the calculation from the final point
               held in X.

          IFAIL= 3
               No further improvement in the approximate solution X is
               possible; XTOL is too small.

          IFAIL= 4
               The iteration is not making good progress. This failure exit
               may indicate that the system does not have a zero or that
               the solution is very close to the origin (see Section 7).
               Otherwise, rerunning C05PBF from a different starting point
               may avoid the region of difficulty.

          7. Accuracy

             ^
          If x is the true solution, C05PBF tries to ensure that

                                  ^              ^
                             ||x-x|| <=XTOL*||x|| .
                                    2            2

                                                     -k
          If this condition is satisfied with XTOL=10  , then the larger
          components of x have k significant decimal digits. There is a
          danger that the smaller components of x may have large relative
          errors, but the fast rate of convergence of C05PBF usually avoids
          the possibility.

          If XTOL is less than machine precision and the above test is
          satisfied with the machine precision in place of XTOL, then the
          routine exits with IFAIL = 3.

          Note: this convergence test is based purely on relative error,
          and may not indicate convergence if the solution is very close to
          the origin.

          The test assumes that the functions and Jacobian are coded
          consistently and that the functions are reasonably well behaved.
          If these conditions are not satisfied then C05PBF may incorrectly
          indicate convergence. The coding of the Jacobian can be checked
          using C05ZAF. If the Jacobian is coded correctly, then the
          validity of the answer can be checked by rerunning C05PBF with a
          tighter tolerance.

          8. Further Comments

          The time required by C05PBF to solve a given problem depends on n
          , the behaviour of the functions, the accuracy requested and the
          starting point. The number of arithmetic operations executed by
                                2
          C05PBF is about 11.5*n  to process each evaluation of the
                                   3
          functions and about 1.3*n  to process each evaluation of the
          Jacobian. Unless FCN can be evaluated quickly, the timing of
          C05PBF will be strongly influenced by the time spent in FCN.

          Ideally the problem should be scaled so that, at the solution,
          the function values are of comparable magnitude.

          9. Example

          To determine the values x ,..., x  which satisfy the tridiagonal
                                   1       9
          equations:

                                  (3-2x )x -2x =-1
                                       1  1   2

                        -x   +(3-2x )x -2x   =-1,  i=2,3,...,8.
                          i-1      i  i   i+1

                                  -x +(3-2x )x =-1.
                                    8      9  9

          The example program is not reproduced here. The source code for
          all example programs is distributed with the NAG Foundation
          Library software and should be available on-line.

\end{verbatim}
\endscroll
\end{page}
\begin{page}{manpageXXc05zaf}{NAG On-line Documentation: c05zaf}
\beginscroll
\begin{verbatim}



     C05ZAF(3NAG)      Foundation Library (12/10/92)      C05ZAF(3NAG)



          C05 -- Roots of One or More Transcendental Equations       C05ZAF
                  C05ZAF -- NAG Foundation Library Routine Document

          Note: Before using this routine, please read the Users' Note for
          your implementation to check implementation-dependent details.
          The symbol (*) after a NAG routine name denotes a routine that is
          not included in the Foundation Library.

          1. Purpose

          C05ZAF checks the user-provided gradients of a set of non-linear
          functions in several variables, for consistency with the
          functions themselves. The routine must be called twice.

          2. Specification

                 SUBROUTINE C05ZAF (M, N, X, FVEC, FJAC, LDFJAC, XP, FVECP,
                1                   MODE, ERR)
                 INTEGER          M, N, LDFJAC, MODE
                 DOUBLE PRECISION X(N), FVEC(M), FJAC(LDFJAC,N), XP(N),
                1                 FVECP(M), ERR(M)

          3. Description

          C05ZAF is based upon the MINPACK routine CHKDER (More et al [1]).
          It checks the ith gradient for consistency with the ith function
          by computing a forward-difference approximation along a suitably
          chosen direction and comparing this approximation with the user-
          supplied gradient along the same direction. The principal
          characteristic of C05ZAF is its invariance under changes in scale
          of the variables or functions.

          4. References

          [1]   More J J, Garbow B S and Hillstrom K E User Guide for
                MINPACK-1. Technical Report ANL-80-74. Argonne National
                Laboratory.

          5. Parameters

           1:  M -- INTEGER                                           Input
               On entry: the number of functions.

           2:  N -- INTEGER                                           Input
               On entry: the number of variables. For use with C05PBF and
               C05PCF(*), M = N.

           3:  X(N) -- DOUBLE PRECISION array                         Input
               On entry: the components of a point x, at which the
               consistency check is to be made. (See Section 8.)

           4:  FVEC(M) -- DOUBLE PRECISION array                      Input
               On entry: when MODE = 2, FVEC must contain the functions
               evaluated at x.

           5:  FJAC(LDFJAC,N) -- DOUBLE PRECISION array               Input
               On entry: when MODE = 2, FJAC must contain the user-
               supplied gradients. (The ith row of FJAC must contain the
               gradient of the ith function evaluated at the point x.)

           6:  LDFJAC -- INTEGER                                      Input
               On entry:
               the first dimension of the array FJAC as declared in the
               (sub)program from which C05ZAF is called.
               Constraint: LDFJAC >= M.

           7:  XP(N) -- DOUBLE PRECISION array                       Output
               On exit: when MODE = 1, XP is set to a neighbouring point
               to X.

           8:  FVECP(M) -- DOUBLE PRECISION array                     Input
               On entry: when MODE = 2, FVECP must contain the functions
               evaluated at XP.

           9:  MODE -- INTEGER                                        Input
               On entry: the value 1 on the first call and the value 2 on
               the second call of C05ZAF.

          10:  ERR(M) -- DOUBLE PRECISION array                      Output
               On exit: when MODE = 2, ERR contains measures of
               correctness of the respective gradients. If there is no loss
               of significance (see Section 8), then if ERR(i) is 1.0 the i
               th user-supplied gradient is correct, whilst if ERR(i) is 0.
               0 the ith gradient is incorrect. For values of ERR(i)
               between 0.0 and 1.0 the categorisation is less certain. In
               general, a value of ERR(i)>0.5 indicates that the ith
               gradient is probably correct.

          6. Error Indicators and Warnings

          None.

          7. Accuracy

          See below.

          8. Further Comments

          The time required by C05ZAF increases with M and N.

          C05ZAF does not perform reliably if cancellation or rounding
          errors cause a severe loss of significance in the evaluation of a
          function. Therefore, none of the components of x should be
          unusually small (in particular, zero) or any other value which
          may cause loss of significance. The relative differences between
          corresponding elements of FVECP and FVEC should be at least two
          orders of magnitude greater than the machine precision.

          9. Example

          This example checks the Jacobian matrix for a problem with 15
          functions of 3 variables. The results indicate that the first 7
          gradients are probably incorrect (this is caused by a deliberate
          error in the code to calculate the Jacobian).

          The example program is not reproduced here. The source code for
          all example programs is distributed with the NAG Foundation
          Library software and should be available on-line.

\end{verbatim}
\endscroll
\end{page}
\begin{page}{manpageXXc06}{NAG On-line Documentation: c06}
\beginscroll
\begin{verbatim}



     C06(3NAG)         Foundation Library (12/10/92)         C06(3NAG)



          C06 -- Summation of Series                    Introduction -- C06
                                    Chapter C06
                                Summation of Series

          1. Scope of the Chapter

          This chapter is concerned with calculating the discrete Fourier
          transform of a sequence of real or complex data values, and
          applying it to calculate convolutions and correlations.

          2. Background to the Problems

          2.1. Discrete Fourier Transforms

          2.1.1.  Complex transforms

          Most of the routines in this chapter calculate the finite
          discrete Fourier transform  (DFT) of a sequence of n complex
          numbers z , for j=0,1,...,n-1. The transform is defined by:
                   j

                                  n-1
                          ^    1  --      (   2(pi)jk)
                          z = --- >  z exp(-i -------)                  (1)
                           k      --  j   (      n   )
                              \/n j=0

          for k=0,1,...,n-1. Note that equation (1) makes sense for all
                                             ^
          integral k and with this extension z  is periodic with period n,
                                              k
               ^  ^                        ^   ^
          i.e. z =z    , and in particular z  =z   .
                k  k+-n                     -k  n-k

                                    ^                                 ^
          If we write z =x +iy  and z =a +ib , then the definition of z
                       j  j   j      k  k   k                          k
          may be written in terms of sines and cosines as:

                            n-1
                         1  -- (     ( 2(pi)jk)      ( 2(pi)jk))
                    a = --- >  (x cos( -------)+y sin( -------))
                     k      -- ( j   (    n   )  j   (    n   ))
                        \/n j=0

                            n-1
                         1  -- (     ( 2(pi)jk)      ( 2(pi)jk))
                    b = --- >  (y cos( -------)-x sin( -------)).
                     k      -- ( j   (    n   )  j   (    n   ))
                        \/n j=0

          The original data values z  may conversely be recovered from the
                                    j
                    ^
          transform z  by an  inverse discrete Fourier transform:
                     k


                                  n-1
                               1  -- ^    (   2(pi)jk)
                          z = --- >  z exp(+i -------)                  (2)
                           j      --  k   (      n   )
                              \/n k=0

          for j=0,1,...,n-1. If we take the complex conjugate of (2), we
                                                               

                                                               ^
          find that the sequence z  is the DFT of the sequence z . Hence
                                  j                             k
                                          ^
          the inverse DFT of the sequence z  may be obtained by: taking the
                                           k
                                    ^
          complex conjugates of the z ; performing a DFT; and taking the
                                     k
          complex conjugates of the result.

          Notes: definitions of the discrete Fourier transform vary.
          Sometimes (2) is used as the definition of the DFT, and (1) as
                                                                      

          the definition of the inverse. Also the scale-factor of 1/\/n may
          be omitted in the definition of the DFT, and replaced by 1/n in
          the definition of the inverse.

          2.1.2. Real transforms

          If the original sequence is purely real valued, i.e. z =x , then
                                                                j  j

                                        n-1
                         ^           1  --      (   2(pi)jk)
                         z =a +ib = --- >  x exp(-i -------)
                          k  k   k      --  j   (      n   )
                                    \/n j=0

              ^                                ^
          and z    is the complex conjugate of z . Thus the DFT of a real
               n-k                              k
          sequence is a particular type of complex sequence, called a
          Hermitian sequence, or  half-complex or  conjugate symmetric with
          the properties:
                         a   =a  b   =-b  b =0 and, if n is even, b   =0.
                          n-k  k  n-k   k  0                       n/2

          Thus a Hermitian sequence of n complex data values can be
          represented by only n, rather than 2n, independent real values.
          This can obviously lead to economies in storage, the following
          scheme being used in this chapter: the real parts a  for
                                                             k
          0<=k<=n/2 are stored in normal order in the first n/2+1 locations
          of an array X of length n; the corresponding non-zero imaginary
          parts are stored in reverse order in the remaining locations of
          X. In other words, if X is declared with bounds (0:n-1) in the
                                                               ^
          user's (sub)program, the real and imaginary parts of z  are
                                                                k
          stored as follows:

                         if n=2s    if n=2s-1

              X(0)       a          a
                          0          0

              X(1)       a          a
                          1          1

              X(2)       a          a
                          2          2

              .          .          .

              .          .          .

              .          .          .

              X(s-1)     a          a
                          s-1        s-1

              X(s)       a          b
                          s          s-1

              X(s+1)     b          b
                          s-1        s-2

              .          .          .

              .          .          .

              .          .          .

              X(n-2)     b          b
                          2          2

              X(n-1)     b          b
                          1          1


                         (     n/2-1                                      )
                       1 (     --   (     ( 2(pi)jk)      ( 2(pi)jk))     )
          Hence   x = ---(a +2 >    (a cos( -------)-b sin( -------))+a   )
                   j     ( 0   --   ( k   (    n   )  k   (    n   ))  n/2)
                      \/n(     k=0                                        )

          where a    = 0 if n is odd.
                 n/2

          2.1.3.  Fourier integral transforms

          The usual application of the discrete Fourier transform is that
          of obtaining an approximation of the  Fourier integral transform


                                +infty
                                /
                          F(s)= |     f(t)exp(-i2(pi)st)dt
                                /
                                -infty

          when f(t) is negligible outside some region (0,c). Dividing the
          region into n equal intervals we have

                                    n-1
                                  c --
                           F(s)~= - >  f exp(-i2(pi)sjc/n)
                                  n --  j
                                    j=0

          and so

                                   n-1
                                 c --
                            F ~= - >  f exp(-i2(pi)jk/n)
                             k   n --  j
                                   j=0

          for k=0,1,...,n-1, where f =f(jc/n) and F =F(k/c).
                                    j              k

          Hence the discrete Fourier transform gives an approximation to
          the Fourier integral transform in the region s=0 to s=n/c.

          If the function f(t) is defined over some more general interval
          (a,b), then the integral transform can still be approximated by
          the discrete transform provided a shift is applied to move the
          point a to the origin.

          2.1.4.  Convolutions and correlations

          One of the most important applications of the discrete Fourier
          transform is to the computation of the discrete convolution or
          correlation of two vectors x and y defined (as in Brigham [1])
          by:

                                n-1
                                --
               convolution: z = >  x y
                             k  --  j k-j
                                j=0

                                n-1
                                -- 
               correlation: w = >  x y
                             k  --  j k+j
                                j=0

          (Here x and y are assumed to be periodic with period n.)

          Under certain circumstances (see Brigham [1]) these can be used
          as approximations to the convolution or correlation integrals
          defined by:

                                    +infty
                                    /
                              z(s)= |     x(t)y(s-t)dt
                                    /
                                    -infty

          and

                          +infty
                          /     
                    w(s)= |     x(t)y(s+t)dt,   -infty<s<+infty.
                          /
                          -infty

          For more general advice on the use of Fourier transforms, see
          Hamming [2]; more detailed information on the fast Fourier
          transform algorithm can be found in Van Loan [3] and Brigham [1].

          2.2. References

          [1]   Brigham E O (1973) The Fast Fourier Transform. Prentice-
                Hall.

          [2]   Hamming R W (1962) Numerical Methods for Scientists and
                Engineers. McGraw-Hill.

          [3]   Van Loan C (1992) Computational Frameworks for the Fast
                Fourier Transform. SIAM Philadelphia.

          3. Recommendations on Choice and Use of Routines

          3.1. One-dimensional Fourier Transforms

          The choice of routine is determined first of all by whether the
          data values constitute a real, Hermitian or general complex
          sequence. It is wasteful of time and storage to use an
          inappropriate routine.

          Two groups, each of three routines, are provided

                                   Group 1      Group 2

              Real sequences       C06EAF       C06FPF

              Hermitian sequences  C06EBF       C06FQF

              General complex      C06ECF       C06FRF
              sequences

          Group 1 routines each compute a single transform of length n,
          without requiring any extra working storage. The Group 1 routines
          impose some restrictions on the value of n, namely that no prime
          factor of n may exceed 19 and the total number of prime factors
          (including repetitions) may not exceed 20 (though the latter
                                                     6
          restriction only becomes relevant when n>10 ).

          Group 2 routines are designed to perform several transforms in a
          single call, all with the same value of n. They do however
          require more working storage. Even on scalar processors, they may
          be somewhat faster than repeated calls to Group 1 routines
          because of reduced overheads and because they pre-compute and
          store the required values of trigonometric functions. Group 2
          routines impose no practical restrictions on the value of n;
          however the fast Fourier transform algorithm ceases to be 'fast'
          if applied to values of n which cannot be expressed as a product
          of small prime factors. All the above routines are particularly
          efficient if the only prime factors of n are 2, 3 or 5.

          If extensive use is to be made of these routines, users who are
          concerned about efficiency are advised to conduct their own
          timing tests.

          To compute inverse discrete Fourier transforms the above routines
          should be used in conjunction with the utility routines C06GBF,
          C06GCF and C06GQF which form the complex conjugate of a Hermitian
          or general sequence of complex data values.

          3.2. Multi-dimensional Fourier Transforms

          C06FUF computes a 2-dimensional discrete Fourier transform of a
          2-dimensional sequence of complex data values. This is defined by

                         n -1 n -1
                          1    2          (   2(pi)j k )   (   2(pi)j k )
          ^         1    --   --          (         1 1)   (         2 2)
          z    = ------- >    >   z    exp(-i ---------)exp(-i ---------).
                         --   --          (      n     )   (      n     )
           k k     /n n  j =0 j =0 j j    (       1    )   (       2    )
            1 2  \/  1 2  1    2    1 2

          3.3. Convolution and Correlation

          C06EKF computes either the discrete convolution or the discrete
          correlation of two real vectors.

          3.4. Index

          Complex conjugate,
               complex sequence                                      C06GCF
               Hermitian sequence                                    C06GBF
               multiple Hermitian sequences                          C06GQF
          Complex sequence from Hermitian sequences                  C06GSF
          Convolution or Correlation
               real vectors                                          C06EKF
          Discrete Fourier Transform
               two-dimensional
                    complex sequence                                 C06FUF
               one-dimensional, multiple transforms
                    complex sequence                                 C06FRF
                    Hermitian sequence                               C06FQF
                    real sequence                                    C06FPF
               one-dimensional, single transforms
                    complex sequence                                 C06ECF
                    Hermitian sequence                               C06EBF
                    real sequence                                    C06EAF



          C06 -- Summation of Series                        Contents -- C06
          Chapter C06

          Summation of Series

          C06EAF  Single 1-D real discrete Fourier transform, no extra
                  workspace

          C06EBF  Single 1-D Hermitian discrete Fourier transform, no extra
                  workspace

          C06ECF  Single 1-D complex discrete Fourier transform, no extra
                  workspace

          C06EKF  Circular convolution or correlation of two real vectors,
                  no extra workspace

          C06FPF  Multiple 1-D real discrete Fourier transforms

          C06FQF  Multiple 1-D Hermitian discrete Fourier transforms

          C06FRF  Multiple 1-D complex discrete Fourier transforms

          C06FUF  2-D complex discrete Fourier transform

          C06GBF  Complex conjugate of Hermitian sequence

          C06GCF  Complex conjugate of complex sequence

          C06GQF  Complex conjugate of multiple Hermitian sequences

          C06GSF  Convert Hermitian sequences to general complex sequences

\end{verbatim}
\endscroll
\end{page}
\begin{page}{manpageXXc06eaf}{NAG On-line Documentation: c06eaf}
\beginscroll
\begin{verbatim}



     C06EAF(3NAG)      Foundation Library (12/10/92)      C06EAF(3NAG)



          C06 -- Summation of Series                                 C06EAF
                  C06EAF -- NAG Foundation Library Routine Document

          Note: Before using this routine, please read the Users' Note for
          your implementation to check implementation-dependent details.
          The symbol (*) after a NAG routine name denotes a routine that is
          not included in the Foundation Library.

          1. Purpose

          C06EAF calculates the discrete Fourier transform of a sequence of
          n real data values. (No extra workspace required.)

          2. Specification

                 SUBROUTINE C06EAF (X, N, IFAIL)
                 INTEGER          N, IFAIL
                 DOUBLE PRECISION X(N)

          3. Description

          Given a sequence of n real data values x , for j=0,1,...,n-1,
                                                  j
          this routine calculates their discrete Fourier transform defined
          by:

                            n-1
                    ^    1  --       (   2(pi)jk)
                    z = --- >  x *exp(-i -------), k=0,1,...,n-1.
                     k      --  j    (      n   )
                        \/n j=0

                                      1
          (Note the scale factor of  --- in this definition.) The
                                       

                                     \/n
                             ^
          transformed values z  are complex, but they form a Hermitian
                              k
                          ^                                ^
          sequence (i.e., z    is the complex conjugate of z ), so they are
                           n-k                              k
          completely determined by n real numbers (see also the Chapter
          Introduction).

          To compute the inverse discrete Fourier transform defined by:

                                   n-1
                           ^    1  --       (   2(pi)jk)
                           w = --- >  x *exp(+i -------),
                            k      --  j    (      n   )
                               \/n j=0

          this routine should be followed by a call of C06GBF to form the
                                    ^
          complex conjugates of the z .
                                     k

          The routine uses the fast Fourier transform (FFT) algorithm
          (Brigham [1]). There are some restrictions on the value of n (see
          Section 5).

          4. References

          [1]   Brigham E O (1973) The Fast Fourier Transform. Prentice-
                Hall.

          5. Parameters

           1:  X(N) -- DOUBLE PRECISION array                  Input/Output
               On entry: if X is declared with bounds (0:N-1) in the (sub)
               program from which C06EAF is called, then X(j) must contain
               x , for j=0,1,...,n-1. On exit: the discrete Fourier
                j
               transform stored in Hermitian form. If the components of the
                         ^
               transform z  are written as a +ib , and if X is declared
                          k                 k   k
               with bounds (0:N-1) in the (sub)program from which C06EAF is
               called, then for 0<=k<=n/2, a  is contained in X(k), and for
                                            k
               1<=k<=(n-1)/2, b  is contained in X(n-k). (See also Section
                               k
               2.1.2 of the Chapter Introduction, and the Example Program.)

           2:  N -- INTEGER                                           Input
               On entry: the number of data values, n. The largest prime
               factor of N must not exceed 19, and the total number of
               prime factors of N, counting repetitions, must not exceed
               20. Constraint: N > 1.

           3:  IFAIL -- INTEGER                                Input/Output
               On entry: IFAIL must be set to 0, -1 or 1. For users not
               familiar with this parameter (described in the Essential
               Introduction) the recommended value is 0.

               On exit: IFAIL = 0 unless the routine detects an error (see
               Section 6).

          6. Error Indicators and Warnings

          Errors detected by the routine:

          IFAIL= 1
               At least one of the prime factors of N is greater than 19.

          IFAIL= 2
               N has more than 20 prime factors.

          IFAIL= 3
               N <= 1.

          7. Accuracy

          Some indication of accuracy can be obtained by performing a
          subsequent inverse transform and comparing the results with the
          original sequence (in exact arithmetic they would be identical).

          8. Further Comments

          The time taken by the routine is approximately proportional to
          n*logn, but also depends on the factorization of n. The routine
          is somewhat faster than average if the only prime factors of n
          are 2, 3 or 5; and fastest of all if n is a power of 2.

          On the other hand, the routine is particularly slow if n has
          several unpaired prime factors, i.e., if the 'square-free' part
          of n has several factors. For such values of n, routine C06FAF(*)
          (which requires an additional n elements of workspace) is
          considerably faster.

          9. Example

          This program reads in a sequence of real data values, and prints
          their discrete Fourier transform (as computed by C06EAF), after
          expanding it from Hermitian form into a full complex sequence.

          It then performs an inverse transform using C06GBF and C06EBF,
          and prints the sequence so obtained alongside the original data
          values.

          The example program is not reproduced here. The source code for
          all example programs is distributed with the NAG Foundation
          Library software and should be available on-line.

\end{verbatim}
\endscroll
\end{page}
\begin{page}{manpageXXc06ebf}{NAG On-line Documentation: c06ebf}
\beginscroll
\begin{verbatim}



     C06EBF(3NAG)      Foundation Library (12/10/92)      C06EBF(3NAG)



          C06 -- Summation of Series                                 C06EBF
                  C06EBF -- NAG Foundation Library Routine Document

          Note: Before using this routine, please read the Users' Note for
          your implementation to check implementation-dependent details.
          The symbol (*) after a NAG routine name denotes a routine that is
          not included in the Foundation Library.

          1. Purpose

          C06EBF calculates the discrete Fourier transform of a Hermitian
          sequence of n complex data values. (No extra workspace required.)

          2. Specification

                 SUBROUTINE C06EBF (X, N, IFAIL)
                 INTEGER          N, IFAIL
                 DOUBLE PRECISION X(N)

          3. Description

          Given a Hermitian sequence of n complex data values z  (i.e., a
                                                               j
          sequence such that z  is real and z    is the complex conjugate
                              0              n-j
          of z , for j=1,2,...,n-1) this routine calculates their discrete
              j
          Fourier transform defined by:

                            n-1
                    ^    1  --       (   2(pi)jk)
                    x = --- >  z *exp(-i -------), k=0,1,...,n-1.
                     k      --  j    (      n   )
                        \/n j=0

                                      1
          (Note the scale factor of  --- in this definition.) The
                                       

                                     \/n
                             ^
          transformed values x  are purely real (see also the the Chapter
                              k
          Introduction).

          To compute the inverse discrete Fourier transform defined by:

                                   n-1
                           ^    1  --       (   2(pi)jk)
                           y = --- >  z *exp(+i -------),
                            k      --  j    (      n   )
                               \/n j=0

          this routine should be preceded by a call of C06GBF to form the
          complex conjugates of the z .
                                     j

          The routine uses the fast Fourier transform (FFT) algorithm
          (Brigham [1]). There are some restrictions on the value of n (see
          Section 5).

          4. References

          [1]   Brigham E O (1973) The Fast Fourier Transform. Prentice-
                Hall.

          5. Parameters

           1:  X(N) -- DOUBLE PRECISION array                  Input/Output
               On entry: the sequence to be transformed stored in
               Hermitian form. If the data values z  are written as x +iy ,
                                                   j                 j   j
               and if X is declared with bounds (0:N-1) in the subroutine
               from which C06EBF is called, then for 0<=j<=n/2, x  is
                                                                 j
               contained in X(j), and for 1<=j<=(n-1)/2, y  is contained in
                                                          j
               X(n-j). (See also Section 2.1.2 of the Chapter Introduction
               and the Example Program.) On exit: the components of the
                                          ^
               discrete Fourier transform x . If X is declared with bounds
                                           k
               (0:N-1) in the (sub)program from which C06EBF is called,
                    ^
               then x  is stored in X(k), for k=0,1,...,n-1.
                     k

           2:  N -- INTEGER                                           Input
               On entry: the number of data values, n. The largest prime
               factor of N must not exceed 19, and the total number of
               prime factors of N, counting repetitions, must not exceed
               20. Constraint: N > 1.

           3:  IFAIL -- INTEGER                                Input/Output
               On entry: IFAIL must be set to 0, -1 or 1. For users not
               familiar with this parameter (described in the Essential
               Introduction) the recommended value is 0.

               On exit: IFAIL = 0 unless the routine detects an error (see
               Section 6).

          6. Error Indicators and Warnings

          Errors detected by the routine:

          IFAIL= 1
               At least one of the prime factors of N is greater than 19.

          IFAIL= 2
               N has more than 20 prime factors.

          IFAIL= 3
               N <= 1.

          7. Accuracy

          Some indication of accuracy can be obtained by performing a
          subsequent inverse transform and comparing the results with the
          original sequence (in exact arithmetic they would be identical).

          8. Further Comments

          The time taken by the routine is approximately proportional to
          n*logn, but also depends on the factorization of n. The routine
          is somewhat faster than average if the only prime factors of n
          are 2, 3 or 5; and fastest of all if n is a power of 2.

          On the other hand, the routine is particularly slow if n has
          several unpaired prime factors, i.e., if the 'square-free' part
          of n has several factors. For such values of n, routine C06FBF(*)
          (which requires an additional n elements of workspace) is
          considerably faster.

          9. Example

          This program reads in a sequence of real data values which is
          assumed to be a Hermitian sequence of complex data values stored
          in Hermitian form. The input sequence is expanded into a full
          complex sequence and printed alongside the original sequence. The
          discrete Fourier transform (as computed by C06EBF) is printed
          out.

          The program then performs an inverse transform using C06EAF and
          C06GBF, and prints the sequence so obtained alongside the
          original data values.

          The example program is not reproduced here. The source code for
          all example programs is distributed with the NAG Foundation
          Library software and should be available on-line.

\end{verbatim}
\endscroll
\end{page}
\begin{page}{manpageXXc06ecf}{NAG On-line Documentation: c06ecf}
\beginscroll
\begin{verbatim}



     C06ECF(3NAG)      Foundation Library (12/10/92)      C06ECF(3NAG)



          C06 -- Summation of Series                                 C06ECF
                  C06ECF -- NAG Foundation Library Routine Document

          Note: Before using this routine, please read the Users' Note for
          your implementation to check implementation-dependent details.
          The symbol (*) after a NAG routine name denotes a routine that is
          not included in the Foundation Library.

          1. Purpose

          C06ECF calculates the discrete Fourier transform of a sequence of
          n complex data values. (No extra workspace required.)

          2. Specification

                 SUBROUTINE C06ECF (X, Y, N, IFAIL)
                 INTEGER          N, IFAIL
                 DOUBLE PRECISION X(N), Y(N)

          3. Description

          Given a sequence of n complex data values z , for j=0,1,...,n-1,
                                                     j
          this routine calculates their discrete Fourier transform defined
          by:

                            n-1
                    ^    1  --       (   2(pi)jk)
                    z = --- >  z *exp(-i -------), k=0,1,...,n-1.
                     k      --  j    (      n   )
                        \/n j=0

                                      1
          (Note the scale factor of  --- in this definition.)
                                       

                                     \/n

          To compute the inverse discrete Fourier transform defined by:

                                   n-1
                           ^    1  --       (   2(pi)jk)
                           w = --- >  z *exp(+i -------),
                            k      --  j    (      n   )
                               \/n j=0

          this routine should be preceded and followed by calls of C06GCF
                                                           ^
          to form the complex conjugates of the z  and the z .
                                                 j          k

          The routine uses the fast Fourier transform (FFT) algorithm
          (Brigham [1]). There are some restrictions on the value of n (see
          Section 5).

          4. References

          [1]   Brigham E O (1973) The Fast Fourier Transform. Prentice-
                Hall.

          5. Parameters

           1:  X(N) -- DOUBLE PRECISION array                  Input/Output
               On entry: if X is declared with bounds (0:N-1) in the (sub)
               program from which C06ECF is called, then X(j) must contain
               x , the real part of z , for j=0,1,...,n-1. On exit: the
                j                    j
               real parts a  of the components of the discrete Fourier
                           k
               transform. If X is declared with bounds (0:N-1) in the (sub)
               program from which C06ECF is called, then a  is contained in
                                                          k
               X(k), for k=0,1,...,n-1.

           2:  Y(N) -- DOUBLE PRECISION array                  Input/Output
               On entry: if Y is declared with bounds (0:N-1) in the (sub)
               program from which C06ECF is called, then Y(j) must contain
               y , the imaginary part of z , for j=0,1,...,n-1. On exit:
                j                         j
               the imaginary parts b  of the components of the discrete
                                    k
               Fourier transform. If Y is declared with bounds (0:N-1) in
               the (sub)program from which C06ECF is called, then b  is
                                                                   k
               contained in Y(k), for k=0,1,...,n-1.

           3:  N -- INTEGER                                           Input
               On entry: the number of data values, n. The largest prime
               factor of N must not exceed 19, and the total number of
               prime factors of N, counting repetitions, must not exceed
               20. Constraint: N > 1.

           4:  IFAIL -- INTEGER                                Input/Output
               On entry: IFAIL must be set to 0, -1 or 1. For users not
               familiar with this parameter (described in the Essential
               Introduction) the recommended value is 0.

               On exit: IFAIL = 0 unless the routine detects an error (see
               Section 6).

          6. Error Indicators and Warnings

          Errors detected by the routine:

          IFAIL= 1
               At least one of the prime factors of N is greater than 19.

          IFAIL= 2
               N has more than 20 prime factors.

          IFAIL= 3
               N <= 1.

          7. Accuracy

          Some indication of accuracy can be obtained by performing a
          subsequent inverse transform and comparing the results with the
          original sequence (in exact arithmetic they would be identical).

          8. Further Comments

          The time taken by the routine is approximately proportional to
          n*logn, but also depends on the factorization of n. The routine
          is somewhat faster than average if the only prime factors of n
          are 2, 3 or 5; and fastest of all if n is a power of 2.

          On the other hand, the routine is particularly slow if n has
          several unpaired prime factors, i.e., if the 'square-free' part
          of n has several factors. For such values of n, routine C06FCF(*)
          (which requires an additional n real elements of workspace) is
          considerably faster.

          9. Example

          This program reads in a sequence of complex data values and
          prints their discrete Fourier transform.

          It then performs an inverse transform using C06GCF and C06ECF,
          and prints the sequence so obtained alongside the original data
          values.

          The example program is not reproduced here. The source code for
          all example programs is distributed with the NAG Foundation
          Library software and should be available on-line.

\end{verbatim}
\endscroll
\end{page}
\begin{page}{manpageXXc06ekf}{NAG On-line Documentation: c06ekf}
\beginscroll
\begin{verbatim}



     C06EKF(3NAG)      Foundation Library (12/10/92)      C06EKF(3NAG)



          C06 -- Summation of Series                                 C06EKF
                  C06EKF -- NAG Foundation Library Routine Document

          Note: Before using this routine, please read the Users' Note for
          your implementation to check implementation-dependent details.
          The symbol (*) after a NAG routine name denotes a routine that is
          not included in the Foundation Library.

          1. Purpose

          C06EKF calculates the circular convolution or correlation of two
          real vectors of period n. No extra workspace is required.

          2. Specification

                 SUBROUTINE C06EKF (JOB, X, Y, N, IFAIL)
                 INTEGER          JOB, N, IFAIL
                 DOUBLE PRECISION X(N), Y(N)

          3. Description

          This routine computes:

          if JOB =1, the discrete convolution of x and y, defined by:

                                  n-1        n-1
                                  --         --
                              z = >  x y   = >  x   y ;
                               k  --  j k-j  --  k-j j
                                  j=0        j=0

          if JOB =2, the discrete correlation of x and y defined by:

                                       n-1
                                       --
                                   w = >  x y   .
                                    k  --  j k+j
                                       j=0

          Here x and y are real vectors, assumed to be periodic, with
          period n, i.e., x =x    =x     =...; z and w are then also
                           j  j+-n  j+-2n
          periodic with period n.

          Note: this usage of the terms 'convolution' and 'correlation' is
          taken from Brigham [1]. The term 'convolution' is sometimes used
          to denote both these computations.

             ^  ^  ^     ^
          If x, y, z and w are the discrete Fourier transforms of these
          sequences,

                        n-1
                ^    1  --       (   2(pi)jk)
          i.e., x = --- >  x *exp(-i -------), etc,
                 k      --  j    (      n   )
                    \/n j=0

                ^      ^ ^
          then  z =\/n.x y
                 k      k k

                       

                ^      ^ ^
          and   w =\/n.x y
                 k      k k

          (the bar denoting complex conjugate).

          This routine calls the same auxiliary routines as C06EAF and
          C06EBF to compute discrete Fourier transforms, and there are some
          restrictions on the value of n.

          4. References

          [1]   Brigham E O (1973) The Fast Fourier Transform. Prentice-
                Hall.

          5. Parameters

           1:  JOB -- INTEGER                                         Input
               On entry: the computation to be performed:
                                    n-1
                                    --
                    if JOB = 1, z = >  x y    (convolution);
                                 k  --  j k-j
                                    j=0

                                    n-1
                                    --
                    if JOB = 2, w = >  x y    (correlation).
                                 k  --  j k+j
                                    j=0
               Constraint: JOB = 1 or 2.

           2:  X(N) -- DOUBLE PRECISION array                  Input/Output
               On entry: the elements of one period of the vector x. If X
               is declared with bounds (0:N-1) in the (sub)program from
               which C06EKF is called, then X(j) must contain x , for
                                                               j
               j=0,1,...,n-1. On exit: the corresponding elements of the
               discrete convolution or correlation.

           3:  Y(N) -- DOUBLE PRECISION array                  Input/Output
               On entry: the elements of one period of the vector y. If Y
               is declared with bounds (0:N-1) in the (sub)program from
               which C06EKF is called, then Y(j) must contain y , for
                                                               j
               j=0,1,...,n-1. On exit: the discrete Fourier transform of
               the convolution or correlation returned in the array X; the
               transform is stored in Hermitian form, exactly as described
               in the document  C06EAF.

           4:  N -- INTEGER                                           Input
               On entry: the number of values, n, in one period of the
               vectors X and Y. The largest prime factor of N must not
               exceed 19, and the total number of prime factors of N,
               counting repetitions, must not exceed 20. Constraint: N > 1.

           5:  IFAIL -- INTEGER                                Input/Output
               On entry: IFAIL must be set to 0, -1 or 1. For users not
               familiar with this parameter (described in the Essential
               Introduction) the recommended value is 0.

               On exit: IFAIL = 0 unless the routine detects an error (see
               Section 6).

          6. Error Indicators and Warnings

          Errors detected by the routine:

          IFAIL= 1
               At least one of the prime factors of N is greater than 19.

          IFAIL= 2
               N has more than 20 prime factors.

          IFAIL= 3
               N <= 1.

          IFAIL= 4
               JOB /= 1 or 2.

          7. Accuracy

          The results should be accurate to within a small multiple of the
          machine precision.

          8. Further Comments

          The time taken by the routine is approximately proportional to
          n*logn, but also depends on the factorization of n. The routine
          is faster than average if the only prime factors are 2, 3 or 5;
          and fastest of all if n is a power of 2.

          The routine is particularly slow if n has several unpaired prime
          factors, i.e., if the 'square free' part of n has several
          factors. For such values of n, routine C06FKF(*) is considerably
          faster (but requires an additional workspace of n elements).

          9. Example

          This program reads in the elements of one period of two real
          vectors x and y and prints their discrete convolution and
          correlation (as computed by C06EKF). In realistic computations
          the number of data values would be much larger.

          The example program is not reproduced here. The source code for
          all example programs is distributed with the NAG Foundation
          Library software and should be available on-line.

\end{verbatim}
\endscroll
\end{page}
\begin{page}{manpageXXc06fpf}{NAG On-line Documentation: c06fpf}
\beginscroll
\begin{verbatim}



     C06FPF(3NAG)      Foundation Library (12/10/92)      C06FPF(3NAG)



          C06 -- Summation of Series                                 C06FPF
                  C06FPF -- NAG Foundation Library Routine Document

          Note: Before using this routine, please read the Users' Note for
          your implementation to check implementation-dependent details.
          The symbol (*) after a NAG routine name denotes a routine that is
          not included in the Foundation Library.

          1. Purpose

          C06FPF computes the discrete Fourier transforms of m sequences,
          each containing n real data values. This routine is designed to
          be particularly efficient on vector processors.

          2. Specification

                 SUBROUTINE C06FPF (M, N, X, INIT, TRIG, WORK, IFAIL)
                 INTEGER          M, N, IFAIL
                 DOUBLE PRECISION X(M*N), TRIG(2*N), WORK(M*N)
                 CHARACTER*1      INIT

          3. Description

                                                   p
          Given m sequences of n real data values x , for j=0,1,...,n-1;
                                                   j
          p=1,2,...,m, this routine simultaneously calculates the Fourier
          transforms of all the sequences defined by:

                     n-1
             ^p   1  --  p    (   2(pi)jk)
             z = --- >  x *exp(-i -------), k=0,1,...,n-1; p=1,2,...,m.
              k      --  j    (      n   )
                 \/n j=0

                                   1
          (Note the scale factor  --- in this definition.)
                                    

                                  \/n

                                 ^p
          The transformed values z  are complex, but for each value of p
                                  k
              ^p                                 ^p
          the z  form a Hermitian sequence (i.e.,z    is the complex
               k                                  n-k
                       ^p
          conjugate of z ), so they are completely determined by mn real
                        k
          numbers (see also the Chapter Introduction).

          The discrete Fourier transform is sometimes defined using a
          positive sign in the exponential term:

                                   n-1
                           ^p   1  --  p    (   2(pi)jk)
                           z = --- >  x *exp(+i -------).
                            k      --  j    (      n   )
                               \/n j=0

          To compute this form, this routine should be followed by a call
                                                          ^p
          to C06GQF to form the complex conjugates of the z .
                                                           k

          The routine uses a variant of the fast Fourier transform (FFT)
          algorithm (Brigham [1]) known as the Stockham self-sorting
          algorithm, which is described in Temperton [2]. Special coding is
          provided for the factors 2, 3, 4, 5 and 6. This routine is
          designed to be particularly efficient on vector processors, and
          it becomes especially fast as M, the number of transforms to be
          computed in parallel, increases.

          4. References

          [1]   Brigham E O (1973) The Fast Fourier Transform. Prentice-
                Hall.

          [2]   Temperton C (1983) Fast Mixed-Radix Real Fourier Transforms.
                J. Comput. Phys. 52 340--350.

          5. Parameters

           1:  M -- INTEGER                                           Input
               On entry: the number of sequences to be transformed, m.
               Constraint: M >= 1.

           2:  N -- INTEGER                                           Input
               On entry: the number of real values in each sequence, n.
               Constraint: N >= 1.

           3:  X(M,N) -- DOUBLE PRECISION array                Input/Output
               On entry: the data must be stored in X as if in a two-
               dimensional array of dimension (1:M,0:N-1); each of the m
               sequences is stored in a row of the array. In other words,
               if the data values of the pth sequence to be transformed are
                           p
               denoted by x , for j=0,1,...,n-1, then the mn elements of
                           j
               the array X must contain the values
                    1  2      m   1  2       m      1    2        m
                   x ,x ,...,x , x ,x ,..., x ,...,x   ,x   ,...,x   .
                    0  0      0   1  1       1      n-1  n-1      n-1
               On exit: the m discrete Fourier transforms stored as if in
               a two-dimensional array of dimension (1:M,0:N-1). Each of
               the m transforms is stored in a row of the array in
               Hermitian form, overwriting the corresponding original
               sequence. If the n components of the discrete Fourier
                         ^p                 p   p                       p
               transform z  are written as a +ib , then for 0<=k<=n/2, a
                          k                 k   k                       k
                                                               p
               is contained in X(p,k), and for 1<=k<=(n-1)/2, b  is
                                                               k
               contained in X(p,n-k). (See also Section 2.1.2 of the
               Chapter Introduction.)

           4:  INIT -- CHARACTER*1                                    Input
               On entry: if the trigonometric coefficients required to
               compute the transforms are to be calculated by the routine
               and stored in the array TRIG, then INIT must be set equal to
               'I' (Initial call).

               If INIT contains 'S' (Subsequent call), then the routine
               assumes that trigonometric coefficients for the specified
               value of n are supplied in the array TRIG, having been
               calculated in a previous call to one of C06FPF, C06FQF or
               C06FRF.

               If INIT contains 'R' (Restart then the routine assumes that
               trigonometric coefficients for the particular value of n are
               supplied in the array TRIG, but does not check that C06FPF,
               C06FQF or C06FRF have previously been called. This option
               allows the TRIG array to be stored in an external file, read
               in and re-used without the need for a call with INIT equal
               to 'I'. The routine carries out a simple test to check that
               the current value of n is consistent with the array TRIG.
               Constraint: INIT = 'I', 'S' or 'R'.

           5:  TRIG(2*N) -- DOUBLE PRECISION array             Input/Output
               On entry: if INIT = 'S' or 'R', TRIG must contain the
               required coefficients calculated in a previous call of the
               routine. Otherwise TRIG need not be set. On exit: TRIG
               contains the required coefficients (computed by the routine
               if INIT = 'I').

           6:  WORK(M*N) -- DOUBLE PRECISION array                Workspace

           7:  IFAIL -- INTEGER                                Input/Output
               On entry: IFAIL must be set to 0, -1 or 1. For users not
               familiar with this parameter (described in the Essential
               Introduction) the recommended value is 0.

               On exit: IFAIL = 0 unless the routine detects an error (see
               Section 6).


          6. Error Indicators and Warnings

          Errors detected by the routine:

          If on entry IFAIL = 0 or -1, explanatory error messages are
          output on the current error message unit (as defined by X04AAF).

          IFAIL= 1
               On entry M < 1.

          IFAIL= 2
               N < 1.

          IFAIL= 3
               INIT is not one of 'I', 'S' or 'R'.

          IFAIL= 4
               INIT = 'S', but none of C06FPF, C06FQF or C06FRF has
               previously been called.

          IFAIL= 5
               INIT = 'S' or 'R', but the array TRIG and the current value
               of N are inconsistent.

          7. Accuracy

          Some indication of accuracy can be obtained by performing a
          subsequent inverse transform and comparing the results with the
          original sequence (in exact arithmetic they would be identical).

          8. Further Comments

          The time taken by the routine is approximately proportional to
          nm*logn, but also depends on the factors of n. The routine is
          fastest if the only prime factors of n are 2, 3 and 5, and is
          particularly slow if n is a large prime, or has large prime
          factors.

          9. Example

          This program reads in sequences of real data values and prints
          their discrete Fourier transforms (as computed by C06FPF). The
          Fourier transforms are expanded into full complex form using
          C06GSF and printed. Inverse transforms are then calculated by
          calling C06GQF followed by C06FQF showing that the original
          sequences are restored.

          The example program is not reproduced here. The source code for
          all example programs is distributed with the NAG Foundation
          Library software and should be available on-line.

\end{verbatim}
\endscroll
\end{page}
\begin{page}{manpageXXc06fqf}{NAG On-line Documentation: c06fqf}
\beginscroll
\begin{verbatim}



     C06FQF(3NAG)      Foundation Library (12/10/92)      C06FQF(3NAG)



          C06 -- Summation of Series                                 C06FQF
                  C06FQF -- NAG Foundation Library Routine Document

          Note: Before using this routine, please read the Users' Note for
          your implementation to check implementation-dependent details.
          The symbol (*) after a NAG routine name denotes a routine that is
          not included in the Foundation Library.

          1. Purpose

          C06FQF computes the discrete Fourier transforms of m Hermitian
          sequences, each containing n complex data values. This routine is
          designed to be particularly efficient on vector processors.

          2. Specification

                 SUBROUTINE C06FQF (M, N, X, INIT, TRIG, WORK, IFAIL)
                 INTEGER          M, N, IFAIL
                 DOUBLE PRECISION X(M*N), TRIG(2*N), WORK(M*N)
                 CHARACTER*1      INIT

          3. Description

                                                                p
          Given m Hermitian sequences of n complex data values z , for
                                                                j
          j=0,1,...,n-1; p=1,2,...,m, this routine simultaneously
          calculates the Fourier transforms of all the sequences defined
          by:

                     n-1
             ^p   1  --  p    (   2(pi)jk)
             x = --- >  z *exp(-i -------), k=0,1,...,n-1; p=1,2,...,m.
              k      --  j    (      n   )
                 \/n j=0

                                   1
          (Note the scale factor  --- in this definition.)
                                    

                                  \/n

          The transformed values are purely real (see also the Chapter
          Introduction).

          The discrete Fourier transform is sometimes defined using a
          positive sign in the exponential term

                                   n-1
                           ^p   1  --  p    (   2(pi)jk)
                           x = --- >  z *exp(+i -------).
                            k      --  j    (      n   )
                               \/n j=0

          To compute this form, this routine should be preceded by a call
                                                          ^p
          to C06GQF to form the complex conjugates of the z .
                                                           j

          The routine uses a variant of the fast Fourier transform (FFT)
          algorithm (Brigham [1]) known as the Stockham self-sorting
          algorithm, which is described in Temperton [2]. Special code is
          included for the factors 2, 3, 4, 5 and 6. This routine is
          designed to be particularly efficient on vector processors, and
          it becomes especially fast as m, the number of transforms to be
          computed in parallel, increases.

          4. References

          [1]   Brigham E O (1973) The Fast Fourier Transform. Prentice-
                Hall.

          [2]   Temperton C (1983) Fast Mixed-Radix Real Fourier Transforms.
                J. Comput. Phys. 52 340--350.

          5. Parameters

           1:  M -- INTEGER                                           Input
               On entry: the number of sequences to be transformed, m.
               Constraint: M >= 1.

           2:  N -- INTEGER                                           Input
               On entry: the number of data values in each sequence, n.
               Constraint: N >= 1.

           3:  X(M,N) -- DOUBLE PRECISION array                Input/Output
               On entry: the data must be stored in X as if in a two-
               dimensional array of dimension (1:M,0:N-1); each of the m
               sequences is stored in a row of the array in Hermitian form.
                                     p                 p   p
               If the n data values z  are written as x +iy , then for
                                     j                 j   j
                           p
               0<=j<=n/2, x  is contained in X(p,j), and for 1<=j<=(n-1)/2,
                           j
                p
               y  is contained in X(p,n-j). (See also Section 2.1.2 of the
                j
               Chapter Introduction.) On exit: the components of the m
               discrete Fourier transforms, stored as if in a two-
               dimensional array of dimension (1:M,0:N-1). Each of the m
               transforms is stored as a row of the array, overwriting the
               corresponding original sequence. If the n components of the
                                                         ^p
               discrete Fourier transform are denoted by x , for
                                                          k
               k=0,1,...,n-1, then the mn elements of the array X contain
               the values
                   ^1 ^2     ^m  ^1 ^2      ^m     ^1   ^2       ^m
                   x ,x ,...,x , x ,x ,..., x ,...,x   ,x   ,...,x   .
                    0  0      0   1  1       1      n-1  n-1      n-1

           4:  INIT -- CHARACTER*1                                    Input
               On entry: if the trigonometric coefficients required to
               compute the transforms are to be calculated by the routine
               and stored in the array TRIG, then INIT must be set equal to
               'I' (Initial call).

               If INIT contains 'S' (Subsequent call), then the routine
               assumes that trigonometric coefficients for the specified
               value of n are supplied in the array TRIG, having been
               calculated in a previous call to one of C06FPF, C06FQF or
               C06FRF.

               If INIT contains 'R' (Restart), then the routine assumes
               that trigonometric coefficients for the particular value of
               N are supplied in the array TRIG, but does not check that
               C06FPF, C06FQF or C06FRF have previously been called. This
               option allows the TRIG array to be stored in an external
               file, read in and re-used without the need for a call with
               INIT equal to 'I'. The routine carries out a simple test to
               check that the current value of n is compatible with the
               array TRIG. Constraint: INIT = 'I', 'S' or 'R'.

           5:  TRIG(2*N) -- DOUBLE PRECISION array             Input/Output
               On entry: if INIT = 'S' or 'R', TRIG must contain the
               required coefficients calculated in a previous call of the
               routine. Otherwise TRIG need not be set. On exit: TRIG
               contains the required coefficients (computed by the routine
               if INIT = 'I').

           6:  WORK(M*N) -- DOUBLE PRECISION array                Workspace

           7:  IFAIL -- INTEGER                                Input/Output
               On entry: IFAIL must be set to 0, -1 or 1. For users not
               familiar with this parameter (described in the Essential
               Introduction) the recommended value is 0.

               On exit: IFAIL = 0 unless the routine detects an error (see
               Section 6).

          6. Error Indicators and Warnings

          Errors detected by the routine:

          If on entry IFAIL = 0 or -1, explanatory error messages are
          output on the current error message unit (as defined by X04AAF).

          IFAIL= 1
               On entry M < 1.

          IFAIL= 2
               On entry N < 1.

          IFAIL= 3
               On entry INIT is not one of 'I', 'S' or 'R'.

          IFAIL= 4
               On entry INIT = 'S', but none of C06FPF, C06FQF and C06FRF
               has previously been called.

          IFAIL= 5
               On entry INIT = 'S' or 'R', but the array TRIG and the
               current value of n are inconsistent.

          7. Accuracy

          Some indication of accuracy can be obtained by performing a
          subsequent inverse transform and comparing the results with the
          original sequence (in exact arithmetic they would be identical).

          8. Further Comments

          The time taken by the routine is approximately proportional to
          nm*logn, but also depends on the factors of n. The routine is
          fastest if the only prime factors of n are 2, 3 and 5, and is
          particularly slow if n is a large prime, or has large prime
          factors.

          9. Example

          This program reads in sequences of real data values which are
          assumed to be Hermitian sequences of complex data stored in
          Hermitian form. The sequences are expanded into full complex form
          using C06GSF and printed. The discrete Fourier transforms are
          then computed (using C06FQF) and printed out. Inverse transforms
          are then calculated by calling C06FPF followed by C06GQF showing
          that the original sequences are restored.

          The example program is not reproduced here. The source code for
          all example programs is distributed with the NAG Foundation
          Library software and should be available on-line.

\end{verbatim}
\endscroll
\end{page}
\begin{page}{manpageXXc06frf}{NAG On-line Documentation: c06frf}
\beginscroll
\begin{verbatim}



     C06FRF(3NAG)      Foundation Library (12/10/92)      C06FRF(3NAG)



          C06 -- Summation of Series                                 C06FRF
                  C06FRF -- NAG Foundation Library Routine Document

          Note: Before using this routine, please read the Users' Note for
          your implementation to check implementation-dependent details.
          The symbol (*) after a NAG routine name denotes a routine that is
          not included in the Foundation Library.

          1. Purpose

          C06FRF computes the discrete Fourier transforms of m sequences,
          each containing n complex data values. This routine is designed
          to be particularly efficient on vector processors.

          2. Specification

                 SUBROUTINE C06FRF (M, N, X, Y, INIT, TRIG, WORK, IFAIL)
                 INTEGER          M, N, IFAIL
                 DOUBLE PRECISION X(M*N), Y(M*N), TRIG(2*N), WORK(2*M*N)
                 CHARACTER*1      INIT

          3. Description

                                                      p
          Given m sequences of n complex data values z , for j=0,1,...,n-1;
                                                      j
          p=1,2,...,m, this routine simultaneously calculates the Fourier
          transforms of all the sequences defined by:

                     n-1
             ^p   1  --  p    (   2(pi)jk)
             z = --- >  z *exp(-i -------), k=0,1,...,n-1; p=1,2,...,m.
              k      --  j    (      n   )
                 \/n j=0

                                   1
          (Note the scale factor  --- in this definition.)
                                    

                                  \/n

          The discrete Fourier transform is sometimes defined using a
          positive sign in the exponential term

                                   n-1
                           ^p   1  --  p    (   2(pi)jk)
                           z = --- >  z *exp(+i -------).
                            k      --  j    (      n   )
                               \/n j=0

          To compute this form, this routine should be preceded and
          followed by a call of C06GCF to form the complex conjugates of
               p         ^p
          the z  and the z .
               j          k

          The routine uses a variant of the fast Fourier transform (FFT)
          algorithm (Brigham [1]) known as the Stockham self-sorting
          algorithm, which is described in Temperton [2]. Special code is
          provided for the factors 2, 3, 4, 5 and 6. This routine is
          designed to be particularly efficient on vector processors, and
          it becomes especially fast as m, the number of transforms to be
          computed in parallel, increases.

          4. References

          [1]   Brigham E O (1973) The Fast Fourier Transform. Prentice-
                Hall.

          [2]   Temperton C (1983) Self-sorting Mixed-radix Fast Fourier
                Transforms. J. Comput. Phys. 52 1--23.

          5. Parameters

           1:  M -- INTEGER                                           Input
               On entry: the number of sequences to be transformed, m.
               Constraint: M >= 1.

           2:  N -- INTEGER                                           Input
               On entry: the number of complex values in each sequence, n.
               Constraint: N >= 1.

           3:  X(M,N) -- DOUBLE PRECISION array                Input/Output

           4:  Y(M,N) -- DOUBLE PRECISION array                Input/Output
               On entry: the real and imaginary parts of the complex data
               must be stored in X and Y respectively as if in a two-
               dimensional array of dimension (1:M,0:N-1); each of the m
               sequences is stored in a row of each array. In other words,
               if the real parts of the pth sequence to be transformed are
                           p
               denoted by x , for j=0,1,...,n-1, then the mn elements of
                           j
               the array X must contain the values
                    1  2      m   1  2      m       1    2        m
                   x ,x ,...,x , x ,x ,...,x ,..., x   ,x   ,...,x   .
                    0  0      0   1  1      1       n-1  n-1      n-1
               On exit: X and Y are overwritten by the real and imaginary
               parts of the complex transforms.

           5:  INIT -- CHARACTER*1                                    Input
               On entry: if the trigonometric coefficients required to
               compute the transforms are to be calculated by the routine
               and stored in the array TRIG, then INIT must be set equal to
               'I' (Initial call).

               If INIT contains 'S' (Subsequent call), then the routine
               assumes that trigonometric coefficients for the specified
               value of n are supplied in the array TRIG, having been
               calculated in a previous call to one of C06FPF, C06FQF or
               C06FRF.

               If INIT contains 'R' (Restart) then the routine assumes that
               trigonometric coefficients for the particular value of n are
               supplied in the array TRIG, but does not check that C06FPF,
               C06FQF or C06FRF have previously been called. This option
               allows the TRIG array to be stored in an external file, read
               in and re-used without the need for a call with INIT equal
               to 'I'. The routine carries out a simple test to check that
               the current value of n is compatible with the array TRIG.
               Constraint: INIT = 'I', 'S' or 'R'.

           6:  TRIG(2*N) -- DOUBLE PRECISION array             Input/Output
               On entry: if INIT = 'S' or 'R', TRIG must contain the
               required coefficients calculated in a previous call of the
               routine. Otherwise TRIG need not be set. On exit: TRIG
               contains the required coefficients (computed by the routine
               if INIT = 'I').

           7:  WORK(2*M*N) -- DOUBLE PRECISION array              Workspace

           8:  IFAIL -- INTEGER                                Input/Output
               On entry: IFAIL must be set to 0, -1 or 1. For users not
               familiar with this parameter (described in the Essential
               Introduction) the recommended value is 0.

               On exit: IFAIL = 0 unless the routine detects an error (see
               Section 6).

          6. Error Indicators and Warnings

          Errors detected by the routine:

          If on entry IFAIL = 0 or -1, explanatory error messages are
          output on the current error message unit (as defined by X04AAF).

          IFAIL= 1
               On entry M < 1.

          IFAIL= 2
               On entry N < 1.

          IFAIL= 3
               On entry INIT is not one of 'I', 'S' or 'R'.

          IFAIL= 4
               On entry INIT = 'S', but none of C06FPF, C06FQF and C06FRF
               has previously been called.

          IFAIL= 5
               On entry INIT = 'S' or 'R', but the array TRIG and the
               current value of n are inconsistent.

          7. Accuracy

          Some indication of accuracy can be obtained by performing a
          subsequent inverse transform and comparing the results with the
          original sequence (in exact arithmetic they would be identical).

          8. Further Comments

          The time taken by the routine is approximately proportional to
          nm*logn, but also depends on the factors of n. The routine is
          fastest if the only prime factors of n are 2, 3 and 5, and is
          particularly slow if n is a large prime, or has large prime
          factors.

          9. Example

          This program reads in sequences of complex data values and prints
          their discrete Fourier transforms (as computed by C06FRF).
          Inverse transforms are then calculated using C06GCF and C06FRF
          and printed out, showing that the original sequences are
          restored.

          The example program is not reproduced here. The source code for
          all example programs is distributed with the NAG Foundation
          Library software and should be available on-line.

\end{verbatim}
\endscroll
\end{page}
\begin{page}{manpageXXc06fuf}{NAG On-line Documentation: c06fuf}
\beginscroll
\begin{verbatim}



     C06FUF(3NAG)      Foundation Library (12/10/92)      C06FUF(3NAG)



          C06 -- Summation of Series                                 C06FUF
                  C06FUF -- NAG Foundation Library Routine Document

          Note: Before using this routine, please read the Users' Note for
          your implementation to check implementation-dependent details.
          The symbol (*) after a NAG routine name denotes a routine that is
          not included in the Foundation Library.

          1. Purpose

          C06FUF computes the two-dimensional discrete Fourier transform of
          a bivariate sequence of complex data values. This routine is
          designed to be particularly efficient on vector processors.

          2. Specification

                 SUBROUTINE C06FUF (M, N, X, Y, INIT, TRIGM, TRIGN, WORK,
                1                   IFAIL)
                 INTEGER          M, N, IFAIL
                 DOUBLE PRECISION X(M*N), Y(M*N), TRIGM(2*M), TRIGN(2*N),
                1                 WORK(2*M*N)
                 CHARACTER*1      INIT

          3. Description

          This routine computes the two-dimensional discrete Fourier
          transform of a bivariate sequence of complex data values z    ,
                                                                    j j
                                                                     1 2
          where j =0,1,...,m-1, j =0,1,...,n-1.
                 1               2

          The discrete Fourier transform is here defined by:

                            m-1  n-1          (       ( j k   j k ))
                ^       1   --   --           (       (  1 1   2 2))
                z    = ---- >    >   z    *exp(-2(pi)i( ----+ ----)),
                            --   --   j j     (       (  m     n  ))
                 k k   \/mn j =0 j =0  1 2
                  1 2        1    2

          where k =0,1,...,m-1, k =0,1,...,n-1.
                 1               2

                                      1
          (Note the scale factor of  ---- in this definition.)
                                       

                                     \/mn

          To compute the inverse discrete Fourier transform, defined with
          exp(+2(pi)i(...)) in the above formula instead of exp(-
          2(pi)i(...)), this routine should be preceded and followed by
          calls of C06GCF to form the complex conjugates of the data values
          and the transform.

          This routine calls C06FRF to perform multiple one-dimensional
          discrete Fourier transforms by the fast Fourier transform (FFT)
          algorithm in Brigham [1]. It is designed to be particularly
          efficient on vector processors.

          4. References

          [1]   Brigham E O (1973) The Fast Fourier Transform. Prentice-
                Hall.

          [2]   Temperton C (1983) Self-sorting Mixed-radix Fast Fourier
                Transforms. J. Comput. Phys. 52 1--23.

          5. Parameters

           1:  M -- INTEGER                                           Input
               On entry: the number of rows, m, of the arrays X and Y.
               Constraint: M >= 1.

           2:  N -- INTEGER                                           Input
               On entry: the number of columns, n, of the arrays X and Y.
               Constraint: N >= 1.

           3:  X(M,N) -- DOUBLE PRECISION array                Input/Output

           4:  Y(M,N) -- DOUBLE PRECISION array                Input/Output
               On entry: the real and imaginary parts of the complex data
               values must be stored in arrays X and Y respectively. If X
               and Y are regarded as two-dimensional arrays of dimension
               (0:M-1,0:N-1), then X(j ,j ) and Y(j ,j ) must contain the
                                      1  2         1  2
               real and imaginary parts of z    . On exit: the real and
                                            j j
                                             1 2
               imaginary parts respectively of the corresponding elements
               of the computed transform.

           5:  INIT -- CHARACTER*1                                    Input
               On entry: if the trigonometric coefficients required to
               compute the transforms are to be calculated by the routine
               and stored in the arrays TRIGM and TRIGN, then INIT must be
               set equal to 'I', (Initial call).

               If INIT contains 'S', (Subsequent call), then the routine
               assumes that trigonometric coefficients for the specified
               values of m and n are supplied in the arrays TRIGM and
               TRIGN, having been calculated in a previous call to the
               routine.

               If INIT contains 'R', (Restart), then the routine assumes
               that trigonometric coefficients for the particular values of
               m and n are supplied in the arrays TRIGM and TRIGN, but does
               not check that the routine has previously been called. This
               option allows the TRIGM and TRIGN arrays to be stored in an
               external file, read in and re-used without the need for a
               call with INIT equal to 'I'. The routine carries out a
               simple test to check that the current values of m and n are
               compatible with the arrays TRIGM and TRIGN. Constraint: INIT
               = 'I', 'S' or 'R'.

           6:  TRIGM(2*M) -- DOUBLE PRECISION array            Input/Output

           7:  TRIGN(2*N) -- DOUBLE PRECISION array            Input/Output
               On entry: if INIT = 'S' or 'R',TRIGM and TRIGN must contain
               the required coefficients calculated in a previous call of
               the routine. Otherwise TRIGM and TRIGN need not be set.

               If m=n the same array may be supplied for TRIGM and TRIGN.
               On exit: TRIGM and TRIGN contain the required coefficients
               (computed by the routine if INIT = 'I').

           8:  WORK(2*M*N) -- DOUBLE PRECISION array              Workspace

           9:  IFAIL -- INTEGER                                Input/Output
               On entry: IFAIL must be set to 0, -1 or 1. For users not
               familiar with this parameter (described in the Essential
               Introduction) the recommended value is 0.

               On exit: IFAIL = 0 unless the routine detects an error (see
               Section 6).

          6. Error Indicators and Warnings

          Errors detected by the routine:

          If on entry IFAIL = 0 or -1, explanatory error messages are
          output on the current error message unit (as defined by X04AAF).

          IFAIL= 1
               On entry M < 1.

          IFAIL= 2
               On entry N < 1.

          IFAIL= 3
               On entry INIT is not one of 'I', 'S' or 'R'.

          IFAIL= 4
               On entry INIT = 'S', but C06FUF has not previously been
               called.

          IFAIL= 5
               On entry INIT = 'S' or 'R', but at least one of the arrays
               TRIGM and TRIGN is inconsistent with the current value of M
               or N.

          7. Accuracy

          Some indication of accuracy can be obtained by performing a
          subsequent inverse transform and comparing the results with the
          original sequence (in exact arithmetic they would be identical).

          8. Further Comments

          The time taken by the routine is approximately proportional to
          mn*log(mn), but also depends on the factorization of the
          individual dimensions m and n. The routine is somewhat faster
          than average if their only prime factors are 2, 3 or 5; and
          fastest of all if they are powers of 2.

          9. Example

          This program reads in a bivariate sequence of complex data values
          and prints the two-dimensional Fourier transform. It then
          performs an inverse transform and prints the sequence so
          obtained, which may be compared to the original data values.

          The example program is not reproduced here. The source code for
          all example programs is distributed with the NAG Foundation
          Library software and should be available on-line.

\end{verbatim}
\endscroll
\end{page}
\begin{page}{manpageXXc06gbf}{NAG On-line Documentation: c06gbf}
\beginscroll
\begin{verbatim}



     C06GBF(3NAG)      Foundation Library (12/10/92)      C06GBF(3NAG)



          C06 -- Summation of Series                                 C06GBF
                  C06GBF -- NAG Foundation Library Routine Document

          Note: Before using this routine, please read the Users' Note for
          your implementation to check implementation-dependent details.
          The symbol (*) after a NAG routine name denotes a routine that is
          not included in the Foundation Library.

          1. Purpose

          C06GBF forms the complex conjugate of a Hermitian sequence of n
          data values.

          2. Specification

                 SUBROUTINE C06GBF (X, N, IFAIL)
                 INTEGER          N, IFAIL
                 DOUBLE PRECISION X(N)

          3. Description

          This is a utility routine for use in conjunction with C06EAF,
          C06EBF, C06FAF(*) or C06FBF(*) to calculate inverse discrete
          Fourier transforms (see the Chapter Introduction).

          4. References

          None.

          5. Parameters

           1:  X(N) -- DOUBLE PRECISION array                  Input/Output
               On entry: if the data values z  are written as x +iy  and
                                             j                 j   j
               if X is declared with bounds (0:N-1) in the (sub)program
               from which C06GBF is called, then for 0<=j<=n/2, X(j) must
               contain x (=x   ), while for n/2<j<=n-1, X(j) must contain
                        j   n-j
               -y (=y   ). In other words, X must contain the Hermitian
                 j   n-j
               sequence in Hermitian form. (See also Section 2.1.2 of the
               Chapter Introduction). On exit: the imaginary parts y  are
                                                                    j
               negated. The real parts x  are not referenced.
                                        j

           2:  N -- INTEGER                                           Input
               On entry: the number of data values, n. Constraint: N >= 1.

           3:  IFAIL -- INTEGER                                Input/Output
               On entry: IFAIL must be set to 0, -1 or 1. For users not
               familiar with this parameter (described in the Essential
               Introduction) the recommended value is 0.

               On exit: IFAIL = 0 unless the routine detects an error (see
               Section 6).

          6. Error Indicators and Warnings

          Errors detected by the routine:

          IFAIL= 1
               N < 1.

          7. Accuracy

          Exact.

          8. Further Comments

          The time taken by the routine is negligible.

          9. Example

          This program reads in a sequence of real data values, calls
          C06EAF followed by C06GBF to compute their inverse discrete
          Fourier transform, and prints this after expanding it from
          Hermitian form into a full complex sequence.

          The example program is not reproduced here. The source code for
          all example programs is distributed with the NAG Foundation
          Library software and should be available on-line.

\end{verbatim}
\endscroll
\end{page}
\begin{page}{manpageXXc06gcf}{NAG On-line Documentation: c06gcf}
\beginscroll
\begin{verbatim}



     C06GCF(3NAG)      Foundation Library (12/10/92)      C06GCF(3NAG)



          C06 -- Summation of Series                                 C06GCF
                  C06GCF -- NAG Foundation Library Routine Document

          Note: Before using this routine, please read the Users' Note for
          your implementation to check implementation-dependent details.
          The symbol (*) after a NAG routine name denotes a routine that is
          not included in the Foundation Library.

          1. Purpose

          C06GCF forms the complex conjugate of a sequence of n data
          values.

          2. Specification

                 SUBROUTINE C06GCF (Y, N, IFAIL)
                 INTEGER          N, IFAIL
                 DOUBLE PRECISION Y(N)

          3. Description

          This is a utility routine for use in conjunction with C06ECF or
          C06FCF(*) to calculate inverse discrete Fourier transforms (see
          the Chapter Introduction).

          4. References

          None.

          5. Parameters

           1:  Y(N) -- DOUBLE PRECISION array                  Input/Output
               On entry: if Y is declared with bounds (0:N-1) in the (sub)
               program which C06GCF is called, then Y(j) must contain the
               imaginary part of the jth data value, for 0<=j<=n-1. On
               exit: these values are negated.

           2:  N -- INTEGER                                           Input
               On entry: the number of data values, n. Constraint: N >= 1.

           3:  IFAIL -- INTEGER                                Input/Output
               On entry: IFAIL must be set to 0, -1 or 1. For users not
               familiar with this parameter (described in the Essential
               Introduction) the recommended value is 0.

               On exit: IFAIL = 0 unless the routine detects an error (see
               Section 6).

          6. Error Indicators and Warnings

          Errors detected by the routine:

          IFAIL= 1
               N < 1.

          7. Accuracy

          Exact.

          8. Further Comments

          The time taken by the routine is negligible.

          9. Example

          This program reads in a sequence of complex data values and
          prints their inverse discrete Fourier transform as computed by
          calling C06GCF, followed by C06ECF and C06GCF again.

          The example program is not reproduced here. The source code for
          all example programs is distributed with the NAG Foundation
          Library software and should be available on-line.

\end{verbatim}
\endscroll
\end{page}
\begin{page}{manpageXXc06gqf}{NAG On-line Documentation: c06gqf}
\beginscroll
\begin{verbatim}



     C06GQF(3NAG)      Foundation Library (12/10/92)      C06GQF(3NAG)



          C06 -- Summation of Series                                 C06GQF
                  C06GQF -- NAG Foundation Library Routine Document

          Note: Before using this routine, please read the Users' Note for
          your implementation to check implementation-dependent details.
          The symbol (*) after a NAG routine name denotes a routine that is
          not included in the Foundation Library.

          1. Purpose

          C06GQF forms the complex conjugates of m Hermitian sequences,
          each containing n data values.

          2. Specification

                 SUBROUTINE C06GQF (M, N, X, IFAIL)
                 INTEGER          M, N, IFAIL
                 DOUBLE PRECISION X(M*N)

          3. Description

          This is a utility routine for use in conjunction with C06FPF and
          C06FQF to calculate inverse discrete Fourier transforms (see the
          Chapter Introduction).

          4. References

          None.

          5. Parameters

           1:  M -- INTEGER                                           Input
               On entry: the number of Hermitian sequences to be
               conjugated, m. Constraint: M >= 1.

           2:  N -- INTEGER                                           Input
               On entry: the number of data values in each Hermitian
               sequence, n. Constraint: N >= 1.

           3:  X(M,N) -- DOUBLE PRECISION array                Input/Output
               On entry: the data must be stored in array X as if in a
               two-dimensional array of dimension (1:M,0:N-1); each of the
               m sequences is stored in a row of the array in Hermitian
                                           p                 p   p
               form. If the n data values z  are written as x +iy , then
                                           j                 j   j
                               p
               for 0<=j<=n/2, x  is contained in X(p,j), and for 1<=j<=(n-
                               j
                      p
               1)/2, y  is contained in X(p,n-j). (See also Section 2.1.2
                      j
               of the Chapter Introduction.) On exit: the imaginary parts
                p                              p
               y  are negated. The real parts x  are not referenced.
                j                              j

           4:  IFAIL -- INTEGER                                Input/Output
               On entry: IFAIL must be set to 0, -1 or 1. For users not
               familiar with this parameter (described in the Essential
               Introduction) the recommended value is 0.

               On exit: IFAIL = 0 unless the routine detects an error (see
               Section 6).

          6. Error Indicators and Warnings

          Errors detected by the routine:

          If on entry IFAIL = 0 or -1, explanatory error messages are
          output on the current error message unit (as defined by X04AAF).

          IFAIL= 1
               On entry M < 1.

          IFAIL= 2
               On entry N < 1.

          7. Accuracy

          Exact.

          8. Further Comments

          None.

          9. Example

          This program reads in sequences of real data values which are
          assumed to be Hermitian sequences of complex data stored in
          Hermitian form. The sequences are expanded into full complex form
          using C06GSF and printed. The sequences are then conjugated
          (using C06GQF) and the conjugated sequences are expanded into
          complex form using C06GSF and printed out.

          The example program is not reproduced here. The source code for
          all example programs is distributed with the NAG Foundation
          Library software and should be available on-line.

\end{verbatim}
\endscroll
\end{page}
\begin{page}{manpageXXc06gsf}{NAG On-line Documentation: c06gsf}
\beginscroll
\begin{verbatim}



     C06GSF(3NAG)      Foundation Library (12/10/92)      C06GSF(3NAG)



          C06 -- Summation of Series                                 C06GSF
                  C06GSF -- NAG Foundation Library Routine Document

          Note: Before using this routine, please read the Users' Note for
          your implementation to check implementation-dependent details.
          The symbol (*) after a NAG routine name denotes a routine that is
          not included in the Foundation Library.

          1. Purpose

          C06GSF takes m Hermitian sequences, each containing n data
          values, and forms the real and imaginary parts of the m
          corresponding complex sequences.

          2. Specification

                 SUBROUTINE C06GSF (M, N, X, U, V, IFAIL)
                 INTEGER          M, N, IFAIL
                 DOUBLE PRECISION X(M*N), U(M*N), V(M*N)

          3. Description

          This is a utility routine for use in conjunction with C06FPF and
          C06FQF (see the Chapter Introduction).

          4. References

          None.

          5. Parameters

           1:  M -- INTEGER                                           Input
               On entry: the number of Hermitian sequences, m, to be
               converted into complex form. Constraint: M >= 1.

           2:  N -- INTEGER                                           Input
               On entry: the number of data values, n, in each sequence.
               Constraint: N >= 1.

           3:  X(M,N) -- DOUBLE PRECISION array                       Input
               On entry: the data must be stored in X as if in a two-
               dimensional array of dimension (1:M,0:N-1); each of the m
               sequences is stored in a row of the array in Hermitian form.
                                     p                 p   p
               If the n data values z  are written as x +iy , then for
                                     j                 j   j
                           p
               0<=j<=n/2, x  is contained in X(p,j), and for 1<=j<=(n-1)/2,
                           j
                p
               y  is contained in X(p,n-j). (See also Section 2.1.2 of the
                j
               Chapter Introduction.)

           4:  U(M,N) -- DOUBLE PRECISION array                      Output

           5:  V(M,N) -- DOUBLE PRECISION array                      Output
               On exit: the real and imaginary parts of the m sequences of
               length n, are stored in U and V respectively, as if in two-
               dimensional arrays of dimension (1:M,0:N-1); each of the m
               sequences is stored as if in a row of each array. In other
               words, if the real parts of the pth sequence are denoted by
                p
               x , for j=0,1,...,n-1 then the mn elements of the array U
                j
               contain the values
                     1  2      m   1  2      m       1    2        m
                    x ,x ,...,x , x ,x ,...,x ,..., x   ,x   ,...,x
                     0  0      0   1  1      1       n-1  n-1      n-1

           6:  IFAIL -- INTEGER                                Input/Output
               On entry: IFAIL must be set to 0, -1 or 1. For users not
               familiar with this parameter (described in the Essential
               Introduction) the recommended value is 0.

               On exit: IFAIL = 0 unless the routine detects an error (see
               Section 6).

          6. Error Indicators and Warnings

          Errors detected by the routine:

          If on entry IFAIL = 0 or -1, explanatory error messages are
          output on the current error message unit (as defined by X04AAF).

          IFAIL= 1
               On entry M < 1.

          IFAIL= 2
               On entry N < 1.

          7. Accuracy

          Exact.

          8. Further Comments

          None.

          9. Example

          This program reads in sequences of real data values which are
          assumed to be Hermitian sequences of complex data stored in
          Hermitian form. The sequences are then expanded into full complex
          form using C06GSF and printed.

          The example program is not reproduced here. The source code for
          all example programs is distributed with the NAG Foundation
          Library software and should be available on-line.

\end{verbatim}
\endscroll
\end{page}