* ======================================================================
* nist guide to available math software.
* fullsource for module dpodi from package cmlib.
* retrieved from camsun on wed jul 19 18:23:00 1995.
* ======================================================================
      subroutine dpodi(a,lda,n,det,job)
c***begin prologue  dpodi
c***date written   780814   (yymmdd)
c***revision date  820801   (yymmdd)
c***category no.  d2b1b,d3b1b
c***keywords  determinant,double precision,factor,inverse,
c             linear algebra,linpack,matrix,positive definite
c***author  moler, c. b., (u. of new mexico)
c***purpose  computes the determinant and inverse of a certain double
c            precision symmetric positive definite matrix (see abstract)
c            using the factors computed by dpoco, dpofa or dqrdc.
c***description

c     dpodi computes the determinant and inverse of a certain
c     double precision symmetric positive definite matrix (see below)
c     using the factors computed by dpoco, dpofa or dqrdc.

c     on entry

c        a       double precision(lda, n)
c                the output  a  from dpoco or dpofa
c                or the output  x  from dqrdc.

c        lda     integer
c                the leading dimension of the array  a .

c        n       integer
c                the order of the matrix  a .

c        job     integer
c                = 11   both determinant and inverse.
c                = 01   inverse only.
c                = 10   determinant only.

c     on return

c        a       if dpoco or dpofa was used to factor  a , then
c                dpodi produces the upper half of inverse(a) .
c                if dqrdc was used to decompose  x , then
c                dpodi produces the upper half of inverse(trans(x)*x)
c                where trans(x) is the transpose.
c                elements of  a  below the diagonal are unchanged.
c                if the units digit of job is zero,  a  is unchanged.

c        det     double precision(2)
c                determinant of  a  or of  trans(x)*x  if requested.
c                otherwise not referenced.
c                determinant = det(1) * 10.0**det(2)
c                with  1.0 .le. det(1) .lt. 10.0
c                or  det(1) .eq. 0.0 .

c     error condition

c        a division by zero will occur if the input factor contains
c        a zero on the diagonal and the inverse is requested.
c        it will not occur if the subroutines are called correctly
c        and if dpoco or dpofa has set info .eq. 0 .

c     linpack.  this version dated 08/14/78 .
c     cleve moler, university of new mexico, argonne national lab.

c     subroutines and functions

c     blas daxpy,dscal
c     fortran mod
c***references  dongarra j.j., bunch j.r., moler c.b., stewart g.w.,
c                 *linpack users  guide*, siam, 1979.
c***routines called  daxpy,dscal
c***end prologue  dpodi
      integer lda,n,job
      double precision a(lda,1)
      double precision det(2)

      double precision t
      double precision s
      integer i,j,jm1,k,kp1

c     compute determinant

c***first executable statement  dpodi
      if (job/10 .eq. 0) go to 70
         det(1) = 1.0d0
         det(2) = 0.0d0
         s = 10.0d0
         do 50 i = 1, n
            det(1) = a(i,i)**2*det(1)
c        ...exit
            if (det(1) .eq. 0.0d0) go to 60
   10       if (det(1) .ge. 1.0d0) go to 20
               det(1) = s*det(1)
               det(2) = det(2) - 1.0d0
            go to 10
   20       continue
   30       if (det(1) .lt. s) go to 40
               det(1) = det(1)/s
               det(2) = det(2) + 1.0d0
            go to 30
   40       continue
   50    continue
   60    continue
   70 continue

c     compute inverse(r)

      if (mod(job,10) .eq. 0) go to 140
         do 100 k = 1, n
            a(k,k) = 1.0d0/a(k,k)
            t = -a(k,k)
            call dscal(k-1,t,a(1,k),1)
            kp1 = k + 1
            if (n .lt. kp1) go to 90
            do 80 j = kp1, n
               t = a(k,j)
               a(k,j) = 0.0d0
               call daxpy(k,t,a(1,k),1,a(1,j),1)
   80       continue
   90       continue
  100    continue

c        form  inverse(r) * trans(inverse(r))

         do 130 j = 1, n
            jm1 = j - 1
            if (jm1 .lt. 1) go to 120
            do 110 k = 1, jm1
               t = a(k,j)
               call daxpy(k,t,a(1,j),1,a(1,k),1)
  110       continue
  120       continue
            t = a(j,j)
            call dscal(j,t,a(1,j),1)
  130    continue
  140 continue
      return
      end

      subroutine dscal(n,da,dx,incx)
c***begin prologue  dscal
c***date written   791001   (yymmdd)
c***revision date  820801   (yymmdd)
c***category no.  d1a6
c***keywords  blas,linear algebra,scale,vector
c***author  lawson, c. l., (jpl)
c           hanson, r. j., (snla)
c           kincaid, d. r., (u. of texas)
c           krogh, f. t., (jpl)
c***purpose  d.p. vector scale x = a*x
c***description

c                b l a s  subprogram
c    description of parameters

c     --input--
c        n  number of elements in input vector(s)
c       da  double precision scale factor
c       dx  double precision vector with n elements
c     incx  storage spacing between elements of dx

c     --output--
c       dx  double precision result (unchanged if n.le.0)

c     replace double precision dx by double precision da*dx.
c     for i = 0 to n-1, replace dx(1+i*incx) with  da * dx(1+i*incx)
c***references  lawson c.l., hanson r.j., kincaid d.r., krogh f.t.,
c                 *basic linear algebra subprograms for fortran usage*,
c                 algorithm no. 539, transactions on mathematical
c                 software, volume 5, number 3, september 1979, 308-323
c***routines called  (none)
c***end prologue  dscal

      double precision da,dx(1)
c***first executable statement  dscal
      if(n.le.0)return
      if(incx.eq.1)goto 20

c        code for increments not equal to 1.

      ns = n*incx
          do 10 i = 1,ns,incx
          dx(i) = da*dx(i)
   10     continue
      return

c        code for increments equal to 1.


c        clean-up loop so remaining vector length is a multiple of 5.

   20 m = mod(n,5)
      if( m .eq. 0 ) go to 40
      do 30 i = 1,m
        dx(i) = da*dx(i)
   30 continue
      if( n .lt. 5 ) return
   40 mp1 = m + 1
      do 50 i = mp1,n,5
        dx(i) = da*dx(i)
        dx(i + 1) = da*dx(i + 1)
        dx(i + 2) = da*dx(i + 2)
        dx(i + 3) = da*dx(i + 3)
        dx(i + 4) = da*dx(i + 4)
   50 continue
      return
      end

      subroutine dpofa(a,lda,n,info)
c***begin prologue  dpofa
c***date written   780814   (yymmdd)
c***revision date  820801   (yymmdd)
c***category no.  d2b1b
c***keywords  double precision,factor,linear algebra,linpack,matrix,
c             positive definite
c***author  moler, c. b., (u. of new mexico)
c***purpose  factors a double precision symmetric positive definite
c            matrix.
c***description

c     dpofa factors a double precision symmetric positive definite
c     matrix.

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     on entry

c        a       double precision(lda, n)
c                the symmetric matrix to be factored.  only the
c                diagonal and upper triangle are used.

c        lda     integer
c                the leading dimension of the array  a .

c        n       integer
c                the order of the matrix  a .

c     on return

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        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     linpack.  this version dated 08/14/78 .
c     cleve moler, university of new mexico, argonne national lab.

c     subroutines and functions

c     blas ddot
c     fortran dsqrt
c***references  dongarra j.j., bunch j.r., moler c.b., stewart g.w.,
c                 *linpack users  guide*, siam, 1979.
c***routines called  ddot
c***end prologue  dpofa
      integer lda,n,info
      double precision a(lda,1)

      double precision ddot,t
      double precision s
      integer j,jm1,k
c     begin block with ...exits to 40


c***first executable statement  dpofa
         do 30 j = 1, n
            info = j
            s = 0.0d0
            jm1 = j - 1
            if (jm1 .lt. 1) go to 20
            do 10 k = 1, jm1
               t = a(k,j) - ddot(k-1,a(1,k),1,a(1,j),1)
               t = t/a(k,k)
               a(k,j) = t
               s = s + t*t
   10       continue
   20       continue
            s = a(j,j) - s
c     ......exit
            if (s .le. 0.0d0) go to 40
            a(j,j) = dsqrt(s)
   30    continue
         info = 0
   40 continue
      return
      end
      double precision function ddot(n,dx,incx,dy,incy)
c***begin prologue  ddot
c***date written   791001   (yymmdd)
c***revision date  820801   (yymmdd)
c***category no.  d1a4
c***keywords  blas,double precision,inner product,linear algebra,vector
c***author  lawson, c. l., (jpl)
c           hanson, r. j., (snla)
c           kincaid, d. r., (u. of texas)
c           krogh, f. t., (jpl)
c***purpose  d.p. inner product of d.p. vectors
c***description

c                b l a s  subprogram
c    description of parameters

c     --input--
c        n  number of elements in input vector(s)
c       dx  double precision vector with n elements
c     incx  storage spacing between elements of dx
c       dy  double precision vector with n elements
c     incy  storage spacing between elements of dy

c     --output--
c     ddot  double precision dot product (zero if n .le. 0)

c     returns the dot product of double precision dx and dy.
c     ddot = sum for i = 0 to n-1 of  dx(lx+i*incx) * dy(ly+i*incy)
c     where lx = 1 if incx .ge. 0, else lx = (-incx)*n, and ly is
c     defined in a similar way using incy.
c***references  lawson c.l., hanson r.j., kincaid d.r., krogh f.t.,
c                 *basic linear algebra subprograms for fortran usage*,
c                 algorithm no. 539, transactions on mathematical
c                 software, volume 5, number 3, september 1979, 308-323
c***routines called  (none)
c***end prologue  ddot

      double precision dx(1),dy(1)
c***first executable statement  ddot
      ddot = 0.d0
      if(n.le.0)return
      if(incx.eq.incy) if(incx-1) 5,20,60
    5 continue

c         code for unequal or nonpositive increments.

      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
         ddot = ddot + dx(ix)*dy(iy)
        ix = ix + incx
        iy = iy + incy
   10 continue
      return

c        code for both increments equal to 1.


c        clean-up loop so remaining vector length is a multiple of 5.

   20 m = mod(n,5)
      if( m .eq. 0 ) go to 40
      do 30 i = 1,m
         ddot = ddot + dx(i)*dy(i)
   30 continue
      if( n .lt. 5 ) return
   40 mp1 = m + 1
      do 50 i = mp1,n,5
         ddot = ddot + dx(i)*dy(i) + dx(i+1)*dy(i+1) +
     1   dx(i + 2)*dy(i + 2) + dx(i + 3)*dy(i + 3) + dx(i + 4)*dy(i + 4)
   50 continue
      return

c         code for positive equal increments .ne.1.

   60 continue
      ns = n*incx
          do 70 i=1,ns,incx
          ddot = ddot + dx(i)*dy(i)
   70     continue
      return
      end


      subroutine daxpy(n,da,dx,incx,dy,incy)
c***begin prologue  daxpy
c***date written   791001   (yymmdd)
c***revision date  820801   (yymmdd)
c***category no.  d1a7
c***keywords  blas,double precision,linear algebra,triad,vector
c***author  lawson, c. l., (jpl)
c           hanson, r. j., (snla)
c           kincaid, d. r., (u. of texas)
c           krogh, f. t., (jpl)
c***purpose  d.p computation y = a*x + y
c***description

c                b l a s  subprogram
c    description of parameters

c     --input--
c        n  number of elements in input vector(s)
c       da  double precision scalar multiplier
c       dx  double precision vector with n elements
c     incx  storage spacing between elements of dx
c       dy  double precision vector with n elements
c     incy  storage spacing between elements of dy

c     --output--
c       dy  double precision result (unchanged if n .le. 0)

c     overwrite double precision dy with double precision da*dx + dy.
c     for i = 0 to n-1, replace  dy(ly+i*incy) with da*dx(lx+i*incx) +
c       dy(ly+i*incy), where lx = 1 if incx .ge. 0, else lx = (-incx)*n
c       and ly is defined in a similar way using incy.
c***references  lawson c.l., hanson r.j., kincaid d.r., krogh f.t.,
c                 *basic linear algebra subprograms for fortran usage*,
c                 algorithm no. 539, transactions on mathematical
c                 software, volume 5, number 3, september 1979, 308-323
c***routines called  (none)
c***end prologue  daxpy

      double precision dx(1),dy(1),da
c***first executable statement  daxpy
      if(n.le.0.or.da.eq.0.d0) return
      if(incx.eq.incy) if(incx-1) 5,20,60
    5 continue

c        code for nonequal or nonpositive increments.

      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        code for both increments equal to 1


c        clean-up loop so remaining vector length is a multiple of 4.

   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

c        code for equal, positive, nonunit increments.

   60 continue
      ns = n*incx
          do 70 i=1,ns,incx
          dy(i) = da*dx(i) + dy(i)
   70     continue
      return
      end