subroutine cslz(a,x,r,m)
      integer m
      double complex a(m),x(m),r(m)
c
c     cslz solves the double complex linear system
c     a * x = b
c     with the c - matrix a .
c
c     on entry
c
c        a       double complex(m)
c                the first row of the c - matrix .
c                on return a is unaltered .
c
c        x       double complex(m)
c                the right hand side vector b .
c
c        r       double complex(m)
c                a work vector .
c
c        m       integer
c                the order of the matrix a .
c
c     on return
c
c        x       the solution vector .
c
c     toeplitz package. this version dated 07/23/82 .
c
c     subroutines and functions
c
c        fortran ... dcmplx,dcos,dfloat,dsin
c
c     internal variables
c
      integer i1,i2
      double precision p,ri,rm,v1,v2
      double complex e,e1,f,f1,t,t1
c
      t1 = x(1)
      x(1) = t1/a(1)
      if (m .eq. 1) go to 50
      rm = dfloat(m)
c
c     compute the inverse discrete fourier transformation
c     of the first row of the matrix and the discrete
c     fourier transformation of the right hand side vector .
c
      t = (0.0d0,0.0d0)
      ri = -1.0d0
      do 20 i1 = 1, m
         ri = ri + 1.0d0
c
c        minimize error in forming multiples of 2*pi .
c
         p = ((201.d0/32.d0)*ri + 1.93530717958647692528d-3*ri)/rm
c
         v1 = dcos(p)
         v2 = dsin(p)
         e = dcmplx(v1,-v2)
         e1 = dcmplx(v1,v2)
         f = a(1)
         f1 = t1
         do 10 i2 = 2, m
            f = e*f + a(i2)
            f1 = e1*f1 + x(i2)
   10    continue
         r(i1) = (e1*f1)/(e*f)
         t = t + r(i1)
   20 continue
c
c     compute the solution of the given system by
c     the inverse discrete fourier transformation .
c
      x(1) = t/rm
      ri = 0.0d0
      do 40 i1 = 2, m
         ri = ri + 1.0d0
c
c        minimize error in forming multiples of 2*pi .
c
         p = ((201.d0/32.d0)*ri + 1.93530717958647692528d-3*ri)/rm
c
         v1 = dcos(p)
         v2 = dsin(p)
         e = dcmplx(v1,-v2)
         f = r(1)
         do 30 i2 = 2, m
            f = e*f + r(i2)
   30    continue
         x(i1) = e*f/rm
   40 continue
   50 continue
      return
      end