c     ..................................................................
c
c        subroutine dsinv
c
c        purpose
c           invert a given symmetric positive definite matrix
c
c        usage
c           call dsinv(a,n,eps,ier)
c
c        description of parameters
c           a      - double precision upper triangular part of given
c                    symmetric positive definite n by n coefficient
c                    matrix.
c                    on return a contains the resultant upper
c                    triangular matrix in double precision.
c           n      - the number of rows (columns) in given matrix.
c           eps    - single precision input constant which is used
c                    as relative tolerance for test on loss of
c                    significance.
c           ier    - resulting error parameter coded as follows
c                    ier=0  - no error
c                    ier=-1 - no result because of wrong input parame-
c                             ter n or because some radicand is non-
c                             positive (matrix a is not positive
c                             definite, possibly due to loss of signi-
c                             ficance)
c                    ier=k  - warning which indicates loss of signifi-
c                             cance. the radicand formed at factoriza-
c                             tion step k+1 was still positive but no
c                             longer greater than abs(eps*a(k+1,k+1)).
c
c        remarks
c           the upper triangular part of given matrix is assumed to be
c           stored columnwise in n*(n+1)/2 successive storage locations.
c           in the same storage locations the resulting upper triangu-
c           lar matrix is stored columnwise too.
c           the procedure gives results if n is greater than 0 and all
c           calculated radicands are positive.
c
c        subroutines and function subprograms required
c           dmfsd
c
c        method
c           solution is done using factorization by subroutine dmfsd.
c
c     ..................................................................
c
      subroutine dsinv (a,n,eps,ier)
c
      dimension a(1)
      double precision a, din, work
c
c        factorize given matrix by means of subroutine dmfsd
c        a = transpose(t) * t
      call dmfsd (a,n,eps,ier)
      if (ier) 90,10,10
c
c        invert upper triangular matrix t
c        prepare inversion-loop
   10 ipiv=n*(n+1)/2
      ind=ipiv
c
c        initialize inversion-loop
      do 60 i=1,n
      din=1.d0/a(ipiv)
      a(ipiv)=din
      min=n
      kend=i-1
      lanf=n-kend
      if (kend) 50,50,20
   20 j=ind
c
c        initialize row-loop
      do 40 k=1,kend
      work=0.d0
      min=min-1
      lhor=ipiv
      lver=j
c
c        start inner loop
      do 30 l=lanf,min
      lver=lver+1
      lhor=lhor+l
   30 work=work+a(lver)*a(lhor)
c        end of inner loop
c
      a(j)=-work*din
   40 j=j-min
c        end of row-loop
c
   50 ipiv=ipiv-min
   60 ind=ind-1
c        end of inversion-loop
c
c        calculate inverse(a) by means of inverse(t)
c        inverse(a) = inverse(t) * transpose(inverse(t))
c        initialize multiplication-loop
      do 80 i=1,n
      ipiv=ipiv+i
      j=ipiv
c
c        initialize row-loop
      do 80 k=i,n
      work=0.d0
      lhor=j
c
c        start inner loop
      do 70 l=k,n
      lver=lhor+k-i
      work=work+a(lhor)*a(lver)
   70 lhor=lhor+l
c        end of inner loop
c
      a(j)=work
   80 j=j+k
c        end of row- and multiplication-loop
c
   90 return
      end