Télécharger la plaquette Cast3M


fr en

Télécharger ntape6.eso

Retour à la liste

Numérotation des lignes : 

ntape6

  1. C NTAPE6    SOURCE    CB215821  26/08/24    21:17:35     12622          
  2.       SUBROUTINE NTAPE6(MCP,MCQ,IVMINU,IVMINL,IVMAXU,IVMAXL,IVLAMB,
  3.      * M,N,NVD,IVFP,IVFQ,MVDU,MVDL,IVB,IVD,IVN,II,KK,IVDR,IDVD,
  4.      * NDR,TERMIN,IVLL,IVUL,IPBASE)
  5.       IMPLICIT INTEGER(I-N)
  6.       IMPLICIT REAL*8(A-H,O-Z)
  7.       LOGICAL TERMIN
  8. -INC TMXMAT
  9. -INC SMLENTI
  10.  
  11. -INC PPARAM
  12. -INC CCOPTIO
  13. -INC SMLREEL
  14.       POINTEUR MLREE4.MLREEL,MLREE5.MLREEL,MLREE6.MLREEL,MLREE7.MLREEL
  15.       IVXL1=0
  16.       IVXU1=0
  17.       IVU1=0
  18.       IVN1=0
  19.       IVD1=0
  20.       IVGM1=0
  21.       IVGE1=0
  22.       IVLAM1=0
  23.       N11 = N + 1
  24.       II=-1
  25.       KK=-1
  26.       MLREEL=IVDR
  27.       JG=N11
  28.       SEGINI MLREE1,MLREE2
  29.       IVPZ=MLREE2
  30.       IVQZ=MLREE1
  31.       MXMAT=MCQ
  32.       CALL MATVE1(XMAT,PROG,M,N11,MLREE1.PROG,1)
  33.       MXMAT=MCP
  34.       CALL MATVE1(XMAT,PROG,M,N11,MLREE2.PROG,1)
  35.       ALPMAX=1.E25
  36.       MLREE1=IVLAMB
  37.       IPB=0
  38.       DO 1 I=1,M
  39.       IF(PROG(I).LT.-1.D-20)   THEN
  40.           IF(ALPMAX+(MLREE1.PROG(I)/PROG(I)).GT.0.D0) THEN
  41.             ALPMAX= -( MLREE1.PROG(I)/PROG(I))
  42.             IPB=I
  43.           ENDIF
  44.       ENDIF
  45.     1 CONTINUE
  46.       IF(IIMPI.EQ.1799)
  47.      *WRITE(IOIMP,FMT='('' ALPHAMAX  IPB = '',E12.5,2X,I3)')ALPMAX,IPB
  48.       NDIS=0
  49.       MLENTI=IDVD
  50.       DO 7 I=1,NVD
  51.          NDIS=NDIS+LECT(I)-1
  52.     7 CONTINUE
  53.       JG=(NDIS-NVD)+2*(N11-NVD)+2
  54.       SEGINI MLREE6,MLENT1,MLENT2
  55.       MLREE6.PROG(1)=0.D0
  56.       MLENT1.LECT(1)=-3
  57.       MLENT2.LECT(1)=0
  58. *
  59. *  CALCUL DES ALPHA CRITIQUES
  60. *  ON COMMENCE PAR LES VARIABLES CONTINUES
  61. *
  62.       IP=1
  63.       MLREEL=IVD
  64.       MLREE1=IVN
  65.       MLREE2=IVMINU
  66.       MLREE3=IVMINL
  67.       MLREE4=IVQZ
  68.       MLREE5=IVPZ
  69.       DO 2 I=NVD+1,N11
  70.          CN=PROG(I)*(MLREE2.PROG(I)**2)
  71.          CN=CN-(MLREE1.PROG(I)*(MLREE3.PROG(I)**2))
  72.          CD=(MLREE3.PROG(I) ** 2) * MLREE5.PROG(I)
  73.          CD=CD-((MLREE2.PROG(I)**2)*MLREE4.PROG(I))
  74.          IP=IP+1
  75.          IF(CD.EQ.0.D0) CD=1.D-35
  76.          MLREE6.PROG(IP)=CN/CD
  77.          MLENT1.LECT(IP)=I-NVD
  78.          MLENT2.LECT(IP)=-1
  79.     2 CONTINUE
  80.       MLREE2=IVMAXU
  81.       MLREE3=IVMAXL
  82.       DO 3 I=NVD+1,N11
  83.          CN=PROG(I)*(MLREE2.PROG(I)**2)
  84.          CN=CN-(MLREE1.PROG(I)*(MLREE3.PROG(I)**2))
  85.          CD=(MLREE3.PROG(I) ** 2) * MLREE5.PROG(I)
  86.          CD=CD-((MLREE2.PROG(I)**2)*MLREE4.PROG(I))
  87.          IP=IP+1
  88.          IF(CD.EQ.0.D0) CD=1.D-35
  89.          MLREE6.PROG(IP)=CN/CD
  90.          MLENT1.LECT(IP)=I-NVD
  91.          MLENT2.LECT(IP)=-1
  92.     3 CONTINUE
  93. *
  94. * POUR LES VARIABLES DISCRETES
  95. *
  96.       IF(NVD.NE.0) THEN
  97.        IF(IIMPI.EQ.1799)
  98.      * WRITE(IOIMP,FMT='('' IL Y A'',I4,'' VARIABLES DISCRETES'')')NVD
  99.         MXMAT=MVDU
  100.         MXMA1=MVDL
  101.         DO 61 I=1,NVD
  102.           DO 4 J = 3,LECT(I)
  103.             CN=PROG(I)*XMAT(I,J)*XMAT(I,J-1)
  104.             CN=CN-(MLREE1.PROG(I)*MXMA1.XMAT(I,J)*MXMA1.XMAT(I,J-1))
  105.             CD=MLREE4.PROG(I)* MXMA1.XMAT(I,J)*MXMA1.XMAT(I,J-1)
  106.             CD=CD-(MLREE5.PROG(I)*XMAT(I,J)*XMAT(I,J-1))
  107.             IP=IP+1
  108.             IF(CD.EQ.0.D0) CD=1.D-35
  109.             MLREE6.PROG(IP)= CN/CD
  110.             IF(ABS(MLREE6.PROG(IP)).LT.1D-10)MLREE6.PROG(IP)=0.D0
  111.             MLENT1.LECT(IP)=I
  112.             MLENT2.LECT(IP)=J
  113.     4     CONTINUE
  114.    61   CONTINUE
  115.       ENDIF
  116.       MAC2=MLENT1
  117.       MAC3=MLENT2
  118.       SEGSUP MLREE4,MLREE5
  119. *
  120. * ON ENTRE ALPHA MAX DANS LE TABLEAU DES ALPHA
  121. *
  122.       IP=IP+1
  123.       MLREE6.PROG(IP)=ALPMAX
  124. *
  125.       MLREEL=MLREE6
  126. *
  127. *  TRI DES ALPHA
  128. *
  129.       JG=PROG(/1)
  130.       SEGINI MLREE1
  131.       SEGINI MLENT1,MLENT2
  132.       MORDRE=MLENT1
  133.       DO 5 I=1,JG
  134.          MLENT1.LECT(I)=I
  135.    5  CONTINUE
  136.       CALL TRIFLO(PROG,MLREE1.PROG,MLENT1.LECT,MLENT2.LECT,JG)
  137.       IF(IIMPI.EQ.1799) WRITE(IOIMP,FMT='('' LISTE DES ALPHA APRES TRI
  138.      * '',/,(1X,5E12.5))')(PROG(I),I=1,IP)
  139.       SEGSUP MLREE1,MLENT2
  140. *
  141. *   ELIMINATION DES ALPHA NON ADMISSIBLES
  142. *
  143.       IDL=0
  144.       IFL=IP+2
  145.       DO 6 I=1,IP
  146.          IF(PROG(I).LT.1.D-20) IDL=I
  147.          IF(PROG(I).GT.ALPMAX)THEN
  148.              IF(IFL.EQ.IP+2) IFL=I-1
  149.              MLENT1.LECT(I)=-2
  150.          ENDIF
  151.    6  CONTINUE
  152.       IFL=MIN(IFL,IP)
  153.       IF(IIMPI.EQ.1799) WRITE(IOIMP,FMT='('' LISTE DES ALPHA APRES TRI
  154.      * ET ELIMINATION '',/,(1X,5E12.5))')(PROG(I),I=IDL,IFL)
  155. *
  156. *   RECHERCHE DU PAS ALPHA OPTIMUM
  157. *
  158.       MLREE1=IVLAMB
  159.       MLREE2=IVDR
  160.       SEGINI,MLREE3=MLREE1
  161.       IVLAM1=MLREE3
  162. *
  163. *  CALCUL DU PREMIER POINT
  164. *
  165.       I=IDL
  166.       INPINF=IDL
  167.       QQZ=PROG(I)
  168.       DO 8 J=1,MLREE1.PROG(/1)
  169.          MLREE3.PROG(J)=QQZ*MLREE2.PROG(J)+MLREE1.PROG(J)
  170.     8 CONTINUE
  171.       IF(IIMPI.EQ.1799) WRITE(IOIMP,FMT='('' VALEUR DU LAMBDA POUR ALPHA
  172.      * =0 '',/,(1X,5E12.5))')(MLREE3.PROG(I),I=1,M)
  173.       CALL NTAPE1(MCP,MCQ,IVFP,IVFQ,IVLAM1,NVD,M,N,MVDU,MVDL,IVMINU,
  174.      * IVMINL,IVMAXU,IVMAXL,IVU1,IVN1,IVD1,IVUL,IVLL,IVXU1,IVXL1)
  175.       CALL NTAPE2(MCP,MCQ,IVXU1,IVXL1,IVB,N,M,IVGE1,
  176.      *IVGM1,IVLAM1,IPBASE)
  177.       XLP1=0.D0
  178.       MLREE4=IVGE1
  179.       DO 9 J=1,MLREE1.PROG(/1)
  180.       XLP1=MLREE2.PROG(J)*MLREE4.PROG(J) +XLP1
  181.     9 CONTINUE
  182.       IF(IIMPI.EQ.1799) WRITE(IOIMP,FMT='(''  1ER ALPHA , PROD SCAL :''
  183.      *  ,1X,E12.5,'' , '',E12.5)')QQZ,XLP1
  184.          MLREE4=IVXU1
  185.          MLREE5=IVU1
  186.          MLREE6=IVN1
  187.          SEGSUP MLREE4,MLREE5,MLREE6
  188.          MLREE4=IVD1
  189.          MLREE5=IVGM1
  190.          MLREE6=IVGE1
  191.          MLREE7=IVXL1
  192.          SEGSUP MLREE4,MLREE5,MLREE6,MLREE7
  193.       IVXL1=0
  194.       IVXU1=0
  195.       IVU1=0
  196.       IVN1=0
  197.       IVD1=0
  198.       IVGM1=0
  199.       IVGE1=0
  200. *
  201. * CALCUL SUR LE DERNIER POINT
  202. *
  203.       INPSUP=IFL
  204.       I=IFL
  205.       QQZ=PROG(I)
  206.       DO 10 J=1,MLREE1.PROG(/1)
  207.          MLREE3.PROG(J)=QQZ*MLREE2.PROG(J)+MLREE1.PROG(J)
  208.    10 CONTINUE
  209.       CALL NTAPE1(MCP,MCQ,IVFP,IVFQ,IVLAM1,NVD,M,N,MVDU,MVDL,IVMINU,
  210.      * IVMINL,IVMAXU,IVMAXL,IVU1,IVN1,IVD1,IVUL,IVLL,IVXU1,IVXL1)
  211.       CALL NTAPE2(MCP,MCQ,IVXU1,IVXL1,IVB,N,M,IVGE1,
  212.      *IVGM1,IVLAM1,IPBASE)
  213.       XLP2=0.D0
  214.       MLREE4=IVGE1
  215.       DO 11 J=1,MLREE1.PROG(/1)
  216.          XLP2=MLREE2.PROG(J)*MLREE4.PROG(J) +XLP2
  217.    11 CONTINUE
  218.       IF(IIMPI.EQ.1799) WRITE(IOIMP,FMT='('' DERNIER ALPHA , PROD SCAL
  219.      * :'',1X,E12.5,'' , '',E12.5)')QQZ,XLP2
  220. *
  221. * RECHERCHE  PAR DICHOTOMIE SUR LES INDICES pour trouver l'interval
  222. * de alpha interessant
  223. *
  224.       IF(XLP1 * XLP2.GE.0.D0) THEN
  225.           ALPHA=ALPMAX
  226.           IF (ALPMAX.LT.1E25) THEN
  227.              KK = -3
  228.              II=IPB
  229.           ENDIF
  230.           GO TO 20
  231.       ENDIF
  232.    12 CONTINUE
  233.       IK=(INPINF+INPSUP)/2
  234.       IF(IK.EQ.INPINF) GO TO 15
  235.          MLREE4=IVXU1
  236.          MLREE5=IVU1
  237.          MLREE6=IVN1
  238.          SEGSUP MLREE4,MLREE5,MLREE6
  239.          MLREE4=IVD1
  240.          MLREE5=IVGM1
  241.          MLREE6=IVGE1
  242.          MLREE7=IVXL1
  243.          SEGSUP MLREE4,MLREE5,MLREE6,MLREE7
  244.       IVXL1=0
  245.       IVXU1=0
  246.       IVU1=0
  247.       IVN1=0
  248.       IVD1=0
  249.       IVGM1=0
  250.       IVGE1=0
  251.       I=IK
  252.       QQZ=PROG(I)
  253.       DO 13 J=1,MLREE1.PROG(/1)
  254.          MLREE3.PROG(J)=QQZ*MLREE2.PROG(J)+MLREE1.PROG(J)
  255.    13 CONTINUE
  256.       CALL NTAPE1(MCP,MCQ,IVFP,IVFQ,IVLAM1,NVD,M,N,MVDU,MVDL,IVMINU,
  257.      * IVMINL,IVMAXU,IVMAXL,IVU1,IVN1,IVD1,IVUL,IVLL,IVXU1,IVXL1)
  258.       CALL NTAPE2(MCP,MCQ,IVXU1,IVXL1,IVB,N,M,IVGE1,
  259.      *IVGM1,IVLAM1,IPBASE)
  260.       XLP3=0.D0
  261.       MLREE4=IVGE1
  262.       DO 14 J=1,MLREE1.PROG(/1)
  263.          XLP3=MLREE2.PROG(J)*MLREE4.PROG(J) +XLP3
  264.    14 CONTINUE
  265.       IF(XLP1*XLP3.LE.0.D0) THEN
  266.          XLP2=XLP3
  267.          INPSUP=I
  268.       ELSE
  269.          XLP1=XLP3
  270.          INPINF=I
  271.       ENDIF
  272.       GO TO 12
  273.    15 CONTINUE
  274. *
  275. * RECHERCHE PAR DICHOTOMIE pour determiner la valeur de alpha donnant le
  276. * min de l(lambda + alpha * grad)
  277. *
  278.       IPP=0
  279.       ALPE=PROG(INPINF)
  280.       ALGR=PROG(INPSUP)
  281.       ALPHA=(ALGR + ALPE)/2.D0
  282.    60 CONTINUE
  283.       IPP=IPP + 1
  284.       IF(IPP.GT.100) THEN
  285.         CALL ERREUR ( 603)
  286.         RETURN
  287.        ENDIF
  288.       ALPHA0=ALPHA
  289.          MLREE4=IVXU1
  290.          MLREE5=IVU1
  291.          MLREE6=IVN1
  292.          SEGSUP MLREE4,MLREE5,MLREE6
  293.          MLREE4=IVD1
  294.          MLREE5=IVGM1
  295.          MLREE6=IVGE1
  296.          MLREE7=IVXL1
  297.          SEGSUP MLREE4,MLREE5,MLREE6,MLREE7
  298.       IVXL1=0
  299.       IVXU1=0
  300.       IVU1=0
  301.       IVN1=0
  302.       IVD1=0
  303.       IVGM1=0
  304.       IVGE1=0
  305. *      IF(IIMPI.EQ.1799) WRITE(IOIMP,FMT='('' ALPHA '',E12.5)')ALPHA
  306.       QQZ=ALPHA
  307.       DO 17 J=1,MLREE1.PROG(/1)
  308.          MLREE3.PROG(J)=QQZ*MLREE2.PROG(J)+MLREE1.PROG(J)
  309.   17  CONTINUE
  310.       CALL NTAPE1(MCP,MCQ,IVFP,IVFQ,IVLAM1,NVD,M,N,MVDU,MVDL,IVMINU,
  311.      *    IVMINL,IVMAXU,IVMAXL,IVU1,IVN1,IVD1,IVUL,IVLL,IVXU1,IVXL1)
  312.       CALL NTAPE2(MCP,MCQ,IVXU1,IVXL1,IVB,N,M,IVGE1,
  313.      *IVGM1,IVLAM1,IPBASE)
  314.       XLP3=0.D0
  315.       MLREE4=IVGE1
  316.       DO 50 J=1,MLREE1.PROG(/1)
  317.          XLP3=MLREE2.PROG(J)*MLREE4.PROG(J) +XLP3
  318.    50 CONTINUE
  319.       IF(XLP1*XLP3.LE.0.D0) THEN
  320.          XLP2=XLP3
  321.          ALGR=ALPHA
  322.       ELSE
  323.          XLP1=XLP3
  324.          ALPE=ALPHA
  325.       ENDIF
  326.       ALPHA =(ALGR + ALPE)/2
  327.       EPSI=1.D-10
  328.       IF(ALPHA.LT.1.D-12.AND.NDR.GE.1)THEN
  329.          TERMIN=.TRUE.
  330.          GO TO 20
  331.       ENDIF
  332.       IF((ABS(ALPHA-ALPHA0)/ALPHA).GE.EPSI)   GO TO 60
  333.    20 CONTINUE
  334.          MLENT1=MAC2
  335.          MLENT2=MAC3
  336.          MLENT3=MORDRE
  337.        IF(IIMPI.EQ.1799)WRITE(IOIMP,FMT='('' ALPHA '',E12.5)')ALPHA
  338.       IF (KK.NE.-3)THEN
  339. *
  340. ****************** TEST SI ARRET SUR PLAN DE DISCONTINUITE ************
  341. *
  342.          ALPHA1=PROG(INPSUP)
  343.          IF((ABS(ALPHA1-ALPHA)/ALPHA).LE.EPSI)THEN
  344.               KK=MLENT2.LECT(MLENT3.LECT(INPSUP))
  345.               II=MLENT1.LECT(MLENT3.LECT(INPSUP))
  346.               ALPHA=ALPHA1
  347.          ENDIF
  348.          ALPHA2=PROG(INPINF)
  349.          IF(((ABS(ALPHA2-ALPHA)/ALPHA).LE.EPSI.AND.ALPHA2.NE.0.D0)
  350.      *     .OR.ALPHA.LT.1.D-12)THEN
  351.               KK=MLENT2.LECT(MLENT3.LECT(INPINF))
  352.               II=MLENT1.LECT(MLENT3.LECT(INPINF))
  353.               ALPHA=ALPHA2
  354.          ENDIF
  355.       ENDIF
  356.       DO 21 J=1,MLREE1.PROG(/1)
  357.          MLREE3.PROG(J)=ALPHA * MLREE2.PROG(J)+MLREE1.PROG(J)
  358.    21 CONTINUE
  359.       IF (KK.EQ.-3)THEN
  360.          MLREE3.PROG(II)=0.
  361.       ENDIF
  362.       ALPHA0=ALPHA
  363.          MLREE4=IVXU1
  364.          MLREE5=IVU1
  365.          MLREE6=IVN1
  366.          SEGSUP MLREE4,MLREE5,MLREE6
  367.          MLREE4=IVD1
  368.          MLREE5=IVGM1
  369.          MLREE6=IVGE1
  370.          MLREE7=IVXL1
  371.          SEGSUP MLREE4,MLREE5,MLREE6,MLREE7
  372.       IVLAMB=IVLAM1
  373.       SEGSUP MLENT1,MLENT2,MLREEL,MLENT3
  374.       RETURN
  375.       END
  376.  
  377.  
  378.  

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