subroutine cbind(x,ldx,m,nx,y,ldy,ny,res,ldres)
c===============================================================c
c     written by Vassilis Hajivassiliou                         c
c     with the aid of Yoosoon Chang.                            c
c     version 1.0, March 22, 1992                               c
c     version 1.1, July 11, 1992                                c
c===============================================================c
c
c-----This routine takes two input matrices x and y with the following
c     dimensions and performs a horizontal concatenation on them.
c
c      x(ldx,m,nx), where ldx is the leading dimension
c      y(ldy,m,ny), where ldy is the leading dimension
c
c     The input matrices x and y must have the same number of rows.
c     The number of columns in the resulting matrix res is the sum
c     of nx and ny.
c     The horizontally concatenated matrix res(ldres,m,nx+ny) is
c     then returned.
c
c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
      implicit double precision (a-h,o-z)
      implicit integer*4 (i-n)
      dimension x(ldx,nx), y(ldy,ny),res(ldres,nx+ny)
c
      do 1 j = 1, nx
         do 2 i = 1,m
            res(i,j) = x(i,j)
 2       continue
 1    continue
c
      do 3 j = nx+1, nx+ny
         do 4 i =1,m
            res(i,j) = y(i,j-nx)
 4       continue
 3    continue
c
      return
      end
      subroutine chol(x,ldx,n,chx,ldchx)
c===============================================================c
c     written by Vassilis Hajivassiliou                         c
c     with the aid of Yoosoon Chang.                            c
c     version 1.0, March 22, 1992                               c
c     version 1.1, July 11, 1992                                c
c===============================================================c
c
c     This routine is based on the (double precision) subroutine DCHDC
c     in DLINP.
c
c-----chol computes the decomposition of a positive matrix x with ldx
c     leading dimension and n effective dimension, and returns the
c     lower triagular cholesky factor of x.
c
c     on entry ------ x(ldx, n), a positive matrix
c     on exit  ------ x:   the upper cholesky factor of x
c                     chx: the lower cholesky factor of x
c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
c
      implicit double precision (a-h, o-z)
      implicit integer*4 (i-n)
      dimension x(ldx,n),chx(ldchx,n)
      dimension work(1,1),jpvt(2000)
c
       call dchdc(x,ldx,n,work,jpvt,job,info)
c
c-----x, returned by the linpack subroutine contains, in its upper
c     half, the cholesky factor of the input matrix x.
c     The following loops removes the retained values of the input
c     matrix x in its lower half.
c
      do 100 j = 1,n
         do 200 i = j+1,n
            x(i,j) = 0.0d0
 200     continue
 100  continue
c
c-----transpose x to get the lower triangular cholesky factor.
c
      call transps(x,ldx,n,n,chx,ldchx)
      return
      end
      subroutine diag(x,ldx,iorder,res,ldres)
c===============================================================c
c     written by Vassilis Hajivassiliou                         c
c     with the aid of Yoosoon Chang.                            c
c     version 1.0, March 22, 1992                               c
c     version 1.1, July 11, 1992                                c
c===============================================================c
c
c-----diag extracts diagnol elements from a (m by n) matrix x, and put
c     it into a vector res, which has iorder elements.
c     iorder = min(m,n)
c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
c
      implicit double precision (a-h,o-z)
      dimension x(ldx,iorder), res(ldres)
c
      do 1 j = 1,iorder
         do 2 i = 1,iorder
            if ( i .eq. j) then
               res(j) = x(i,j)
            endif
 2       continue
 1    continue
      return
      end
      subroutine diagrv(x,ldx,iorder,nx,vec,ldvec)
c===============================================================c
c     written by Vassilis Hajivassiliou                         c
c     with the aid of Yoosoon Chang.                            c
c     version 1.0, March 22, 1992                               c
c     version 1.1, July 11, 1992                                c
c===============================================================c
c
c-----diagrv replaces the diagnol elements of a (m by n) matrix x with
c     leading dimension ldx with a vector, vec(ldvec).
c     iorder = min(m,n)
c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
c
      implicit double precision (a-h, o-z)
      dimension x(ldx,iorder), vec(ldvec)
c
      do 1 j = 1,nx
         do 2 i = 1,iorder
            if (i .eq. j) then
               x(i,j) = vec(i)
            endif
 2       continue
 1    continue
      return
      end
      subroutine eye(eyemat,ldx,iorder)
c===============================================================c
c     written by Vassilis Hajivassiliou                         c
c     with the aid of Yoosoon Chang.                            c
c     version 1.0, March 22, 1992                               c
c     version 1.1, July 11, 1992                                c
c===============================================================c
c
c-----eye creates "iorder" dimensional identity matrix with a leading
c     dimension ldx.
c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
c
      implicit double precision (a-h,o-z)
      dimension eyemat(ldx,iorder)
c
      do 1 j = 1,iorder
         do 2 i = 1,iorder
            if (i .eq. j) then
               eyemat(i,j) = 1.0d0
            else
               eyemat(i,j) = 0.0d0
            endif
 2       continue
 1    continue
      return
      end
      subroutine hadcomp(x,ldx,y,ldy,res,ldres,mx,nx,coper)
c===============================================================c
c     written by Vassilis Hajivassiliou                         c
c     with the aid of Yoosoon Chang.                            c
c     version 1.0, March 22, 1992                               c
c     version 1.1, July 11, 1992                                c
c===============================================================c
c
c-----hadcomp does element by element comparision of the matrices x and
c     y with the same dimensions (mx by nx) according to the coper code,
c     and stores the resulting values into the res matrix of the same
c     dimesion.
c     ldx,ldy,and ldres are the leading dimensions for x,y, and res
c     respectively.
c
c          coper         logical operation
c          -----         ------------------
c          'gt'                >
c          'ge'                >=
c          'eq'                =
c          'lt'                <
c          'le'                <=
c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
c
      implicit double precision (a-h,o-z)
      implicit integer*4 (i-n)
      dimension x(ldx,nx), y(ldy,nx), res(ldres,nx)
      character*2 coper
c
      call lc(coper)
c
      if (coper .eq. 'eq') then
         do 101 j = 1,nx
            do 102 i = 1,mx
               if (x(i,j) .eq. y(i,j)) then
                  res(i,j) = 1.0d0
               else
                  res(i,j) = 0.0d0
               endif
 102        continue
 101     continue
      else if (coper .eq. 'lt') then
         do 201 j = 1,nx
            do 202 i = 1,mx
               if (x(i,j) .lt. y(i,j)) then
                  res(i,j) = 1.0d0
               else
                  res(i,j) = 0.0d0
               endif
 202        continue
 201     continue
      else if (coper .eq. 'gt') then
         do 301 j = 1,nx
            do 302 i = 1,mx
               if (x(i,j) .gt. y(i,j)) then
                  res(i,j) = 1.0d0
               else
                  res(i,j) = 0.0d0
               endif
 302        continue
 301     continue
      else if (coper .eq. 'le') then
         do 401 j = 1,nx
            do 402 i = 1,mx
               if (x(i,j) .le. y(i,j)) then
                  res(i,j) = 1.0d0
               else
                  res(i,j) = 0.0d0
               endif
 402        continue
 401     continue
      else if (coper .eq. 'gt') then
         do 501 j = 1,nx
            do 502 i = 1,mx
               if (x(i,j) .ge. y(i,j)) then
                  res(i,j) = 1.0d0
               else
                  res(i,j) = 0.0d0
               endif
 502        continue
 501     continue
      else
         print *, coper,'is an invalid logical operator'
      endif
      return
      end
      subroutine hadoper(x,ldx,y,ldy,res,ldres,mx,nx,ioper)
c===============================================================c
c     written by Vassilis Hajivassiliou                         c
c     with the aid of Yoosoon Chang.                            c
c     version 1.0, March 22, 1992                               c
c     version 1.1, July 11, 1992                                c
c===============================================================c
c
c-----hadoper does element by element arithmetic operations of the
c     matrices x and y with the same dimensions (mx by nx) according to
c     the "ioper" code,and stores the resulting values into the res
c     matrix of the same dimesion.
c     ldx,ldy,and ldres are the leading dimensions for x,y, and res
c     respectively.
c
c     ioper         operation
c     ------------------------
c       1              +
c       2              -
c       3              *
c       4              /
c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
c
      implicit double precision (a-h,o-z)
      implicit integer*4 (i-n)
      dimension x(ldx,nx), y(ldy,nx), res(ldres,nx)
c
      if (ioper .eq. 1) then
         do 101 j = 1,nx
            do 102 i = 1,mx
               res(i,j) = x(i,j) + y(i,j)
 102        continue
 101     continue
      else if (ioper .eq. 2) then
         do 201 j = 1,nx
            do 202 i = 1,mx
               res(i,j) = x(i,j) - y(i,j)
 202        continue
 201     continue
      else if (ioper .eq. 3) then
         do 301 j = 1,nx
            do 302 i = 1,mx
               res(i,j) = x(i,j) * y(i,j)
 302        continue
 301     continue
      else if (ioper .eq. 4) then
         do 401 j = 1,nx
            do 402 i = 1,mx
               res(i,j) = x(i,j) / y(i,j)
 402        continue
 401     continue
      else
         print *, ioper,'is an invalid code.'
      endif
      return
      end
      function intfloor(x)
      integer*4 intfloor
      double precision x
c===============================================================c
c     written by Vassilis Hajivassiliou                         c
c     with the aid of Yoosoon Chang.                            c
c     version 1.0, March 22, 1992                               c
c     version 1.1, July 11, 1992                                c
c===============================================================c
c
c     intfloor takes double precision number x and truncates fractional
c     part down toward negative infinity. So, on exit, x has a
c     integer value.
c
c     e.g.) intfloor(-14.10d0) = -15
c           intfloor(4.99d0) = 4
c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
c
      if (x .ge. 0.0d0) then
ccc   ixfloor = idint(x)
        ixfloor = x
      else
ccc   ixfloor = idint(x) - 1
        ixfloor = x - 1
      endif
      x=ixfloor
      intfloor=ixfloor
      return
      end
      subroutine matcopy(y,ldy,x,ldx,mstart,mend,nstart,nend)
c===============================================================c
c     written by Vassilis Hajivassiliou                         c
c     with the aid of Yoosoon Chang.                            c
c     version 1.0, March 22, 1992                               c
c     version 1.1, July 11, 1992                                c
c===============================================================c
c
c-----matcopy copies a matrix x[mstart=mend,nstart=nend] with ldx
c     leading dimesion to a matrix y with the dimensions (mend-mstart+1)
c     by (nend-nstart+1) and leading dimension ldy.
c
c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
c
      implicit double precision (a-h, o-z)
      implicit integer*4 (i-n)
      dimension x(ldx,nend), y(ldy,nend)
c
      do 1 j = nstart,nend
         do 2 i = mstart,mend
            y(i-mstart+1,j-nstart+1) = x(i,j)
 2       continue
 1    continue
c
      return
      end
      subroutine maxc(x,ldx,mx,nx,res,ldres)
c===============================================================c
c     written by Vassilis Hajivassiliou                         c
c     with the aid of Yoosoon Chang.                            c
c     version 1.0, March 22, 1992                               c
c     version 1.1, July 11, 1992                                c
c===============================================================c
c
c     maxc is based on the double precision function FMAX in VFLIB.
c
c-----maxc takes the maximal value of each column of a matirx x(mx,nx)
c     with a leading dimension ldx, and stores them in res(nx).
c     The resulting vector res(nx) with a leading dimension ldres
c     is then returned.
c
c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
c
      implicit double precision (a-h, o-z)
      dimension x(ldx,nx), res(ldres)
c
      do 30 j = 1,nx
         res(j) = fmax(x(1,j),mx)
 30   continue
      return
      end
      subroutine mcdfn(x,ldx,mx,nx,res,ldres)
c===============================================================c
c     written by Vassilis Hajivassiliou                         c
c     with the aid of Yoosoon Chang.                            c
c     version 1.0, March 22, 1992                               c
c     version 1.1, July 11, 1992                                c
c===============================================================c
c
c-----mcdfn computes the cumulative distribution function(cdf) at each
c     element of the input matrix x(mx,nx) with a leading dimension ldx,
c     using the double precision function PHEEV from VFLIB.
c     Res(ldres,nx) stores the resulting values and is returned.
c
c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
c
      implicit double precision (a-h, o-z)
      implicit integer*4 (i-n)
      dimension x(ldx,nx), res(ldres,nx)
c
      do 1 j = 1,nx
         do 2 i = 1,mx
            res(i,j) = pheev(x(i,j))
 2       continue
 1    continue
c
      return
      end
      subroutine mcdfni(x,ldx,mx,nx,res,ldres)
c===============================================================c
c     written by Vassilis Hajivassiliou                         c
c     with the aid of Yoosoon Chang.                            c
c     version 1.0, March 22, 1992                               c
c     version 1.1, July 11, 1992                                c
c===============================================================c
c
c     This routine is based on the double precision function PHINV
c     in the VFLIB.
c
c-----mcdfn takes a (mx by nx) probability matrix x with a leading
c     dimension ldx and calculates, for each element of x,
c     the normal inverse cdf(quantile). Res(ldres,nx) stores the
c     resulting values and is returned.
c
c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
c
      implicit double precision (a-h, o-z)
      implicit integer*4 (i-n)
      dimension x(ldx,nx), res(ldres,nx)
c
      do 1 j = 1,nx
         do 2 i=1,mx
            res(i,j) = phinv(x(i,j))
 2       continue
 1    continue
c
      return
      end
      subroutine mdscal(x,ldx,res,ldres,mx,nx,scale)
c===============================================================c
c     written by Vassilis Hajivassiliou                         c
c     with the aid of Yoosoon Chang.                            c
c     version 1.0, March 22, 1992                               c
c     version 1.1, July 11, 1992                                c
c===============================================================c
c
c-----mdscal scales a matrix x(mx,nx) with a leading dimesion ldx,
c     by the value of the given scale.
c     Res(ldres,nx) stores the scaled values of x, and is returned.
c
c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
c
      implicit double precision (a-h, o-z)
      implicit integer*4 (i-n)
      dimension x(ldx,nx), res(ldres,nx)
c
      do 1 j = 1,nx
         do 2 i=1,mx
            res(i,j) = x(i,j)*scale
 2       continue
 1    continue
      return
      end
      subroutine meanc(x,ldx,mx,nx,res,ldres)
c===============================================================c
c     written by Vassilis Hajivassiliou                         c
c     with the aid of Yoosoon Chang.                            c
c     version 1.0, March 22, 1992                               c
c     version 1.1, July 11, 1992                                c
c===============================================================c
c
c-----meanc takes a matrix x(mx,nx) with a leading dimension ldx, and
c     computes the mean of every column of the matrix.
c     Res(nx) with a leading dimensin ldres stores the mean of each
c     column and is returned.
c
c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
c
      implicit double precision (a-h,o-z)
      dimension x(ldx,nx), res(ldres)
c
cvahfix      dm = m
      dmx = mx
      do 10 j = 1,nx
         sum = 0.0d0
         do 20 i = 1,mx
            sum = sum + x(i,j)
 20      continue
cha         temp = sum
cha         res(j) = temp/dble(float(mx))
cvahfix         res(j) = sum/dm
         res(j) = sum/dmx
 10   continue
      return
      end
      subroutine minc(x,ldx,mx,nx,res,ldres)
c===============================================================c
c     written by Vassilis Hajivassiliou                         c
c     with the aid of Yoosoon Chang.                            c
c     version 1.0, March 22, 1992                               c
c     version 1.1, July 11, 1992                                c
c===============================================================c
c
c     minc is based on the double precision function FMIN from VFLIB.
c
c-----minc takes the minimum value of each column of a matrix x(mx,nx),
c     with a leading dimension ldx, and stores them in res(nx) with a
c     leading dimension ldres. The resulting vector res(ldres) is then
c     returned.
c
c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
c
      implicit double precision (a-h,o-z)
      implicit integer*4 (i-n)
      dimension x(ldx,nx), res(ldres)
c
      do 10 j = 1,nx
         res(j) = fmin(x(1,j),mx)
 10   continue
      return
      end
      subroutine mrnorm(x,ldx,mx,nx,iseed)
c===============================================================c
c     written by Vassilis Hajivassiliou                         c
c     with the aid of Yoosoon Chang.                            c
c     version 1.0, March 22, 1992                               c
c     version 1.1, July 11, 1992                                c
c===============================================================c
c
c     This routine is based on the double precision function RANORM in
c     the VFLIB.
c
c-----mrnorm creates and returns a (mx by nx) matrix of standard normal
c     random numbers with a leading dimension ldx, using a seed given
c     by the argument iseed.
c
c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
c
      implicit double precision (a-h,o-z)
      implicit integer*4 (i-n)
      dimension x(ldx,nx)
c
      do 1 j = 1,nx
         do 2 i = 1,mx
            x(i,j) = ranorm(0.0d0,1.0d0,iseed)
 2       continue
 1    continue
      return
      end
      subroutine mrunif(x,ldx,mx,nx,iseed)
c===============================================================c
c     written by Vassilis Hajivassiliou                         c
c     with the aid of Yoosoon Chang.                            c
c     version 1.0, March 22, 1992                               c
c     version 1.1, July 11, 1992                                c
c===============================================================c
c
c-----mrunif generates and returns a matrix x(mx,nx) with a leading
c     dimension ldx, of uniform random numbers on the interval [0,1].
c     The argument iseed is given to the function RANDU from VFLIB.
c     The argument iseed is the current random number generator seed
c     and RANDU updates it.
c
c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
c
      implicit double precision (a-h,o-z)
      implicit integer*4 (i-n)
      dimension x(ldx,nx)
c
      do 1 j = 1,nx
         do 2 i = 1,mx
cvahfix            x(i,j) = unifrm(iseed)
            x(i,j) = randu(iseed)
 2       continue
 1    continue
      return
      end
      subroutine packr(x,ldx,m,n,codemis,icount)
c===============================================================c
c     written by Vassilis Hajivassiliou                         c
c     with the aid of Yoosoon Chang.                            c
c     version 1.0, March 22, 1992                               c
c     version 1.1, July 11, 1992                                c
c===============================================================c
c
c     packr takes m by n matrix with ldx leading dimension and checks
c     each element of x for codemis. If any element of a row matches
c     with the codemis, the row is deleted.
c     On exit, the matrix x is free of elements=codemis)
c
c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
      implicit double precision (a-h,o-z)
      dimension x(ldx,n)
c
      icount = 0
      i = 0
      iflag = 1
 999  continue
      if ( i .ne. m ) then
         i = i + 1
         do 10 j = 1,n
            if ( x(i,j) .eq. codemis ) then
               iflag = 0
               goto 999
           endif
           iflag = 1
 10      continue
         if (iflag .eq. 1) then
            icount = icount + 1
            do 20 j = 1,n
               x(icount,j) = x(i,j)
 20         continue
            go to 999
         endif
      endif
c
      do 30 j=1,n
         do 40 i = icount+1, m
            x(i,j) = 0.0d0
 40      continue
 30   continue
c
      return
      end
      subroutine prodc(x,ldx,mx,nx,res,ldres)
c===============================================================c
c     written by Vassilis Hajivassiliou                         c
c     with the aid of Yoosoon Chang.                            c
c     version 1.0, March 22, 1992                               c
c     version 1.1, July 11, 1992                                c
c===============================================================c
c
c returns products of columns.
c
c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
c
      implicit double precision (a-h,o-z)
      implicit integer*4 (i-n)
      dimension x(ldx,nx), res(ldres)
c
      do 10 j = 1,nx
         prod = 1.0d0
         do 20 i=1,mx
            prod = prod * x(i,j)
 20      continue
         res(j) = prod
 10   continue
      return
      end
      subroutine rbind(x,ldx,mx,n,y,ldy,my,res,ldres)
c===============================================================c
c     written by Vassilis Hajivassiliou                         c
c     with the aid of Yoosoon Chang.                            c
c     version 1.0, March 22, 1992                               c
c     version 1.1, July 11, 1992                                c
c===============================================================c
c
c     rbind takes two input matrices of the following dimension and
c     performs a vertical concatenation on them.
c
c      x(ldx,mx,n), where ldx is the leading dimension
c      y(ldy,my,n), where ldy is the leading dimension
c
c     The input matrices x and y must have the same number of columns.
c     The number of rows in the resulting matrix res is the sum
c     of mx and my.
c     The vertically concatenated matrix res(ldres,mx+my,n)is then
c     returned.
c
c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
c
      implicit double precision (a-h,o-z)
      implicit integer*4 (i-n)
      dimension x(ldx,n), y(ldy,n),res(ldres,n)
c
      do 1 j = 1, n
         do 2 i = 1,mx
            res(i,j) = x(i,j)
 2       continue
         do 3 i = mx+1, mx+my
            res(i,j) = y(i-mx,j)
 3       continue
 1    continue
c
      return
      end
      subroutine stdc(x,ldx,mx,nx,res,ldres)
c===============================================================c
c     written by Vassilis Hajivassiliou                         c
c     with the aid of Yoosoon Chang.                            c
c     version 1.0, March 22, 1992                               c
c     version 1.1, July 11, 1992                                c
c===============================================================c
c
c-----stdc takes a matrix x(ldx,mx,nx) where ldx is the leading
c     dimension of x and computes the standard deviation of the
c     elements in each column of x.
c     The vector res(nx) with a leading dimension ldres stores the
c     resulting standard deviations and is returned.
c
c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
      implicit double precision (a-h,o-z)
      dimension x(ldx,nx), res(ldres)
c
cvahfix      dm = m
      dmx = mx
      do 3 j = 1,nx
         sum = 0.0d0
         do 4 i = 1,mx
            sum = sum + x(i,j)
 4       continue
         temp = sum
cha         cmean = temp/dble(float(mx))
cvahfix         cmean = temp/dm
         cmean = temp/dmx
         var = 0.0d0
         do 5 k = 1,mx
            var = var + (x(k,j) - cmean)**2.0d0
 5       continue
         ttemp = var
         res(j) = ttemp **(0.5d0)
 3    continue
      return
      end
      subroutine sumc(x,ldx,mx,nx,res,ldres)
c===============================================================c
c     written by Vassilis Hajivassiliou                         c
c     with the aid of Yoosoon Chang.                            c
c     version 1.0, March 22, 1992                               c
c     version 1.1, July 11, 1992                                c
c===============================================================c
c
c-----sumc takes a matrix x(mx,nx) with a leading dimension ldx, and
c     computes the sum of every column of the matrix.
c     The resulting vector res(nx) with a leading dimension ldres
c     stores the sums of each column and is returned.
c
c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
c
      implicit double precision (a-h,o-z)
      implicit integer*4 (i-n)
      dimension x(ldx,nx), res(ldres)
c
      do 10 j = 1,nx
         sum = 0.0d0
         do 20 i=1,mx
            sum = sum + x(i,j)
 20      continue
         res(j) = sum
 10   continue
      return
      end
      subroutine transps(x,ldx,mx,nx,res,ldres)
c===============================================================c
c     written by Vassilis Hajivassiliou                         c
c     with the aid of Yoosoon Chang.                            c
c     version 1.0, March 22, 1992                               c
c     version 1.1, July 11, 1992                                c
c===============================================================c
c
c-----transps transposes a (mx by nx) matrix with leading dimension ldx.
c     It returns the transposed(nx by mx) matrix res with leading
c     dimension ldres.
c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
c
      implicit double precision (a-h,o-z)
      implicit integer*4 (i-n)
      dimension x(ldx,nx),res(ldres,mx)
c
      do 3 j = 1,nx
         do 4 i = 1,mx
            res(j,i) = x(i,j)
 4       continue
 3    continue
      return
      end