      subroutine cexist(cfile,nrc)
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     vah 3:29 am, monday, october 27, 1986
c
c     query the existence of a file using the inquire statement
c     nrc=0 if file exists, = -1 if not
c
c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
      character*(*) cfile
      logical lexist
      integer*4 nrc
      nrc=-1
      inquire(file=cfile,exist=lexist)
c     print *,' cmscmd,nrc :',cmscmd,nrc
      if(lexist) nrc=0
      return
      end
      double precision function chicdf(cs,df,ier)
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   constructed from the
c
c   imsl routine name   - mdch
c
c   by vah, 16:28:44.31, tue, jul 30, 1991
c
c   changed into full double precision.  calls to mdnor replaced by
c   calls to pheev
c
c-----------------------------------------------------------------------
c
c   computer            - ibm77: double
c
c   latest revision     -  16:30:02.36, tue, jul 30, 1991
c
c   purpose             - chi-squared probability distribution function
c
c   usage               - p=chicdf(cs,df,ier)
c
c   arguments    cs     - input value for which the probability is
c                           computed. cs must be greater than or equal
c                           to zero.
c                df     - input number of degrees of freedom of the
c                           chi-squared distribution. df must be greater
c                           than or equal to .5 and less than or equal
c                           to 200,000.
c                ier    - error parameter. (output)
c                         terminal error
c                           ier = 129 indicates that cs or df was
c                             specified incorrectly.
c                         warning error
c                           ier = 34 indicates that the normal pdf
c                             would have produced an underflow.
c
c   returns:     p      - output probability that a random variable
c                           which follows the chi-squared distribution
c                           with df degrees of freedom is less than or
c                           equal to cs.
c   precision/hardware  - double/all
c
c   reqd. routines      - pheev, dgamma
c
c-----------------------------------------------------------------------
c
c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
c                                  specifications for arguments
      implicit double precision (a-h,o-z)
ccc      real               cs,df,p
      double precision   cs,df,p
c                                  specifications for local variables
ccc      real               pt2
      double precision   pt2
      double precision   a,z,dgam,eps,w,w1,b,z1,half,one,thrten,thrd
      double precision   dgamma
ccc      real               rinfm,x,c
ccc      notice that c is declared but never used, so removed
      double precision   rinfm,x
      data               eps/1.0d-6/,half/5.d-1/,thrten/13.d0/,one/1.d0/
      data               thrd/.3333333333333333d0/
      data               pt2/.2222222e0/
ccc      data               rinfm/-0.7237e+76/
      data               rinfm/-0.7237d+76/
      func(w,a,z)=w*dexp(a*dlog(z)-z)
c                                  first executable statement
c                                  test for invalid input values
      if (df .ge. .5 .and. df .le. 2.e5 .and. cs .ge. 0.0) go to 5
      ier=129
      p=rinfm
      go to 9000
    5 ier=0
c                                  set p=0. if cs is less than or
c                                  equal to 10.**(-12)
      if (cs .gt. 1.e-12) go to 15
   10 p=0.0
      go to 9005
   15 if(df.le.100.) go to 20
c                                  use normal distribution approximation
c                                  for large degrees of freedom
      if(cs.lt.2.0) go to 10
      x=((cs/df)**thrd-(one-pt2/df))/sqrt(pt2/df)
      if (x .gt. 5.0) go to 50
      if (x .lt. -18.8055) go to 55
ccc      call mdnor (x,p)
      p=pheev(x)
      go to 9005
c                                  initialization for calculation using
c                                  incomplete gamma function
   20 if (cs .gt. 200.) go to 50
      a=half*df
      z=half*cs
      dgam = dgamma(a)
      w=dmax1(half*a,thrten)
      if (z .ge. w) go to 35
      if (df .gt. 25. .and. cs .lt. 2.) go to 10
c                                  calculate using equation no. 6.5.29
      w=one/(dgam*a)
      w1=w
         do 25 i=1,50
         b=i
         w1=w1*z/(a+b)
         if (w1 .le. eps*w) go to 30
         w=w+w1
   25    continue
   30 p=func(w,a,z)
      go to 9005
c                                  calculate using equation no. 6.5.32
   35 z1=one/z
      b=a-one
      w1=b*z1
      w=one+w1
         do 40 i=2,50
         b=b-one
         w1=w1*b*z1
         if (w1 .le. eps*w) go to 45
         w=w+w1
   40    continue
   45 w=z1*func(w,a,z)
      p=one-w/dgam
      go to 9005
   50 p=1.0
      go to 9005
c                                  warning error - underflow would have
c                                  occurred
   55 p=0.0
      ier=34
 9000 continue
ccc      call uertst (ier,'mdch  ')
ccc 9005 return
 9005 continue
      chicdf=p
      return
      end
      subroutine daxpy(n,da,dx,incx,dy,incy)
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     constant times a vector plus a vector.
c     uses unrolled loops for increments equal to one.
c     jack dongarra, linpack, 3/11/78.
c
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
      subroutine dchdc(a,lda,p,work,jpvt,job,info)
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
      integer lda,p,jpvt(1),job,info
      double precision a(lda,1),work(1)
c
c     dchdc computes the cholesky decomposition of a positive definite
c     matrix.  a pivoting option allows the user to estimate the
c     condition of a positive definite matrix or determine the rank
c     of a positive semidefinite matrix.
c
c     on entry
c
c         a      double precision(lda,p).
c                a contains the matrix whose decomposition is to
c                be computed.  onlt the upper half of a need be stored.
c                the lower part of the array a is not referenced.
c
c         lda    integer.
c                lda is the leading dimension of the array a.
c
c         p      integer.
c                p is the order of the matrix.
c
c         work   double precision.
c                work is a work array.
c
c         jpvt   integer(p).
c                jpvt contains integers that control the selection
c                of the pivot elements, if pivoting has been requested.
c                each diagonal element a(k,k)
c                is placed in one of three classes according to the
c                value of jpvt(k).
c
c                   if jpvt(k) .gt. 0, then x(k) is an initial
c                                      element.
c
c                   if jpvt(k) .eq. 0, then x(k) is a free element.
c
c                   if jpvt(k) .lt. 0, then x(k) is a final element.
c
c                before the decomposition is computed, initial elements
c                are moved by symmetric row and column interchanges to
c                the beginning of the array a and final
c                elements to the end.  both initial and final elements
c                are frozen in place during the computation and only
c                free elements are moved.  at the k-th stage of the
c                reduction, if a(k,k) is occupied by a free element
c                it is interchanged with the largest free element
c                a(l,l) with l .ge. k.  jpvt is not referenced if
c                job .eq. 0.
c
c        job     integer.
c                job is an integer that initiates column pivoting.
c                if job .eq. 0, no pivoting is done.
c                if job .ne. 0, pivoting is done.
c
c     on return
c
c         a      a contains in its upper half the cholesky factor
c                of the matrix a as it has been permuted by pivoting.
c
c         jpvt   jpvt(j) contains the index of the diagonal element
c                of a that was moved into the j-th position,
c                provided pivoting was requested.
c
c         info   contains the index of the last positive diagonal
c                element of the cholesky factor.
c
c     for positive definite matrices info = p is the normal return.
c     for pivoting with positive semidefinite matrices info will
c     in general be less than p.  however, info may be greater than
c     the rank of a, since rounding error can cause an otherwise zero
c     element to be positive. indefinite systems will always cause
c     info to be less than p.
c
c     linpack. this version dated 08/14/78 .
c     j.j. dongarra and g.w. stewart, argonne national laboratory and
c     university of maryland.
c
c
c     blas daxpy,dswap
c     fortran dsqrt
c
c     internal variables
c
      integer pu,pl,plp1,i,j,jp,jt,k,kb,km1,kp1,l,maxl
      double precision temp
      double precision maxdia
      logical swapk,negk
c
      pl = 1
      pu = 0
      info = p
      if (job .eq. 0) go to 160
c
c        pivoting has been requested. rearrange the
c        the elements according to jpvt.
c
         do 70 k = 1, p
            swapk = jpvt(k) .gt. 0
            negk = jpvt(k) .lt. 0
            jpvt(k) = k
            if (negk) jpvt(k) = -jpvt(k)
            if (.not.swapk) go to 60
               if (k .eq. pl) go to 50
                  call dswap(pl-1,a(1,k),1,a(1,pl),1)
                  temp = a(k,k)
                  a(k,k) = a(pl,pl)
                  a(pl,pl) = temp
                  plp1 = pl + 1
                  if (p .lt. plp1) go to 40
                  do 30 j = plp1, p
                     if (j .ge. k) go to 10
                        temp = a(pl,j)
                        a(pl,j) = a(j,k)
                        a(j,k) = temp
                     go to 20
   10                continue
                     if (j .eq. k) go to 20
                        temp = a(k,j)
                        a(k,j) = a(pl,j)
                        a(pl,j) = temp
   20                continue
   30             continue
   40             continue
                  jpvt(k) = jpvt(pl)
                  jpvt(pl) = k
   50          continue
               pl = pl + 1
   60       continue
   70    continue
         pu = p
         if (p .lt. pl) go to 150
         do 140 kb = pl, p
            k = p - kb + pl
            if (jpvt(k) .ge. 0) go to 130
               jpvt(k) = -jpvt(k)
               if (pu .eq. k) go to 120
                  call dswap(k-1,a(1,k),1,a(1,pu),1)
                  temp = a(k,k)
                  a(k,k) = a(pu,pu)
                  a(pu,pu) = temp
                  kp1 = k + 1
                  if (p .lt. kp1) go to 110
                  do 100 j = kp1, p
                     if (j .ge. pu) go to 80
                        temp = a(k,j)
                        a(k,j) = a(j,pu)
                        a(j,pu) = temp
                     go to 90
   80                continue
                     if (j .eq. pu) go to 90
                        temp = a(k,j)
                        a(k,j) = a(pu,j)
                        a(pu,j) = temp
   90                continue
  100             continue
  110             continue
                  jt = jpvt(k)
                  jpvt(k) = jpvt(pu)
                  jpvt(pu) = jt
  120          continue
               pu = pu - 1
  130       continue
  140    continue
  150    continue
  160 continue
      do 270 k = 1, p
c
c        reduction loop.
c
         maxdia = a(k,k)
         kp1 = k + 1
         maxl = k
c
c        determine the pivot element.
c
         if (k .lt. pl .or. k .ge. pu) go to 190
            do 180 l = kp1, pu
               if (a(l,l) .gt. maxdia) go to 170
                  maxdia = a(l,l)
                  maxl = l
  170          continue
  180       continue
  190    continue
c
c        quit if the pivot element is not positive.
c
         if (maxdia .gt. 0.0d0) go to 200
            info = k - 1
c     ......exit
            go to 280
  200    continue
         if (k .eq. maxl) go to 210
c
c           start the pivoting and update jpvt.
c
            km1 = k - 1
            call dswap(km1,a(1,k),1,a(1,maxl),1)
            a(maxl,maxl) = a(k,k)
            a(k,k) = maxdia
            jp = jpvt(maxl)
            jpvt(maxl) = jpvt(k)
            jpvt(k) = jp
  210    continue
c
c        reduction step. pivoting is contained across the rows.
c
         work(k) = dsqrt(a(k,k))
         a(k,k) = work(k)
         if (p .lt. kp1) go to 260
         do 250 j = kp1, p
            if (k .eq. maxl) go to 240
               if (j .ge. maxl) go to 220
                  temp = a(k,j)
                  a(k,j) = a(j,maxl)
                  a(j,maxl) = temp
               go to 230
  220          continue
               if (j .eq. maxl) go to 230
                  temp = a(k,j)
                  a(k,j) = a(maxl,j)
                  a(maxl,j) = temp
  230          continue
  240       continue
            a(k,j) = a(k,j)/work(k)
            work(j) = a(k,j)
            temp = -a(k,j)
            call daxpy(j-k,temp,work(kp1),1,a(kp1,j),1)
  250    continue
  260    continue
  270 continue
  280 continue
      return
      end
      subroutine  dcopy(n,dx,incx,dy,incy)
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     copies a vector, x, to a vector, y.
c     uses unrolled loops for increments equal to one.
c     jack dongarra, linpack, 3/11/78.
c
c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
      double precision dx(1),dy(1)
      integer i,incx,incy,ix,iy,m,mp1,n
c
      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
        dy(iy) = 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,7)
      if( m .eq. 0 ) go to 40
      do 30 i = 1,m
        dy(i) = dx(i)
   30 continue
      if( n .lt. 7 ) return
   40 mp1 = m + 1
      do 50 i = mp1,n,7
        dy(i) = dx(i)
        dy(i + 1) = dx(i + 1)
        dy(i + 2) = dx(i + 2)
        dy(i + 3) = dx(i + 3)
        dy(i + 4) = dx(i + 4)
        dy(i + 5) = dx(i + 5)
        dy(i + 6) = dx(i + 6)
   50 continue
      return
      end
      double precision function ddot(n,dx,incx,dy,incy)
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     forms the dot product of two vectors.
c     uses unrolled loops for increments equal to one.
c     jack dongarra, linpack, 3/11/78.
c
c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
      double precision dx(1),dy(1),dtemp
      integer i,incx,incy,ix,iy,m,mp1,n
c
      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
      function dgamma (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                                  specifications for arguments
      double precision   x,dgamma
c                                  specifications for local variables
      dimension          p(9),q(8),p4(7)
      double precision   a,b,den,p,p4,q,pi,t,top,big1,
     1                   xinf,xmin,y,sign,r
      integer            i,ier,j
      logical            mflag
c                                  coefficients for minimax
c                                  approximation to gamma(x),
c                                  2.0 .le. x .le. 3.0
      data               p(1)/-5.966047488753637d01/,
     *                   p(2)/5.864023793062003d01/,
     *                   p(3)/-1.364106217165365d03/,
     *                   p(4)/-8.117569271425580d02/,
     *                   p(5)/-1.569414683149179d04/,
     *                   p(6)/-1.525979925758372d04/,
     *                   p(7)/-7.264059615964330d04/,
     *                   p(8)/-8.972275718101010d-01/,
     *                   p(9)/3.349618189847578d00/
      data               q(1)/4.103991474182904d02/,
     *                   q(2)/-2.262590291514875d03/,
     *                   q(3)/2.494325576714903d03/,
     *                   q(4)/2.362106244383048d04/,
     *                   q(5)/-5.741873227396418d04/,
     *                   q(6)/-7.257239715408240d04/,
     *                   q(7)/-9.491399521686949d00/,
     *                   q(8)/-3.255006939455704d01/
c                                  coefficients for minimax
c                                  approximation to ln(gamma(x)),
c                                  12.0 .le. x
      data               p4(1)/8.40596949829d-04/,
     1                   p4(2)/-5.9523334141881d-04/,
     2                   p4(3)/7.9365078409227d-04/,
     *                   p4(4)/-2.777777777769526d-03/,
     *                   p4(5)/8.333333333333333d-02/,
     *                   p4(6)/9.189385332046727d-01/,
     6                   p4(7)/-1.7816839846d-03/
      data               pi/3.141592653589793d0/
      data               xinf/1.7d+38/
c                                  dgamma(xmin) .approx. xinf
c                                  dgamma(big1) .approx. xinf
      data               xmin/5.8775d-39/
      data               big1/34.844d0/
c                                  first executable statement
      ier = 0
      mflag = .false.
      t = x
      if (dabs(t).gt.xmin) go to 5
      ier = 130
      dgamma = xinf
      if (t.le.0.0d0) dgamma = -xinf
      go to 9000
    5 if (dabs(t).lt.big1) go to 10
      ier = 129
      dgamma = xinf
      go to 9000
   10 if (t.gt.0.0d0) go to 25
c                                  argument is negative
      mflag = .true.
      t = -t
      r = dint(t)
      sign = 1.0d0
      if (dmod(t,2.0d0).eq.0.0d0) sign = -1.0d0
      r = t-r
      if (r.ne.0.0d0) go to 20
      ier = 130
      dgamma = xinf
      if (sign.eq.-1.0d0) dgamma = -xinf
      go to 9000
c                                  argument is not a negative integer
   20 r = pi/dsin(r*pi)*sign
      t = t+1.0d0
c                                  evaluate approximation for dgamma(t)
c                                    t .gt. xmin
   25 if (t.gt.12.0d0) go to 60
      i = t
      a = 1.0d0
      if (i.gt.2) go to 40
      i = i+1
      go to (30,35,50),i
c                                  0.0 .lt. t .lt. 1.0
   30 a = a/(t*(t+1.0d0))
      t = t+2.0d0
      go to 50
c                                  1.0 .le. t .lt. 2.0
   35 a = a/t
      t = t+1.0d0
      go to 50
c                                  3.0 .le. t .le. 12.0
   40 do 45 j=3,i
         t = t-1.0d0
         a = a*t
   45 continue
c                                  2.0 .le. t .le. 3.0
   50 top = p(8)*t+p(9)
      den = t+q(8)
      do 55 j=1,7
         top = top*t+p(j)
         den = den*t+q(j)
   55 continue
      y = (top/den)*a
      if (mflag) y = r/y
      dgamma = y
      go to 9005
c                                  t .gt. 12.0
   60 top = dlog(t)
      top = t*(top-1.0d0)-.5d0*top
      t = 1.0d0/t
      b = t*t
      a = p4(7)
      do 65 j = 1,5
   65 a = a*b+p4(j)
      y = a*t+p4(6)+top
      y = dexp(y)
      if (mflag) y = r/y
      dgamma = y
      go to 9005
 9000 continue
c      call uertst(-ier,6hmgamad)
c      call uertst(ier,6hdgamma)
 9005 return
      end
      subroutine dsidi(a,lda,n,kpvt,det,inert,work,job)
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
      integer lda,n,job
      double precision a(lda,1),work(1)
      double precision det(2)
      integer kpvt(1),inert(3)
c
c     dsidi computes the determinant, inertia and inverse
c     of a double precision symmetric matrix using the factors from
c     dsifa.
c
c     on entry
c
c        a       double precision(lda,n)
c                the output from dsifa.
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        kpvt    integer(n)
c                the pivot vector from dsifa.
c
c        work    double precision(n)
c                work vector.  contents destroyed.
c
c        job     integer
c                job has the decimal expansion  abc  where
c                   if  c .ne. 0, the inverse is computed,
c                   if  b .ne. 0, the determinant is computed,
c                   if  a .ne. 0, the inertia is computed.
c
c                for example, job = 111  gives all three.
c
c     on return
c
c        variables not requested by job are not used.
c
c        a      contains the upper triangle of the inverse of
c               the original matrix.  the strict lower triangle
c               is never referenced.
c
c        det    double precision(2)
c               determinant of original matrix.
c               determinant = det(1) * 10.0**det(2)
c               with 1.0 .le. dabs(det(1)) .lt. 10.0
c               or det(1) = 0.0.
c
c        inert  integer(3)
c               the inertia of the original matrix.
c               inert(1)  =  number of positive eigenvalues.
c               inert(2)  =  number of negative eigenvalues.
c               inert(3)  =  number of zero eigenvalues.
c
c     error condition
c
c        a division by zero may occur if the inverse is requested
c        and  dsico  has set rcond .eq. 0.0
c        or  dsifa  has set  info .ne. 0 .
c
c     linpack. this version dated 08/14/78 .
c     james bunch, univ. calif. san diego, argonne nat. lab
c
c     subroutines and functions
c
c     blas daxpy,dcopy,ddot,dswap
c     fortran dabs,iabs,mod
c
c     internal variables.
c
c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
      double precision akkp1,ddot,temp
      double precision ten,d,t,ak,akp1
      integer j,jb,k,km1,ks,kstep
      logical noinv,nodet,noert
c
      noinv = mod(job,10) .eq. 0
      nodet = mod(job,100)/10 .eq. 0
      noert = mod(job,1000)/100 .eq. 0
c
      if (nodet .and. noert) go to 140
         if (noert) go to 10
            inert(1) = 0
            inert(2) = 0
            inert(3) = 0
   10    continue
         if (nodet) go to 20
            det(1) = 1.0d0
            det(2) = 0.0d0
            ten = 10.0d0
   20    continue
         t = 0.0d0
         do 130 k = 1, n
            d = a(k,k)
c
c           check if 1 by 1
c
            if (kpvt(k) .gt. 0) go to 50
c
c              2 by 2 block
c              use det (d  s)  =  (d/t * c - t) * t  ,  t = dabs(s)
c                      (s  c)
c              to avoid underflow/overflow troubles.
c              take two passes through scaling.  use  t  for flag.
c
               if (t .ne. 0.0d0) go to 30
                  t = dabs(a(k,k+1))
                  d = (d/t)*a(k+1,k+1) - t
               go to 40
   30          continue
                  d = t
                  t = 0.0d0
   40          continue
   50       continue
c
            if (noert) go to 60
               if (d .gt. 0.0d0) inert(1) = inert(1) + 1
               if (d .lt. 0.0d0) inert(2) = inert(2) + 1
               if (d .eq. 0.0d0) inert(3) = inert(3) + 1
   60       continue
c
            if (nodet) go to 120
               det(1) = d*det(1)
               if (det(1) .eq. 0.0d0) go to 110
   70             if (dabs(det(1)) .ge. 1.0d0) go to 80
                     det(1) = ten*det(1)
                     det(2) = det(2) - 1.0d0
                  go to 70
   80             continue
   90             if (dabs(det(1)) .lt. ten) go to 100
                     det(1) = det(1)/ten
                     det(2) = det(2) + 1.0d0
                  go to 90
  100             continue
  110          continue
  120       continue
  130    continue
  140 continue
c
c     compute inverse(a)
c
      if (noinv) go to 270
         k = 1
  150    if (k .gt. n) go to 260
            km1 = k - 1
            if (kpvt(k) .lt. 0) go to 180
c
c              1 by 1
c
               a(k,k) = 1.0d0/a(k,k)
               if (km1 .lt. 1) go to 170
                  call dcopy(km1,a(1,k),1,work,1)
                  do 160 j = 1, km1
                     a(j,k) = ddot(j,a(1,j),1,work,1)
                     call daxpy(j-1,work(j),a(1,j),1,a(1,k),1)
  160             continue
                  a(k,k) = a(k,k) + ddot(km1,work,1,a(1,k),1)
  170          continue
               kstep = 1
            go to 220
  180       continue
c
c              2 by 2
c
               t = dabs(a(k,k+1))
               ak = a(k,k)/t
               akp1 = a(k+1,k+1)/t
               akkp1 = a(k,k+1)/t
               d = t*(ak*akp1 - 1.0d0)
               a(k,k) = akp1/d
               a(k+1,k+1) = ak/d
               a(k,k+1) = -akkp1/d
               if (km1 .lt. 1) go to 210
                  call dcopy(km1,a(1,k+1),1,work,1)
                  do 190 j = 1, km1
                     a(j,k+1) = ddot(j,a(1,j),1,work,1)
                     call daxpy(j-1,work(j),a(1,j),1,a(1,k+1),1)
  190             continue
                  a(k+1,k+1) = a(k+1,k+1) + ddot(km1,work,1,a(1,k+1),1)
                  a(k,k+1) = a(k,k+1) + ddot(km1,a(1,k),1,a(1,k+1),1)
                  call dcopy(km1,a(1,k),1,work,1)
                  do 200 j = 1, km1
                     a(j,k) = ddot(j,a(1,j),1,work,1)
                     call daxpy(j-1,work(j),a(1,j),1,a(1,k),1)
  200             continue
                  a(k,k) = a(k,k) + ddot(km1,work,1,a(1,k),1)
  210          continue
               kstep = 2
  220       continue
c
c           swap
c
            ks = iabs(kpvt(k))
            if (ks .eq. k) go to 250
               call dswap(ks,a(1,ks),1,a(1,k),1)
               do 230 jb = ks, k
                  j = k + ks - jb
                  temp = a(j,k)
                  a(j,k) = a(ks,j)
                  a(ks,j) = temp
  230          continue
               if (kstep .eq. 1) go to 240
                  temp = a(ks,k+1)
                  a(ks,k+1) = a(k,k+1)
                  a(k,k+1) = temp
  240          continue
  250       continue
            k = k + kstep
         go to 150
  260    continue
  270 continue
      return
      end
       subroutine dsifa(a,lda,n,kpvt,info)
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
      integer lda,n,kpvt(1),info
      double precision a(lda,1)
c
c     dsifa factors a double precision symmetric matrix by elimination
c     with symmetric pivoting.
c
c     to solve  a*x = b , follow dsifa by dsisl.
c     to compute  inverse(a)*c , follow dsifa by dsisl.
c     to compute  determinant(a) , follow dsifa by dsidi.
c     to compute  inertia(a) , follow dsifa by dsidi.
c     to compute  inverse(a) , follow dsifa by dsidi.
c
c     on entry
c
c        a       double precision(lda,n)
c                the symmetric matrix to be factored.
c                only the diagonal and upper triangle are used.
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       a block diagonal matrix and the multipliers which
c                were used to obtain it.
c                the factorization can be written  a = u*d*trans(u)
c                where  u  is a product of permutation and unit
c                upper triangular matrices , trans(u) is the
c                transpose of  u , and  d  is block diagonal
c                with 1 by 1 and 2 by 2 blocks.
c
c        kpvt    integer(n)
c                an integer vector of pivot indices.
c
c        info    integer
c                = 0  normal value.
c                = k  if the k-th pivot block is singular. this is
c                     not an error condition for this subroutine,
c                     but it does indicate that dsisl or dsidi may
c                     divide by zero if called.
c
c     linpack. this version dated 08/14/78 .
c     james bunch, univ. calif. san diego, argonne nat. lab.
c
c     subroutines and functions
c
c     blas daxpy,dswap,idamax
c     fortran dabs,dmax1,dsqrt
c
c     internal variables
c
c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
      double precision ak,akm1,bk,bkm1,denom,mulk,mulkm1,t
      double precision absakk,alpha,colmax,rowmax
      integer imax,imaxp1,j,jj,jmax,k,km1,km2,kstep,idamax
      logical swap
c
c
c     initialize
c
c     alpha is used in choosing pivot block size.
      alpha = (1.0d0 + dsqrt(17.0d0))/8.0d0
c
      info = 0
c
c     main loop on k, which goes from n to 1.
c
      k = n
   10 continue
c
c        leave the loop if k=0 or k=1.
c
c     ...exit
         if (k .eq. 0) go to 200
         if (k .gt. 1) go to 20
            kpvt(1) = 1
            if (a(1,1) .eq. 0.0d0) info = 1
c     ......exit
            go to 200
   20    continue
c
c        this section of code determines the kind of
c        elimination to be performed.  when it is completed,
c        kstep will be set to the size of the pivot block, and
c        swap will be set to .true. if an interchange is
c        required.
c
         km1 = k - 1
         absakk = dabs(a(k,k))
c
c        determine the largest off-diagonal element in
c        column k.
c
         imax = idamax(k-1,a(1,k),1)
         colmax = dabs(a(imax,k))
         if (absakk .lt. alpha*colmax) go to 30
            kstep = 1
            swap = .false.
         go to 90
   30    continue
c
c           determine the largest off-diagonal element in
c           row imax.
c
            rowmax = 0.0d0
            imaxp1 = imax + 1
            do 40 j = imaxp1, k
               rowmax = dmax1(rowmax,dabs(a(imax,j)))
   40       continue
            if (imax .eq. 1) go to 50
               jmax = idamax(imax-1,a(1,imax),1)
               rowmax = dmax1(rowmax,dabs(a(jmax,imax)))
   50       continue
            if (dabs(a(imax,imax)) .lt. alpha*rowmax) go to 60
               kstep = 1
               swap = .true.
            go to 80
   60       continue
            if (absakk .lt. alpha*colmax*(colmax/rowmax)) go to 70
               kstep = 1
               swap = .false.
            go to 80
   70       continue
               kstep = 2
               swap = imax .ne. km1
   80       continue
   90    continue
         if (dmax1(absakk,colmax) .ne. 0.0d0) go to 100
c
c           column k is zero.  set info and iterate the loop.
c
            kpvt(k) = k
            info = k
         go to 190
  100    continue
         if (kstep .eq. 2) go to 140
c
c           1 x 1 pivot block.
c
            if (.not.swap) go to 120
c
c              perform an interchange.
c
               call dswap(imax,a(1,imax),1,a(1,k),1)
               do 110 jj = imax, k
                  j = k + imax - jj
                  t = a(j,k)
                  a(j,k) = a(imax,j)
                  a(imax,j) = t
  110          continue
  120       continue
c
c           perform the elimination.
c
            do 130 jj = 1, km1
               j = k - jj
               mulk = -a(j,k)/a(k,k)
               t = mulk
               call daxpy(j,t,a(1,k),1,a(1,j),1)
               a(j,k) = mulk
  130       continue
c
c           set the pivot array.
c
            kpvt(k) = k
            if (swap) kpvt(k) = imax
         go to 190
  140    continue
c
c           2 x 2 pivot block.
c
            if (.not.swap) go to 160
c
c              perform an interchange.
c
               call dswap(imax,a(1,imax),1,a(1,k-1),1)
               do 150 jj = imax, km1
                  j = km1 + imax - jj
                  t = a(j,k-1)
                  a(j,k-1) = a(imax,j)
                  a(imax,j) = t
  150          continue
               t = a(k-1,k)
               a(k-1,k) = a(imax,k)
               a(imax,k) = t
  160       continue
c
c           perform the elimination.
c
            km2 = k - 2
            if (km2 .eq. 0) go to 180
               ak = a(k,k)/a(k-1,k)
               akm1 = a(k-1,k-1)/a(k-1,k)
               denom = 1.0d0 - ak*akm1
               do 170 jj = 1, km2
                  j = km1 - jj
                  bk = a(j,k)/a(k-1,k)
                  bkm1 = a(j,k-1)/a(k-1,k)
                  mulk = (akm1*bk - bkm1)/denom
                  mulkm1 = (ak*bkm1 - bk)/denom
                  t = mulk
                  call daxpy(j,t,a(1,k),1,a(1,j),1)
                  t = mulkm1
                  call daxpy(j,t,a(1,k-1),1,a(1,j),1)
                  a(j,k) = mulk
                  a(j,k-1) = mulkm1
  170          continue
  180       continue
c
c           set the pivot array.
c
            kpvt(k) = 1 - k
            if (swap) kpvt(k) = -imax
            kpvt(k-1) = kpvt(k)
  190    continue
         k = k - kstep
      go to 10
  200 continue
      return
      end
      subroutine  dswap (n,dx,incx,dy,incy)
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     interchanges two vectors.
c     uses unrolled loops for increments equal one.
c     jack dongarra, linpack, 3/11/78.
c
c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
      double precision dx(1),dy(1),dtemp
      integer i,incx,incy,ix,iy,m,mp1,n
c
      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 not equal
c         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 = dx(ix)
        dx(ix) = dy(iy)
        dy(iy) = dtemp
        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,3)
      if( m .eq. 0 ) go to 40
      do 30 i = 1,m
        dtemp = dx(i)
        dx(i) = dy(i)
        dy(i) = dtemp
   30 continue
      if( n .lt. 3 ) return
   40 mp1 = m + 1
      do 50 i = mp1,n,3
        dtemp = dx(i)
        dx(i) = dy(i)
        dy(i) = dtemp
        dtemp = dx(i + 1)
        dx(i + 1) = dy(i + 1)
        dy(i + 1) = dtemp
        dtemp = dx(i + 2)
        dx(i + 2) = dy(i + 2)
        dy(i + 2) = dtemp
   50 continue
      return
      end
      subroutine f3swap(x,y)
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     required for pheev3
c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
      implicit double precision(a-h,o-z)
      z=x
      x=y
      y=z
      return
      end
      subroutine fillm(vec,npi,npj,val)
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+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
      double precision vec(npi,npj),val
      do 100 j=1,npj
         do 101 i=1,npi
            vec(i,j)=val
  101    continue
  100 continue
      return
      end
      subroutine fillv(a,np,val)
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+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
      double precision a(np),val
      do 1 i=1,np
         a(i)=val
 1    continue
      return
      end
      double precision function fmax(vec,n)
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+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
      implicit double precision(a-h,o-z)
      dimension vec(n)
      tempmax=-1.d20
      do 1 iii=1,n
            temp=vec(iii)
            if(temp .gt. tempmax) tempmax=temp
 1    continue
      fmax=tempmax
      return
      end
      double precision function fmin(vec,n)
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+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
      implicit double precision(a-h,o-z)
      dimension vec(n)
      tempmin=1.d20
      do 1 iii=1,n
            temp=vec(iii)
            if(temp .lt. tempmin) tempmin=temp
 1    continue
      fmin=tempmin
      return
      end
      double precision function hsec()
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+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
      implicit double precision (a-h,o-z)
      external sec
 
      hsec = 100.d0 * sec()
      return
      end
      double precision function hsec1990()
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 vah 17:25:23.43, sat, nov 24, 1990
c
c returns number of hundredth's of a second since jan 1, 1990
c
c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
      implicit double precision(a-h,o-z)
      dimension idaysinm(12),idate(3),itime(5)
      external hsec,timdat
 
      idaysinm(1)=31
      idaysinm(2)=28
      idaysinm(3)=31
      idaysinm(4)=30
      idaysinm(5)=31
      idaysinm(6)=30
      idaysinm(7)=31
      idaysinm(8)=31
      idaysinm(9)=30
      idaysinm(10)=31
      idaysinm(11)=30
      idaysinm(12)=31
 
      call timdat(itime,idate)
      idom=idate(1)
      month=idate(2)
      year=idate(3)
      if(year/100.d0 .lt. 1.d0) year=year+1900.d0
      if(dmod(year,4.d0) .eq. 0.d0) idaysinm(2)=29
 
      doy=idom
      do 1 i=1,month-1
         doy=doy+idaysinm(i)
    1 continue
 
      hsec1990=hsec()+100.*3600.*24.*((year-1990.)*365.+doy)
      return
      end
      integer function icoccfst(needle,haystack)
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
c
c     vah 3/10/1988,10:51:56:23
c         20:40:48.97, Sat, Sep 1, 1990
c
c     Returns integer location where needle occurs first in haystack.
c     Returns 0 if no match.
c
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
      character*(*) needle,haystack
      lenn=len(needle)
      lenh=len(haystack)
      nnn=0
      maxlen=lenh-lenn+1
      do 1 i=1,maxlen
cdbg         print *,i,needle(1:lenn),haystack(i:(i+lenn-1)),
cdbg     * (needle(1:lenn).eq.haystack(i:(i+lenn-1)))
         if(needle(1:lenn).eq.haystack(i:(i+lenn-1))) then
            nnn=i
            goto 2
         endif
    1 continue
    2 continue
      icoccfst=nnn
cdbg      print *,'icoccfst returns',nnn,needle,haystack
cdbg      pause
      return
      end
      integer function icocclst(needle,haystack)
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
c
c     vah 14:30:53.47, Sat, Sep 1, 1990
c
c     Returns integer location where needle occurs last in haystack.
c     Returns 0 if no match.
c
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
      character*(*) needle,haystack
      lenn=len(needle)
      lenh=len(haystack)
      nnn=0
      maxlen=lenh-lenn+1
      do 1 i=1,maxlen
cdbg         print *,i,needle(1:lenn),haystack(lenh-lenn+2-i:lenh-i+1),
cdbg     *     (needle(1:lenn).eq.haystack(lenh-lenn+2-i:lenh-i+1))
         if(needle(1:lenn).eq.haystack(lenh-lenn+2-i:lenh-i+1)) then
            nnn=lenh-lenn+2-i
            goto 2
         endif
    1 continue
    2 continue
      icocclst=nnn
cdbg      print *,'icocclst returns',nnn,needle,haystack
cdbg      pause
      return
      end
      integer function idamax(n,dx,incx)
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     finds the index of element having max. absolute value.
c     jack dongarra, linpack, 3/11/78.
c
c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
      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
c=======================================================================
c     dlinpack routines for inversion end
c=======================================================================
      end
      subroutine invv(a,b,ia,ib,c,n1,npr,nzr,nnr,rcond,*)
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     vah : 6:14 pm, friday, july 11, 1986
c***********************************************************************
c
c     nb: assumes a is symmetric
c                      *********
c
c     on entry,
c        *****
c     a contains the matrix to be inverted,
c     ia is the lead dim of a, ib the lead dim of b,
c     n1 the effective dim.
c
c     on exit,
c        ****
c     a is the inverse, b contains the eigenvectors,
c     c is the vector of eigenvalues in descending order,
c     npr=# of +ve roots, nnr=# of -ve roots, nzr=# of 0 roots
c
c     error-returns
c     ************
c     case 1: eigen value calculations did not converge(matev2 int=2)
c     case 2: a root is smaller or equal
c     to machine absolute precision (eps2), without attempting to
c     invert the matrix. in both cases, a will be wrong in general
c
c     written by vassilis hajivassiliou 30mar85,
c     using matev2(a,b,n1,ia,ib,int)
c     from gqopt.
c
c     eps1= machine relative accuracy (max # of significant digits)
c     eps2= machine absolute accuracy (smallest epsilon > 0)
c***********************************************************************
c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
       implicit double precision(a-h,o-z)
       dimension a(ia,n1),b(ib,n1),c(n1)
ccc       real*4 x
       logical novec
      data eps1/1.d-12/
      data eps2/1.d-50/
ccc      common/bprec/eps1,eps2
 
      int=1
c     since we'll need the eigen vectors
 
       novec = int .eq. 0
 1     n = n1
       xn = n
       if(n-2) 2,13,3
 2     if(.not. novec) b(1,1) = 1.d0
       return
 3     nc = n - 2
       do 12 k = 1,nc
       kp = k + 1
       hq=0.d0
       do 9 ml1=kp,n
 9     hq=hq+a(ml1,k)*a(ml1,k)
       s=dsign(dsqrt(hq),a(k+1,k))
       if(dabs(s) .lt. 1.d-30) go to 12
       a(k+1,k) = a(k+1,k) + s
       alpha = 1.d0/dsqrt(s * a(k+1,k))
       do 4 kr = kp,n
 4     a(kr,k) = alpha*a(kr, k)
       do 5 kr = kp, n
       hq=0.d0
       do 10 ml2=kp,n
 10    hq=hq+a(kr,ml2)*a(ml2,k)
 5     a(k,kr)=hq
       hq=0.d0
       do 11 ml3=kp,n
 11    hq=hq+a(k,ml3)*a(ml3,k)
       beta=.5d0*hq
       do 6 kr = kp,n
 6     a(k,kr) = a(k,kr) - beta*a(kr,k)
       do 8 kr = kp,n
       do 7 j = kr,n
       a(kr,j) = a(kr, j)-a(k,j)*a(kr,k) - a(k,kr)*a(j,k)
 7     a(j,kr) = a(kr, j)
 8     continue
 12    a(k,k+1) = -s
 13    continue
       hq=0.d0
       do 14 ml4=1,n
 14    hq=hq+a(ml4,ml4)*a(ml4,ml4)
       nn=n-1
       do 15 ml9=1,nn
 15    hq=hq+2.d0*a(ml9,ml9+1)*a(ml9,ml9+1)
       s=dsqrt(hq)
       eps = dmax1(s*eps1/dsqrt(xn),eps2)
       if (novec) go to 280
       do 210 k = 1,n
       do 208 j = 1,n
 208   b(j,k) = 0.d0
 210   b(k,k) = 1.d0
       if (n .eq. 2) go to 280
       do 230 k = 1,nc
        kp = k + 1
       do 225 kr =  2,n
       hq=0.d0
       do 228 ml5=kp,n
 228   hq=hq+b(kr,ml5)*a(ml5,k)
       s=hq
       do 225 j = kp,n
 225   b(kr,j) = b(kr,j) - s*a(j,k)
 230   continue
 280   continue
       nl = n
 304   nf = 1
 305   if(nf .ne. nl) go to 308
 306   if(nf .eq. 1) go to 999
       nl = nf - 1
       go to 304
 308   nd = nl - nf
       nda = nl - 1
       itcnt = 0
       itlim = 30*(nl+1-nf)
 309   do 310 kk = 1,nd
       k = nl - kk
 3100  if(dabs(a(k,k+1)).le. eps) go to 315
 310   continue
       go to 322
 315   nf = k + 1
       go to 305
 322   q = a(nf, nf+1)
       p = a(nf,nf) - 0.5d0*(a(nl,nl) + a(nl-1,nl-1)) - dsign(
     x dsqrt(.25d0*(a(nl,nl)-a(nl-1,nl-1))**2 + a(nl-1,nl)**2) ,
     y a(nl,nl) + a(nl-1,nl-1) )
       itcnt = itcnt + 1
       do 350 k = nf,nda
       if(k .eq. nf) go to 325
       p = a(k-1,k)
       q = r
 325   s = p**2 + q**2
       pp = p**2/s
       qq = q**2/s
       pq = p*q/s
       t = 2.d0*pq*a(k, k+1)
       a(k,k+1) = (pp-qq)*a(k,k+1) + pq*(a(k+1,k+1) - a(k,k) )
       u = a(k,k)
       a(k,k) = pp*u + qq*a(k+1,k+1) + t
       a(k+1,k+1) = qq*u + pp*a(k+1,k+1) - t
       s = dsqrt(s)
       p = p/s
       q = q/s
       if(k .ne. nf) a(k-1,k) = p*a(k-1,k) + q*r
       if(k .eq. nda) goto 340
       r = q*a(k+1,k+2)
       a(k+1,k+2) = p*a(k+1,k+2)
 340   if(novec) go to 350
       do 348 j = 1,n
       t = p*b(j,k) + q*b(j,k+1)
       b(j,k+1) =-q*b(j,k) + p*b(j,k+1)
 348   b(j,k) = t
 350   continue
       if(itcnt .le. itlim) go to 309
       int = 2
       go to 306
 999   continue
       nd = n-1
       do 800 k = 1,nd
       s = a(k,k)
       kr = k
       kp = kr+1
       do 710 j = kp,n
       if(a(j,j).le.s)  go to 710
       kr = j
       s = a(kr,kr)
 710   continue
       if(kr .eq.k)  go to 800
       a(kr,kr) = a(k,k)
       a(k,k) = s
       if (novec) go to 800
       do 725  j=1,n
       t = b(j,k)
       b(j,k) = b(j,kr)
 725   b(j,kr) = t
 800   continue
 
c     check convergence of roots:
      if(int.eq.2) then
         write(*,*) ' eigen values calculations did not converge'
         write(*,*) ' inverse maybe unreliable'
         return 1
      endif
 
c     copy eigenvalues:
      do 161 k=1,n1
 161  c(k)=a(k,k)
 
      npr=0
      nzr=0
      nnr=0
 
      do 171 i=1,n1
       if(c(i).gt.0.d0) then
          npr=npr+1
       elseif(c(i).lt.0.d0) then
          nnr=nnr+1
       else
          nzr=nzr+1
       endif
 171  continue
 
      rcond=-999.d0
 
      do 160 i=1,n1
      do 160 j=1,i
      a(i,j)=0.d0
      do 155 k=1,n1
      if(dabs(c(k)).le.eps2) goto 979
 155  a(i,j)=a(i,j)+b(i,k)*b(j,k)/c(k)
 160  a(j,i)=a(i,j)
 
c     now we know that roots are fine:
      rcond=c(1)/c(n1)
 
      return
 
 979  write(*,998) k,c(k)
 998  format('   at least one root smaller than machine accuracy:'/
     *       '   root(',i2,')=',d15.7/'   no inversion attempted')
       return 1
 
      end
      function iseednew(iseed0,icycle)
      integer*4 iseednew
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     vah 14:38:37.59, Mon, Mar 16, 1992
c
c     Starts from iseed0, and calls randu icycle times.
c     For the same iseed0 value, it will thus return a different seed
c     for every different value of icycle.
c
c     iseed0 and icycle are NOT changed.
c
c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
      double precision randu,temp
      integer*4 itemp,iseed0
 
      itemp=iseed0
      do 1 i=1,icycle
         temp=randu(itemp)
    1 continue
      iseednew=itemp
      return
      end
      subroutine lc(ch)
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     changes the character variable ch to lower case
c     vah 3:26 am, monday, october 27, 1986
c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
      integer*4 length
      character*(*) ch
      character*1 c1
      length=len(ch)
      do 1 i=1,length
      c1=ch(i:i)
      ic=ichar(c1)
c      print *,'ic:',ic
      if(ic.gt.64.and.ic.lt.91) ch(i:i)=char(ic+32)
 1    continue
      return
      end
      integer function ltrim(string)
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     vah 6:47 pm, saturday, february 27, 1988
      character*(*) string
c
c    return the blank-stripped length
c
c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
      do 10 l=len(string),1,-1
      ltrim=l
      if (string(l:l) .ne. ' ') return
10    continue
      return
      end
      subroutine machine(machnam)
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     vah 8:41:23.91, Tue, Mar 24, 1992
c
c     returns the value of the environment variable 'machine'
c     requires machine-dependent subprogram GETEVAR
c
c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
      character*(*) machnam
      call getevar('machine',machnam)
      return
      end
      subroutine matptv(a,lda,b,ldb,c,ldc,n,k,l)
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     c=a-transpose*b
c
c     a(lda,k) b(ldb,l) c(ldc,l)
c       nxk      nxl      kxl
c
c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
      double precision a(lda,k),b(ldb,l),c(ldc,l),hq
      do 1 i=1,k
         do 2 j=1,l
            hq=0.d00
            do 3 ijk=1,n
               hq=hq+a(ijk,i)*b(ijk,j)
    3       continue
            c(i,j)=hq
    2    continue
    1 continue
      return
      end
      subroutine matpv(a,lda,b,ldb,c,ldc,n,k,l)
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     c=a*b
c
c     a(lda,k) b(ldb,l) c(ldc,l)
c       nxk      kxl      nxl
c
c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
      double precision a(lda,k),b(ldb,l),c(ldc,l),hq
      do 1 i=1,n
         do 2 j=1,l
            hq=0.d00
            do 3 ijk=1,k
               hq=hq+a(i,ijk)*b(ijk,j)
    3       continue
            c(i,j)=hq
    2    continue
    1 continue
      return
      end
      subroutine mprint(a,lda,n,m,aname,nfile,ndig)
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     vah : note that aname is character*(*)
c
c      dbpm   trace print of matrix
c      prints the matrix a in table of n rows and m columns
c      lda is the no. of rows in dimension of a
c      aname is the name of array
c      nfile is the file unit number where the output should go
c      ndig=6,7 or 8 for  6 per line in d16.8
c      ndig=3,4 or 5 for  8 per line in d12.4
c      ndig=0,1 or 2 for 10 per line in d9.1
c
c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
      implicit double precision (a-h,o-z)
      character*(*) aname
      dimension a(lda,m),outst(10)
      write(nfile,1) aname,n,m
    1 format(/5x,' array ',a8,' (',i4,' by ',i4,' )'/)
      j1=1
c      max no. of columns
      jmax=max0(10-(ndig/3)*2,6)
      jj=min0(jmax,m)
  100 do 20 i=1,n
         do 21 j=j1,jj
   21       outst(j-j1+1)=a(i,j)
         if(jmax.ge.10)
     *   write(nfile,10)(outst(j-j1+1),j=j1,jj)
         if(jmax.ge.8.and.jmax.lt.10)
     *   write(nfile,8)(outst(j-j1+1),j=j1,jj)
         if(jmax.lt.8)
     *   write(nfile,6)(outst(j-j1+1),j=j1,jj)
    6    format(6 (1pd16.8))
    8    format(8 (1pd12.4))
   10    format(10(1pd 9.1))
c
   20 continue
      if(jj.ge.m)write(nfile,2)
      if(jj.ge.m)return
      j1=jj+1
      jj=min0(jj+jmax,m)
      if(n.gt.1)write(nfile,2)
    2 format(/)
      goto 100
      end
      double precision function nrmpdf(x,xm,sd)
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
c     vah  2:48 pm, friday, april 18, 1986                    c
c         11:38 pm, wednesday, september 24, 1986             c
c     msl  4:20 pm, thursday, october 6, 1988                 c
c     vah+msl, 5:40 pm, wed, dec 7, 1988 (pareto pdf)         c
c ----------------------------------------------------------- c
c ----------------------------------------------------------- c
c                                                             c
c     1.  double precision function nrmpdf(x,xm,sd)           c
c     2.  double precision function unfpdf(x,a,b)             c
c     3.  double precision function betpdf(arg,a,d)           c
c     4.  double precision function binpdf(x,t,p)             c
c     5.  double precision function caupdf(x,xloc,scale)      c
c     6.  double precision function combnr(xn,xr)             c
c     7.  double precision function exppdf(x,theta)           c
c     8.  double precision function gampdf(x,theta,t)         c
c     9.  double precision function parpdf(x,theta)           c
c    10.  double precision function permnr(xn,xr)             c
c    11.  double precision function dgamcp(x)                 c
c    12.  double precision function xfact(x)                  c
c    13.  double precision function logpdf(arg,xm,std)        c
c    14.  double precision function dexpdf(arg,xm,std)        c
c                                                             c
c ----------------------------------------------------------- c
      implicit double precision (a-h,o-z)
      implicit integer*4 (i-n)
      external sfee
c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
      arg=(x-xm)/sd
      nrmpdf=sfee(arg)/sd
      return
      end
      double precision function ph3del(x,y)
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     required for pheev3
c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
      implicit double precision(a-h,o-z)
      if(x*y)1,3,2
    1 ph3del=1.
      return
    2 ph3del=0.
      return
    3 if(x+y)1,2,2
      end
      double precision function ph3s(x,a,b,phix,txa)
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     required for pheev3
      implicit double precision(a-h,o-z)
c     computes steck's (1958 ann math stat) s-fn with max error e-9,
c     though computer rounding errors may reduce this accuracy. calls
c     function subprograms ph3stk(x,a,b) in which dabs(b) < 1.1,
c     pheev(x), and ph3t(x,a,phix) note that phix = pheev(x)
c     and txa = ph3t(x,a,phix).
      ba = dabs(b)
      if(ba.gt.1.1)go to 1
      ph3s=ph3stk(x,a,b)
      return
    1 ta=dabs(txa)
      aa=dabs(a)
      biv=1./ba
      xf=dabs(phix-.5)
      xa=dabs(x)
      xab=xa*aa*ba
      phixab=pheev(xab)
      xfa=phixab-.5
      if(aa.le.1.)go to 2
      ph3s=
     * xfa*ta-xf*ph3t(xab,biv,phixab)+.25*(xf-xfa)+ph3stk(xab,biv,1./aa)
      go to 3
    2 ab=aa*ba
      ph3s=.25*(xf+.5)+xfa*ta-ph3stk(xab,1./ab,aa)-ph3stk(xa,ab,biv)
    3 if(x.ge.0.)go to 4
      ph3s=.1591549*datan(ba/dsqrt(1.+a*a*(1.+b*b)))-s
    4 if(b.ge.0.)return
      ph3s=-ph3s
      return
      end
      double precision function ph3stk(x,a,b)
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     required for pheev3
      implicit double precision(a-h,o-z)
c        accurate to .0001
      z=1.+a*a
      z1=dsqrt(z+(a*.5*b)**2)
      z2=b/dsqrt(z+(a*b)**2)
      ph3stk=.1591549*pheev(x*z1)*datan(z2)
      return
      end
      double precision function ph3t(z,a,phix)
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     required for pheev3
      implicit double precision(a-h,o-z)
c        owens t double precision function accurate to .00005
c        phix is a dummy argument
      x1=z
      a1=a
      ff1=0.
      f1=1.
      ph3t=0.
      if(a .ge. 0.)go to 10
      a1=-a1
      f1=-1.
   10 if(a1 .le. 1.)go to 15
      x1=dabs(x1)
      gg1=pheev(x1)
      x1=a1*x1
      a1=1./a1
      gg2=pheev(x1)
      ff1=.5*(gg1+gg2)-gg1*gg2
      ff1=ff1*f1
      f1=-f1
   15 theta=datan(a1)
      x2a=x1*x1
      if(x2a .gt. 150.)go to 20
      ph3t=f1*theta*.15915494*dexp(-.5*x2a*a1/theta)
      ph3t=ph3t*(1.+.00868*(x1*a1)**4)
      ph3t=ph3t+ff1
   20 return
      end
      double precision function pheev(y)
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     this essentially uses the imsl real*8 function derf, provided
c     by clint cummins for use in quail5.1 . the modification using
c     the fact that phee(z)=0.5+0.5*erf(z/sqr2) was done by vah .
c
c                                  specifications for arguments
c       sqr2, pt5 added by vah
c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
      double precision   derf,y,sqr2,pt5
c                                  specifications for local variables
      dimension          p(5),q(4),p1(9),q1(8),p2(6),q2(5)
      double precision   p,q,p1,q1,p2,q2,xmin,xlarge,sqrpi,x,
     *                   res,xsq,xnum,xden,xi
      integer            isw,i
c     added by vah
      data sqr2   /1.41421356237309d00/,pt5/0.5d0/
c                                  coefficients for 0.0 .le. y .lt.
c                                  .477
      data               p(1)/113.8641541510502d0/,
     *                   p(2)/377.4852376853020d0/,
     *                   p(3)/3209.377589138469d0/,
     *                   p(4)/.1857777061846032d0/,
     *                   p(5)/3.161123743870566d0/
      data               q(1)/244.0246379344442d0/,
     *                   q(2)/1282.616526077372d0/,
     *                   q(3)/2844.236833439171d0/,
     *                   q(4)/23.60129095234412d0/
c                                  coefficients for .477 .le. y
c                                  .le. 4.0
      data               p1(1)/8.883149794388376d0/,
     *                   p1(2)/66.11919063714163d0/,
     *                   p1(3)/298.6351381974001d0/,
     *                   p1(4)/881.9522212417691d0/,
     *                   p1(5)/1712.047612634071d0/,
     *                   p1(6)/2051.078377826071d0/,
     *                   p1(7)/1230.339354797997d0/,
     *                   p1(8)/2.153115354744038d-8/,
     *                   p1(9)/.5641884969886701d0/
      data               q1(1)/117.6939508913125d0/,
     *                   q1(2)/537.1811018620099d0/,
     *                   q1(3)/1621.389574566690d0/,
     *                   q1(4)/3290.799235733460d0/,
     *                   q1(5)/4362.619090143247d0/,
     *                   q1(6)/3439.367674143722d0/,
     *                   q1(7)/1230.339354803749d0/,
     *                   q1(8)/15.74492611070983d0/
c                                  coefficients for 4.0 .lt. y
      data               p2(1)/-3.603448999498044d-01/,
     *                   p2(2)/-1.257817261112292d-01/,
     *                   p2(3)/-1.608378514874228d-02/,
     *                   p2(4)/-6.587491615298378d-04/,
     *                   p2(5)/-1.631538713730210d-02/,
     *                   p2(6)/-3.053266349612323d-01/
      data               q2(1)/1.872952849923460d0/,
     *                   q2(2)/5.279051029514284d-01/,
     *                   q2(3)/6.051834131244132d-02/,
     *                   q2(4)/2.335204976268692d-03/,
     *                   q2(5)/2.568520192289822d0/
c                                  constants
      data               xmin/1.0d-10/,xlarge/6.375d0/
      data               sqrpi/.5641895835477563d0/
c                                  first executable statement
c       modified by vah
      x = y/sqr2
      isw = 1
      if (x.ge.0.0d0) go to 5
      isw = -1
      x = -x
    5 if (x.lt..477d0) go to 10
      if (x.le.4.0d0) go to 25
      if (x.lt.xlarge) go to 35
      res = 1.0d0
      go to 50
c                                  abs(y) .lt. .477, evaluate
c                                  approximation for derf
   10 if (x.lt.xmin) go to 20
      xsq = x*x
      xnum = p(4)*xsq+p(5)
      xden = xsq+q(4)
      do 15 i=1,3
         xnum = xnum*xsq+p(i)
         xden = xden*xsq+q(i)
   15 continue
      res = x*xnum/xden
      go to 50
   20 res = x*p(3)/q(3)
      go to 50
c                                  .477 .le. abs(y) .le. 4.0
c                                  evaluate approximation for derf
   25 xsq = x*x
      xnum = p1(8)*x+p1(9)
      xden = x+q1(8)
      do 30 i=1,7
         xnum = xnum*x+p1(i)
         xden = xden*x+q1(i)
   30 continue
      res = xnum/xden
      go to 45
c                                  4.0 .lt. abs(y), evaluate
c                                  minimax approximation for derf
   35 xsq = x*x
      xi = 1.0d0/xsq
      xnum = p2(5)*xi+p2(6)
      xden = xi+q2(5)
      do 40 i=1,4
         xnum = xnum*xi+p2(i)
         xden = xden*xi+q2(i)
   40 continue
      res = (sqrpi+xi*xnum/xden)/x
   45 res = res*dexp(-xsq)
      res = 1.0d0-res
   50 if (isw.eq.-1) res = -res
      derf = res
c       added by vah
      pheev=pt5+pt5*derf
      return
      end
c
c
c
      double precision function pheev2(z1,z2,r)
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  program to compute bivariate normal cdf
c  cacm algorith number 462
c  based on algorithm in owen, annals of math stat, 1956:1075-90
c  accuracy is controlled by idig--number of sig digits to right of
c  decimal point in answer.  more accurate than old hausman routine
c  with idig=15 time is approximately .05 seconds on pc w/8087
c  time varies with values.  can be faster
c
c  modified by vah 6:49 pm, tuesday, december 9, 1986
c
c  routines used:
c      pheev(z) is standard cumulative normal at z
c
c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
      implicit double precision(a-h,o-z)
      implicit integer*2(i-n)
      external pheev
ccc      data twopi /6.2831853071795865770d00/
ccc      data pt5,pt25,zero,one,two /.5d0,.25d0,0.d0,1.d0,2.d0/
 
      if(r .eq. 0.d0)then
            pheev2=pheev(z1)*pheev(z2)
            return
      endif
      sqr=1.d0-r**2
      if(sqr .le. 0.d0)then
            write(*,801)
  801       format(' rho on boundary or outside range: rho=',d13.5)
            stop
      endif
      sqr=dsqrt(sqr)
      ah=-z1
      ak=-z2
      twopi=6.283185307179587d0
      b=0.d0
ccc      idig=15
      gh=.5d0*pheev(-ah)
      gk=.5d0*pheev(-ak)
ccc      if(idig .ne. 15)then
ccc            c0n=.5d0*twopi/(10.d0**idig)
ccc      else
ccc            c0n=twopi*.5d-15
            c0n=twopi*.5d-13
ccc      endif
      if(ah) 170,150,170
  150 if(ak) 190,160,190
  160 b=datan(r/sqr)/twopi+.25d0
      go to 350
  170 b=gh
      if(ah*ak) 180,200,190
  180 b=b-.5d0
  190 b=b+gk
      if(ah)200,340,200
  200 wh=-ah
      wk=(ak/ah-r)/sqr
      gw=2.d0*gh
      is=-1
  210 sgn=-1.d0
      t=0.d0
      if(wk) 220,320,220
  220 if(dabs(wk)-1.d0)270,230,240
  230 t=wk*gw*(1.d0-gw)*.5d0
      go to 310
  240 sgn=-sgn
      wh=wh*wk
      g2=pheev(wh)
      wk=1.d0/wk
      if(wk)250,260,260
  250 b=b+.5d0
  260 b=b-(gw+g2)*.5d0+gw*g2
  270 h2=wh*wh
      a2=wk*wk
      h4=h2*.5d0
      ex=dexp(-h4)
      w2=h4*ex
      ap=1.d0
      s2=ap-ex
      sp=ap
      s1=0.d0
      sn=s1
      c0nex=dabs(c0n/wk)
      go to 290
  280 sn=sp
      sp=sp+1.d0
      s2=s2-w2
      w2=w2*h4/sp
      ap=-ap*a2
  290 cn=ap*s2/(sn+sp)
      s1=s1+cn
      if(dabs(cn)-c0nex)300,300,280
  300 t=(datan(wk)-wk*s1)/twopi
  310 b=b+sgn*t
  320 if(is)330,350,350
  330 if(ak)340,350,340
  340 wh=-ak
      wk=(ah/ak-r)/sqr
      gw=2.d0*gk
      is=1
      go to 210
  350 if(b .lt. 0.d0)then
            pheev2=0.d0
      elseif(b .gt. 1.d0)then
            pheev2=1.d0
      else
            pheev2=b
      endif
      return
      end
      subroutine pheev3(z,x1,x2,x3,rho12,rho31,rho23,*)
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     vah 12:55 am, sunday, may 22, 1988
c         based on farber's program tnorm
c
c     for the trivariate normal r.v. (x,y,w) with zero means, unit
c     variance and corr(x,y)=rho12,corr(y,w)rho23,corr(w,x)=rho31,tnor
c     uses the method of section 2 of steck (1958) to compute z= pr
c     (x<x1,y<x2,w<x3) accuracy of the result depends on the
c     subprograms ph3t and ph3stk. for more detail re accuracy see daley
c     (1973) anu stat dept (ias) rep. ref.: g.p. steck (1958) ann math
c     statis vol 29,pp 780-800.
c        subroutine ph3t and ph3stk use closed form approximations
c        accur to about .00025,  time is about .003
c
c     requires routines:
c      1. f3swap
c      2. ph3del
c      3. ph3t
c      4. ph3s
c      5. ph3stk
c      6. pheev
c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
      implicit double precision(a-h,o-z)
      implicit integer*4(i-n)
      dimension x(6),xiv(3),xr(6),phix(6),r(6)
      dimension y(6),yr(3),yrt(6)
      dimension a(6),b(6),txa(6),sxab(6)
      x(1)=x1
      x(2)=x2
      x(3)=x3
      do 8 i=1,3
      if(dabs(x(i)).ge.1.e-15)go to 8
      x(i)=1.e-15
    8 continue
      r(1)=rho12
      r(2)=rho23
      r(3)=rho31
      ind=0
c     ind=0 for all x'x of same sign, =1 otherwise. if x's are of
c     different signs, then permute to have x(3) of diff't sign
c     and change x(3),r(2), and r(3). signs.
      if(x(1))2,1,2
    1 if(x(2)*x(3))7,10,10
    2 if(x(1)*x(2))4,3,3
    3 if(x(1)*x(3))7,10,10
    4 if(x(1)*x(3))5,6,6
    5 call f3swap(x(1),x(3))
      call f3swap(r(1),r(2))
      go to 7
    6 call f3swap(x(2),x(3))
      call f3swap(r(1),r(3))
    7 ind=1
      x(3)=-x(3)
      r(2)=-r(2)
      r(3)=-r(3)
   10 do 11 i=1,3
   11 y(i)=1.-r(i)*r(i)
      det=y(1)+y(2)+y(3)-2.+2.*r(1)*r(2)*r(3)
      if(det.gt.1.e-15)go to 13
c     if(det.lt.-1.e-4)go to 31
      if(det .lt. -1.e-4)return1
      rr=.999999
   13 det=1./dsqrt(det)
      z=0.
      yrt(1)=1./dsqrt(y(1))
      yrt(2)=1./dsqrt(y(2))
      yrt(3)=1./dsqrt(y(3))
      do 16 i=1,3
      j=i+3
      x(j)=x(i)
      xiv(i)=1./x(i)
      phix(i)=pheev(x(i))
      phix(j)=phix(i)
      y(j)=y(i)
      yrt(j)=yrt(i)
      r(j)=r(i)
      yr(i)=r(i)*r(i+2)-r(i+1)
      xr(i)=x(i+1)-x(i)*r(i)
      if(dabs(xr(i)).gt.1.e-25)go to 14
      xr(i)=1.e-25
   14 xr(j)=x(i+2)-x(i)*r(i+2)
      if(dabs(xr(j)).ge.1.e-25)go to 15
      xr(j)=1.e-25
   15 a(i)=xr(i)*xiv(i)*yrt(i)
      a(j)=xr(j)*xiv(i)*yrt(i+2)
      b(i)=det*(yr(i)+y(i)*xr(j)/xr(i))
      b(j)=det*(yr(i)+y(i+2)*xr(i)/xr(j))
      z=z+(1.-ph3del(a(i),a(j)))*phix(i)
   16 continue
      z=.5*z
      do 17 i=1,6
      txa(i)=ph3t(x(i),a(i),phix(i))
      sxab(i)=ph3s(x(i),a(i),b(i),phix(i),txa(i))
      z=z-.5*txa(i)-sxab(i)
   17 continue
      if(ind.eq.0)return
      z=.5*(phix(1)+phix(5)-ph3del(x(1),x(2)))-txa(1)-txa(5)-z
      return
      end
      double precision function phinv(ppp)
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     28 mar 86
c     12:33 pm, monday, june 30, 1986
c*****inverse normal c.d.f. (pheev(ppp)=phinv) using
c*****formula 26.2.23 of abramovitz and stegun p.933 10th edition.
c*****ppp unchanged on exit. accuracy<4.5e-4
c*****
c*****written by v.a. hajivassiliou 24sep84
c*****
c*****uses inverse interpolation to 3rd order so as to
c*****raise accuracy-see a&s example 5 p.954
c*****
c*****modification done 14apr85
c*****
c*****subprogram(s) used: pheev
c*****
c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
      implicit double precision (a-h,o-z)
      external pheev
      data pt5neg/-0.5d0/,pt5/0.5d0/,
     *     zero  / 0.d0 /,one/1.d0 /,two/2.d0/,six/6.d0/
      data c0,c1,c2,d1,d2,d3/2.515517d0,.802853d0,.010328d0,
     *                       1.432788d0,.189269d0,.001308d0/
      data osr2pi/0.39894228040143267793d00/
      iflip=0
      p=ppp
      if(ppp.eq.one) go to 998
      if(ppp.eq.zero) go to 999
      if(ppp.gt.pt5) then
       p=one-ppp
       iflip=1
      end if
c*****
c     p<0.5 as required by inversion:
c*****
      t=dsqrt(dlog(one/p**2))
      x0=t-(c0+c1*t+c2*t**2)/(one+d1*t+d2*t**2+d3*t**3)
c*****
c     up to now x0<0:
c*****
      x0=-x0
c*****
c     now x0>0 as required by inverse interpolation:
c*****
      sfeex0=osr2pi*dexp(pt5neg*x0**2)
      t=(p-pheev(x0))/sfeex0
      ph=x0+t+x0*t**2/two+(two*x0**2+one)*t**3/six
c*****
c     back to original x0<0, ph<0:
c*****
      ph=-ph
      if(iflip.eq.0) ph=-ph
      phinv=ph
      return
 998  phinv=18.d0
c*****a crude way of dealing with ppp=1
      return
 999  phinv=-18.d0
c*****a crude way of dealing with ppp=0
      return
      end
c
c
c
      double precision function randu(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     vah : 28 mar 86
c     1:07 pm, monday, june 30, 1986
c*****
c*****this is a pseudo random number generator,
c*****randu being the realization of a uniform on [0,1]
c*****pseudo r.v.. requires a seed (iseed) which is updated
c*****and returned. written by v.a. hajivassiliou, 24sep84,
c*****to exactly reproduce the pr1me uniform pseudo random
c*****generator rand$a(iseed).
c*****ref:
c*****lewis,goodman & miller,ibm systems journal,v.8#2 1969,pp.136-145
c*****
c*****uses multiplier number 1 for fishman and moore, "a statistical
c*****evaluation of multiplicative congruential random number generators
c*****with modulus 2**31-1", jasa 77, march 1982.  this was found to
c*****be quite acceptable in the tests, and the fastest.
c
c*****
c*****note: the imsl function ggubs is exactly the same algorithm,
c*****      except that it divides by 2147483648 , potentially
c*****      speeding the process up on 32 bit machines
c*****
ccc can't re-mention randu in rm profort:
ccc      double precision randu,xm,xb,seed,xnews
c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
      implicit integer*4 (i-n)
      double precision xm,xb,seed,xnews
      integer*4 iseed
      data xb,xm/16807.d0,2147483647.d0/
      seed=dble(iseed)
      xnews=dmod(xb*seed,xm)
      randu=xnews/xm
ccc      iseed=intl(xnews)
ccc      for rmpf:
      iseed=int(xnews)
ccc      iseed=jfix(xnews)
      return
      end
      double precision function ranorm(xm,sd,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*****           random number generators package
c*****           vassilis hajivassiliou october 1984
c***** last modified : 5:45 pm, wednesday, september 24, 1986
c**********************************************************************
c     this set contains:
c     *****************
c     1. ranorm(xm,sd,iseed)         : (tsp)
c     2. binrmb(x,y,xm,ym,xs,ys,rho,iseed)
c     3. birnrm(x,y,xm,ym,xs,ys,rho,iseed)
c     4. phinv(p)
c     5. randb(prob,iseed)
c     6. randc(xm,scale,iseed)
c     7. randcv(xm,scale,iseed)
c     8. rande(xm,iseed)
c     9. randl(xm,sd,iseed)
c    10. randp(xlambda,iseed)
c    11. randpa(xtheta,iseed)
c    12. randu(iseed)
c    13. ranorv(xm,sd,iseed)        : (using phinv)
c    14. ranrmb(xm,sd,iseed)
c    15. randde(xm,sd,iseed)
c
c**********************************************************************
      implicit double precision(a-h,o-z)
      integer*4 iseed
c**********************************************************************
c
c     normal(xm,sd) random number generator . vassilis hajivassiliou
c
c     october 1984
c
c     this is from cacm algorithm #488, originally implemented
c     for use in tsp4.0 by  bronwyn h. hall(oct82).
c
c**********************************************************************
c     subprogram(s) used: randu
c**********************************************************************
c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
      dimension d(60),d1(30),d2(30)
      equivalence(d1(1),d(1)),(d2(1),d(31))
c
      data d1/
     *    0.674489750d0,0.475859630d0,0.383771164d0,0.328611323d0,
     *    0.291142827d0,0.263684322d0,0.242508452d0,0.225567444d0,
     *    0.211634166d0,0.199924267d0,0.189910758d0,0.181225181d0,
     *    0.173601400d0,0.166841909d0,0.160796729d0,0.155349717d0,
     *    0.150409384d0,0.145902577d0,0.141770033d0,0.137963174d0,
     *    0.134441762d0,0.131172150d0,0.128125965d0,0.125279090d0,
     *    0.122610883d0,0.120103560d0,0.117741707d0,0.115511892d0,
     *    0.113402349d0,0.111402720d0/
      data d2/
     *    0.109503852d0,0.107697617d0,0.105976772d0,0.104334841d0,
     *    0.102766012d0,0.101265052d0,0.099827234d0,0.098448282d0,
     *    0.097124309d0,0.095851778d0,0.094627461d0,0.093448407d0,
     *    0.092311909d0,0.091215482d0,0.090156838d0,0.089133867d0,
     *    0.088144619d0,0.087187293d0,0.086260215d0,0.085361834d0,
     *    0.084490706d0,0.083645487d0,0.082824924d0,0.082027847d0,
     *    0.081253162d0,0.080499844d0,0.079766932d0,0.079053527d0,
     *    0.078358781d0,0.077681899d0/
c
      u = randu(iseed)
        a = 0.0d0
        i = 0
   10   continue
        if (u.eq.1.0d0) u = randu(iseed)
           u = u+u
        if (u.lt.1.0d0) go to 20
           u = u-1.0d0
           i = i+1
           a = a-d(i)
           go to 10
   20    continue
         w = d(i+1)*u
         v = w*(0.5d0*w-a)
   30    continue
         u = randu(iseed)
         if (v.le.u) go to 40
         v = randu(iseed)
         if (u.gt.v) go to 30
         u = (v-u)/(1.0d0-u)
         go to 20
   40    continue
         u = (u-v)/(1.0d0-v)
         u = u+u
         if (u.lt.1.0d0) go to 50
            u = u-1.0d0
            rn = w-a
            go to 100
   50    continue
         rn = a-w
  100 continue
      ranorm=xm+sd*rn
      return
      end
      double precision function sfee(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+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
      implicit double precision(a-h,o-z)
      data osr2pi/0.39894228040143267793d00/,pt5neg/-0.5d0/
      sfee=0.d0
      xsq=x*x
      if(xsq.gt.1000.d0) return
      sfee=osr2pi*dexp(pt5neg*xsq)
      return
      end
      subroutine sorts(y,n)
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
ccc      subroutine sorts(y,n,ierr)
c
ccc      integer*4 n,ierr
      integer*4 n
      double precision y(n)
c
c     shell sort n values in t() from smallest to largest.
c     *****
c
c     note that local system sort utilities are likely to be
c     more efficient, and should be substituted whenever possible.
c
c     "abc's of eda" manual 1981 p.313
c     written 11/03/85 17:17:29
c
c     local variables
c
ccc      integer*2 i, j, j1, igap, nmg
c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
      integer*4 i, j, j1, igap, nmg
      double precision temp
c
      if(n .ge. 1) goto 10
ccc         ierr=1
         goto 999
 10   if(n .eq.1) goto 999
c
c     one element is always sorted
c
      igap=n
 20   igap=igap/2
      nmg=n-igap
      do 40 j1=1,nmg
         i=j1+igap
c
c     do j=j1,1,-igap
c
         j=j1
 30   if(y(j) .le. y(i)) goto 40
c
c     swap out-of-order pair
c
      temp=y(i)
      y(i)=y(j)
      y(j)=temp
c
c     keep old pointer for next time through
c
      i=j
      j=j-igap
      if(j .ge. 1) goto 30
 40   continue
      if(igap .gt. 1) goto 20
 999  return
      end
c
c
c
      subroutine srtrnk(yvec,n,ivec)
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     vah 8:11 pm, tuesday, october 28, 1986
c         16:50:51.97, sun, apr 9, 1989
c
      implicit double precision(a-h,o-z)
      implicit integer (i-n)
c
c     shell sort n values in yvec() from smallest to largest.
c     *****
c     vah 7:29 pm, tuesday, october 28, 1986
c         modified to return also array ivec.
c     on entry, ivec(1) controls the operations to be performed:
c     if ivec(1) ge 0: yvec is sorted from min to max (call this xsort)
c                      ivec will contain the original locations of the
c                      elements of yvec (call this iloc)
c     if ivec(1) lt 0: yvec is sorted from min to max in place, i.e.
c                      the original yvec is returned
c                      ivec will contain the ranks of the yvec
c                      (call this irnk)
c-----------------------------------------------------------------------
c     the full relations are:
c
c            information returned with       information returned with
c  original         ivec(1) ge 0                 with ivec(1) lt 0
c  --------  -------------------------       -------------------------
c  yvec      yvec      ivec                  yvec      ivec
c  ----      ----      ----                  ----      ----
c  x         xsort     iloc                  x         irnk
c  -         -----     ----                  -         ----
c  2         1         4                     2         2
c  4         2         1                     4         3
c  7         4         2                     7         5
c  1         5         5                     1         1
c  5         7         3                     5         4
c
c            to obtain x and irnk            to obtain xsort and iloc
c            from xsort, iloc                from x, irnk
c            --------------------            ------------------------
c            x(iloc(j)=xsort(j)              xsort(irnk(j))=x(j)
c            irnk(iloc(j))=j                 iloc(irnk(j))=j
c
c-----------------------------------------------------------------------
c
c     note that local system sort utilities are likely to be
c     more efficient, and should be substituted whenever possible.
c
c     "abc's of eda" manual 1981 p.313
c     written 11/03/85 17:17:29
c
c     local variables
c
ccc      integer*2 i, j, j1, igap, nmg
c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
      dimension yvec(n),ivec(n)
 
      if(ivec(1).lt.0) then
            isort=-1
      else
            isort=1
      endif
c
      if(n .ge. 1) goto 10
ccc         ierr=-1
         ivec(1)=1
         goto 999
 10   if(n .eq.1) goto 999
c
c     one element is always sorted
c
      do 555 i=1,n
 555  ivec(i)=i
      igap=n
 20   igap=igap/2
      nmg=n-igap
      do 40 j1=1,nmg
         i=j1+igap
c
c     do j=j1,1,-igap
c
         j=j1
 30   if(yvec(j) .le. yvec(i)) goto 40
c
c     swap out-of-order pair
c
      temp=yvec(i)
      itemp=ivec(i)
      yvec(i)=yvec(j)
      ivec(i)=ivec(j)
      yvec(j)=temp
      ivec(j)=itemp
c
c     keep old pointer for next time through
c
      i=j
      j=j-igap
      if(j .ge. 1) goto 30
 40   continue
      if(igap .gt. 1) goto 20
ccc      do 777 j=1,n
ccc      irnk(ivec(j))=j
ccc 777  continue
 
      if(isort.le.0) then
            next0=1
            inext=ivec(1)
            temp0=yvec(1)
 887        continue
            do 888 jjj=1,n
 
c      print *,'switching',inext
                  temp1=yvec(inext)
                  yvec(inext)=temp0
                  temp0=temp1
 
                  next=ivec(inext)
                  ivec(inext)=-next0
                  next0=inext
                  inext=next
 
                  if(inext.lt.0) then
c      print *,'inext lt 0 next0',inext,next0
                        do 886 iii=1,n
                        if(ivec(iii).gt.0) then
                               inext=ivec(iii)
                               next0=iii
c      print *,'restarting with inext,next0:',inext,next0
                               goto 887
                        endif
 886                    continue
                        goto 889
                  endif
 888        continue
 889        continue
            do 884 iii=1,n
                  ivec(iii)=-ivec(iii)
 884        continue
      endif
 
c      print *,'ivec',(ivec(kkk),kkk=1,n)
c      print *,'yvec',(yvec(kkk),kkk=1,n)
 999  return
      end
      subroutine uc(ch)
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     changes the character variable ch to upper case
c     vah 5:13 pm, thursday, april 3, 1986
c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
      integer*4 length
      character*(*) ch
      character*1 c1
      length=len(ch)
      do 1 i=1,length
      c1=ch(i:i)
      ic=ichar(c1)
c      print *,'ic:',ic
      if(ic.gt.96.and.ic.lt.123) ch(i:i)=char(ic-32)
 1    continue
      return
      end
      subroutine when(nfunit)
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+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
      integer*4 nfunit
      integer*4 itim(5),idat(3)
      character*3 cm(12)
      cm(1)='jan'
      cm(2)='feb'
      cm(3)='mar'
      cm(4)='apr'
      cm(5)='may'
      cm(6)='jun'
      cm(7)='jul'
      cm(8)='aug'
      cm(9)='sep'
      cm(10)='oct'
      cm(11)='nov'
      cm(12)='dec'
      call timdat(itim,idat)
      if(idat(3).lt.100) idat(3)=idat(3)+1900
      write(nfunit,100) idat(1),cm(idat(2)),idat(3),(itim(i),i=1,4)
 100  format('0','date:',i4,' ',a3,' ',i4,2x,
     *              'time:',2x,i2.2,':',i2.2,':',i2.2,'.',i2.2/
     *           1x,'----', 16x,'----')
      return
      end
