Télécharger la plaquette Cast3M


fr en

Télécharger ottob3.eso

Retour à la liste

Numérotation des lignes : 

ottob3

  1. C OTTOB3    SOURCE    CHAT      05/01/13    02:07:00     5004
  2.       SUBROUTINE OTTOB3(A,B,LADIM,XX,KERRE)
  3. C======================================================================
  4. C
  5. C     ENTREES
  6. C     -------
  7. C     A(3,3) = MATRICE SYMETRIQUE
  8. C     B(3,3) = MATRICE SYMETRIQUE
  9. C     LADIM   = 2 OU 3  SI 2 ON NE S OCCUPE QUE DE A(2,2)
  10. C                       SI 3                    DE A(3,3)
  11. C     SORTIES
  12. C     -------
  13. C
  14. C       ON SORT XX(1) < XX(2) < XX(3)
  15. C
  16. C===============================================================
  17.       IMPLICIT INTEGER(I-N)
  18.       IMPLICIT REAL*8 (A-H,O-Z)
  19.  
  20. -INC PPARAM
  21. -INC CCOPTIO
  22.       DIMENSION A(3,3),B(3,3),XX(3)
  23. *
  24.       KERRE=0
  25. *
  26.       IF(LADIM.EQ.3) GOTO 700
  27. *
  28. *     CAS DIMENSION 2
  29. *
  30. *     A0 + A1*X + A2* X2
  31. *
  32.       A0=A(1,1)*A(2,2)-A(1,2)*A(1,2)
  33.       A1=A(1,1)*B(2,2)+A(2,2)*B(1,1)-2.D0*A(1,2)*B(1,2)
  34.       A2=B(1,1)*B(2,2)-B(1,2)*B(1,2)
  35.  
  36. *
  37.       IF(IIMPI.EQ.42) THEN
  38.          WRITE(IOIMP,70615) A0,A1,A2
  39. 70615    FORMAT(2X,'OTTOB3 -  A0 A1 A2 ',3(1X,1PE12.5)/)
  40.       ENDIF
  41. *
  42.       IF(A2.EQ.0.D0) THEN
  43.          IF(A1.EQ.0.D0) THEN
  44.             WRITE(IOIMP,76000)
  45. 76000       FORMAT(2X,'OTTOB3 - A2 ET A1 SONT NULS ' )
  46.             KERRE=70
  47.             RETURN
  48.          ENDIF
  49.          XX1=-A0/A1
  50.          XX2=1.D4
  51.       ELSE IF(A1.EQ.0.D0) THEN
  52.          DIS=-A0/A2
  53.          IF(DIS.LT.0.D0) THEN
  54.             WRITE(IOIMP,76001) DIS
  55. 76001       FORMAT(2X,'OTTOB3 - RAPPORT DIS NEGATIF = ',1PE12.5/)
  56.             KERRE=70
  57.             RETURN
  58.          ENDIF
  59.          XX1=-SQRT(DIS)
  60.          XX2= SQRT(DIS)
  61.       ELSE
  62.          AUX=4.D0*A2*A0/(A1*A1)
  63.          IF(ABS(AUX).LT.1.D-8) THEN
  64.             XX1=-A0/A1 - A0*A0*A2/(A1**3)
  65.             XX2=-A1/A2 -XX1
  66.          ELSE
  67.             DIS = A1*A1-4.D0*A2*A0
  68. *
  69.             IF(DIS.LT.0.D0) THEN
  70.                AUX=1.D-12*A1*A1
  71.                IF(ABS(DIS).LT.AUX ) THEN
  72.                   DIS= MAX(DIS,0.D0)
  73.                ELSE
  74.                  WRITE(IOIMP,76002) DIS , AUX
  75. 76002            FORMAT(2X,'OTTOB3 - DISCRIMINANT NEGATIF = ',
  76.      &                  1PE12.5,2X,' AUX=',1PE12.5/)
  77.                  KERRE=70
  78.                  RETURN
  79.               ENDIF
  80.             ENDIF
  81.             DIS=SQRT(DIS)
  82.             XX1=(-A1-DIS)/(2.D0*A2)
  83.             XX2=(-A1+DIS)/(2.D0*A2)
  84.          ENDIF
  85.       ENDIF
  86.       XX(1)=MIN(XX1,XX2)
  87.       XX(2)=MAX(XX1,XX2)
  88.       RETURN
  89. *
  90. *     CAS DIMENSION 3
  91. *
  92. *     A0 + A1*X + A2* X2 + A3* X3
  93. *
  94.  700  CONTINUE
  95. *
  96.       A0=A(1,1)*(A(2,2)*A(3,3)-A(2,3)*A(2,3))
  97.      &  -A(1,2)*(A(3,3)*A(1,2)-A(1,3)*A(2,3))
  98.      &  +A(1,3)*(A(1,2)*A(2,3)-A(2,2)*A(1,3))
  99. *
  100.       A1=A(1,1)*(A(3,3)*B(2,2)+A(2,2)*B(3,3)-2.D0*A(2,3)*B(2,3))
  101.      &  -A(1,2)*(A(3,3)*B(1,2)+A(1,2)*B(3,3)
  102.      &          -A(1,3)*B(2,3)-A(2,3)*B(1,3))
  103.      &  +A(1,3)*(A(2,3)*B(1,2)+A(1,2)*B(2,3)
  104.      &          -A(1,3)*B(2,2)-A(2,2)*B(1,3))
  105.      &  +B(1,1)*(A(2,2)*A(3,3)-A(2,3)*A(2,3))
  106.      &  -B(1,2)*(A(3,3)*A(1,2)-A(2,3)*A(1,3))
  107.      &  +B(1,3)*(A(1,2)*A(2,3)-A(2,2)*A(1,3))
  108. *
  109.       A2=B(1,1)*(A(3,3)*B(2,2)+A(2,2)*B(3,3)-2.D0*A(2,3)*B(2,3))
  110.      &  -B(1,2)*(A(3,3)*B(1,2)+A(1,2)*B(3,3)
  111.      &          -A(1,3)*B(2,3)-A(2,3)*B(1,3))
  112.      &  +B(1,3)*(A(2,3)*B(1,2)+A(1,2)*B(2,3)
  113.      &          -A(1,3)*B(2,2)-A(2,2)*B(1,3))
  114.      &  +A(1,1)*(B(2,2)*B(3,3)-B(2,3)*B(2,3))
  115.      &  -A(1,2)*(B(3,3)*B(1,2)-B(2,3)*B(1,3))
  116.      &  +A(1,3)*(B(1,2)*B(2,3)-B(2,2)*B(1,3))
  117. *
  118.       A3=B(1,1)*(B(2,2)*B(3,3)-B(2,3)*B(2,3))
  119.      &  -B(1,2)*(B(3,3)*B(1,2)-B(1,3)*B(2,3))
  120.      &  +B(1,3)*(B(1,2)*B(2,3)-B(2,2)*B(1,3))
  121. *
  122.          IF(IIMPI.EQ.42) THEN
  123.             WRITE(IOIMP,70616) A0,A1,A2,A3
  124. 70616       FORMAT(2X,'OTTOB3 -  A0 A1 A2 A3 ',4(1X,1PE12.5)/)
  125.          ENDIF
  126. ******
  127.       IF(ABS(A3).LT.1.D-6) THEN
  128.          A3=0.D0
  129.       ENDIF
  130. ********
  131.  
  132. *
  133.       IF(A3.EQ.0.D0) THEN
  134. *
  135.          IF(A2.EQ.0.D0) THEN
  136.             IF(A1.EQ.0.D0) THEN
  137.                WRITE(IOIMP,76003)
  138. 76003          FORMAT(2X,'OTTOB3 -  A3, A2 ET A1 SONT NULS ' )
  139.                KERRE=70
  140.                RETURN
  141.             ENDIF
  142.             XX(1)=-A0/A1
  143.             XX(2)=1.D4
  144.             XX(3)=1.D4
  145.          ELSE
  146.             DIS = A1*A1-4.D0*A2*A0
  147. ***********
  148.             DIS = MAX ( DIS, 0.D0)
  149. ***********
  150.  
  151.  
  152.             IF(DIS.LT.0.D0) THEN
  153.                WRITE(IOIMP,76001) DIS
  154.                KERRE=70
  155.                RETURN
  156.             ENDIF
  157.             DIS=SQRT(DIS)
  158.             XX(1)=(-A1-DIS)/(2.D0*A2)
  159.             XX(2)=(-A1+DIS)/(2.D0*A2)
  160.             XX(3)=1.D4
  161.          ENDIF
  162.          RETURN
  163. *
  164.       ELSE
  165. *
  166. *        MISE SOUS FORME POUR DEGRE3
  167. *
  168.          AS0=A0/A3
  169.          AS1=A1/A3
  170.          AS2=A2/A3
  171. *
  172.          CALL DEGRE3(AS0,AS1,AS2,XR1,XI1,XR2,XI2,XR3,XI3)
  173. *
  174.          IF(IIMPI.EQ.42) THEN
  175.             WRITE(IOIMP,76005) XR1,XI1,XR2,XI2,XR3,XI3
  176. 76005       FORMAT(2X,'OTTOB3 -  XR1 = ',1PE12.5,2X,'XI1=',1PE12.5/
  177.      &             2X,'          XR2 = ',1PE12.5,2X,'XI2=',1PE12.5/
  178.      &             2X,'          XR3 = ',1PE12.5,2X,'XI3=',1PE12.5/)
  179.          ENDIF
  180. *
  181. *        xi1 est toujours nul
  182. *
  183.          IF(XI2.NE.0.D0) THEN
  184.             XX(1)=XR1
  185.             XX(2)=1.D4
  186.             XX(3)=1.D4
  187.          ELSE
  188.             IF (XR1.LE.XR2.AND.XR1.LE.XR3)THEN
  189.                 XX(1)=XR1
  190.                 XX(2)=XR2
  191.                 XX(3)=XR3
  192.                 IF (XR3.LE.XR2) THEN
  193.                    XX(2)=XR3
  194.                    XX(3)=XR2
  195.                 ENDIF
  196. *
  197.             ELSE IF (XR2.LE.XR3.AND.XR2.LE.XR1)THEN
  198.                 XX(1)=XR2
  199.                 XX(2)=XR3
  200.                 XX(3)=XR1
  201.                 IF (XR1.LE.XR3) THEN
  202.                    XX(2)=XR1
  203.                    XX(3)=XR3
  204.                 ENDIF
  205. *
  206.             ELSE IF (XR3.LE.XR1.AND.XR3.LE.XR2)THEN
  207.                 XX(1)=XR3
  208.                 XX(2)=XR1
  209.                 XX(3)=XR2
  210.                 IF (XR2.LE.XR1) THEN
  211.                    XX(2)=XR2
  212.                    XX(3)=XR1
  213.                 ENDIF
  214.             ENDIF
  215.          ENDIF
  216.       ENDIF
  217.       RETURN
  218.       END
  219.  
  220.  
  221.  
  222.  

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