Télécharger la plaquette Cast3M


fr en

Télécharger detq3d.eso

Retour à la liste

Numérotation des lignes : 

detq3d

  1. C DETQ3D    SOURCE    PV090527  23/02/13    21:15:05     11592          
  2.       subroutine detq3d (a,a2,ngf,n,d,ifault)
  3.  
  4. c*********************************************************************72
  5. c
  6. cc detq computes the determinant of an orthogonal matrix.
  7. c
  8. c  Licensing:
  9. c
  10. c    This code is distributed under the GNU LGPL license.
  11. c
  12. c  Modified:
  13. c
  14. c    16 January 2020
  15. c
  16. c  Author:
  17. c
  18. c    Original FORTRAN77 version by J C Gower.
  19. c    This version by John Burkardt
  20. c
  21. c  Reference:
  22. c
  23. c    J C Gower,
  24. c    Algorithm AS 82:
  25. c    The determinant of an orthogonal matrix,
  26. c    Applied Statistics,
  27. c    Volume 24, Number 1, 1975, page 150-153.
  28. c
  29. c  Input:
  30. c
  31. c    real A(N,N), the orthogonal matrix stored by rows or columns.
  32. c
  33. c    integer N, the order of the matrix.
  34. c
  35. c  Output:
  36. c
  37. c    real D, the determinant of A.
  38. c
  39. c    integer IFAULT, 
  40. c    0, no error occurred.
  41. c    1, an error was detected.
  42.  
  43.  
  44.       integer n
  45.       integer ngf
  46.       real*8 a(ngf,ngf+1)
  47.       real*8 a2(ngf*ngf)
  48.       real*8 d
  49.       integer i
  50.       integer ifault
  51.       integer j
  52.       integer k
  53.       integer p
  54.       integer q
  55.       integer r
  56.       integer s
  57.       real*8 tol
  58.       real*8 x
  59.       real*8 y
  60.  
  61.       ifault = 0
  62.       tol = 1.0d-8
  63.  
  64.       if ( n .le. 0 ) then
  65.         write ( *, '(a)' ) ' '
  66.         write ( *, '(a)' ) 'detq - Fatal error!'
  67.         write ( *, '(a)' ) '  n <= 0'
  68.         ifault = 1
  69.         return
  70.       end if
  71.  
  72.       k = 0
  73.       do i = 1, n
  74.         do j = 1, n
  75.           k = k + 1
  76.           a2(k) = a(i,j)
  77.         end do
  78.       end do
  79.  
  80.       d = 1.0D+00
  81.       r = 1
  82.  
  83.       do k = 2, n + 1
  84.  
  85.         q = r
  86.         x = a2(r)
  87.         y = sign ( 1.0D+00, x )
  88.         d = d * y
  89.         y = - 1.0D+00 / ( x + y )
  90.         x = abs ( x ) - 1.0D+00
  91.  
  92.         if ( tol .lt. abs ( x ) ) then
  93.  
  94.           if ( 0.0D+00 .lt. x ) then
  95.             write ( *, '(a)' ) ' '
  96.             write ( *, '(a)' ) 'detq - Fatal error!'
  97.             write ( *, '(a)' ) '  x < 0.0'
  98.             write ( *, '(a,g14.6)' ) '  x = ', x
  99.             ifault = 1
  100.             return
  101.           end if
  102.  
  103.           if ( k .eq. (n + 1) ) then
  104.             write ( *, '(a)' ) ' '
  105.             write ( *, '(a)' ) 'detq - Fatal error!'
  106.             write ( *, '(a)' ) '  k == n + 1'
  107.             ifault = 1
  108.             return
  109.           end if
  110.  
  111.           do i = k, n
  112.             q = q + n
  113.             x = a2(q) * y
  114.             p = r
  115.             s = q
  116.             do j = k, n
  117.               p = p + 1
  118.               s = s + 1
  119.               a2(s) = a2(s) + x * a2(p)
  120.             end do
  121.  
  122.           end do
  123.  
  124.         end if
  125.  
  126.         r = r + n + 1
  127.  
  128.       end do
  129.  
  130.       return
  131.       end
  132.  
  133.  
  134.  
  135.  

© Cast3M 2003 - Tous droits réservés.
Mentions légales