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