subroutine cccslz(a,x,r,m,l,k,lda)

c*********************************************************************72
c
cc CCCSLZ solves a linear system with a CCC matrix.
c
c  Discussion:
c
c     cccslz solves the double complex linear system
c     a * x = b
c     with the ccc - matrix a .
c
c     on entry
c
c        a       double complex(m*l,k)
c                the first row of outer blocks of the ccc - matrix .
c                each outer block is represented by its first row
c                of inner blocks. each inner block is represented
c                by its first row. on return a has been destroyed .
c
c        x       double complex(m*l*k)
c                the right hand side vector b .
c
c        r       double complex(max(m,2*l,2*k))
c                a work vector .
c
c        m       integer
c                the order of the inner blocks of the matrix a .
c
c        l       integer
c                the number of inner blocks in a row or column
c                of an outer block of the matrix a .
c
c        k       integer
c                the number of outer blocks in a row or column
c                of the matrix a .
c
c        lda     integer
c                the leading dimension of the array a .
c
c     on return
c
c        x       the solution vector .
c
c     toeplitz package. this version dated 07/23/82 .
c
      integer m,l,k,lda
      double complex a(lda,k),x(m,l,k),r(1)
      integer i1,i2,i3,ml
      double precision rk
c
      rk = dfloat(k)
      ml = m*l
c
c     reduce the ccc - matrix to a block-diagonal matrix
c     by the inverse discrete fourier transformation .
c
      call salwz(a,r,r(k+1),ml,k,lda,-1)
c
c     compute the discrete fourier transformation of
c     the right hand side vector .
c
      call salwz(x,r,r(k+1),ml,k,ml,1)
c
c     solve the block-diagonal system, blocks of which
c     are cc - matrices .
c
      do 10 i3 = 1, k
         call ccslz(a(1,i3),x(1,1,i3),r,m,l,m)
   10 continue
c
c     compute the solution of the given system by
c     the inverse discrete fourier transformation .
c
      call salwz(x,r,r(k+1),ml,k,ml,-1)
c
      do 40 i3 = 1, k
         do 30 i2 = 1, l
            do 20 i1 = 1, m
               x(i1,i2,i3) = x(i1,i2,i3)/rk
   20       continue
   30    continue
   40 continue
      return
      end
      subroutine ccgslz(a,x,r,m,l,k,lda)

c*********************************************************************72
c
cc CCGSLZ solves a linear system with a CCG matrix.
c
c  Discussion:
c
c     ccgslz solves the double complex linear system
c     a * x = b
c     with the ccg - matrix a .
c
c     on entry
c
c        a       double complex(m**2*l,k)
c                the first row of outer blocks of the ccg - matrix .
c                each outer block is represented by its first row
c                of inner blocks. each inner block is represented
c                by columns. on return a has been destroyed .
c
c        x       double complex(m*l*k)
c                the right hand side vector b .
c
c        r       double complex(max(m,2*l,2*k))
c                a work vector .
c
c        m       integer
c                the order of the inner blocks of the matrix a .
c
c        l       integer
c                the number of inner blocks in a row or column
c                of an outer block of the matrix a .
c
c        k       integer
c                the number of outer blocks in a row or column
c                of the matrix a .
c
c        lda     integer
c                the leading dimension of the array a .
c
c     on return
c
c        x       the solution vector .
c
c     toeplitz package. this version dated 07/23/82 .
c
      integer m,l,k,lda
      double complex a(lda,k),x(m,l,k),r(1)
      integer i1,i2,i3,ml,mm
      double precision rk
c
      rk = dfloat(k)
      mm = m**2
      ml = m*l
c
c     reduce the ccg - matrix to a block-diagonal matrix
c     by the inverse discrete fourier transformation .
c
      call salwz(a,r,r(k+1),mm*l,k,lda,-1)
c
c     compute the discrete fourier transformation of
c     the right hand side vector .
c
      call salwz(x,r,r(k+1),ml,k,ml,1)
c
c     solve the block-diagonal system, blocks of which
c     are cg - matrices .
c
      do 10 i3 = 1, k
         call cgslz(a(1,i3),x(1,1,i3),r,m,l,mm)
   10 continue
c
c     compute the solution of the given system by
c     the inverse discrete fourier transformation .
c
      call salwz(x,r,r(k+1),ml,k,ml,-1)
c
      do 40 i3 = 1, k
         do 30 i2 = 1, l
            do 20 i1 = 1, m
               x(i1,i2,i3) = x(i1,i2,i3)/rk
   20       continue
   30    continue
   40 continue
      return
      end
      subroutine ccslz(a,x,r,m,l,lda)

c*********************************************************************72
c
cc CCSLZ solves a linear system with a CC matrix.
c
c  Discussion:
c
c     ccslz solves the double complex linear system
c     a * x = b
c     with the cc - matrix a .
c
c     on entry
c
c        a       double complex(m,l)
c                the first row of blocks of the cc - matrix .
c                each block is represented by its first row .
c                on return a has been destroyed .
c
c        x       double complex(m*l)
c                the right hand side vector b .
c
c        r       double complex(max(m,2*l))
c                a work vector .
c
c        m       integer
c                the order of the blocks of the matrix a .
c
c        l       integer
c                the number of blocks in a row or column
c                of the matrix a .
c
c        lda     integer
c                the leading dimension of the array a .
c
c     on return
c
c        x       the solution vector .
c
c     toeplitz package. this version dated 07/23/82 .
c
      integer m,l,lda
      double complex a(lda,l),x(m,l),r(1)
      integer i1,i2
      double precision rl
c
      rl = dfloat(l)
c
c     reduce the cc - matrix to a block-diagonal matrix
c     by the inverse discrete fourier transformation .
c
      call salwz(a,r,r(l+1),m,l,lda,-1)
c
c     compute the discrete fourier transformation of
c     the right hand side vector .
c
      call salwz(x,r,r(l+1),m,l,m,1)
c
c     solve the block-diagonal system, blocks of which
c     are c - matrices .
c
      do 10 i2 = 1, l
         call cslz(a(1,i2),x(1,i2),r,m)
   10 continue
c
c     compute the solution of the given system by
c     the inverse discrete fourier transformation .
c
      call salwz(x,r,r(l+1),m,l,m,-1)
c
      do 30 i2 = 1, l
         do 20 i1 = 1, m
            x(i1,i2) = x(i1,i2)/rl
   20    continue
   30 continue
      return
      end
      subroutine cctslz(a,x,r,m,l,k,lda)

c*********************************************************************72
c
cc CCTSLZ solves a linear system with a CCT matrix.
c
c  Discussion:
c
c     cctslz solves the double complex linear system
c     a * x = b
c     with the cct - matrix a .
c
c     on entry
c
c        a       double complex((2*m - 1)*l,k)
c                the first row of outer blocks of the cct - matrix .
c                each outer block is represented by its first row
c                of inner blocks. each inner block is represented
c                by its first row followed by its first column
c                beginning with the second element .
c                on return a has been destroyed .
c
c        x       double complex(m*l*k)
c                the right hand side vector b .
c
c        r       double complex(max(2*m - 2,2*l,2*k))
c                a work vector .
c
c        m       integer
c                the order of the inner blocks of the matrix a .
c
c        l       integer
c                the number of inner blocks in a row or column
c                of an outer block of the matrix a .
c
c        k       integer
c                the number of outer blocks in a row or column
c                of the matrix a .
c
c        lda     integer
c                the leading dimension of the array a .
c
c     on return
c
c        x       the solution vector .
c
c     toeplitz package. this version dated 07/23/82 .
c
      integer m,l,k,lda
      double complex a(lda,k),x(m,l,k),r(1)
      integer i1,i2,i3,m2,ml
      double precision rk
c
      rk = dfloat(k)
      m2 = 2*m - 1
      ml = m*l
c
c     reduce the cct - matrix to a block-diagonal matrix
c     by the inverse discrete fourier transformation .
c
      call salwz(a,r,r(k+1),m2*l,k,lda,-1)
c
c     compute the discrete fourier transformation of
c     the right hand side vector .
c
      call salwz(x,r,r(k+1),ml,k,ml,1)
c
c     solve the block-diagonal system, blocks of which
c     are ct - matrices .
c
      do 10 i3 = 1, k
         call ctslz(a(1,i3),x(1,1,i3),r,m,l,m2)
   10 continue
c
c     compute the solution of the given system by
c     the inverse discrete fourier transformation .
c
      call salwz(x,r,r(k+1),ml,k,ml,-1)
c
      do 40 i3 = 1, k
         do 30 i2 = 1, l
            do 20 i1 = 1, m
               x(i1,i2,i3) = x(i1,i2,i3)/rk
   20       continue
   30    continue
   40 continue
      return
      end
      subroutine cgslz(a,x,r,m,l,lda)

c*********************************************************************72
c
cc CGSLZ solves a linear system with a CG matrix.
c
c  Discussion:
c
c     cgslz solves the double complex linear system
c     a * x = b
c     with the cg - matrix a .
c
c     on entry
c
c        a       double complex(m**2,l)
c                the first row of blocks of the cg - matrix .
c                each block is represented by columns .
c                on return a has been destroyed .
c
c        x       double complex(m*l)
c                the right hand side vector b .
c
c        r       double complex(max(m,2*l))
c                a work vector .
c
c        m       integer
c                the order of the blocks of the matrix a .
c
c        l       integer
c                the number of blocks in a row or column
c                of the matrix a .
c
c        lda     integer
c                the leading dimension of the array a .
c
c     on return
c
c        x       the solution vector .
c
c     toeplitz package. this version dated 07/23/82 .
c
      integer m,l,lda
      double complex a(lda,l),x(m,l),r(1)
      integer i1,i2,ii
      double precision rl
c
      rl = dfloat(l)
c
c     reduce the cg - matrix to a block-diagonal matrix
c     by the inverse discrete fourier transformation .
c
      call salwz(a,r,r(l+1),m**2,l,lda,-1)
c
c     compute the discrete fourier transformation of
c     the right hand side vector .
c
      call salwz(x,r,r(l+1),m,l,m,1)
c
c     solve the block-diagonal system, blocks of which
c     are g - matrices .
c
      do 10 i2 = 1, l
         call zgefa(a(1,i2),m,m,r,ii)
         call zgesl(a(1,i2),m,m,r,x(1,i2),0)
   10 continue
c
c     compute the solution of the given system by
c     the inverse discrete fourier transformation .
c
      call salwz(x,r,r(l+1),m,l,m,-1)
c
      do 30 i2 = 1, l
         do 20 i1 = 1, m
            x(i1,i2) = x(i1,i2)/rl
   20    continue
   30 continue
      return
      end
      subroutine cqrd(a,q,s,m,l,ldq,lds)

c*********************************************************************72
c
cc CQRD computes the QR factorization of a CC matrix.
c
c  Discussion:
c
c     cqrd computes the qr factorization in the form
c     a * r(inverse) = q
c     of the double precision column-circulant matrix a .
c
c     on entry
c
c        a       double precision(m)
c                the first column of the column-circulant matrix .
c                on return a is unaltered .
c
c        m       integer
c                the number of rows of the matrices a and q .
c                m must be at least as large as l .
c
c        l       integer
c                the number of columns of the matrices a and q
c                and the order of the upper triangular matrix s .
c
c        ldq     integer
c                the leading dimension of the array q .
c
c        lds     integer
c                the leading dimension of the array s .
c
c     on return
c
c        q       double precision(m,l)
c                the q matrix of the factorization .
c                the columns of q are orthonormal .
c
c        s       double precision(l,l)
c                the inverse of the r matrix of the factorization .
c                elements below the main diagonal are not accessed .
c
c     toeplitz package. this version dated 07/23/82 .
c
      integer m,l,ldq,lds
      double precision a(m),q(ldq,l),s(lds,l)
      integer i,j,j1,ji
      double precision scale,dnrm2
      double precision c,ddot
c
c     initialization (last column of q used as work vector) .
c
      do 10 i = 1, m
         q(i,1) = a(i)
         q(i,l) = a(i)
   10 continue
c
c     recurrent process for the lattice algorithm with normalization .
c
      do 70 j1 = 1, l
         j = j1 + 1
         scale = 1.0d0/dnrm2(m,q(1,j1),1)
         if (j1 .eq. l) go to 60
            c = -scale*(q(m,j1)*q(1,l) +
     *           ddot(m-1,q(1,j1),1,q(2,l),1))/dnrm2(m,q(1,l),1)
            q(1,j) = q(m,j1) + c*q(1,l)
            do 20 i = 2, m
               q(i,j) = q(i-1,j1) + c*q(i,l)
   20       continue
            if (j .eq. l) go to 30
               q(1,l) = q(1,l) + c*q(m,j1)
               call daxpy(m-1,c,q(1,j1),1,q(2,l),1)
   30       continue
            s(1,j) = c
            if (j .eq. 2) go to 50
               do 40 i = 2, j1
                  ji = j - i
                  s(i,j) = s(i-1,j1) + c*s(ji,j1)
   40          continue
   50       continue
   60    continue
         call dscal(m,scale,q(1,j1),1)
         s(j1,j1) = 1.0d0
         call dscal(j1,scale,s(1,j1),1)
   70 continue
      return
      end
      subroutine cqrz(a,q,s,m,l,ldq,lds)

c*********************************************************************72
c
cc CQRZ computes the QR factorization of a CCC matrix.
c
c  Discussion:
c
c     cqrz computes the qr factorization in the form
c     a * r(inverse) = q
c     of the double complex column-circulant matrix a .
c
c     on entry
c
c        a       double complex(m)
c                the first column of the column-circulant matrix .
c                on return a is unaltered .
c
c        m       integer
c                the number of rows of the matrices a and q .
c                m must be at least as large as l .
c
c        l       integer
c                the number of columns of the matrices a and q
c                and the order of the upper triangular matrix s .
c
c        ldq     integer
c                the leading dimension of the array q .
c
c        lds     integer
c                the leading dimension of the array s .
c
c     on return
c
c        q       double complex(m,l)
c                the q matrix of the factorization .
c                the columns of q are orthonormal .
c
c        s       double complex(l,l)
c                the inverse of the r matrix of the factorization .
c                elements below the main diagonal are not accessed .
c
c     toeplitz package. this version dated 07/23/82 .
c
      integer m,l,ldq,lds
      double complex a(m),q(ldq,l),s(lds,l)
      integer i,j,j1,ji
      double precision scale,dznrm2
      double complex c,zdotc
c
c     initialization (last column of q used as work vector) .
c
      do 10 i = 1, m
         q(i,1) = a(i)
         q(i,l) = a(i)
   10 continue
c
c     recurrent process for the lattice algorithm with normalization .
c
      do 70 j1 = 1, l
         j = j1 + 1
         scale = 1.0d0/dznrm2(m,q(1,j1),1)
         if (j1 .eq. l) go to 60
            c = -scale*(dconjg(q(m,j1))*q(1,l) +
     *           zdotc(m-1,q(1,j1),1,q(2,l),1))/dznrm2(m,q(1,l),1)
            q(1,j) = q(m,j1) + c*q(1,l)
            do 20 i = 2, m
               q(i,j) = q(i-1,j1) + c*q(i,l)
   20       continue
            if (j .eq. l) go to 30
               q(1,l) = q(1,l) + c*q(m,j1)
               call zaxpy(m-1,c,q(1,j1),1,q(2,l),1)
   30       continue
            s(1,j) = c
            if (j .eq. 2) go to 50
               do 40 i = 2, j1
                  ji = j - i
                  s(i,j) = s(i-1,j1) + c*s(ji,j1)
   40          continue
   50       continue
   60    continue
         call zdscal(m,scale,q(1,j1),1)
         s(j1,j1) = (1.0d0,0.0d0)
         call zdscal(j1,scale,s(1,j1),1)
   70 continue
      return
      end
      subroutine cslz ( a, x, r, m )

c*********************************************************************72
c
cc CSLZ solves a linear system with a C matrix.
c
c  Discussion:
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
      integer m
      double complex a(m),x(m),r(m)

      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

         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
      subroutine ctgslz(a,x,r,m,l,k,lda)

c*********************************************************************72
c
cc CTGSLZ solves a linear system with a CTG matrix.
c
c  Discussion:
c
c     ctgslz solves the double complex linear system
c     a * x = b
c     with the ctg - matrix a .
c
c     on entry
c
c        a       double complex(m**2*(2*l - 1),k)
c                the first row of outer blocks of the ctg - matrix .
c                each outer block is represented by its first row
c                of inner blocks followed by its first column
c                of inner blocks beginning with the second block .
c                each inner block is represented by columns .
c                on return a has been destroyed .
c
c        x       double complex(m*l*k)
c                the right hand side vector b .
c
c        r       double complex(max(m**2*(2*l + 3) + m,2*k))
c                a work vector .
c
c        m       integer
c                the order of the inner blocks of the matrix a .
c
c        l       integer
c                the number of inner blocks in a row or column
c                of an outer block of the matrix a .
c
c        k       integer
c                the number of outer blocks in a row or column
c                of the matrix a .
c
c        lda     integer
c                the leading dimension of the array a .
c
c     on return
c
c        x       the solution vector .
c
c     toeplitz package. this version dated 07/23/82 .
c
      integer m,l,k,lda
      double complex a(lda,k),x(m,l,k),r(1)
      integer i1,i2,i3,ml,mm
      double precision rk
c
      rk = dfloat(k)
      mm = m**2
      ml = m*l
c
c     reduce the ctg - matrix to a block-diagonal matrix
c     by the inverse discrete fourier transformation .
c
      call salwz(a,r,r(k+1),mm*(2*l - 1),k,lda,-1)
c
c     compute the discrete fourier transformation of
c     the right hand side vector .
c
      call salwz(x,r,r(k+1),ml,k,ml,1)
c
c     solve the block-diagonal system, blocks of which
c     are tg - matrices .
c
      do 10 i3 = 1, k
         call tgslz(a(1,i3),x(1,1,i3),r,m,l,mm)
   10 continue
c
c     compute the solution of the given system by
c     the inverse discrete fourier transformation .
c
      call salwz(x,r,r(k+1),ml,k,ml,-1)
c
      do 40 i3 = 1, k
         do 30 i2 = 1, l
            do 20 i1 = 1, m
               x(i1,i2,i3) = x(i1,i2,i3)/rk
   20       continue
   30    continue
   40 continue
      return
      end
      subroutine ctslz(a,x,r,m,l,lda)

c*********************************************************************72
c
cc CTSLZ solves a linear system with a CT matrix.
c
c  Discussion:
c
c     ctslz solves the double complex linear system
c     a * x = b
c     with the ct - matrix a .
c
c     on entry
c
c        a       double complex(2*m - 1,l)
c                the first row of blocks of the ct - matrix .
c                each block is represented by its first row
c                followed by its first column beginning with the
c                second element. on return a has been destroyed .
c
c        x       double complex(m*l)
c                the right hand side vector b .
c
c        r       double complex(max(2*m - 2,2*l))
c                a work vector .
c
c        m       integer
c                the order of the blocks of the matrix a .
c
c        l       integer
c                the number of blocks in a row or column
c                of the matrix a .
c
c        lda     integer
c                the leading dimension of the array a .
c
c     on return
c
c        x       the solution vector .
c
c     toeplitz package. this version dated 07/23/82 .
c
      integer m,l,lda
      double complex a(lda,l),x(m,l),r(1)
      integer i1,i2
      double precision rl
c
      rl = dfloat(l)
c
c     reduce the ct - matrix to a block-diagonal matrix
c     by the inverse discrete fourier transformation .
c
      call salwz(a,r,r(l+1),2*m - 1,l,lda,-1)
c
c     compute the discrete fourier transformation of
c     the right hand side vector .
c
      call salwz(x,r,r(l+1),m,l,m,1)
c
c     solve the block-diagonal system, blocks of which
c     are t - matrices .
c
      do 10 i2 = 1, l
         call tslz(a(1,i2),x(1,i2),r,m)
   10 continue
c
c     compute the solution of the given system by
c     the inverse discrete fourier transformation .
c
      call salwz(x,r,r(l+1),m,l,m,-1)
c
      do 30 i2 = 1, l
         do 20 i1 = 1, m
            x(i1,i2) = x(i1,i2)/rl
   20    continue
   30 continue
      return
      end
      subroutine daxpy(n,da,dx,incx,dy,incy)

c*********************************************************************72
c
cc DAXPY adds a multiple of one vector to another.
c
c  Discussion:
c
c     constant times a vector plus a vector.
c     uses unrolled loops for increments equal to one.
c     jack dongarra, linpack, 3/11/78.
c
      double precision dx(1),dy(1),da
      integer i,incx,incy,ixiy,m,mp1,n
c
      if(n.le.0)return
      if (da .eq. 0.0d0) return
      if(incx.eq.1.and.incy.eq.1)go to 20
c
c        code for unequal increments or equal increments
c          not equal to 1
c
      ix = 1
      iy = 1
      if(incx.lt.0)ix = (-n+1)*incx + 1
      if(incy.lt.0)iy = (-n+1)*incy + 1
      do 10 i = 1,n
        dy(iy) = dy(iy) + da*dx(ix)
        ix = ix + incx
        iy = iy + incy
   10 continue
      return
c
c        code for both increments equal to 1
c
c
c        clean-up loop
c
   20 m = mod(n,4)
      if( m .eq. 0 ) go to 40
      do 30 i = 1,m
        dy(i) = dy(i) + da*dx(i)
   30 continue
      if( n .lt. 4 ) return
   40 mp1 = m + 1
      do 50 i = mp1,n,4
        dy(i) = dy(i) + da*dx(i)
        dy(i + 1) = dy(i + 1) + da*dx(i + 1)
        dy(i + 2) = dy(i + 2) + da*dx(i + 2)
        dy(i + 3) = dy(i + 3) + da*dx(i + 3)
   50 continue
      return
      end
      function dcabs1(z)

c*********************************************************************72
c
cc DCABS1 computes the L1 norm of a complex value.
c
c  Discussion:
c
      double precision dcabs1
      double complex z,zz
      double precision t(2)
      equivalence (zz,t(1))

      zz = z
      dcabs1 = dabs(t(1)) + dabs(t(2))

      return
      end
      function ddot(n,dx,incx,dy,incy)

c*********************************************************************72
c
cc DDOT forms the dot product of two vectors.
c
c  Discussion:
c
c     forms the dot product of two vectors.
c     uses unrolled loops for increments equal to one.
c     jack dongarra, linpack, 3/11/78.
c
      double precision ddot
      double precision dx(1),dy(1),dtemp
      integer i,incx,incy,ix,iy,m,mp1,n

      ddot = 0.0d0
      dtemp = 0.0d0
      if(n.le.0)return
      if(incx.eq.1.and.incy.eq.1)go to 20
c
c        code for unequal increments or equal increments
c          not equal to 1
c
      ix = 1
      iy = 1
      if(incx.lt.0)ix = (-n+1)*incx + 1
      if(incy.lt.0)iy = (-n+1)*incy + 1
      do 10 i = 1,n
        dtemp = dtemp + dx(ix)*dy(iy)
        ix = ix + incx
        iy = iy + incy
   10 continue
      ddot = dtemp
      return
c
c        code for both increments equal to 1
c
c
c        clean-up loop
c
   20 m = mod(n,5)
      if( m .eq. 0 ) go to 40
      do 30 i = 1,m
        dtemp = dtemp + dx(i)*dy(i)
   30 continue
      if( n .lt. 5 ) go to 60
   40 mp1 = m + 1
      do 50 i = mp1,n,5
        dtemp = dtemp + dx(i)*dy(i) + dx(i + 1)*dy(i + 1) +
     *   dx(i + 2)*dy(i + 2) + dx(i + 3)*dy(i + 3) + dx(i + 4)*dy(i + 4)
   50 continue
   60 ddot = dtemp
      return
      end
      subroutine dgefa(a,lda,n,ipvt,info)

c*********************************************************************72
c
cc DGEFA factors a general matrix.
c
c  Discussion:
c
c     dgefa factors a double precision matrix by gaussian elimination.
c
c     dgefa is usually called by dgeco, but it can be called
c     directly with a saving in time if  rcond  is not needed.
c     (time for dgeco) = (1 + 9/n)*(time for dgefa) .
c
c     on entry
c
c        a       double precision(lda, n)
c                the matrix to be factored.
c
c        lda     integer
c                the leading dimension of the array  a .
c
c        n       integer
c                the order of the matrix  a .
c
c     on return
c
c        a       an upper triangular matrix and the multipliers
c                which were used to obtain it.
c                the factorization can be written  a = l*u  where
c                l  is a product of permutation and unit lower
c                triangular matrices and  u  is upper triangular.
c
c        ipvt    integer(n)
c                an integer vector of pivot indices.
c
c        info    integer
c                = 0  normal value.
c                = k  if  u(k,k) .eq. 0.0 .  this is not an error
c                     condition for this subroutine, but it does
c                     indicate that dgesl or dgedi will divide by zero
c                     if called.  use  rcond  in dgeco for a reliable
c                     indication of singularity.
c
c     linpack. this version dated 08/14/78 .
c     cleve moler, university of new mexico, argonne national lab.
c
      integer lda,n,ipvt(1),info
      double precision a(lda,1)
      double precision t
      integer idamax,j,k,kp1,l,nm1
c
c     gaussian elimination with partial pivoting
c
      info = 0
      nm1 = n - 1
      if (nm1 .lt. 1) go to 70
      do 60 k = 1, nm1
         kp1 = k + 1
c
c        find l = pivot index
c
         l = idamax(n-k+1,a(k,k),1) + k - 1
         ipvt(k) = l
c
c        zero pivot implies this column already triangularized
c
         if (a(l,k) .eq. 0.0d0) go to 40
c
c           interchange if necessary
c
            if (l .eq. k) go to 10
               t = a(l,k)
               a(l,k) = a(k,k)
               a(k,k) = t
   10       continue
c
c           compute multipliers
c
            t = -1.0d0/a(k,k)
            call dscal(n-k,t,a(k+1,k),1)
c
c           row elimination with column indexing
c
            do 30 j = kp1, n
               t = a(l,j)
               if (l .eq. k) go to 20
                  a(l,j) = a(k,j)
                  a(k,j) = t
   20          continue
               call daxpy(n-k,t,a(k+1,k),1,a(k+1,j),1)
   30       continue
         go to 50
   40    continue
            info = k
   50    continue
   60 continue
   70 continue
      ipvt(n) = n
      if (a(n,n) .eq. 0.0d0) info = n
      return
      end
      subroutine dgesl(a,lda,n,ipvt,b,job)

c*********************************************************************72
c
cc DGESL solves a linear system with a factored general matrix.
c
c  Discussion:
c
c     dgesl solves the double precision system
c     a * x = b  or  trans(a) * x = b
c     using the factors computed by dgeco or dgefa.
c
c     on entry
c
c        a       double precision(lda, n)
c                the output from dgeco or dgefa.
c
c        lda     integer
c                the leading dimension of the array  a .
c
c        n       integer
c                the order of the matrix  a .
c
c        ipvt    integer(n)
c                the pivot vector from dgeco or dgefa.
c
c        b       double precision(n)
c                the right hand side vector.
c
c        job     integer
c                = 0         to solve  a*x = b ,
c                = nonzero   to solve  trans(a)*x = b  where
c                            trans(a)  is the transpose.
c
c     on return
c
c        b       the solution vector  x .
c
c     error condition
c
c        a division by zero will occur if the input factor contains a
c        zero on the diagonal.  technically this indicates singularity
c        but it is often caused by improper arguments or improper
c        setting of lda .  it will not occur if the subroutines are
c        called correctly and if dgeco has set rcond .gt. 0.0
c        or dgefa has set info .eq. 0 .
c
c     to compute  inverse(a) * c  where  c  is a matrix
c     with  p  columns
c           call dgeco(a,lda,n,ipvt,rcond,z)
c           if (rcond is too small) go to ...
c           do 10 j = 1, p
c              call dgesl(a,lda,n,ipvt,c(1,j),0)
c        10 continue
c
c     linpack. this version dated 08/14/78 .
c     cleve moler, university of new mexico, argonne national lab.
c
      integer lda,n,ipvt(1),job
      double precision a(lda,1),b(1)
      double precision ddot,t
      integer k,kb,l,nm1
c
      nm1 = n - 1
      if (job .ne. 0) go to 50
c
c        job = 0 , solve  a * x = b
c        first solve  l*y = b
c
         if (nm1 .lt. 1) go to 30
         do 20 k = 1, nm1
            l = ipvt(k)
            t = b(l)
            if (l .eq. k) go to 10
               b(l) = b(k)
               b(k) = t
   10       continue
            call daxpy(n-k,t,a(k+1,k),1,b(k+1),1)
   20    continue
   30    continue
c
c        now solve  u*x = y
c
         do 40 kb = 1, n
            k = n + 1 - kb
            b(k) = b(k)/a(k,k)
            t = -b(k)
            call daxpy(k-1,t,a(1,k),1,b(1),1)
   40    continue
      go to 100
   50 continue
c
c        job = nonzero, solve  trans(a) * x = b
c        first solve  trans(u)*y = b
c
         do 60 k = 1, n
            t = ddot(k-1,a(1,k),1,b(1),1)
            b(k) = (b(k) - t)/a(k,k)
   60    continue
c
c        now solve trans(l)*x = y
c
         if (nm1 .lt. 1) go to 90
         do 80 kb = 1, nm1
            k = n - kb
            b(k) = b(k) + ddot(n-k,a(k+1,k),1,b(k+1),1)
            l = ipvt(k)
            if (l .eq. k) go to 70
               t = b(l)
               b(l) = b(k)
               b(k) = t
   70       continue
   80    continue
   90    continue
  100 continue
      return
      end
      function dnrm2 ( n, dx, incx)

c*********************************************************************72
c
cc DNRM2 computes the L2 norm of a vector.
c
c  Discussion:
c
c     euclidean norm of the n-vector stored in dx() with storage
c     increment incx .
c     if    n .le. 0 return with result = 0.
c     if n .ge. 1 then incx must be .ge. 1
c
c           c.l.lawson, 1978 jan 08
c
c     four phase method     using two built-in constants that are
c     hopefully applicable to all machines.
c         cutlo = maximum of  dsqrt(u/eps)  over all known machines.
c         cuthi = minimum of  dsqrt(v)      over all known machines.
c     where
c         eps = smallest no. such that eps + 1. .gt. 1.
c         u   = smallest positive no.   (underflow limit)
c         v   = largest  no.            (overflow  limit)
c
c     brief outline of algorithm..
c
c     phase 1    scans zero components.
c     move to phase 2 when a component is nonzero and .le. cutlo
c     move to phase 3 when a component is .gt. cutlo
c     move to phase 4 when a component is .ge. cuthi/m
c     where m = n for x() real and m = 2*n for complex.
c
      double precision dnrm2
      integer          next
      double precision   dx(1), cutlo, cuthi, hitest, sum, xmax,zero,one
      data   zero, one /0.0d0, 1.0d0/
c
c     values for cutlo and cuthi..
c     from the environmental parameters listed in the imsl converter
c     document the limiting values are as follows..
c     cutlo, s.p.   u/eps = 2**(-102) for  honeywell.  close seconds are
c                   univac and dec at 2**(-103)
c                   thus cutlo = 2**(-51) = 4.44089e-16
c     cuthi, s.p.   v = 2**127 for univac, honeywell, and dec.
c                   thus cuthi = 2**(63.5) = 1.30438e19
c     cutlo, d.p.   u/eps = 2**(-67) for honeywell and dec.
c                   thus cutlo = 2**(-33.5) = 8.23181d-11
c     cuthi, d.p.   same as s.p.  cuthi = 1.30438d19
c     data cutlo, cuthi / 8.232d-11,  1.304d19 /
c     data cutlo, cuthi / 4.441e-16,  1.304e19 /
      data cutlo, cuthi / 8.232d-11,  1.304d19 /
c
      if(n .gt. 0) go to 10
         dnrm2  = zero
         go to 300
c
   10 assign 30 to next
      sum = zero
      nn = n * incx
c                                                 begin main loop
      i = 1
   20    go to next,(30, 50, 70, 110)
   30 if( dabs(dx(i)) .gt. cutlo) go to 85
      assign 50 to next
      xmax = zero
c
c                        phase 1.  sum is zero
c
   50 if( dx(i) .eq. zero) go to 200
      if( dabs(dx(i)) .gt. cutlo) go to 85
c
c                                prepare for phase 2.
      assign 70 to next
      go to 105
c
c                                prepare for phase 4.
c
  100 i = j
      assign 110 to next
      sum = (sum / dx(i)) / dx(i)
  105 xmax = dabs(dx(i))
      go to 115
c
c                   phase 2.  sum is small.
c                             scale to avoid destructive underflow.
c
   70 if( dabs(dx(i)) .gt. cutlo ) go to 75
c
c                     common code for phases 2 and 4.
c                     in phase 4 sum is large.  scale to avoid overflow.
c
  110 if( dabs(dx(i)) .le. xmax ) go to 115
         sum = one + sum * (xmax / dx(i))**2
         xmax = dabs(dx(i))
         go to 200
c
  115 sum = sum + (dx(i)/xmax)**2
      go to 200
c
c
c                  prepare for phase 3.
c
   75 sum = (sum * xmax) * xmax
c
c
c     for real or d.p. set hitest = cuthi/n
c     for complex      set hitest = cuthi/(2*n)
c
   85 hitest = cuthi/float( n )
c
c                   phase 3.  sum is mid-range.  no scaling.
c
      do 95 j =i,nn,incx
      if(dabs(dx(j)) .ge. hitest) go to 100
   95    sum = sum + dx(j)**2
      dnrm2 = dsqrt( sum )
      go to 300
c
  200 continue
      i = i + incx
      if ( i .le. nn ) go to 20
c
c              end of main loop.
c
c              compute square root and adjust for scaling.
c
      dnrm2 = xmax * dsqrt(sum)
  300 continue
      return
      end
      subroutine dscal(n,da,dx,incx)

c*********************************************************************72
c
cc DSCAL scales a vector by a constant.
c
c  Discussion:
c
c     scales a vector by a constant.
c     uses unrolled loops for increment equal to one.
c     jack dongarra, linpack, 3/11/78.
c
      double precision da,dx(1)
      integer i,incx,m,mp1,n,nincx
c
      if(n.le.0)return
      if(incx.eq.1)go to 20
c
c        code for increment not equal to 1
c
      nincx = n*incx
      do 10 i = 1,nincx,incx
        dx(i) = da*dx(i)
   10 continue
      return
c
c        code for increment equal to 1
c
c
c        clean-up loop
c
   20 m = mod(n,5)
      if( m .eq. 0 ) go to 40
      do 30 i = 1,m
        dx(i) = da*dx(i)
   30 continue
      if( n .lt. 5 ) return
   40 mp1 = m + 1
      do 50 i = mp1,n,5
        dx(i) = da*dx(i)
        dx(i + 1) = da*dx(i + 1)
        dx(i + 2) = da*dx(i + 2)
        dx(i + 3) = da*dx(i + 3)
        dx(i + 4) = da*dx(i + 4)
   50 continue
      return
      end
      subroutine dtrdi(t,ldt,n,det,job,info)

c*********************************************************************72
c
cc DTRDI computes the determinant and inverse of a triangular matrix.
c
c  Discussion:
c
c     dtrdi computes the determinant and inverse of a double precision
c     triangular matrix.
c
c     on entry
c
c        t       double precision(ldt,n)
c                t contains the triangular matrix. the zero
c                elements of the matrix are not referenced, and
c                the corresponding elements of the array can be
c                used to store other information.
c
c        ldt     integer
c                ldt is the leading dimension of the array t.
c
c        n       integer
c                n is the order of the system.
c
c        job     integer
c                = 010       no det, inverse of lower triangular.
c                = 011       no det, inverse of upper triangular.
c                = 100       det, no inverse.
c                = 110       det, inverse of lower triangular.
c                = 111       det, inverse of upper triangular.
c
c     on return
c
c        t       inverse of original matrix if requested.
c                otherwise unchanged.
c
c        det     double precision(2)
c                determinant of original matrix if requested.
c                otherwise not referenced.
c                determinant = det(1) * 10.0**det(2)
c                with  1.0 .le. dabs(det(1)) .lt. 10.0
c                or  det(1) .eq. 0.0 .
c
c        info    integer
c                info contains zero if the system is nonsingular
c                and the inverse is requested.
c                otherwise info contains the index of
c                a zero diagonal element of t.
c
c
c     linpack. this version dated 08/14/78 .
c     cleve moler, university of new mexico, argonne national lab.
c
      integer ldt,n,job,info
      double precision t(ldt,1),det(2)
      double precision temp
      double precision ten
      integer i,j,k,kb,km1,kp1
c
c     begin block permitting ...exits to 180
c
c        compute determinant
c
         if (job/100 .eq. 0) go to 70
            det(1) = 1.0d0
            det(2) = 0.0d0
            ten = 10.0d0
            do 50 i = 1, n
               det(1) = t(i,i)*det(1)
c           ...exit
               if (det(1) .eq. 0.0d0) go to 60
   10          if (dabs(det(1)) .ge. 1.0d0) go to 20
                  det(1) = ten*det(1)
                  det(2) = det(2) - 1.0d0
               go to 10
   20          continue
   30          if (dabs(det(1)) .lt. ten) go to 40
                  det(1) = det(1)/ten
                  det(2) = det(2) + 1.0d0
               go to 30
   40          continue
   50       continue
   60       continue
   70    continue
c
c        compute inverse of upper triangular
c
         if (mod(job/10,10) .eq. 0) go to 170
            if (mod(job,10) .eq. 0) go to 120
c              begin block permitting ...exits to 110
                  do 100 k = 1, n
                     info = k
c              ......exit
                     if (t(k,k) .eq. 0.0d0) go to 110
                     t(k,k) = 1.0d0/t(k,k)
                     temp = -t(k,k)
                     call dscal(k-1,temp,t(1,k),1)
                     kp1 = k + 1
                     if (n .lt. kp1) go to 90
                     do 80 j = kp1, n
                        temp = t(k,j)
                        t(k,j) = 0.0d0
                        call daxpy(k,temp,t(1,k),1,t(1,j),1)
   80                continue
   90                continue
  100             continue
                  info = 0
  110          continue
            go to 160
  120       continue
c
c              compute inverse of lower triangular
c
               do 150 kb = 1, n
                  k = n + 1 - kb
                  info = k
c     ............exit
                  if (t(k,k) .eq. 0.0d0) go to 180
                  t(k,k) = 1.0d0/t(k,k)
                  temp = -t(k,k)
                  if (k .ne. n) call dscal(n-k,temp,t(k+1,k),1)
                  km1 = k - 1
                  if (km1 .lt. 1) go to 140
                  do 130 j = 1, km1
                     temp = t(k,j)
                     t(k,j) = 0.0d0
                     call daxpy(n-k+1,temp,t(k,k),1,t(k,j),1)
  130             continue
  140             continue
  150          continue
               info = 0
  160       continue
  170    continue
  180 continue
      return
      end
      function dznrm2( n, zx, incx)

c*********************************************************************72
c
cc DZNRM2 computes the L2 norm of a vector.
c
c  Discussion:
c
c     unitary norm of the complex n-vector stored in zx() with storage
c     increment incx .
c     if    n .le. 0 return with result = 0.
c     if n .ge. 1 then incx must be .ge. 1
c
c           c.l.lawson , 1978 jan 08
c
c     four phase method     using two built-in constants that are
c     hopefully applicable to all machines.
c         cutlo = maximum of  sqrt(u/eps)  over all known machines.
c         cuthi = minimum of  sqrt(v)      over all known machines.
c     where
c         eps = smallest no. such that eps + 1. .gt. 1.
c         u   = smallest positive no.   (underflow limit)
c         v   = largest  no.            (overflow  limit)
c
c     brief outline of algorithm..
c
c     phase 1    scans zero components.
c     move to phase 2 when a component is nonzero and .le. cutlo
c     move to phase 3 when a component is .gt. cutlo
c     move to phase 4 when a component is .ge. cuthi/m
c     where m = n for x() real and m = 2*n for complex.
c
      double precision dznrm2
      logical imag, scale
      integer          next
      double precision cutlo, cuthi, hitest, sum, xmax, absx, zero, one
      double complex      zx(1)
      double precision dreal,dimag
      double complex zdumr,zdumi
      dreal(zdumr) = zdumr
      dimag(zdumi) = (0.0d0,-1.0d0)*zdumi
      data         zero, one /0.0d0, 1.0d0/
c
c     values for cutlo and cuthi..
c     from the environmental parameters listed in the imsl converter
c     document the limiting values are as follows..
c     cutlo, s.p.   u/eps = 2**(-102) for  honeywell.  close seconds are
c                   univac and dec at 2**(-103)
c                   thus cutlo = 2**(-51) = 4.44089e-16
c     cuthi, s.p.   v = 2**127 for univac, honeywell, and dec.
c                   thus cuthi = 2**(63.5) = 1.30438e19
c     cutlo, d.p.   u/eps = 2**(-67) for honeywell and dec.
c                   thus cutlo = 2**(-33.5) = 8.23181d-11
c     cuthi, d.p.   same as s.p.  cuthi = 1.30438d19
c     data cutlo, cuthi / 8.232d-11,  1.304d19 /
c     data cutlo, cuthi / 4.441e-16,  1.304e19 /
      data cutlo, cuthi / 8.232d-11,  1.304d19 /
c
      if(n .gt. 0) go to 10
         dznrm2  = zero
         go to 300
c
   10 assign 30 to next
      sum = zero
      nn = n * incx
c                                                 begin main loop
      do 210 i=1,nn,incx
         absx = dabs(dreal(zx(i)))
         imag = .false.
         go to next,(30, 50, 70, 90, 110)
   30 if( absx .gt. cutlo) go to 85
      assign 50 to next
      scale = .false.
c
c                        phase 1.  sum is zero
c
   50 if( absx .eq. zero) go to 200
      if( absx .gt. cutlo) go to 85
c
c                                prepare for phase 2.
      assign 70 to next
      go to 105
c
c                                prepare for phase 4.
c
  100 assign 110 to next
      sum = (sum / absx) / absx
  105 scale = .true.
      xmax = absx
      go to 115
c
c                   phase 2.  sum is small.
c                             scale to avoid destructive underflow.
c
   70 if( absx .gt. cutlo ) go to 75
c
c                     common code for phases 2 and 4.
c                     in phase 4 sum is large.  scale to avoid overflow.
c
  110 if( absx .le. xmax ) go to 115
         sum = one + sum * (xmax / absx)**2
         xmax = absx
         go to 200
c
  115 sum = sum + (absx/xmax)**2
      go to 200
c
c
c                  prepare for phase 3.
c
   75 sum = (sum * xmax) * xmax
c
   85 assign 90 to next
      scale = .false.
c
c     for real or d.p. set hitest = cuthi/n
c     for complex      set hitest = cuthi/(2*n)
c
      hitest = cuthi/float( n )
c
c                   phase 3.  sum is mid-range.  no scaling.
c
   90 if(absx .ge. hitest) go to 100
         sum = sum + absx**2
  200 continue
c                  control selection of real and imaginary parts.
c
      if(imag) go to 210
         absx = dabs(dimag(zx(i)))
         imag = .true.
      go to next,(  50, 70, 90, 110 )
c
  210 continue
c
c              end of main loop.
c              compute square root and adjust for scaling.
c
      dznrm2 = dsqrt(sum)
      if(scale) dznrm2 = dznrm2 * xmax
  300 continue
      return
      end
      function idamax(n,dx,incx)

c*********************************************************************72
c
cc IDAMAX locates the element of maximum absolute value in a vector.
c
c  Discussion:
c
c     finds the index of element having max. absolute value.
c     jack dongarra, linpack, 3/11/78.
c
      integer idamax
      double precision dx(1),dmax
      integer i,incx,ix,n
c
      idamax = 0
      if( n .lt. 1 ) return
      idamax = 1
      if(n.eq.1)return
      if(incx.eq.1)go to 20
c
c        code for increment not equal to 1
c
      ix = 1
      dmax = dabs(dx(1))
      ix = ix + incx
      do 10 i = 2,n
         if(dabs(dx(ix)).le.dmax) go to 5
         idamax = i
         dmax = dabs(dx(ix))
    5    ix = ix + incx
   10 continue
      return
c
c        code for increment equal to 1
c
   20 dmax = dabs(dx(1))
      do 30 i = 2,n
         if(dabs(dx(i)).le.dmax) go to 30
         idamax = i
         dmax = dabs(dx(i))
   30 continue
      return
      end
      function izamax(n,zx,incx)

c*********************************************************************72
c
cc IZAMAX locates the element of maximum absolute value in a vector.
c
c  Discussion:
c
c     finds the index of element having max. absolute value.
c     jack dongarra, 3/11/78.
c
      integer izamax
      double complex zx(1)
      double precision smax
      double complex zdum
      double precision dcabs1
c
      izamax = 1
      if(n.le.1)return
      if(incx.eq.1)go to 20
c
c        code for increments not equal to 1
c
      ix = 1
      if(incx.lt.0)ix = (-n+1)*incx + 1
      smax = dcabs1(zx(1))
      ix = ix + incx
      do 10 i = 2,n
         if(dcabs1(zx(ix)).le.smax) go to 5
         izamax = ix
         smax = dcabs1(zx(ix))
    5    ix = ix + incx
   10 continue
      return
c
c        code for increments equal to 1
c
   20 smax = dcabs1(zx(1))
      do 30 i = 2,n
         if(dcabs1(zx(i)).le.smax) go to 30
         izamax = i
         smax = dcabs1(zx(i))
   30 continue
      return
      end
      subroutine salwz(a,r1,r2,m,l,lda,job)

c*********************************************************************72
c
cc SALWZ: direct or inverse discrete Fourier transform of rectangular matrix.
c
c  Discussion:
c
c     salwz computes the direct or inverse discrete fourier
c     transformation for rows of a double complex rectangular matrix .
c
c     on entry
c
c        a       double complex(m,l)
c                the input matrix .
c
c        r1      double complex(l)
c                a work vector .
c
c        r2      double complex(l)
c                a work vector .
c
c        m       integer
c                the number of rows of the matrix a .
c
c        l       integer
c                the number of columns of the matrix a .
c
c        lda     integer
c                the leading dimension of the array a .
c
c        job     integer
c                = 1 for direct fourier transformation .
c                = -1 for inverse fourier transformation .
c
c     on return
c
c        a       the transformed rows of the matrix .
c
c     toeplitz package. this version dated 07/23/82 .
c
      integer m,l,lda,job
      double complex a(lda,l),r1(l),r2(l)
      integer i,i1,i2
      double precision p,ri,rl,v1,v2
      double complex e,f
c
      if (l .eq. 1) go to 60
      rl = dfloat(l)
c
      r1(1) = (1.0d0,0.0d0)
      ri = 0.0d0
      do 10 i1 = 2, l
         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)/rl
c
         v1 = dcos(p)
         v2 = dsin(p)
         if (job .eq. (-1)) v2 = -v2
         r1(i1) = dcmplx(v1,v2)
   10 continue
      do 50 i = 1, m
         do 30 i1 = 1, l
            e = r1(i1)
            f = a(i,1)
            do 20 i2 = 2, l
               f = e*f + a(i,i2)
   20       continue
            r2(i1) = e*f
   30    continue
         do 40 i1 = 1, l
            a(i,i1) = r2(i1)
   40    continue
   50 continue
   60 continue
      return
      end
      subroutine tgsld1(a1,a2,b,x,c1,c2,r1,r2,r3,r5,r6,r,m,l,lda)

c*********************************************************************72
c
cc TGSLD1 solves a linear system with a TG matrix.
c
c  Discussion:
c
c     tgsld1 solves the double precision linear system
c     a * x = b
c     with the tg - matrix a .
c
c     on entry
c
c        a1      double precision(m**2,l)
c                the first row of blocks of the tg - matrix a .
c                each block is represented by columns .
c                on return a1 is unaltered .
c
c        a2      double precision(m**2,l - 1)
c                the first column of blocks of the tg - matrix a
c                beginning with the second block. each block is
c                represented by columns. on return a2 is unaltered .
c
c        b       double precision(m*l)
c                the right hand side vector .
c                on return b is unaltered .
c
c        c1      double precision(m,m,l - 1)
c                a work array .
c
c        c2      double precision(m,m,l - 1)
c                a work array .
c
c        r1      double precision(m,m)
c                a work array .
c
c        r2      double precision(m,m)
c                a work array .
c
c        r3      double precision(m,m)
c                a work array .
c
c        r5      double precision(m,m)
c                a work array .
c
c        r6      double precision(m,m)
c                a work array .
c
c        r       double precision(m)
c                a work vector .
c
c        m       integer
c                the order of the blocks of the matrix a .
c
c        l       integer
c                the number of blocks in a row or column
c                of the matrix a .
c
c        lda     integer
c                the leading dimension of the array a .
c
c     on return
c
c        x       double precision(m*l)
c                the solution vector. x may coincide with b .
c
c     toeplitz package. this version dated 07/23/82 .
c
      integer m,l,lda
      double precision a1(lda,l),a2(lda,1),b(m,l),x(m,l),c1(m,m,1),
     *     c2(m,m,1),r1(m,m),r2(m,m),r3(m,m),r5(m,m),r6(m,m),r(m)
      integer i,i1,i2,i3,ii,j,n,n1,n2
c
c     solve the system with the principal minor of order m .
c
      i3 = 1
      do 20 j = 1, m
         do 10 i = 1, m
            c1(i,j,1) = a1(i3,1)
            r1(i,j) = a1(i3,1)
            r3(i,j) = r1(i,j)
            i3 = i3 + 1
   10    continue
         x(j,1) = b(j,1)
   20 continue
      call dgefa(r3,m,m,r,ii)
      call dgesl(r3,m,m,r,x(1,1),0)
      if (l .eq. 1) go to 420
c
c     recurrent process for solving the system
c     with the tg - matrix for n = 2, l .
c
      do 410 n = 2, l
c
c        compute multiples of the first and last block columns of
c        the inverse of the principal minor of order m*n .
c
         n1 = n - 1
         n2 = n - 2
         i3 = 1
         do 40 j = 1, m
            do 30 i = 1, m
               r5(i,j) = a2(i3,n1)
               r6(i,j) = a1(i3,n)
               i3 = i3 + 1
   30       continue
   40    continue
         if (n .eq. 2) go to 100
            do 60 j = 1, m
               do 50 i = 1, m
                  c1(i,j,n1) = r2(i,j)
   50          continue
   60       continue
            do 90 i1 = 1, n2
               i2 = n - i1
               do 80 j = 1, m
                  i3 = 1
                  do 70 i = 1, m
c for in-line daxpy, activate next 5 lines and deactivate following 3 .
c                    do 65 ii = 1, m
c                       r5(ii,j) = r5(ii,j) + c1(i,j,i2)*a2(i3,i1)
c                       r6(ii,j) = r6(ii,j) + c2(i,j,i1)*a1(i3,i1+1)
c                       i3 = i3 + 1
c  65                continue
                     call daxpy(m,c1(i,j,i2),a2(i3,i1),1,r5(1,j),1)
                     call daxpy(m,c2(i,j,i1),a1(i3,i1+1),1,r6(1,j),1)
                     i3 = i3 + m
   70             continue
   80          continue
   90       continue
  100    continue
         do 120 j = 1, m
            do 110 i = 1, m
               r2(i,j) = -r5(i,j)
  110       continue
            call dgesl(r3,m,m,r,r2(1,j),0)
  120    continue
         do 140 j = 1, m
            do 130 i = 1, m
               r3(i,j) = r6(i,j)
               r6(i,j) = -c1(i,j,1)
  130       continue
  140    continue
         do 160 j = 1, m
            do 150 i = 1, m
c for in-line daxpy, activate next 3 lines and deactivate following 1 .
c              do 145 ii = 1, m
c                 c1(ii,j,1) = c1(ii,j,1) + r2(i,j)*r3(ii,i)
c 145          continue
               call daxpy(m,r2(i,j),r3(1,i),1,c1(1,j,1),1)
  150       continue
  160    continue
         call dgefa(r6,m,m,r,ii)
         do 180 j = 1, m
            call dgesl(r6,m,m,r,r3(1,j),0)
            do 170 i = 1, m
c for in-line daxpy, activate next 3 lines and deactivate following 1 .
c              do 165 ii = 1, m
c                 r1(ii,j) = r1(ii,j) + r3(i,j)*r5(ii,i)
c 165          continue
               call daxpy(m,r3(i,j),r5(1,i),1,r1(1,j),1)
  170       continue
  180    continue
         if (n .eq. 2) go to 320
            do 200 j = 1, m
               do 190 i = 1, m
                  r6(i,j) = c2(i,j,1)
  190          continue
  200       continue
            do 310 i1 = 2, n1
               if (i1 .eq. n1) go to 230
                  do 220 j = 1, m
                     do 210 i = 1, m
                        r5(i,j) = c2(i,j,i1)
  210                continue
  220             continue
  230          continue
               do 260 j = 1, m
                  do 240 i = 1, m
                     c2(i,j,i1) = r6(i,j)
  240             continue
                  do 250 i = 1, m
c for in-line daxpy, activate next 3 lines and deactivate following 1 .
c                    do 245 ii = 1, m
c                       c2(ii,j,i1) = c2(ii,j,i1) + r3(i,j)*c1(ii,i,i1)
c 245                continue
                     call daxpy(m,r3(i,j),c1(1,i,i1),1,c2(1,j,i1),1)
  250             continue
  260          continue
               do 280 j = 1, m
                  do 270 i = 1, m
c for in-line daxpy, activate next 3 lines and deactivate following 1 .
c                    do 265 ii = 1, m
c                       c1(ii,j,i1) = c1(ii,j,i1) + r2(i,j)*r6(ii,i)
c 265                continue
                     call daxpy(m,r2(i,j),r6(1,i),1,c1(1,j,i1),1)
  270             continue
  280          continue
               do 300 j = 1, m
                  do 290 i = 1, m
                     r6(i,j) = r5(i,j)
  290             continue
  300          continue
  310       continue
  320    continue
         do 340 j = 1, m
            do 330 i = 1, m
               c2(i,j,1) = r3(i,j)
  330       continue
  340    continue
c
c        compute the solution of the system with the
c        principal minor of order m*n .
c
         do 360 j = 1, m
            do 350 i = 1, m
               r3(i,j) = r1(i,j)
  350       continue
            x(j,n) = b(j,n)
  360    continue
         do 380 i1 = 1, n1
            i2 = n - i1
            i3 = 1
            do 370 i = 1, m
c for in-line daxpy, activate next 4 lines and deactivate following 2 .
c              do 365 ii = 1, m
c                 x(ii,n) = x(ii,n) - x(i,i2)*a2(i3,i1)
c                 i3 = i3 + 1
c 365          continue
               call daxpy(m,-x(i,i2),a2(i3,i1),1,x(1,n),1)
               i3 = i3 + m
  370       continue
  380    continue
         call dgefa(r3,m,m,r,ii)
         call dgesl(r3,m,m,r,x(1,n),0)
         do 400 i1 = 1, n1
            do 390 i = 1, m
c for in-line daxpy, activate next 3 lines and deactivate following 1 .
c              do 385 ii = 1, m
c                 x(ii,i1) = x(ii,i1) + x(i,n)*c2(ii,i,i1)
c 385          continue
               call daxpy(m,x(i,n),c2(1,i,i1),1,x(1,i1),1)
  390       continue
  400    continue
  410 continue
  420 continue
      return
      end
      subroutine tgsld(a,x,r,m,l,lda)

c*********************************************************************72
c
cc TGSLD solves a linear system with a TG matrix.
c
c  Discussion:
c
c     tgsld calls tgsld1 to solve the double precision linear system
c     a * x = b
c     with the tg - matrix a .
c
c     on entry
c
c        a       double precision(m**2,2*l - 1)
c                the first row of blocks of the tg - matrix
c                followed by its first column of blocks beginning
c                with the second block. each block is represented
c                by columns. on return a is unaltered .
c
c        x       double precision(m*l)
c                the right hand side vector b .
c
c        r       double precision(m**2*(2*l + 3) + m)
c                a work vector .
c
c        m       integer
c                the order of the blocks of the matrix a .
c
c        l       integer
c                the number of blocks in a row or column
c                of the matrix a .
c
c        lda     integer
c                the leading dimension of the array a .
c
c     on return
c
c        x       the solution vector .
c
c     toeplitz package. this version dated 07/23/82 .
c
      integer m,l,lda
      double precision a(lda,1),x(m,l),r(1)
      integer mm,mml,mml1,mml2,mml3,mml4,mml5,mml6
c
c     call subroutine tgsld1
c
      mm = m**2
      mml = mm*(l - 1) + 1
      mml1 = 2*mml - 1
      mml2 = mml1 + mm
      mml3 = mml2 + mm
      mml4 = mml3 + mm
      mml5 = mml4 + mm
      mml6 = mml5 + mm
c
      call tgsld1(a,a(1,l+1),x,x,r,r(mml),r(mml1),r(mml2),
     *            r(mml3),r(mml4),r(mml5),r(mml6),m,l,lda)
c
      return
      end
      subroutine tgslz1(a1,a2,b,x,c1,c2,r1,r2,r3,r5,r6,r,m,l,lda)

c*********************************************************************72
c
cc TGSLZ1 solves a linear system with a TG matrix.
c
      integer m,l,lda
      double complex a1(lda,l),a2(lda,1),b(m,l),x(m,l),c1(m,m,1),
     *     c2(m,m,1),r1(m,m),r2(m,m),r3(m,m),r5(m,m),r6(m,m),r(m)
c
c     tgslz1 solves the double complex linear system
c     a * x = b
c     with the tg - matrix a .
c
c     on entry
c
c        a1      double complex(m**2,l)
c                the first row of blocks of the tg - matrix a .
c                each block is represented by columns .
c                on return a1 is unaltered .
c
c        a2      double complex(m**2,l - 1)
c                the first column of blocks of the tg - matrix a
c                beginning with the second block. each block is
c                represented by columns. on return a2 is unaltered .
c
c        b       double complex(m*l)
c                the right hand side vector .
c                on return b is unaltered .
c
c        c1      double complex(m,m,l - 1)
c                a work array .
c
c        c2      double complex(m,m,l - 1)
c                a work array .
c
c        r1      double complex(m,m)
c                a work array .
c
c        r2      double complex(m,m)
c                a work array .
c
c        r3      double complex(m,m)
c                a work array .
c
c        r5      double complex(m,m)
c                a work array .
c
c        r6      double complex(m,m)
c                a work array .
c
c        r       double complex(m)
c                a work vector .
c
c        m       integer
c                the order of the blocks of the matrix a .
c
c        l       integer
c                the number of blocks in a row or column
c                of the matrix a .
c
c        lda     integer
c                the leading dimension of the array a .
c
c     on return
c
c        x       double complex(m*l)
c                the solution vector. x may coincide with b .
c
c     toeplitz package. this version dated 07/23/82 .
c
c     subroutines and functions
c
c        linpack ... zaxpy,zgefa,zgesl
c                ... (for in-line zaxpy, see directions in comments)
c
c     internal variables
c
      integer i,i1,i2,i3,ii,j,n,n1,n2
c
c     solve the system with the principal minor of order m .
c
      i3 = 1
      do 20 j = 1, m
         do 10 i = 1, m
            c1(i,j,1) = a1(i3,1)
            r1(i,j) = a1(i3,1)
            r3(i,j) = r1(i,j)
            i3 = i3 + 1
   10    continue
         x(j,1) = b(j,1)
   20 continue
      call zgefa(r3,m,m,r,ii)
      call zgesl(r3,m,m,r,x(1,1),0)
      if (l .eq. 1) go to 420
c
c     recurrent process for solving the system
c     with the tg - matrix for n = 2, l .
c
      do 410 n = 2, l
c
c        compute multiples of the first and last block columns of
c        the inverse of the principal minor of order m*n .
c
         n1 = n - 1
         n2 = n - 2
         i3 = 1
         do 40 j = 1, m
            do 30 i = 1, m
               r5(i,j) = a2(i3,n1)
               r6(i,j) = a1(i3,n)
               i3 = i3 + 1
   30       continue
   40    continue
         if (n .eq. 2) go to 100
            do 60 j = 1, m
               do 50 i = 1, m
                  c1(i,j,n1) = r2(i,j)
   50          continue
   60       continue
            do 90 i1 = 1, n2
               i2 = n - i1
               do 80 j = 1, m
                  i3 = 1
                  do 70 i = 1, m
c for in-line zaxpy, activate next 5 lines and deactivate following 3 .
c                    do 65 ii = 1, m
c                       r5(ii,j) = r5(ii,j) + c1(i,j,i2)*a2(i3,i1)
c                       r6(ii,j) = r6(ii,j) + c2(i,j,i1)*a1(i3,i1+1)
c                       i3 = i3 + 1
c  65                continue
                     call zaxpy(m,c1(i,j,i2),a2(i3,i1),1,r5(1,j),1)
                     call zaxpy(m,c2(i,j,i1),a1(i3,i1+1),1,r6(1,j),1)
                     i3 = i3 + m
   70             continue
   80          continue
   90       continue
  100    continue
         do 120 j = 1, m
            do 110 i = 1, m
               r2(i,j) = -r5(i,j)
  110       continue
            call zgesl(r3,m,m,r,r2(1,j),0)
  120    continue
         do 140 j = 1, m
            do 130 i = 1, m
               r3(i,j) = r6(i,j)
               r6(i,j) = -c1(i,j,1)
  130       continue
  140    continue
         do 160 j = 1, m
            do 150 i = 1, m
c for in-line zaxpy, activate next 3 lines and deactivate following 1 .
c              do 145 ii = 1, m
c                 c1(ii,j,1) = c1(ii,j,1) + r2(i,j)*r3(ii,i)
c 145          continue
               call zaxpy(m,r2(i,j),r3(1,i),1,c1(1,j,1),1)
  150       continue
  160    continue
         call zgefa(r6,m,m,r,ii)
         do 180 j = 1, m
            call zgesl(r6,m,m,r,r3(1,j),0)
            do 170 i = 1, m
c for in-line zaxpy, activate next 3 lines and deactivate following 1 .
c              do 165 ii = 1, m
c                 r1(ii,j) = r1(ii,j) + r3(i,j)*r5(ii,i)
c 165          continue
               call zaxpy(m,r3(i,j),r5(1,i),1,r1(1,j),1)
  170       continue
  180    continue
         if (n .eq. 2) go to 320
            do 200 j = 1, m
               do 190 i = 1, m
                  r6(i,j) = c2(i,j,1)
  190          continue
  200       continue
            do 310 i1 = 2, n1
               if (i1 .eq. n1) go to 230
                  do 220 j = 1, m
                     do 210 i = 1, m
                        r5(i,j) = c2(i,j,i1)
  210                continue
  220             continue
  230          continue
               do 260 j = 1, m
                  do 240 i = 1, m
                     c2(i,j,i1) = r6(i,j)
  240             continue
                  do 250 i = 1, m
c for in-line zaxpy, activate next 3 lines and deactivate following 1 .
c                    do 245 ii = 1, m
c                       c2(ii,j,i1) = c2(ii,j,i1) + r3(i,j)*c1(ii,i,i1)
c 245                continue
                     call zaxpy(m,r3(i,j),c1(1,i,i1),1,c2(1,j,i1),1)
  250             continue
  260          continue
               do 280 j = 1, m
                  do 270 i = 1, m
c for in-line zaxpy, activate next 3 lines and deactivate following 1 .
c                    do 265 ii = 1, m
c                       c1(ii,j,i1) = c1(ii,j,i1) + r2(i,j)*r6(ii,i)
c 265                continue
                     call zaxpy(m,r2(i,j),r6(1,i),1,c1(1,j,i1),1)
  270             continue
  280          continue
               do 300 j = 1, m
                  do 290 i = 1, m
                     r6(i,j) = r5(i,j)
  290             continue
  300          continue
  310       continue
  320    continue
         do 340 j = 1, m
            do 330 i = 1, m
               c2(i,j,1) = r3(i,j)
  330       continue
  340    continue
c
c        compute the solution of the system with the
c        principal minor of order m*n .
c
         do 360 j = 1, m
            do 350 i = 1, m
               r3(i,j) = r1(i,j)
  350       continue
            x(j,n) = b(j,n)
  360    continue
         do 380 i1 = 1, n1
            i2 = n - i1
            i3 = 1
            do 370 i = 1, m
c for in-line zaxpy, activate next 4 lines and deactivate following 2 .
c              do 365 ii = 1, m
c                 x(ii,n) = x(ii,n) - x(i,i2)*a2(i3,i1)
c                 i3 = i3 + 1
c 365          continue
               call zaxpy(m,-x(i,i2),a2(i3,i1),1,x(1,n),1)
               i3 = i3 + m
  370       continue
  380    continue
         call zgefa(r3,m,m,r,ii)
         call zgesl(r3,m,m,r,x(1,n),0)
         do 400 i1 = 1, n1
            do 390 i = 1, m
c for in-line zaxpy, activate next 3 lines and deactivate following 1 .
c              do 385 ii = 1, m
c                 x(ii,i1) = x(ii,i1) + x(i,n)*c2(ii,i,i1)
c 385          continue
               call zaxpy(m,x(i,n),c2(1,i,i1),1,x(1,i1),1)
  390       continue
  400    continue
  410 continue
  420 continue
      return
      end
      subroutine tgslz(a,x,r,m,l,lda)

c*********************************************************************72
c
cc TGSLZ solves a linear system with a TG matrix.
c
      integer m,l,lda
      double complex a(lda,1),x(m,l),r(1)
c
c     tgslz calls tgslz1 to solve the double complex linear system
c     a * x = b
c     with the tg - matrix a .
c
c     on entry
c
c        a       double complex(m**2,2*l - 1)
c                the first row of blocks of the tg - matrix
c                followed by its first column of blocks beginning
c                with the second block. each block is represented
c                by columns. on return a is unaltered .
c
c        x       double complex(m*l)
c                the right hand side vector b .
c
c        r       double complex(m**2*(2*l + 3) + m)
c                a work vector .
c
c        m       integer
c                the order of the blocks of the matrix a .
c
c        l       integer
c                the number of blocks in a row or column
c                of the matrix a .
c
c        lda     integer
c                the leading dimension of the array 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        toeplitz package ... tgslz1
c
c     internal variables
c
      integer mm,mml,mml1,mml2,mml3,mml4,mml5,mml6
c
c     call subroutine tgslz1
c
      mm = m**2
      mml = mm*(l - 1) + 1
      mml1 = 2*mml - 1
      mml2 = mml1 + mm
      mml3 = mml2 + mm
      mml4 = mml3 + mm
      mml5 = mml4 + mm
      mml6 = mml5 + mm
c
      call tgslz1(a,a(1,l+1),x,x,r,r(mml),r(mml1),r(mml2),
     *            r(mml3),r(mml4),r(mml5),r(mml6),m,l,lda)
c
      return
      end
      subroutine timestamp ( )

c*********************************************************************72
c
cc TIMESTAMP prints out the current YMDHMS date as a timestamp.
c
c  Licensing:
c
c    This code is distributed under the GNU LGPL license.
c
c  Modified:
c
c    12 January 2007
c
c  Author:
c
c    John Burkardt
c
c  Parameters:
c
c    None
c
      implicit none

      character * ( 8 ) ampm
      integer d
      character * ( 8 ) date
      integer h
      integer m
      integer mm
      character * ( 9 ) month(12)
      integer n
      integer s
      character * ( 10 ) time
      integer y

      save month

      data month /
     &  'January  ', 'February ', 'March    ', 'April    ', 
     &  'May      ', 'June     ', 'July     ', 'August   ', 
     &  'September', 'October  ', 'November ', 'December ' /

      call date_and_time ( date, time )

      read ( date, '(i4,i2,i2)' ) y, m, d
      read ( time, '(i2,i2,i2,1x,i3)' ) h, n, s, mm

      if ( h .lt. 12 ) then
        ampm = 'AM'
      else if ( h .eq. 12 ) then
        if ( n .eq. 0 .and. s .eq. 0 ) then
          ampm = 'Noon'
        else
          ampm = 'PM'
        end if
      else
        h = h - 12
        if ( h .lt. 12 ) then
          ampm = 'PM'
        else if ( h .eq. 12 ) then
          if ( n .eq. 0 .and. s .eq. 0 ) then
            ampm = 'Midnight'
          else
            ampm = 'AM'
          end if
        end if
      end if

      write ( *, 
     &  '(i2,1x,a,1x,i4,2x,i2,a1,i2.2,a1,i2.2,a1,i3.3,1x,a)' ) 
     &  d, month(m), y, h, ':', n, ':', s, '.', mm, ampm

      return
      end
      subroutine tsld1(a1,a2,b,x,c1,c2,m)

c*********************************************************************72
c
cc TSLD1 solves a linear system with a T matrix.
c
      integer m
      double precision a1(m),a2(1),b(m),x(m),c1(1),c2(1)
c
c     tsld1 solves the double precision linear system
c     a * x = b
c     with the t - matrix a .
c
c     on entry
c
c        a1      double precision(m)
c                the first row of the t - matrix a .
c                on return a1 is unaltered .
c
c        a2      double precision(m - 1)
c                the first column of the t - matrix a
c                beginning with the second element .
c                on return a2 is unaltered .
c
c        b       double precision(m)
c                the right hand side vector .
c                on return b is unaltered .
c
c        c1      double precision(m - 1)
c                a work vector .
c
c        c2      double precision(m - 1)
c                a work vector .
c
c        m       integer
c                the order of the matrix a .
c
c     on return
c
c        x       double precision(m)
c                the solution vector. x may coincide with b .
c
c     toeplitz package. this version dated 07/23/82 .
c
c     internal variables
c
      integer i1,i2,n,n1,n2
      double precision r1,r2,r3,r5,r6
c
c     solve the system with the principal minor of order 1 .
c
      r1 = a1(1)
      x(1) = b(1)/r1
      if (m .eq. 1) go to 80
c
c     recurrent process for solving the system
c     with the t - matrix for n = 2, m .
c
      do 70 n = 2, m
c
c        compute multiples of the first and last columns of
c        the inverse of the principal minor of order n .
c
         n1 = n - 1
         n2 = n - 2
         r5 = a2(n1)
         r6 = a1(n)
         if (n .eq. 2) go to 20
            c1(n1) = r2
            do 10 i1 = 1, n2
               i2 = n - i1
               r5 = r5 + a2(i1)*c1(i2)
               r6 = r6 + a1(i1+1)*c2(i1)
   10       continue
   20    continue
         r2 = -r5/r1
         r3 = -r6/r1
         r1 = r1 + r5*r3
         if (n .eq. 2) go to 40
            r6 = c2(1)
            c2(n1) = 0.0d0
            do 30 i1 = 2, n1
               r5 = c2(i1)
               c2(i1) = c1(i1)*r3 + r6
               c1(i1) = c1(i1) + r6*r2
               r6 = r5
   30       continue
   40    continue
         c2(1) = r3
c
c        compute the solution of the system with the
c        principal minor of order n .
c
         r5 = 0.0d0
         do 50 i1 = 1, n1
            i2 = n - i1
            r5 = r5 + a2(i1)*x(i2)
   50    continue
         r6 = (b(n) - r5)/r1
         do 60 i1 = 1, n1
            x(i1) = x(i1) + c2(i1)*r6
   60    continue
         x(n) = r6
   70 continue
   80 continue
      return
      end
      subroutine tsld(a,x,r,m)

c*********************************************************************72
c
cc TSLD solves a linear system with a T matrix.
c
c  Discussion:
c
c     tsld calls tsld1 to solve the double precision linear system
c     a * x = b
c     with the t - matrix a .
c
c     on entry
c
c        a       double precision(2*m - 1)
c                the first row of the t - matrix followed by its
c                first column beginning with the second element .
c                on return a is unaltered .
c
c        x       double precision(m)
c                the right hand side vector b .
c
c        r       double precision(2*m - 2)
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
      integer m

      double precision a(1),x(m),r(1)

      call tsld1 ( a, a(m+1), x, x, r, r(m), m )

      return
      end
      subroutine tslz1(a1,a2,b,x,c1,c2,m)

c*********************************************************************72
c
cc TSLZ1 solves a linear system with a T matrix.
c
      integer m
      double complex a1(m),a2(1),b(m),x(m),c1(1),c2(1)
c
c     tslz1 solves the double complex linear system
c     a * x = b
c     with the t - matrix a .
c
c     on entry
c
c        a1      double complex(m)
c                the first row of the t - matrix a .
c                on return a1 is unaltered .
c
c        a2      double complex(m - 1)
c                the first column of the t - matrix a
c                beginning with the second element .
c                on return a2 is unaltered .
c
c        b       double complex(m)
c                the right hand side vector .
c                on return b is unaltered .
c
c        c1      double complex(m - 1)
c                a work vector .
c
c        c2      double complex(m - 1)
c                a work vector .
c
c        m       integer
c                the order of the matrix a .
c
c     on return
c
c        x       double complex(m)
c                the solution vector. x may coincide with b .
c
c     toeplitz package. this version dated 07/23/82 .
c
c     internal variables
c
      integer i1,i2,n,n1,n2
      double complex r1,r2,r3,r5,r6
c
c     solve the system with the principal minor of order 1 .
c
      r1 = a1(1)
      x(1) = b(1)/r1
      if (m .eq. 1) go to 80
c
c     recurrent process for solving the system
c     with the t - matrix for n = 2, m .
c
      do 70 n = 2, m
c
c        compute multiples of the first and last columns of
c        the inverse of the principal minor of order n .
c
         n1 = n - 1
         n2 = n - 2
         r5 = a2(n1)
         r6 = a1(n)
         if (n .eq. 2) go to 20
            c1(n1) = r2
            do 10 i1 = 1, n2
               i2 = n - i1
               r5 = r5 + a2(i1)*c1(i2)
               r6 = r6 + a1(i1+1)*c2(i1)
   10       continue
   20    continue
         r2 = -r5/r1
         r3 = -r6/r1
         r1 = r1 + r5*r3
         if (n .eq. 2) go to 40
            r6 = c2(1)
            c2(n1) = (0.0d0,0.0d0)
            do 30 i1 = 2, n1
               r5 = c2(i1)
               c2(i1) = c1(i1)*r3 + r6
               c1(i1) = c1(i1) + r6*r2
               r6 = r5
   30       continue
   40    continue
         c2(1) = r3
c
c        compute the solution of the system with the
c        principal minor of order n .
c
         r5 = (0.0d0,0.0d0)
         do 50 i1 = 1, n1
            i2 = n - i1
            r5 = r5 + a2(i1)*x(i2)
   50    continue
         r6 = (b(n) - r5)/r1
         do 60 i1 = 1, n1
            x(i1) = x(i1) + c2(i1)*r6
   60    continue
         x(n) = r6
   70 continue
   80 continue
      return
      end
      subroutine tslz(a,x,r,m)

c*********************************************************************72
c
cc TSLZ solves a linear system with a T matrix.
c
c  Discussion:
c
c    tslz calls tslz1 to solve the double complex linear system
c      a * x = b
c    with the t - matrix a.
c
c     on entry
c
c        a       double complex(2*m - 1)
c                the first row of the t - matrix followed by its
c                first column beginning with the second element .
c                on return a is unaltered .
c
c        x       double complex(m)
c                the right hand side vector b .
c
c        r       double complex(2*m - 2)
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        toeplitz package ... tslz1
c
c     call subroutine tslz1
c
      integer m
      double complex a(1),x(m),r(1)

      call tslz1(a,a(m+1),x,x,r,r(m),m)

      return
      end
      subroutine zaxpy(n,za,zx,incx,zy,incy)

c*********************************************************************72
c
cc ZAXPY adds a multiple of one vector to another.
c
c     constant times a vector plus a vector.
c     jack dongarra, 3/11/78.
c
      double complex zx(1),zy(1),za
      double precision dcabs1
      if(n.le.0)return
      if (dcabs1(za) .eq. 0.0d0) return
      if (incx.eq.1.and.incy.eq.1)go to 20
c
c        code for unequal increments or equal increments
c          not equal to 1
c
      ix = 1
      iy = 1
      if(incx.lt.0)ix = (-n+1)*incx + 1
      if(incy.lt.0)iy = (-n+1)*incy + 1
      do 10 i = 1,n
        zy(iy) = zy(iy) + za*zx(ix)
        ix = ix + incx
        iy = iy + incy
   10 continue
      return
c
c        code for both increments equal to 1
c
   20 do 30 i = 1,n
        zy(i) = zy(i) + za*zx(i)
   30 continue
      return
      end
      function zdotc(n,zx,incx,zy,incy)

c*********************************************************************72
c
cc ZDOTC computes the dot product of two vectors.
c
c     forms the dot product of a vector.
c     jack dongarra, 3/11/78.
c
      double complex zdotc
      double complex zx(1),zy(1),ztemp
      ztemp = (0.0d0,0.0d0)
      zdotc = (0.0d0,0.0d0)
      if(n.le.0)return
      if(incx.eq.1.and.incy.eq.1)go to 20
c
c        code for unequal increments or equal increments
c          not equal to 1
c
      ix = 1
      iy = 1
      if(incx.lt.0)ix = (-n+1)*incx + 1
      if(incy.lt.0)iy = (-n+1)*incy + 1
      do 10 i = 1,n
        ztemp = ztemp + dconjg(zx(ix))*zy(iy)
        ix = ix + incx
        iy = iy + incy
   10 continue
      zdotc = ztemp
      return
c
c        code for both increments equal to 1
c
   20 do 30 i = 1,n
        ztemp = ztemp + dconjg(zx(i))*zy(i)
   30 continue
      zdotc = ztemp
      return
      end
      subroutine zdscal(n,da,zx,incx)

c*********************************************************************72
c
cc ZDSCAL scales a vector by a constant.
c
c     scales a vector by a constant.
c     jack dongarra, 3/11/78.
c
      double complex zx(1)
      double precision da
      if(n.le.0)return
      if(incx.eq.1)go to 20
c
c        code for increments not equal to 1
c
      ix = 1
      if(incx.lt.0)ix = (-n+1)*incx + 1
      do 10 i = 1,n
        zx(ix) = dcmplx(da,0.0d0)*zx(ix)
        ix = ix + incx
   10 continue
      return
c
c        code for increments equal to 1
c
   20 do 30 i = 1,n
        zx(i) = dcmplx(da,0.0d0)*zx(i)
   30 continue
      return
      end
      subroutine zgefa(a,lda,n,ipvt,info)

c*********************************************************************72
c
cc ZGEFA factors a general matrix.
c
      integer lda,n,ipvt(1),info
      complex*16 a(lda,1)
c
c     zgefa factors a complex*16 matrix by gaussian elimination.
c
c     zgefa is usually called by zgeco, but it can be called
c     directly with a saving in time if  rcond  is not needed.
c     (time for zgeco) = (1 + 9/n)*(time for zgefa) .
c
c     on entry
c
c        a       complex*16(lda, n)
c                the matrix to be factored.
c
c        lda     integer
c                the leading dimension of the array  a .
c
c        n       integer
c                the order of the matrix  a .
c
c     on return
c
c        a       an upper triangular matrix and the multipliers
c                which were used to obtain it.
c                the factorization can be written  a = l*u  where
c                l  is a product of permutation and unit lower
c                triangular matrices and  u  is upper triangular.
c
c        ipvt    integer(n)
c                an integer vector of pivot indices.
c
c        info    integer
c                = 0  normal value.
c                = k  if  u(k,k) .eq. 0.0 .  this is not an error
c                     condition for this subroutine, but it does
c                     indicate that zgesl or zgedi will divide by zero
c                     if called.  use  rcond  in zgeco for a reliable
c                     indication of singularity.
c
c     linpack. this version dated 08/14/78 .
c     cleve moler, university of new mexico, argonne national lab.
c
c     subroutines and functions
c
c     blas zaxpy,zscal,izamax
c     fortran dabs
c
c     internal variables
c
      complex*16 t
      integer izamax,j,k,kp1,l,nm1
c
      complex*16 zdum
      double precision cabs1
      double precision dreal,dimag
      complex*16 zdumr,zdumi
      dreal(zdumr) = zdumr
      dimag(zdumi) = (0.0d0,-1.0d0)*zdumi
      cabs1(zdum) = dabs(dreal(zdum)) + dabs(dimag(zdum))
c
c     gaussian elimination with partial pivoting
c
      info = 0
      nm1 = n - 1
      if (nm1 .lt. 1) go to 70
      do 60 k = 1, nm1
         kp1 = k + 1
c
c        find l = pivot index
c
         l = izamax(n-k+1,a(k,k),1) + k - 1
         ipvt(k) = l
c
c        zero pivot implies this column already triangularized
c
         if (cabs1(a(l,k)) .eq. 0.0d0) go to 40
c
c           interchange if necessary
c
            if (l .eq. k) go to 10
               t = a(l,k)
               a(l,k) = a(k,k)
               a(k,k) = t
   10       continue
c
c           compute multipliers
c
            t = -(1.0d0,0.0d0)/a(k,k)
            call zscal(n-k,t,a(k+1,k),1)
c
c           row elimination with column indexing
c
            do 30 j = kp1, n
               t = a(l,j)
               if (l .eq. k) go to 20
                  a(l,j) = a(k,j)
                  a(k,j) = t
   20          continue
               call zaxpy(n-k,t,a(k+1,k),1,a(k+1,j),1)
   30       continue
         go to 50
   40    continue
            info = k
   50    continue
   60 continue
   70 continue
      ipvt(n) = n
      if (cabs1(a(n,n)) .eq. 0.0d0) info = n
      return
      end
      subroutine zgesl(a,lda,n,ipvt,b,job)

c*********************************************************************72
c
cc ZGESL solves a linear system involving a factored general matrix.
c
      integer lda,n,ipvt(1),job
      complex*16 a(lda,1),b(1)
c
c     zgesl solves the complex*16 system
c     a * x = b  or  ctrans(a) * x = b
c     using the factors computed by zgeco or zgefa.
c
c     on entry
c
c        a       complex*16(lda, n)
c                the output from zgeco or zgefa.
c
c        lda     integer
c                the leading dimension of the array  a .
c
c        n       integer
c                the order of the matrix  a .
c
c        ipvt    integer(n)
c                the pivot vector from zgeco or zgefa.
c
c        b       complex*16(n)
c                the right hand side vector.
c
c        job     integer
c                = 0         to solve  a*x = b ,
c                = nonzero   to solve  ctrans(a)*x = b  where
c                            ctrans(a)  is the conjugate transpose.
c
c     on return
c
c        b       the solution vector  x .
c
c     error condition
c
c        a division by zero will occur if the input factor contains a
c        zero on the diagonal.  technically this indicates singularity
c        but it is often caused by improper arguments or improper
c        setting of lda .  it will not occur if the subroutines are
c        called correctly and if zgeco has set rcond .gt. 0.0
c        or zgefa has set info .eq. 0 .
c
c     to compute  inverse(a) * c  where  c  is a matrix
c     with  p  columns
c           call zgeco(a,lda,n,ipvt,rcond,z)
c           if (rcond is too small) go to ...
c           do 10 j = 1, p
c              call zgesl(a,lda,n,ipvt,c(1,j),0)
c        10 continue
c
c     linpack. this version dated 08/14/78 .
c     cleve moler, university of new mexico, argonne national lab.
c
c     subroutines and functions
c
c     blas zaxpy,zdotc
c     fortran dconjg
c
c     internal variables
c
      complex*16 zdotc,t
      integer k,kb,l,nm1
      double precision dreal,dimag
      complex*16 zdumr,zdumi
      dreal(zdumr) = zdumr
      dimag(zdumi) = (0.0d0,-1.0d0)*zdumi
c
      nm1 = n - 1
      if (job .ne. 0) go to 50
c
c        job = 0 , solve  a * x = b
c        first solve  l*y = b
c
         if (nm1 .lt. 1) go to 30
         do 20 k = 1, nm1
            l = ipvt(k)
            t = b(l)
            if (l .eq. k) go to 10
               b(l) = b(k)
               b(k) = t
   10       continue
            call zaxpy(n-k,t,a(k+1,k),1,b(k+1),1)
   20    continue
   30    continue
c
c        now solve  u*x = y
c
         do 40 kb = 1, n
            k = n + 1 - kb
            b(k) = b(k)/a(k,k)
            t = -b(k)
            call zaxpy(k-1,t,a(1,k),1,b(1),1)
   40    continue
      go to 100
   50 continue
c
c        job = nonzero, solve  ctrans(a) * x = b
c        first solve  ctrans(u)*y = b
c
         do 60 k = 1, n
            t = zdotc(k-1,a(1,k),1,b(1),1)
            b(k) = (b(k) - t)/dconjg(a(k,k))
   60    continue
c
c        now solve ctrans(l)*x = y
c
         if (nm1 .lt. 1) go to 90
         do 80 kb = 1, nm1
            k = n - kb
            b(k) = b(k) + zdotc(n-k,a(k+1,k),1,b(k+1),1)
            l = ipvt(k)
            if (l .eq. k) go to 70
               t = b(l)
               b(l) = b(k)
               b(k) = t
   70       continue
   80    continue
   90    continue
  100 continue
      return
      end
      subroutine zscal(n,za,zx,incx)

c*********************************************************************72
c
cc ZSCAL scales a vector by a constant.
c
c    scales a vector by a constant.
c    jack dongarra, 3/11/78.
c
      double complex za,zx(1)
      if(n.le.0)return
      if(incx.eq.1)go to 20
c
c        code for increments not equal to 1
c
      ix = 1
      if(incx.lt.0)ix = (-n+1)*incx + 1
      do 10 i = 1,n
        zx(ix) = za*zx(ix)
        ix = ix + incx
   10 continue
      return
c
c        code for increments equal to 1
c
   20 do 30 i = 1,n
        zx(i) = za*zx(i)
   30 continue
      return
      end
      subroutine ztrdi(t,ldt,n,det,job,info)

c*********************************************************************72
c
cc ZTRDI computes the determinant and inverse of a triangular matrix.
c
      integer ldt,n,job,info
      complex*16 t(ldt,1),det(2)
c
c     ztrdi computes the determinant and inverse of a complex*16
c     triangular matrix.
c
c     on entry
c
c        t       complex*16(ldt,n)
c                t contains the triangular matrix. the zero
c                elements of the matrix are not referenced, and
c                the corresponding elements of the array can be
c                used to store other information.
c
c        ldt     integer
c                ldt is the leading dimension of the array t.
c
c        n       integer
c                n is the order of the system.
c
c        job     integer
c                = 010       no det, inverse of lower triangular.
c                = 011       no det, inverse of upper triangular.
c                = 100       det, no inverse.
c                = 110       det, inverse of lower triangular.
c                = 111       det, inverse of upper triangular.
c
c     on return
c
c        t       inverse of original matrix if requested.
c                otherwise unchanged.
c
c        det     complex*16(2)
c                determinant of original matrix if requested.
c                otherwise not referenced.
c                determinant = det(1) * 10.0**det(2)
c                with  1.0 .le. cabs1(det(1)) .lt. 10.0
c                or  det(1) .eq. 0.0 .
c
c        info    integer
c                info contains zero if the system is nonsingular
c                and the inverse is requested.
c                otherwise info contains the index of
c                a zero diagonal element of t.
c
c
c     linpack. this version dated 08/14/78 .
c     cleve moler, university of new mexico, argonne national lab.
c
c     subroutines and functions
c
c     blas zaxpy,zscal
c     fortran dabs,dcmplx,mod
c
c     internal variables
c
      complex*16 temp
      double precision ten
      integer i,j,k,kb,km1,kp1
c
      complex*16 zdum
      double precision cabs1
      double precision dreal,dimag
      complex*16 zdumr,zdumi
      dreal(zdumr) = zdumr
      dimag(zdumi) = (0.0d0,-1.0d0)*zdumi
      cabs1(zdum) = dabs(dreal(zdum)) + dabs(dimag(zdum))
c     begin block permitting ...exits to 180
c
c        compute determinant
c
         if (job/100 .eq. 0) go to 70
            det(1) = (1.0d0,0.0d0)
            det(2) = (0.0d0,0.0d0)
            ten = 10.0d0
            do 50 i = 1, n
               det(1) = t(i,i)*det(1)
c           ...exit
               if (cabs1(det(1)) .eq. 0.0d0) go to 60
   10          if (cabs1(det(1)) .ge. 1.0d0) go to 20
                  det(1) = dcmplx(ten,0.0d0)*det(1)
                  det(2) = det(2) - (1.0d0,0.0d0)
               go to 10
   20          continue
   30          if (cabs1(det(1)) .lt. ten) go to 40
                  det(1) = det(1)/dcmplx(ten,0.0d0)
                  det(2) = det(2) + (1.0d0,0.0d0)
               go to 30
   40          continue
   50       continue
   60       continue
   70    continue
c
c        compute inverse of upper triangular
c
         if (mod(job/10,10) .eq. 0) go to 170
            if (mod(job,10) .eq. 0) go to 120
c              begin block permitting ...exits to 110
                  do 100 k = 1, n
                     info = k
c              ......exit
                     if (cabs1(t(k,k)) .eq. 0.0d0) go to 110
                     t(k,k) = (1.0d0,0.0d0)/t(k,k)
                     temp = -t(k,k)
                     call zscal(k-1,temp,t(1,k),1)
                     kp1 = k + 1
                     if (n .lt. kp1) go to 90
                     do 80 j = kp1, n
                        temp = t(k,j)
                        t(k,j) = (0.0d0,0.0d0)
                        call zaxpy(k,temp,t(1,k),1,t(1,j),1)
   80                continue
   90                continue
  100             continue
                  info = 0
  110          continue
            go to 160
  120       continue
c
c              compute inverse of lower triangular
c
               do 150 kb = 1, n
                  k = n + 1 - kb
                  info = k
c     ............exit
                  if (cabs1(t(k,k)) .eq. 0.0d0) go to 180
                  t(k,k) = (1.0d0,0.0d0)/t(k,k)
                  temp = -t(k,k)
                  if (k .ne. n) call zscal(n-k,temp,t(k+1,k),1)
                  km1 = k - 1
                  if (km1 .lt. 1) go to 140
                  do 130 j = 1, km1
                     temp = t(k,j)
                     t(k,j) = (0.0d0,0.0d0)
                     call zaxpy(n-k+1,temp,t(k,k),1,t(k,j),1)
  130             continue
  140             continue
  150          continue
               info = 0
  160       continue
  170    continue
  180 continue
      return
      end