anchor.diff:  c
   c     Subprograms called:
   c
  -c     Changed to LAPACK and BLAS ROUTINES ... dcopy, dpotrf, 
   c       Linpack ... dcopy, dpofa, dtrsl.
   c
   c
  
positive.code:          p(i2) = v(i2) + sum
    20  continue  
  c     Solve the triangular system
  c      call dtrsl(wt,m,col,p(col+1),11,info)
        call dpotrs('L',col,1,wt,m,p(col+1),col,info)
        if (info .ne. 0) return
   
  c     	solve D^(1/2)p1=v1.
  
negative.code:    70 return
  
        end
  
  c====================== The end of dpmeps ==============================
  
        subroutine dpofa(a,lda,n,info)
        integer lda,n,info
        double precision a(lda,1)
  c
  c     dpofa factors a double precision symmetric positive definite
  c     matrix.
  c
  c     dpofa is usually called by dpoco, but it can be called
  c     directly with a saving in time if  rcond  is not needed.
  c     (time for dpoco) = (1 + 18/n)*(time for dpofa) .
  c
  c     on entry
  c
  c        a       double precision(lda, n)
  c                the symmetric matrix to be factored.  only the
  c                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       an upper triangular matrix  r  so that  a = trans(r)*r
  c                where  trans(r)  is the transpose.
  c                the strict lower triangle is unaltered.
  c                if  info .ne. 0 , the factorization is not complete.
  c
  c        info    integer
  c                = 0  for normal return.
  c                = k  signals an error condition.  the leading minor
  c                     of order  k  is not positive definite.
  c
  c     This is just a wrapper that calls LAPACK, but with the LINPACK
  c     calling convention.
  
        call dpotrf('U', n, a, lda, info)
        end
  
  c====================== The end of dpofa ===============================
  
        subroutine dtrsl(t, ldt, n, b, job, info)
        integer ldt, n, job, info
        double precision t(ldt,1), b(1)
  c
  c
  c     dtrsl solves systems of the form
  c
  c                   t * x = b
  c     or
  c                   trans(t) * x = b
  c
  c     where t is a triangular matrix of order n. here trans(t)
  c     denotes the transpose of the matrix t.
  c
  c     on entry
  c
  c         t         double precision(ldt,n)
  c                   t contains the matrix of the system. the zero
  c                   elements of the matrix are not referenced, and
  c                   the corresponding elements of the array can be
  c                   used to store other information.
  c
  c         ldt       integer
  c                   ldt is the leading dimension of the array t.
  c
  c         n         integer
  c                   n is the order of the system.
  c
  c         b         double precision(n).
  c                   b contains the right hand side of the system.
  c
  c         job       integer
  c                   job specifies what kind of system is to be solved.
  c                   if job is
  c
  c                        00   solve t*x=b, t lower triangular,
  c                        01   solve t*x=b, t upper triangular,
  c                        10   solve trans(t)*x=b, t lower triangular,
  c                        11   solve trans(t)*x=b, t upper triangular.
  c
  c     on return
  c
  c         b         b contains the solution, if info .eq. 0.
  c                   otherwise b is unaltered.
  c
  c         info      integer
  c                   info contains zero if the system is nonsingular.
  c                   otherwise info contains the index of
  c                   the first zero diagonal element of t.
  c
  c     This is just a wrapper that calls LAPACK, but with the LINPACK
  c     calling convention.
  
        character*1 uplo, trans
  
        if (job .eq. 00) then
            uplo = 'L'
            trans = 'N'
        else if (job .eq. 01) then
            uplo = 'U'
            trans = 'N'
        else if (job .eq. 10) then
            uplo = 'L'
            trans = 'T'
        else if (job .eq. 11) then
            uplo = 'U'
            trans = 'T'
        endif
        call dtrtrs(uplo, trans, 'N', n, 1, t, ldt, b, n, info)
        end
  c====================== The end of dtrsl ==============================
  