subroutine aptdet4 (s, tol, det)

ccbeg.
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c                             SUBROUTINE APTDET4
c
c     call aptdet4 (s, tol, det)
c
c     Version:  aptdet4  Updated    2003 November 21 14:30.
c               aptdet4  Originated 2003 November 20 17:00.
c
c     Author:   Arthur L. Edwards, LLNL, L-298, Telephone (925) 422-4123.
c
c
c     Purpose:  To find the determinant det of the 4x4 matrix s.
c               If the value of det is less than the estimated error in its
c               calculation, based on tol, det will be truncated to zero.
c
c               If det is zero, one or more of the following are true:
c               one or more rows are all zero;
c               one or more columns are all zero;
c               two rows are equal or differ only by a common multiple;
c               two columns are equal or differ only by a common multiple;
c               a row is a linear combination of two or more other rows;
c               a column is a linear combination of two or more other column;
c               the four row vectors are all coplanar in four-space;
c               the four column vectors are all coplanar in four-space.
c
c     Input:    s.
c
c     Output:   det.
c
c     Glossary:
c
c     det       Output   The determinant of the 4x4 matrix s.
c
c     s         Input    The 4x4 square matrix s.  Double-subscripted, size 4x4.
c
c     tol       Input    Numerical tolerance limit, if > 0.
c                          For 64-bit floating point arithmetic, recommend
c                          1.e-5 to 1.e-11.
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
ccend.

c.... Dimensioned arguments.

      real s(4,4)                     ! The 4x4 matrix s.

c.... Local variables.

      real deterr                     ! The estimated error in det.

cbugc***DEBUG begins.
cbug      integer i                       ! Row index in s.
cbug      integer j                       ! Column index in s.
cbug 9901 format (/ 'aptdet4 finding determinant of 4x4 matrix.  tol=',
cbug     &  1pe22.14 / (1p4e20.12))
cbug      write ( 3, 9901) tol, ((s(i,j), j = 1, 4), i = 1, 4)
cbugc***DEBUG ends.

c.... Find the determinant det of the matrix s.

      det = (s(2,2) * (s(3,3) * s(4,4) - s(3,4) * s(4,3))  +
     &       s(2,3) * (s(3,4) * s(4,2) - s(3,2) * s(4,4))  +
     &       s(2,4) * (s(3,2) * s(4,3) - s(3,3) * s(4,2))) * s(1,1) -
     &      (s(2,1) * (s(3,3) * s(4,4) - s(3,4) * s(4,3))  +
     &       s(2,3) * (s(3,4) * s(4,1) - s(3,1) * s(4,4))  +
     &       s(2,4) * (s(3,1) * s(4,3) - s(3,3) * s(4,1))) * s(1,2) +
     &      (s(2,1) * (s(3,2) * s(4,4) - s(3,4) * s(4,2))  +
     &       s(2,2) * (s(3,4) * s(4,1) - s(3,1) * s(4,4))  +
     &       s(2,4) * (s(3,1) * s(4,2) - s(3,2) * s(4,1))) * s(1,3) -
     &      (s(2,1) * (s(3,2) * s(4,3) - s(3,3) * s(4,2))  +
     &       s(2,2) * (s(3,3) * s(4,1) - s(3,1) * s(4,3))  +
     &       s(2,3) * (s(3,1) * s(4,2) - s(3,2) * s(4,1))) * s(1,4)

c.... See if the result should be tested for truncation to zero.

      deterr = 0.0

      if (tol .gt. 0.0) then          ! Truncate small result to zero.

      deterr = 4.0 * tol * (
     &         (abs (s(2,2)) * (abs (s(3,3) * s(4,4))   +
     &                          abs (s(3,4) * s(4,3)))  +
     &          abs (s(2,3)) * (abs (s(3,4) * s(4,2))   +
     &                          abs (s(3,2) * s(4,4)))  +
     &          abs (s(2,4)) * (abs (s(3,2) * s(4,3))   +
     &                          abs (s(3,3) * s(4,2)))) * abs (s(1,1)) +
     &         (abs (s(2,1)) * (abs (s(3,3) * s(4,4))   +
     &                          abs (s(3,4) * s(4,3)))  +
     &          abs (s(2,3)) * (abs (s(3,4) * s(4,1))   +
     &                          abs (s(3,1) * s(4,4)))  +
     &          abs (s(2,4)) * (abs (s(3,1) * s(4,3))   +
     &                          abs (s(3,3) * s(4,1)))) * abs (s(1,2)) +
     &         (abs (s(2,1)) * (abs (s(3,2) * s(4,4))   +
     &                          abs (s(3,4) * s(4,2)))  +
     &          abs (s(2,2)) * (abs (s(3,4) * s(4,1))   +
     &                          abs (s(3,1) * s(4,4)))  +
     &          abs (s(2,4)) * (abs (s(3,1) * s(4,2))   +
     &                          abs (s(3,2) * s(4,1)))) * abs (s(1,3)) +
     &         (abs (s(2,1)) * (abs (s(3,2) * s(4,3))   +
     &                          abs (s(3,3) * s(4,2)))  +
     &          abs (s(2,2)) * (abs (s(3,3) * s(4,1))   +
     &                          abs (s(3,1) * s(4,3)))  +
     &          abs (s(2,3)) * (abs (s(3,1) * s(4,2))   +
     &                          abs (s(3,2) * s(4,1)))) * abs (s(1,4)))

        if (abs (det) .le. deterr) det = 0.0

      endif                           ! Tested tol.

cbugc***DEBUG begins.
cbug 9902 format (/ 'aptdet4 results:' /
cbug     &  ' det =',1pe22.14,'  deterr =',1pe22.14)
cbug      write ( 3, 9902) det, deterr
cbugc***DEBUG ends.

      return

c.... End of subroutine aptdet4.      (+1 line.)
      end

UCRL-WEB-209832