c     ..................................................................
c
c        subroutine dmfsd
c
c        purpose
c           factor a given symmetric positive definite matrix
c
c        usage
c           call dmfsd(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           the product of returned diagonal terms is equal to the
c           square-root of the determinant of the given matrix.
c
c        subroutines and function subprograms required
c           none
c
c        method
c           solution is done using the square-root method of cholesky.
c           the given matrix is represented as product of two triangular
c           matrices, where the left hand factor is the transpose of
c           the returned right hand factor.
c
c     ..................................................................
c
      subroutine dmfsd (a,n,eps,ier)
c
      dimension a(1)
      double precision dpiv, dsum, a
c
c        test on wrong input parameter n
      if (n-1) 120,10,10
   10 ier=0
c
c        initialize diagonal-loop
      kpiv=0
      do 110 k=1,n
      kpiv=kpiv+k
      ind=kpiv
      lend=k-1
c
c        calculate tolerance
      tol=abs(eps*sngl(a(kpiv)))
c
c        start factorization-loop over k-th row
      do 110 i=k,n
      dsum=0.d0
      if (lend) 20,40,20
c
c        start inner loop
   20 do 30 l=1,lend
      lanf=kpiv-l
      lind=ind-l
   30 dsum=dsum+a(lanf)*a(lind)
c        end of inner loop
c
c        transform element a(ind)
   40 dsum=a(ind)-dsum
      if (i-k) 100,50,100
c
c        test for negative pivot element and for loss of significance
   50 if (sngl(dsum)-tol) 60,60,90
   60 if (dsum) 120,120,70
   70 if (ier) 80,80,90
   80 ier=k-1
c
c        compute pivot element
   90 dpiv=dsqrt(dsum)
      a(kpiv)=dpiv
      dpiv=1.d0/dpiv
      go to 110
c
c        calculate terms in row
  100 a(ind)=dsum*dpiv
  110 ind=ind+i
c        end of diagonal-loop
c
      return
  120 ier=-1
      return
      end