C NTAPE6 SOURCE CHAT 05/01/13 02:02:22 5004 SUBROUTINE NTAPE6(MCP,MCQ,IVMINU,IVMINL,IVMAXU,IVMAXL,IVLAMB, * M,N,NVD,IVFP,IVFQ,MVDU,MVDL,IVB,IVD,IVN,II,KK,IVDR,IDVD, * NDR,TERMIN,IVLL,IVUL,IPBASE) IMPLICIT INTEGER(I-N) IMPLICIT REAL*8(A-H,O-Z) LOGICAL TERMIN -INC TMXMAT -INC SMLENTI -INC CCOPTIO -INC SMLREEL POINTEUR MLREE4.MLREEL,MLREE5.MLREEL,MLREE6.MLREEL,MLREE7.MLREEL IVXL1=0 IVXU1=0 IVU1=0 IVN1=0 IVD1=0 IVGM1=0 IVGE1=0 IVLAM1=0 N11 = N + 1 II=-1 KK=-1 MLREEL=IVDR JG=N11 SEGINI MLREE1,MLREE2 IVPZ=MLREE2 IVQZ=MLREE1 MXMAT=MCQ CALL MATVE1(XMAT,PROG,M,N11,MLREE1.PROG,1) MXMAT=MCP CALL MATVE1(XMAT,PROG,M,N11,MLREE2.PROG,1) ALPMAX=1.E25 MLREE1=IVLAMB IPB=0 DO 1 I=1,M IF(PROG(I).LT.-1.D-20) THEN IF(ALPMAX+(MLREE1.PROG(I)/PROG(I)).GT.0.D0) THEN ALPMAX= -( MLREE1.PROG(I)/PROG(I)) IPB=I ENDIF ENDIF 1 CONTINUE IF(IIMPI.EQ.1799) *WRITE(IOIMP,FMT='('' ALPHAMAX IPB = '',E12.5,2X,I3)')ALPMAX,IPB NDIS=0 MLENTI=IDVD DO 7 I=1,NVD NDIS=NDIS+LECT(I)-1 7 CONTINUE JG=(NDIS-NVD)+2*(N11-NVD)+2 SEGINI MLREE6,MLENT1,MLENT2 MLREE6.PROG(1)=0.D0 MLENT1.LECT(1)=-3 MLENT2.LECT(1)=0 * * CALCUL DES ALPHA CRITIQUES * ON COMMENCE PAR LES VARIABLES CONTINUES * IP=1 MLREEL=IVD MLREE1=IVN MLREE2=IVMINU MLREE3=IVMINL MLREE4=IVQZ MLREE5=IVPZ DO 2 I=NVD+1,N11 CN=PROG(I)*(MLREE2.PROG(I)**2) CN=CN-(MLREE1.PROG(I)*(MLREE3.PROG(I)**2)) CD=(MLREE3.PROG(I) ** 2) * MLREE5.PROG(I) CD=CD-((MLREE2.PROG(I)**2)*MLREE4.PROG(I)) IP=IP+1 IF(CD.EQ.0.D0) CD=1.D-35 MLREE6.PROG(IP)=CN/CD MLENT1.LECT(IP)=I-NVD MLENT2.LECT(IP)=-1 2 CONTINUE MLREE2=IVMAXU MLREE3=IVMAXL DO 3 I=NVD+1,N11 CN=PROG(I)*(MLREE2.PROG(I)**2) CN=CN-(MLREE1.PROG(I)*(MLREE3.PROG(I)**2)) CD=(MLREE3.PROG(I) ** 2) * MLREE5.PROG(I) CD=CD-((MLREE2.PROG(I)**2)*MLREE4.PROG(I)) IP=IP+1 IF(CD.EQ.0.D0) CD=1.D-35 MLREE6.PROG(IP)=CN/CD MLENT1.LECT(IP)=I-NVD MLENT2.LECT(IP)=-1 3 CONTINUE * * POUR LES VARIABLES DISCRETES * IF(NVD.NE.0) THEN IF(IIMPI.EQ.1799) * WRITE(IOIMP,FMT='('' IL Y A'',I4,'' VARIABLES DISCRETES'')')NVD MXMAT=MVDU MXMA1=MVDL DO 4 I=1,NVD DO 4 J = 3,LECT(I) CN=PROG(I)*XMAT(I,J)*XMAT(I,J-1) CN=CN-(MLREE1.PROG(I)*MXMA1.XMAT(I,J)*MXMA1.XMAT(I,J-1)) CD=MLREE4.PROG(I)* MXMA1.XMAT(I,J)*MXMA1.XMAT(I,J-1) CD=CD-(MLREE5.PROG(I)*XMAT(I,J)*XMAT(I,J-1)) IP=IP+1 IF(CD.EQ.0.D0) CD=1.D-35 MLREE6.PROG(IP)= CN/CD IF(ABS(MLREE6.PROG(IP)).LT.1D-10)MLREE6.PROG(IP)=0.D0 MLENT1.LECT(IP)=I MLENT2.LECT(IP)=J 4 CONTINUE ENDIF MAC2=MLENT1 MAC3=MLENT2 SEGSUP MLREE4,MLREE5 * * ON ENTRE ALPHA MAX DANS LE TABLEAU DES ALPHA * IP=IP+1 MLREE6.PROG(IP)=ALPMAX * MLREEL=MLREE6 * * TRI DES ALPHA * JG=PROG(/1) SEGINI MLREE1 SEGINI MLENT1,MLENT2 MORDRE=MLENT1 DO 5 I=1,JG MLENT1.LECT(I)=I 5 CONTINUE CALL TRIFLO(PROG,MLREE1.PROG,MLENT1.LECT,MLENT2.LECT,JG) IF(IIMPI.EQ.1799) WRITE(IOIMP,FMT='('' LISTE DES ALPHA APRES TRI * '',/,(1X,5E12.5))')(PROG(I),I=1,IP) SEGSUP MLREE1,MLENT2 * * ELIMINATION DES ALPHA NON ADMISSIBLES * IDL=0 IFL=IP+2 DO 6 I=1,IP IF(PROG(I).LT.1.D-20) IDL=I IF(PROG(I).GT.ALPMAX)THEN IF(IFL.EQ.IP+2) IFL=I-1 MLENT1.LECT(I)=-2 ENDIF 6 CONTINUE IFL=MIN(IFL,IP) IF(IIMPI.EQ.1799) WRITE(IOIMP,FMT='('' LISTE DES ALPHA APRES TRI * ET ELIMINATION '',/,(1X,5E12.5))')(PROG(I),I=IDL,IFL) * * RECHERCHE DU PAS ALPHA OPTIMUM * MLREE1=IVLAMB MLREE2=IVDR SEGINI,MLREE3=MLREE1 IVLAM1=MLREE3 * * CALCUL DU PREMIER POINT * I=IDL INPINF=IDL QQZ=PROG(I) DO 8 J=1,MLREE1.PROG(/1) MLREE3.PROG(J)=QQZ*MLREE2.PROG(J)+MLREE1.PROG(J) 8 CONTINUE IF(IIMPI.EQ.1799) WRITE(IOIMP,FMT='('' VALEUR DU LAMBDA POUR ALPHA * =0 '',/,(1X,5E12.5))')(MLREE3.PROG(I),I=1,M) CALL NTAPE1(MCP,MCQ,IVFP,IVFQ,IVLAM1,NVD,M,N,MVDU,MVDL,IVMINU, * IVMINL,IVMAXU,IVMAXL,IVU1,IVN1,IVD1,IVUL,IVLL,IVXU1,IVXL1) CALL NTAPE2(MCP,MCQ,IVXU1,IVXL1,IVB,N,M,IVGE1, *IVGM1,IVLAM1,IPBASE) XLP1=0.D0 MLREE4=IVGE1 DO 9 J=1,MLREE1.PROG(/1) XLP1=MLREE2.PROG(J)*MLREE4.PROG(J) +XLP1 9 CONTINUE IF(IIMPI.EQ.1799) WRITE(IOIMP,FMT='('' 1ER ALPHA , PROD SCAL :'' * ,1X,E12.5,'' , '',E12.5)')QQZ,XLP1 MLREE4=IVXU1 MLREE5=IVU1 MLREE6=IVN1 SEGSUP MLREE4,MLREE5,MLREE6 MLREE4=IVD1 MLREE5=IVGM1 MLREE6=IVGE1 MLREE7=IVXL1 SEGSUP MLREE4,MLREE5,MLREE6,MLREE7 IVXL1=0 IVXU1=0 IVU1=0 IVN1=0 IVD1=0 IVGM1=0 IVGE1=0 * * CALCUL SUR LE DERNIER POINT * INPSUP=IFL I=IFL QQZ=PROG(I) DO 10 J=1,MLREE1.PROG(/1) MLREE3.PROG(J)=QQZ*MLREE2.PROG(J)+MLREE1.PROG(J) 10 CONTINUE CALL NTAPE1(MCP,MCQ,IVFP,IVFQ,IVLAM1,NVD,M,N,MVDU,MVDL,IVMINU, * IVMINL,IVMAXU,IVMAXL,IVU1,IVN1,IVD1,IVUL,IVLL,IVXU1,IVXL1) CALL NTAPE2(MCP,MCQ,IVXU1,IVXL1,IVB,N,M,IVGE1, *IVGM1,IVLAM1,IPBASE) XLP2=0.D0 MLREE4=IVGE1 DO 11 J=1,MLREE1.PROG(/1) XLP2=MLREE2.PROG(J)*MLREE4.PROG(J) +XLP2 11 CONTINUE IF(IIMPI.EQ.1799) WRITE(IOIMP,FMT='('' DERNIER ALPHA , PROD SCAL * :'',1X,E12.5,'' , '',E12.5)')QQZ,XLP2 * * RECHERCHE PAR DICHOTOMIE SUR LES INDICES pour trouver l'interval * de alpha interessant * IF(XLP1 * XLP2.GE.0.D0) THEN ALPHA=ALPMAX IF (ALPMAX.LT.1E25) THEN KK = -3 II=IPB ENDIF GO TO 20 ENDIF 12 CONTINUE IK=(INPINF+INPSUP)/2 IF(IK.EQ.INPINF) GO TO 15 MLREE4=IVXU1 MLREE5=IVU1 MLREE6=IVN1 SEGSUP MLREE4,MLREE5,MLREE6 MLREE4=IVD1 MLREE5=IVGM1 MLREE6=IVGE1 MLREE7=IVXL1 SEGSUP MLREE4,MLREE5,MLREE6,MLREE7 IVXL1=0 IVXU1=0 IVU1=0 IVN1=0 IVD1=0 IVGM1=0 IVGE1=0 I=IK QQZ=PROG(I) DO 13 J=1,MLREE1.PROG(/1) MLREE3.PROG(J)=QQZ*MLREE2.PROG(J)+MLREE1.PROG(J) 13 CONTINUE CALL NTAPE1(MCP,MCQ,IVFP,IVFQ,IVLAM1,NVD,M,N,MVDU,MVDL,IVMINU, * IVMINL,IVMAXU,IVMAXL,IVU1,IVN1,IVD1,IVUL,IVLL,IVXU1,IVXL1) CALL NTAPE2(MCP,MCQ,IVXU1,IVXL1,IVB,N,M,IVGE1, *IVGM1,IVLAM1,IPBASE) XLP3=0.D0 MLREE4=IVGE1 DO 14 J=1,MLREE1.PROG(/1) XLP3=MLREE2.PROG(J)*MLREE4.PROG(J) +XLP3 14 CONTINUE IF(XLP1*XLP3.LE.0.D0) THEN XLP2=XLP3 INPSUP=I ELSE XLP1=XLP3 INPINF=I ENDIF GO TO 12 15 CONTINUE * * RECHERCHE PAR DICHOTOMIE pour determiner la valeur de alpha donnant le * min de l(lambda + alpha * grad) * IPP=0 ALPE=PROG(INPINF) ALGR=PROG(INPSUP) ALPHA=(ALGR + ALPE)/2.D0 60 CONTINUE IPP=IPP + 1 IF(IPP.GT.100) THEN CALL ERREUR ( 603) RETURN ENDIF ALPHA0=ALPHA MLREE4=IVXU1 MLREE5=IVU1 MLREE6=IVN1 SEGSUP MLREE4,MLREE5,MLREE6 MLREE4=IVD1 MLREE5=IVGM1 MLREE6=IVGE1 MLREE7=IVXL1 SEGSUP MLREE4,MLREE5,MLREE6,MLREE7 IVXL1=0 IVXU1=0 IVU1=0 IVN1=0 IVD1=0 IVGM1=0 IVGE1=0 * IF(IIMPI.EQ.1799) WRITE(IOIMP,FMT='('' ALPHA '',E12.5)')ALPHA QQZ=ALPHA DO 17 J=1,MLREE1.PROG(/1) MLREE3.PROG(J)=QQZ*MLREE2.PROG(J)+MLREE1.PROG(J) 17 CONTINUE CALL NTAPE1(MCP,MCQ,IVFP,IVFQ,IVLAM1,NVD,M,N,MVDU,MVDL,IVMINU, * IVMINL,IVMAXU,IVMAXL,IVU1,IVN1,IVD1,IVUL,IVLL,IVXU1,IVXL1) CALL NTAPE2(MCP,MCQ,IVXU1,IVXL1,IVB,N,M,IVGE1, *IVGM1,IVLAM1,IPBASE) XLP3=0.D0 MLREE4=IVGE1 DO 50 J=1,MLREE1.PROG(/1) XLP3=MLREE2.PROG(J)*MLREE4.PROG(J) +XLP3 50 CONTINUE IF(XLP1*XLP3.LE.0.D0) THEN XLP2=XLP3 ALGR=ALPHA ELSE XLP1=XLP3 ALPE=ALPHA ENDIF ALPHA =(ALGR + ALPE)/2 EPSI=1.D-10 IF(ALPHA.LT.1.D-12.AND.NDR.GE.1)THEN TERMIN=.TRUE. GO TO 20 ENDIF IF((ABS(ALPHA-ALPHA0)/ALPHA).GE.EPSI) GO TO 60 20 CONTINUE MLENT1=MAC2 MLENT2=MAC3 MLENT3=MORDRE IF(IIMPI.EQ.1799)WRITE(IOIMP,FMT='('' ALPHA '',E12.5)')ALPHA IF (KK.NE.-3)THEN * ****************** TEST SI ARRET SUR PLAN DE DISCONTINUITE ************ * ALPHA1=PROG(INPSUP) IF((ABS(ALPHA1-ALPHA)/ALPHA).LE.EPSI)THEN KK=MLENT2.LECT(MLENT3.LECT(INPSUP)) II=MLENT1.LECT(MLENT3.LECT(INPSUP)) ALPHA=ALPHA1 ENDIF ALPHA2=PROG(INPINF) IF(((ABS(ALPHA2-ALPHA)/ALPHA).LE.EPSI.AND.ALPHA2.NE.0.D0) * .OR.ALPHA.LT.1.D-12)THEN KK=MLENT2.LECT(MLENT3.LECT(INPINF)) II=MLENT1.LECT(MLENT3.LECT(INPINF)) ALPHA=ALPHA2 ENDIF ENDIF DO 21 J=1,MLREE1.PROG(/1) MLREE3.PROG(J)=ALPHA * MLREE2.PROG(J)+MLREE1.PROG(J) 21 CONTINUE IF (KK.EQ.-3)THEN MLREE3.PROG(II)=0. ENDIF ALPHA0=ALPHA MLREE4=IVXU1 MLREE5=IVU1 MLREE6=IVN1 SEGSUP MLREE4,MLREE5,MLREE6 MLREE4=IVD1 MLREE5=IVGM1 MLREE6=IVGE1 MLREE7=IVXL1 SEGSUP MLREE4,MLREE5,MLREE6,MLREE7 IVLAMB=IVLAM1 SEGSUP MLENT1,MLENT2,MLREEL,MLENT3 RETURN END