PROGRAM CDCDEV EXTERNAL SIMUL,PROSCA INTEGER N INTEGER IMP,IO,MODE,ITER,MEMAX,IZ(10),IZS(1) DIMENSION X(10),G(10),DX(10),RZ(52),RZS(1) DOUBLE PRECISION DZ(60),DZS(1) N=4 DO 50 I=1,N X(I)=5. 50 CONTINUE CALL SIMUL(2,N,X,F,G,IZ,RZ,DZ) DX(1)=1.E-4 DF1=1. EPS=1.E-4 ZERO=1.E-13 IMP=4 IO=6 ITER=100 NSIM=200 MEMAX=4 CALL M1FC1(SIMUL,PROSCA,N,X,F,G,DX,DF1,EPS,ZERO,IMP,IO,MODE,ITER 1,NSIM,MEMAX,IZ,RZ,DZ,IZS,RZS,DZS) END SUBROUTINE SIMUL(INDIC,N,X,F,G,IS,RS,DS) INTEGER INDIC,N,IS(1) REAL X(1),F,G(1),RS(1) DOUBLE PRECISION DS(1) IF(INDIC.GT.1) GO TO 100 PRINT 10 10 FORMAT(5X,'INDIC=1'/) RETURN 100 F=0.0 DO 50 I=1,N F=F+I*X(I)**2 G(I)=I*X(I) 50 CONTINUE F=F/2 F1=0.5*X(1)**2+2.5*X(2)**2+2.5*X(3)**2+5.0*X(4)**2 F1=F1/2. IF(F.GT.F1) RETURN G(1)=0.5*X(1) G(2)=2.5*X(2) G(3)=2.5*X(3) G(4)=5.0*X(4) RETURN END SUBROUTINE PROSCA(N,X,Y,PS,IZS,RZS,DZS) DIMENSION X(1),Y(1),IZS(1),RZS(1) DOUBLE PRECISION PS,DZS(1) PS=0.D0 DO 10 I=1,N 10 PS=PS+DBLE(X(I))*DBLE(Y(I)) RETURN END SUBROUTINE M1FC1 (SIMUL,PROSCA,N,XN,FN,G,DXMIN,DF1,EPSF, 1 ZERO,IMP,IO,MODE,ITER,NSIM,MEMAX,IZ,RZ,DZ,IZS,RZS,DZS) C DIMENSION IZ=2*(MEMAX+1) C DIMENSION RZ=5*N+(N+4)*MEMAX C DIMENSION DZ=(MEMAX+9)*MEMAX+8 EXTERNAL SIMUL,PROSCA DIMENSION IZ(*),RZ(*),XN(N),G(N),IZS(*),RZS(*) DOUBLE PRECISION DZS(*),DZ(*) IF (N.GT.0 .AND. DF1.GT.0. 1 .AND. EPSF.GE.0. .AND. ZERO.GE.0. 1 .AND. ITER.GE.0 .AND. NSIM.GE.0 1 .AND. MEMAX.GE.1 .AND. DXMIN.GT.0.) GO TO 10 MODE=2 WRITE (IO,1001) 1001 FORMAT (25H M1FC1 APPEL INCOHERENT) GO TO 999 10 NS=1 NGD=NS+N NX=NGD+N NSA=NX+N NGG=NSA+N NAL=NGG+N NAPS=NAL+MEMAX NANC=NAPS+MEMAX NPOIDS=NANC+MEMAX NQ=NPOIDS+MEMAX NJC=1 NIC=NJC+MEMAX+1 NR=1 NA=NR+(MEMAX+1)*(MEMAX+1) NE=NA+MEMAX+1 NRR=NE+MEMAX+1 NXGA=NRR+MEMAX+1 NY=NXGA+MEMAX+1 NW1=NY+MEMAX+1 NW2=NW1+MEMAX+1 C NIZ=2*(MEMAX+1) NRZ=NQ+N*MEMAX-1 NDZ=NW2+MEMAX IF (IMP.GT.0) WRITE (IO,1000) N,MEMAX,NIZ,NRZ,NDZ 1000 FORMAT (22H ENTREE DANS M1FC1. N=,I4,8H MEMAX=,I3/ 122H DIMENSIONS MINIMALES,2X,3HIZ(,I4,8H) RZ(,I6,8H) DZ(, 1I6,1H)/) DO 110 I=1,NIZ 110 IZ(I)=0 DO 120 I=1,NRZ 120 RZ(I)=0. DO 130 I=1,NDZ 130 DZ(I)=0.D0 CALL M1FC1A(SIMUL,PROSCA,N,MODE,XN,FN,G,DF1,EPSF,DXMIN,IMP,ZERO, 1 IO,NTOT,ITER,NSIM,MEMAX,RZ(NS),RZ(NGD),RZ(NX),RZ(NSA), 1 RZ(NGG),RZ(NAL),RZ(NAPS),RZ(NANC),RZ(NPOIDS),RZ(NQ), 1 IZ(NJC),IZ(NIC),DZ(NR),DZ(NA),DZ(NE),DZ(NRR),DZ(NXGA),DZ(NY), 1 DZ(NW1),DZ(NW2),IZS,RZS,DZS) IZ(1)=NTOT 999 RETURN END SUBROUTINE M1FC1A(SIMUL,PROSCA,N,MODE,XN,FN,G,DF0,EPS0,DX,IMP, 1 ZERO,IO,NTOT,ITER,NSIM,MEMAX,S,GD,X,SA,GG,AL,APS,ANC,POIDS, 1 Q,JC,IC,R,A,E,RR,XGA,Y,W1,W2,IZS,RZS,DZS) C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C MINIMISATION OF A FUNCTION HEMIDERIVABLE BY A BUNDLE METHOD . C THE DIRECTION IS OBTAINED BY PROJECTION OF ORIGIN C ON A POLYHEDRE GENERATED BY AN ENSEMBLE OF GRADIENTS ALREADY CALCULATED C AND THE LINE SEARCH YIELdS EITHER A DESCENT STEP OR A NUL STEP. C THE ALGORITHM MINIMIZES F A EPS0 PRES (IF CONVEXITE) C WHERE EPS0 IS A TOLERANCE PROVIDE BY THE USER. C C MODE C <=0 MODE=INDIC DE SIMUL C 1 FIN NORMALE C 2 APPEL INCOHERENT C 3 REDUIRE L'ECHELLE DES X C 4 MAX ITERATIONS C 5 MAX SIMULATIONS C 6 IMPOSSIBLE D'ALLER AU DELA DE DX C 7 EPRF2 MIS EN ECHEC C 8 ON COMMENCE A BOUCLER C IMP C <0 INDIC=1 TOUTES LES -IMP ITERATIONS C 0 PAS D'IMPRESSIONS C 1 IMPRESSIONS INITIALES ET FINALES C 2 IMPRESSIONS A CHAQUE CONVERGENCE C 3 UNE IMPRESSION PAR ITERATION C 4 INFORMATIONS M1FC1 ET MLIS2 C >4 DEBUG C 5 TOLERANCES DIVERSES C 6 POIDS C >6 EPRF2 C -------------------------------------------------- C FAIT APPEL AUX SUBROUTINES SUIVANTES: C --------SUBROUTINE EPRF2 (CALCUL DE LA DIRECTION) C --------SUBROUTINES EREMF1 ET EFINF1 (ESCLAVES DE EPRF2) C --------SUBROUTINE ERDF1 (REDUCTION DU FAISCEAU) C --------SUBROUTINE MLIS2 (RECHERCHE LINEAIRE) C --------SUBROUTINE SIMUL (MODULE DE SIMULATION) C --------SUBROUTINE PROSCA (PRODUIT DE DUALITE DONNANT LE GRADIENT) C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC EXTERNAL SIMUL,PROSCA DIMENSION XN(N),G(N),IZS(*),RZS(*),X(N),GD(N),GG(N) DIMENSION S(N),SA(N),JC(*),IC(*) DIMENSION Q(*),AL(MEMAX),APS(MEMAX),ANC(MEMAX),POIDS(MEMAX) DOUBLE PRECISION DZS(*),R(*),A(*),E(*),RR(*),XGA(*),Y(*),W1(*), 1W2(*),PS 1000 FORMAT (19H M1FC1 ITER NSIM,6X,2HFN,11X,3HEPS,7X,2HS2, 19X,1HU,5X,2HNV) 1001 FORMAT (25H M1FC1 APPEL INCOHERENT) 1002 FORMAT(/6H M1FC1,21H TABLEAU DES POIDS/ 1 (6H M1FC1,3X,7E10.3)) 1003 FORMAT(/6H M1FC1,17H FIN PAR INDIC=,I2,11H DANS SIMUL) 1004 FORMAT(6H M1FC1,I7,I5,E16.7,16H CONVERGENCE A,E10.3,5H PRES, 13H (,E9.2,1H)) 1005 FORMAT(6H M1FC1,I7,I5,E16.7,20H FAISCEAU REDUIT A, 1 I3,10H GRADIENTS) 1006 FORMAT (/22H M1FC1 FIN SUR NSIM) 1007 FORMAT(6H M1FC1,3X,I4,I5,2X,E14.7,3E10.2,I3) 1008 FORMAT(55H M1FC1 LA NORME INITIALE DU GRADIENT (D'APRES PROSCA), 1 10H EST NULLE) 1009 FORMAT (6H M1FC1,10X,6HLOGIC=,I2,4X,3HRO=,E10.3, 1 4X,4HTPS=,E10.3,4X,4HTNC=,E10.3) 1010 FORMAT (6H M1FC1,12X,6HDIAM2=,E10.3,4X,5HETA2=,E10.3,4X, 1 3HAP=,E10.3) 1011 FORMAT(/37H M1FC1 LA DIRECTION NE PIVOTE PLUS) 1012 FORMAT(/24H M1FC1 FIN SUR ITER =,I4) 1013 FORMAT(/31H M1FC1 FIN ANORMALE DE EPRF2) 1014 FORMAT(/23H M1FC1 FIN SUR DXMIN) 1015 FORMAT(/53H M1FC1 ATTENTION ON BUTE SUR TMAX, REDUIRE L'ECHELLE) 1016 FORMAT(/21H M1FC1 FIN NORMALE) 1017 FORMAT (1X) 1018 FORMAT (/25H M1FC1 FIN SUR INDIC=0) 1019 FORMAT (/6H M1FC1,30X,12HETA2 FORCE A,E10.3) C C INITIALISATIONS C ITMAX=ITER ITER=0 ITIMP=0 NAPMAX=NSIM NSIM=1 LOGIC=1 LOGIC2=0 TMAX=1.E20 EPS=DF0 EPSM=EPS DF=DF0 MODE=1 NTOT=0 IFLAG=0 C C INITIALISATION DU FAISCEAU C CALCUL DU DIAMETRE DE L'EPURE ET DU TEST D'ARRET C APS(1)=0. ANC(1)=0. POIDS(1)=0. NTA=0 KGRAD=1 MEMAX1=MEMAX+1 DO 50 I=1,N 50 Q(I)=-G(I) CALL PROSCA (N,G,G,PS,IZS,RZS,DZS) IF (PS.GT.0.D0) GO TO 60 MODE=2 IF (IMP.NE.0) WRITE(IO,1008) GO TO 900 60 DIAM2=100.*DF0*DF0/SNGL(PS) ETA2=1.E-2*EPS0*EPS0/DIAM2 AP=ZERO*DF0/DIAM2 IF(IMP.GT.2) WRITE (IO,1000) C C BOUCLE C 100 ITER=ITER+1 ITIMP=ITIMP+1 IF(ITER.LT.ITMAX) GO TO 110 IF(IMP.GT.0) WRITE (IO,1012) ITER MODE=4 GO TO 900 110 NTOT=NTOT+1 IF(LOGIC.EQ.3) RO=RO*SQRT(S2) IF (ITIMP.NE.-IMP) GO TO 200 ITIMP=0 INDIC=1 CALL SIMUL(INDIC,N,XN,F,G,IZS,RZS,DZS) C C CALCUL DE LA DIRECTION C 200 EPS=AMIN1(EPS,EPSM) EPS=AMAX1(EPS,EPS0) CALL EREMF1 (PROSCA,IFLAG,N,NTOT,NTA,MEMAX1,Q,POIDS,E,A,R, 1 IZS,RZS,DZS) CALL EPRF2 (IFLAG,NTOT,NV,IO,ZERO,S2,EPS,AL, 1 IMP,U,ETA2,MEMAX1,JC,IC,R,A,E,RR,XGA,Y,W1,W2) C C FIN ANORMALE DE EPRF2 C IF(IFLAG.EQ.0) GO TO 250 IF(IMP.GT.0) WRITE (IO,1013) MODE=7 GO TO 900 250 NTA=NTOT CALL EFINF1 (N,NV,JC,XGA,Q,S) U=AMAX1(U,0.) S2=AMAX1(S2,0.) C C TESTS D'ARRET (NB. SI NR G EST INTERIEUR A G C ALORS NR G EST "NUL") C IF (NV.LT.N+2) GO TO 260 ETA2=AMAX1(ETA2,10.*S2) IF (IMP.GE.2) WRITE(IO,1019) ETA2 260 IF (S2.GT.ETA2) GO TO 300 C C CALCUL DE LA PRECISION Z=0. DO 270 K=1,NV J=JC(K)-1 IF (J.GT.0) Z=Z+SNGL(XGA(K))*POIDS(J) 270 CONTINUE EPSM=AMIN1(EPSM,Z) IF (IMP.GE.2) WRITE (IO,1004) ITER,NSIM,FN,EPSM,S2 IF (EPSM.GT.EPS0) GO TO 280 MODE=1 IF (IMP.GT.0) WRITE (IO,1016) GO TO 900 C C DIMINUTION DE EPSILON 280 EPSM=AMAX1(0.1*EPSM,EPS0) EPS=EPSM IF (LOGIC.EQ.3) TOL=0.01*EPS IFLAG=2 GO TO 200 C C SUITE DES ITERATIONS C IMPRESSIONS C 300 IF (IMP.GT.3) WRITE (IO,1017) IF (IMP.GT.2) WRITE (IO,1007) ITER,NSIM,FN,EPS,S2,U,NV IF (IMP.GE.6) WRITE (IO,1002) (POIDS(I),I=1,NTOT) C TEST DE NON-PIVOTAGE IF (LOGIC.NE.3) GO TO 350 Z=0. DO 310 I=1,N Z1=S(I)-SA(I) 310 Z=Z+Z1*Z1 IF(Z.GT.10.*ZERO*ZERO*S2) GO TO 350 IF(IMP.GT.0) WRITE (IO,1011) MODE=8 GO TO 900 C C RECHERCHE LINEAIRE C 350 IFLAG=3 S3=S2+U*EPS IF (LOGIC.EQ.3) GO TO 365 RO=2.*DF/S3 TOL=0.01*EPS GO TO 370 365 RO=RO/SQRT(S2) TOL=AMAX1(0.6*TOL,0.01*EPS0) 370 FA=FN ALFA=0.2 BETA=0.1 FPN=-S3 IF (MEMAX.EQ.1) TOL=0. C CALCUL DE LA RESOLUTION MINIMALE, FONCTION DE DX TMIN=0. DO 372 I=1,N 372 TMIN=AMAX1(TMIN,ABS(S(I)/DX)) TMIN=1./TMIN IF (ITER.EQ.1) ROA=RO CALL MLIS2(SIMUL,PROSCA,N,XN,FN,FPN,RO,TMIN,TMAX,S,S2,G,GD,ALFA, 1 BETA,IMP,IO,LOGIC,NSIM,NAPMAX,X,TOL,AP,TPS,TNC,GG,IZS,RZS,DZS) IF (LOGIC.EQ.0 .OR. LOGIC.EQ.2 .OR. LOGIC.EQ.3) GO TO 380 C SORTIE PAR ANOMALIE DANS MLIS2 IF (IMP.LE.0) GO TO 375 IF (LOGIC.EQ.6 .OR. LOGIC.LT.0) WRITE (IO,1014) IF (LOGIC.EQ.4) WRITE (IO,1006) IF (LOGIC.EQ.5) WRITE (IO,1018) IF (LOGIC.EQ.1) WRITE (IO,1015) 375 IF (LOGIC.EQ.1) MODE=3 IF (LOGIC.EQ.4) MODE=5 IF (LOGIC.EQ.5) MODE=0 IF (LOGIC.EQ.6) MODE=6 IF (LOGIC.LT.0) MODE=LOGIC GO TO 900 380 IF (LOGIC.NE.3) GO TO 385 DO 382 I=1,N 382 SA(I)=S(I) 385 IF (ITER.GT.1) GO TO 390 C C 1ERE ITERATION, AJUSTEMENT DE AP, DIAM ET ETA IF (LOGIC.EQ.0) TPS=(FN-FA)-RO*FPN AP=ZERO*ZERO*ABS(TPS)/(S2*RO*RO) AJUST=RO/ROA IF (LOGIC.NE.3) DIAM2=DIAM2*AJUST*AJUST IF (LOGIC.NE.3) ETA2=ETA2/(AJUST*AJUST) IF (IMP.GE.2) WRITE (IO,1010) DIAM2,ETA2,AP 390 MM=MEMAX-1 IF (LOGIC.EQ.2) MM=MEMAX-2 IF (NTOT.LE.MM) GO TO 400 C C REDUCTION DU FAISCEAU POUR ENTRER LE NOUVEL ELEMENT C CALL ERDF1(PROSCA,N,NTOT,MM,KGRAD,AL,Q,S,POIDS,APS,ANC,MEMAX1, 1 R,E,IC,IZS,RZS,DZS) IFLAG=1 NTA=NTOT IF (IMP.GE.2) WRITE (IO,1005) ITER,NSIM,FN,NTOT C 400 IF (IMP.GE.5) WRITE (IO,1009) LOGIC,RO,TPS,TNC IF (LOGIC.EQ.3) GO TO 500 C C ITERATION DE DESCENTE C IFLAG=MIN0(IFLAG,2) DF=FA-FN IF (NTOT.EQ.0) GO TO 500 C C ACTUALISATION DES POIDS C S3N=RO*SQRT(S2) DO 430 K=1,NTOT NK=(K-1)*N Z=0. DO 420 I=1,N NKI=NK+I C Y(K) CONTIENT -(GK,S) (ATTENTION A L'ARTIFICE) 420 Z=Z+Q(NKI)*S(I) Y(K)=Z Z1=ABS(APS(K)+(-DF+RO*Z)) Z2=ANC(K)+S3N POIDS(K)=AMAX1(Z1,AP*Z2*Z2) APS(K)=Z1 ANC(K)=Z2 430 CONTINUE C C ACTUALISATION DE EPS C EPS=RO*S3 KGRAD=NTOT+1 C C NOUVEL ELEMENT DU FAISCEAU (POUR LES TROIS TYPES DE PAS) C 500 NT1=NTOT+1 IF (LOGIC.EQ.3) GO TO 510 APS(NT1)=0. ANC(NT1)=0. POIDS(NT1)=0. GO TO 520 510 APS(NT1)=TPS ANC(NT1)=SQRT(TNC) POIDS(NT1)=AMAX1(TPS,AP*TNC) 520 NK=NTOT*N DO 530 I=1,N NKI=NK+I 530 Q(NKI)=-G(I) C C TRAITEMENT POUR LOGIC=2 (ON AJOUTE ENCORE UN GRADIENT) IF(LOGIC.NE.2) GO TO 550 NTOT=NTOT+1 LOGIC=3 LOGIC2=1 DO 540 I=1,N 540 G(I)=GD(I) GO TO 390 550 LOGIC=LOGIC-LOGIC2 LOGIC2=0 GO TO 100 C C EPILOGUE C 900 IF (ITER.LE.1) GO TO 990 DO 910 I=1,N 910 G(I)=-S(I) 990 RETURN END SUBROUTINE MLIS2 (SIMUL,PROSCA,N,XN,FN,FPN,T,TMIN,TMAX,D,D2,G,GD, 1 AMD,AMF,IMP,IO,LOGIC,NAP,NAPMAX,X,TOL,A,TPS,TNC,GG,IZS,RZS,DZS) C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C SUBROUTINE EFFECTUANT UNE RECHERCHE LINEAIRE SUR 0 TMAX C PARTANT DU POINT XN DANS LA DIRECTION D. C SOUS L'HYPOTHESE D'HEMIDERIVABILITE, DONNE C UN PAS SERIEUX, BLOQUE, NUL OU SEMI SERIEUX-NUL (2 GRADIENTS). C NECESSITE FPN < 0 ESTIMANT LA DERIVEE A L'ORIGINE. C APPELLE SIMUL SYSTEMATIQUEMENT AVEC INDIC = 4 C C LOGIC C 0 DESCENTE SERIEUSE C 1 DESCENTE BLOQUEE C 2 PAS SEMISERIEUX-NUL C 3 PAS NUL, ENRICHISSEMENT DU FAISEAU C 4 NAP > NAPMAX C 5 RETOUR A L'UTILISATEUR C 6 NON HEMI-DERIVABLE (AU-DELA DE DX) C < 0 CONTRAINTE IMPLICITE ACTIVE C C IMP C =0 PAS D'IMPRESSIONS C >0 MESSAGE EN CAS DE FIN ANORMALE C >3 INFORMATIONS POUR CHAQUE ESSAI DE T C ---------------------------------------- C FAIT APPEL AUX SUBROUTINES: C -------SIMUL(INDIC,N,X,F,G,IZS,RZS,DZS) C -------PROSCA(N,U,V,PS,IZS,RZS,DZS) C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC EXTERNAL SIMUL,PROSCA DIMENSION XN(N),D(N),G(N),X(N),IZS(*),RZS(*),GG(N),GD(N) DOUBLE PRECISION PS,Z1,DZS(*) 1000 FORMAT (/4X,9H MLIS2 ,4X,4HFPN=,E10.3,4H D2=,E9.2, 1 7H TMIN=,E9.2,6H TMAX=,E9.2) 1001 FORMAT (/4X,6H MLIS2,10X,17HTMIN FORCE A TMAX) 1002 FORMAT (4X,6H MLIS2,36X,1HI,E10.3,2E11.3) 1003 FORMAT (4X,6H MLIS2,E13.3,2E11.3,2H I) 1004 FORMAT (4X,6H MLIS2,36X,1HI,E10.3,7H INDIC=,I3) 1006 FORMAT (4X,6H MLIS2,3X,20HCONTRAINTE IMPLICITE,I4,7H ACTIVE) 1007 FORMAT (/4X,6H MLIS2,3X,12HFIN SUR TMIN) 1010 FORMAT (/4X,6H MLIS2,3X,I5,22H SIMULATIONS ATTEINTES) 1011 FORMAT (/4X,6H MLIS2,3X,31HARRET DEMANDE PAR L'UTILISATEUR) C C INITIALISATIONS C TESF=AMF*FPN TESD=AMD*FPN TD=0. TG=0. FG=FN FPG=FPN TA=0. FA=FN FPA=FPN INDICA=1 LOGIC=0 C ELIMINATION D'UN T INITIAL RIDICULEMENT PETIT IF (T.GT.TMIN) GO TO 20 T=TMIN IF (T.LE.TMAX) GO TO 20 IF (IMP.GT.0) WRITE (IO,1001) TMIN=TMAX 20 IF (FN+T*FPN.LT.FN+0.9*T*FPN) GO TO 30 T=2.*T GO TO 20 C 30 IF(T.LT.TMAX) GO TO 40 T=TMAX LOGIC=1 40 IF (IMP.GE.4) WRITE (IO,1000) FPN,D2,TMIN,TMAX DO 50 I=1,N 50 X(I)=XN(I)+T*D(I) C C BOUCLE C 100 NAP=NAP+1 IF(NAP.LE.NAPMAX) GO TO 150 C SORTIE PAR MAXIMUM DE SIMULATIONS LOGIC=4 IF(IMP.GE.4) WRITE(IO,1010) NAP IF (TG.EQ.0.) GO TO 999 FN=FG DO 120 I=1,N G(I)=GG(I) 120 XN(I)=XN(I)+TG*D(I) GO TO 999 150 INDIC=4 CALL SIMUL(INDIC,N,X,F,G,IZS,RZS,DZS) IF(INDIC.NE.0) GO TO 200 C C ARRET DEMANDE PAR L'UTILISATEUR LOGIC=5 FN=F DO 170 I=1,N 170 XN(I)=X(I) IF(IMP.GE.4) WRITE(IO,1011) GO TO 999 C C LES TESTS ELEMENTAIRES SONT FAITS, ON Y VA C TOUT D'ABORD, OU EN SOMMES NOUS ? C 200 IF(INDIC.GT.0) GO TO 210 TD=T INDICD=INDIC LOGIC=0 IF (IMP.GE.4) WRITE(IO,1004) T,INDIC T=TG+0.1*(TD-TG) GO TO 905 C C CALCUL DE LA DERIVEE DIRECTIONNELLE H'(T) 210 CALL PROSCA(N,G,D,PS,IZS,RZS,DZS) FP=SNGL(PS) C C TEST DE DESCENTE (PREMIERE INEGALITE POUR UN PAS SERIEUX) FFN=F-FN IF(FFN.LT.T*TESF) GO TO 300 TD=T FD=F FPD=FP DO 230 I=1,N 230 GD(I)=G(I) INDICD=INDIC LOGIC=0 IF(IMP.GE.4) WRITE(IO,1002) T,FFN,FP IF(TG.NE.0.) GO TO 500 C TESTS POUR UN PAS NUL (SI TG=0) IF(FPD.LT.TESD) GO TO 500 TPS=(FN-F)+TD*FPD TNC=D2*TD*TD P=AMAX1(A*TNC,TPS) IF(P.GT.TOL) GO TO 500 LOGIC=3 GO TO 999 C C DESCENTE 300 IF(IMP.GE.4) WRITE(IO,1003) T,FFN,FP C C TEST DE DERIVEE (DEUXIEME INEGALITE POUR UN PAS SERIEUX) IF(FP.LT.TESD) GO TO 320 C C SORTIE, LE PAS EST SERIEUX LOGIC=0 FN=F FPN=FP DO 310 I=1,N 310 XN(I)=X(I) GO TO 999 C 320 IF (LOGIC.EQ.0) GO TO 350 C C SORTIE PAR DESCENTE BLOQUEE FN=F FPN=FP DO 330 I=1,N 330 XN(I)=X(I) GO TO 999 C C ON A UNE DESCENTE 350 TG=T FG=F FPG=FP DO 360 I=1,N 360 GG(I)=G(I) C IF(TD.NE.0.) GO TO 500 C EXTRAPOLATION TA=T T=9.*TG Z=FPN+3.*FP-4.*FFN/TG IF(Z.GT.0.) T=AMIN1(T,TG*AMAX1(1.,-FP/Z)) T=TG+T IF(T.LT.TMAX) GO TO 900 LOGIC=1 T=TMAX GO TO 900 C C INTERPOLATION C 500 IF(INDICA.GT.0 .AND. INDICD.GT.0) GO TO 510 TA=T T=0.9*TG+0.1*TD GO TO 900 510 TEST=0.1*(TD-TG) C APPROXIMATION CUBIQUE PS=DBLE(FP)+DBLE(FPA)-3.D0*DBLE(FA-F)/DBLE(TA-T) Z1=PS*PS-DBLE(FP)*DBLE(FPA) IF (Z1.GE.0.D0) GO TO 520 IF (FP.LT.0.) TC=TD IF (FP.GE.0.) TC=TG GO TO 600 520 Z1=DSQRT(Z1) IF (T-TA.LT.0.) Z1=-Z1 SIGN=(T-TA)/ABS(T-TA) IF ((SNGL(PS)+FP)*SIGN.GT.0.) GO TO 550 DEN=2.D0*PS+DBLE(FP)+DBLE(FPA) ANUM=Z1-DBLE(FP)-PS IF (ABS((T-TA)*ANUM).GE.(TD-TG)*ABS(DEN)) GO TO 530 TC=T+ANUM*(TA-T)/DEN GO TO 600 530 TC=TD GO TO 600 550 TC=T+FP*(TA-T)/SNGL(PS+DBLE(FP)+Z1) 600 MC=0 IF (TC.LT.TG) MC=-1 IF (TC.GT.TD) MC=1 TC=AMAX1(TC,TG+TEST) TC=AMIN1(TC,TD-TEST) C APPROXIMATION POLYHEDRIQUE PS=DBLE(FPD)-DBLE(FPG) IF (PS.NE.0.D0) GO TO 620 TP=0.5*(TD+TG) GO TO 650 620 TP=((DBLE(FG)-DBLE(FPG)*DBLE(TG))- 1 (DBLE(FD)-DBLE(FPD)*DBLE(TD)))/PS 650 MP=0 IF (TP.LT.TG) MP=-1 IF (TP.GT.TD) MP=1 TP=AMAX1(TP,TG+TEST) TP=AMIN1(TP,TD-TEST) C NOUVEAU T PAR APPROXIMATION CP COMPLETE SECURISEE TA=T IF (MC.EQ.0 .AND. MP.EQ.0) T=AMIN1(TC,TP) IF (MC.EQ.0 .AND. MP.NE.0) T=TC IF (MC.NE.0 .AND. MP.EQ.0) T=TP IF (MC.EQ.1 .AND. MP.EQ.1) T=TD-TEST IF (MC.EQ.-1 .AND. MP.EQ.-1) T=TG+TEST IF (MC*MP.EQ.-1) T=0.5*(TG+TD) C C FIN DE BOUCLE C 900 FA=F FPA=FP 905 INDICA=INDIC C PEUT-ON FAIRE LOGIC=2 ? IF (TD.EQ.0.) GO TO 920 IF (INDICD.LT.0) GO TO 920 IF (TD-TG.GT.10.*TMIN) GO TO 920 IF (FPD.LT.TESD) GO TO 920 TPS=(FG-FD)+(TD-TG)*FPD TNC=D2*(TD-TG)*(TD-TG) P=AMAX1(A*TNC,TPS) IF(P.GT.TOL) GO TO 920 C SORTIE PAR PAS SEMISERIEUX-NUL LOGIC=2 FN=FG FPN=FPG T=TG DO 910 I=1,N XN(I)=XN(I)+TG*D(I) 910 G(I)=GG(I) GO TO 999 C C TEST D'ARRET SUR LA PROXIMITE DE TG ET TD C 920 IF (TD.EQ.0.) GO TO 990 IF (TD-TG.LE.TMIN) GO TO 950 DO 930 I=1,N Z=XN(I)+T*D(I) IF (Z.NE.X(I) .AND. Z.NE.XN(I)) GO TO 990 930 CONTINUE C ARRET SUR DX OU DE SECOURS 950 LOGIC=6 IF (INDICD.LT.0) LOGIC=INDICD IF (TG.EQ.0.) GO TO 970 FN=FG DO 960 I=1,N XN(I)=XN(I)+TG*D(I) 960 G(I)=GG(I) 970 IF (IMP.LE.0) GO TO 999 IF (LOGIC.LT.0) WRITE(IO,1006) LOGIC IF (LOGIC.EQ.6) WRITE(IO,1007) GO TO 999 C C RECOPIAGE DE X ET BOUCLE 990 DO 995 I=1,N 995 X(I)=XN(I)+T*D(I) GO TO 100 C 999 RETURN END SUBROUTINE ERDF1(PROSCA,N,NTOT,NINF,KGRAD, 1 AL,Q,S,POIDS,APS,ANC,MM1,R,E,IC,IZS,RZS,DZS) DIMENSION AL(NTOT),Q(*),POIDS(NTOT),APS(NTOT),ANC(NTOT),IC(MM1), 1 S(N),IZS(*),RZS(*) DOUBLE PRECISION E(MM1),R(*),PS,DZS(*) C C THIS SUBROUTINE REDUCES A BUNDLE C OF SIZE NTOT IN RN C TO A SIZE NO GREATER THAN NINF C IF(NTOT.LE.NINF) GO TO 900 IF (NINF.GT.0) GO TO 100 C C PURE GRADIENT METHOD C NTOT=0 KGRAD=0 GO TO 900 C C REDUCTION TO THE CORRAL 100 NT1=0 DO 150 J=1,NTOT IF(AL(J).EQ.0. .AND. POIDS(J).NE.0.) GO TO 150 NT1=NT1+1 IC(NT1)=J IF(J.EQ.NT1) GO TO 130 NJ=N*(J-1) NN=N*(NT1-1) DO 110 I=1,N NN=NN+1 NJ=NJ+1 110 Q(NN)=Q(NJ) AL(NT1)=AL(J) POIDS(NT1)=POIDS(J) APS(NT1)=APS(J) ANC(NT1)=ANC(J) E(NT1+1)=E(J+1) 130 IF (POIDS(J).EQ.0.) KGRAD=NT1 NN=NT1*MM1+1 NJ=J*MM1+1 DO 140 K=1,NT1 NJK=NJ+IC(K) NN=NN+1 140 R(NN)=R(NJK) 150 CONTINUE NTOT=NT1 IF(NTOT.LE.NINF) GO TO 900 C C CORRAL STILL TOO LARGE C SAVE THE NEAR C CALL PROSCA (N,S,S,PS,IZS,RZS,DZS) E(2)=1.D0 Z=0. Z1=0. Z2=0. DO 370 K=1,NTOT Z1=Z1+AL(K)*APS(K) Z2=Z2+AL(K)*ANC(K) 370 Z=Z+AL(K)*POIDS(K) APS(1)=Z1 ANC(1)=Z2 POIDS(1)=Z R(MM1+2)=PS IF (NINF.GT.1) GO TO 400 NTOT=1 KGRAD=0 DO 380 I=1,N 380 Q(I)=S(I) GO TO 900 C SAVE THE GRADIENT 400 NN=(KGRAD-1)*N DO 470 I=1,N NJ=N+I NN=NN+1 Q(NJ)=Q(NN) 470 Q(I)=S(I) CALL PROSCA(N,Q(N+1),S,PS,IZS,RZS,DZS) E(3)=1.D0 NN=(MM1+1)*KGRAD+1 R(2*MM1+3)=R(NN) R(2*MM1+2)=PS APS(2)=0. ANC(2)=0. POIDS(2)=0. KGRAD=2 NTOT=2 900 RETURN END SUBROUTINE EPRF2 (IFLAG,NTOT,NV,IO,ZERO,S2, 1 EPS,AL,IMP,U,ETA,MM1,JC,IC,R,A,E,RR,XPR,Y,W1,W2) COMMON /EPRF2C/ U1,NC C THE DIMENSION IS MM1*MM1 FOR R DIMENSION AL(NTOT),JC(MM1),IC(MM1) DOUBLE PRECISION A(MM1),E(MM1),R(*),RR(MM1),XPR(MM1), 1 Y(MM1),W1(MM1),W2(MM1) C C ***** ON ENTRY ***** C C IFLAG=0-1 INITIALIZE ON ONE SUBGRADIENT (MU IN) C C IFLAG=2 " " " " " " " C AND STRIVE TO ENTER BY PRIORITY THE C POINTS OF THE PREVIOUS CORRAL AT THE C BEGINNING OF THE ITERATIONS. C C IFLAG=3 INITIALIZE ON THE PREVIOUS PROJECTION C (WITH ITS CORRESPONDING CORRAL) C C C ***** ON EXIT ***** C C IFLAG=0 NORMAL END C C 1 OLD SOLUTION IS ALREADY OPTIMAL C C 2 CONSTRAINTS NON FEASIBLE C C 3 TRYING TO ENTER A VARIABLE C THAT IS ALREADY IN THE CORRAL C C 4 STARTING TO LOOP C C C C C IMP > 5 ONE PRINTS FINAL INFORMATION C C C IMP > 6 ONE PRINTS INFORMATION AT EACH ITERATION C C C IMP > 7 ONE PRINTS ALSO C C - AT EACH ITERATION THE CHOLESKI MATRIX C - AND THE INITIAL INFORMATION SUCH AS (PI,PJ) ... C C C DOUBLE PRECISION PS0,PS,PS1,PS2,PS12,DET,TOL,V1,V2,U1,U2, 1DEPS,DZERO,TETA,DMU,W1S,W2S,W12S 1001 FORMAT(27H EPSILON SMALLER THAN A) 1003 FORMAT(3H A=,10D10.3,/(6X,10D10.3)) 1004 FORMAT(7H (G,G)=,10D10.3,/(7X,10D10.3)) 1005 FORMAT(27H START WITH VARIABLES 1 AND,I4) 1006 FORMAT(7H (S,S)=,E12.4,10H VARIABLE,I4, 12H (,E12.4,12H) COMING IN.) 1007 FORMAT(9H VARIABLE,I4,2H (,I4,3H) =,D11.3,11H GOING OUT., 1 17H FEASIBLE (S,S)=,D11.4,12H UNFEASIBLE=,D11.4) 1008 FORMAT(15H INITIAL CORRAL/(20I6)) 1010 FORMAT(12H EPSILON =,D10.3) 1011 FORMAT(3H X=,10D11.3,/(3X,10D11.3)) 1012 FORMAT(10H CHOLESKI,,10D11.3,/(10X,10D11.3)) 1013 FORMAT(22H DUPLICATE VARIABLE ,I3) 1014 FORMAT(14H FINISHED WITH,I3,10H GRADIENTS,I3, 111H VARIABLES./ 1 7H (S,S)=,E11.4,6H TEST=,E11.4/ 1 32H COST OF THE EXTRA CONSTRAINT U=,D12.5) 1015 FORMAT(20I6) 1016 FORMAT(28H EPRF2 IS APPARENTLY LOOPING) 1018 FORMAT(//) 1019 FORMAT(47H ERROR FROM EPRF2. OLD SOLUTION ALREADY OPTIMAL) 1020 FORMAT(7H (S,S)=,E12.4,5H U1=,D12.3,23H VARIABLE 1 COMING IN.) C C C **** BEGIN **** C C PREPARE VARIOUS DATA C C NITER=0 NT1=NTOT+1 ITMAX=10*NTOT DEPS=DBLE(EPS) INCR=0 K00=1 W1S=0.D0 W2S=0.D0 W12S=0.D0 GAMA=.99 DZERO=10.D0*DBLE(ZERO*ZERO) C INITIAL PRINTOUTS IF(IMP.LE.7) GO TO 100 WRITE(IO,1003) (A(J),J=1,NT1) WRITE(IO,1010) DEPS DO 95 J=1,NT1 MEJ=(J-1)*MM1 WRITE(IO,1004) (R(MEJ+JJ),JJ=1,J) 95 CONTINUE C C INITIAL POINT C 100 IF(IFLAG.NE.3) GO TO 110 IF(IMP.GT.6) WRITE(IO,1008) (JC(K),K=1,NV) J0=NT1 PS=U1*(A(NT1)-DEPS) MENT=(NT1-1)*MM1 DO 103 K=1,NV JK=MENT+JC(K) 103 PS=PS+XPR(K)*R(JK) IF(SNGL(PS).LT.S2) GO TO 107 IF (IMP.GT.0) WRITE(IO,1019) IFLAG=1 RETURN 107 NV=NV+1 NC=NC+1 JC(NV)=J0 NITER=1 GO TO 300 110 IF(IFLAG.LE.1) GO TO 140 C SAVE THE CORRAL OF PREVIOUS CALL DO 120 I=1,NT1 120 IC(I)=0 DO 130 K=1,NV JK=JC(K) 130 IC(JK)=1 IC(NT1)=1 C INITIALIZE WITH ONE FEASIBLE GRADIENT 140 JC(1)=1 NV=2 NC=1 JC(2)=0 DO 150 J=2,NT1 IF(A(J).GT.DEPS) GO TO 150 JC(2)=J 150 CONTINUE IF(JC(2).GT.0) GO TO 160 IF (IMP.GT.0) WRITE(IO,1001) IFLAG=2 RETURN 160 J=JC(2) RR(1)=1.D0 JJ=(J-1)*MM1+J PS=1.D0+R(JJ) IF (PS.GT.0.D0) GO TO 170 IFLAG=3 RETURN 170 RR(2)=DSQRT(PS) R(2)=A(J) DO 180 I=1,NT1 180 XPR(I)=0.D0 XPR(1)=DEPS-A(J) XPR(2)=1.D0 U1=0.D0 U2=-R(JJ) IF(IMP.GT.6) WRITE(IO,1005) J C C STOPPING CRITERION C 200 NITER=NITER+1 IF(IMP.GT.6) WRITE(IO,1011) (XPR(I),I=1,NV) IF(NITER.LE.ITMAX) GO TO 205 IF (IMP.GT.0) WRITE(IO,1016) IFLAG=4 RETURN 205 S2=-SNGL(DEPS*U1+U2) IF(S2.LE.ETA) GO TO 900 SP=GAMA*S2 C FIRST COMPUTE ALL THE TESTS, C AND TEST WITH THE CORRAL OF PREVIOUS CALL J0=0 DO 220 J=2,NT1 PS=U1*(A(J)-DEPS) DO 210 K=1,NV JJ=JC(K) IF(JJ.EQ.1) GO TO 210 J1=MAX0(J,JJ) J2=MIN0(J,JJ) JJ=(J1-1)*MM1+J2 PS=PS+XPR(K)*R(JJ) 210 CONTINUE Y(J)=PS IF(IFLAG.NE.2) GO TO 220 IF(IC(J).NE.1) GO TO 220 IF(SNGL(PS).GE.SP) GO TO 220 J0=J SP=SNGL(PS) 220 CONTINUE IF(J0.EQ.0) GO TO 240 IF(SP.GE.GAMA*S2) GO TO 240 PS1=DABS(U1*(DEPS-A(J0))) DO 230 K=1,NV J=JC(K) IF(J.EQ.J0) GO TO 240 IF(J.EQ.1) GO TO 230 J1=MAX0(J0,J) J2=MIN0(J0,J) JJ=(J1-1)*MM1+J2 PS1=PS1+XPR(K)*DABS(U1*(2.D0*DEPS-A(J))+2.D0*Y(J)-R(JJ)) 230 CONTINUE PS1=PS1*1000.D0*DZERO IF(SP.GT.S2-SNGL(PS1)) GO TO 240 IC(J0)=0 GO TO 280 C NOW THE REMAINING ONES 240 J0=0 SP=GAMA*S2 DO 260 J=2,NT1 IF(IFLAG.EQ.2.AND.IC(J).EQ.1) GO TO 260 IF(SNGL(Y(J)).GE.SP) GO TO 260 SP=SNGL(Y(J)) J0=J 260 CONTINUE IF(J0.EQ.0) GO TO 290 PS1=DABS(U1*(DEPS-A(J0))) DO 270 K=1,NV J=JC(K) IF(J.EQ.1) GO TO 270 J1=MAX0(J0,J) J2=MIN0(J0,J) JJ=(J1-1)*MM1+J2 PS1=PS1+XPR(K)*DABS(U1*(2.D0*DEPS-A(J))+2.D0*Y(J)-R(JJ)) 270 CONTINUE PS1=PS1*1000.D0*DZERO IF(SP.GT.S2-SNGL(PS1)) GO TO 290 280 NC=NC+1 NV=NV+1 JC(NV)=J0 IF(IMP.GT.6) WRITE(IO,1006) S2,J0,SP GO TO 300 C FIRST SET OF OPTIMALITY CONDITIONS SATISFIED 290 IF(U1.GE.-DBLE(FLOAT(NV))*DZERO) GO TO 900 J0=1 NV=NV+1 JC(NV)=1 IF(IMP.GT.6) WRITE(IO,1020) S2,U1 C C AUGMENTING R C 300 NV1=NV-1 DO 305 K=1,NV1 IF(JC(K).NE.J0) GO TO 305 IF (IMP.GT.0) WRITE(IO,1013) J0 IFLAG=3 RETURN 305 CONTINUE J=JC(1) J1=MAX0(J,J0) J2=MIN0(J,J0) JJ=(J1-1)*MM1+J2 R(NV)=(A(J)*A(J0)+E(J)*E(J0)+R(JJ))/RR(1) PS0=R(NV)*R(NV) IF(NV1.EQ.1) GO TO 330 DO 320 K=2,NV1 J=JC(K) J1=MAX0(J,J0) J2=MIN0(J,J0) JJ=(J1-1)*MM1+J2 PS=A(J)*A(J0)+E(J)*E(J0)+R(JJ) K1=K-1 DO 310 KK=1,K1 J1=(KK-1)*MM1+K J2=(KK-1)*MM1+NV 310 PS=PS-R(J1)*R(J2) MEK=K1*MM1+NV R(MEK)=PS/RR(K) 320 PS0=PS0+R(MEK)*R(MEK) JJ=(J0-1)*MM1+J0 PS0=A(J0)*A(J0)+E(J0)*E(J0)+R(JJ)-PS0 IF (PS0.GT.0.D0) GO TO 330 IFLAG=3 RETURN 330 RR(NV)=DSQRT(PS0) IF(NITER.LE.1) GO TO 400 INCR=1 K00=NV C C SOLVING THE CORRAL-SYSTEM C 400 K=K00 IF(K.GT.NV) GO TO 430 IF(IMP.LE.7) GO TO 410 WRITE(IO,1012) RR(1) IF(NV.EQ.1) GO TO 410 DO 404 L=2,NV K1=L-1 WRITE(IO,1012) (R((KK-1)*MM1+L),KK=1,K1),RR(L) 404 CONTINUE 410 J=JC(K) PS1=A(J) PS2=E(J) IF(K.EQ.1) GO TO 420 K1=K-1 DO 415 KK=1,K1 JJ=(KK-1)*MM1+K PS0=R(JJ) PS1=PS1-PS0*W1(KK) 415 PS2=PS2-PS0*W2(KK) 420 PS0=RR(K) W1(K)=PS1/PS0 W2(K)=PS2/PS0 K=K+1 IF(K.LE.NV) GO TO 410 C TWO-TWO SYSTEM 430 K=1 IF(INCR.EQ.1) K=NV 440 W1S=W1S+W1(K)*W1(K) W2S=W2S+W2(K)*W2(K) W12S=W12S+W1(K)*W2(K) K=K+1 IF(K.LE.NV) GO TO 440 DET=W1S*W2S-W12S*W12S PS2=W2S*DEPS-W12S PS1=W1S-W12S*DEPS 450 V1=PS2/DET V2=PS1/DET 460 U1=DEPS-V1 U2=1.D0-V2 IF(NV.EQ.NC+1) U1=0.D0 C BACKWARD Y(NV)=(V1*W1(NV)+V2*W2(NV))/RR(NV) IF(NV.EQ.1) GO TO 500 DO 480 L=2,NV K=NV-L+1 K1=K+1 PS=V1*W1(K)+V2*W2(K) MEK=(K-1)*MM1 DO 470 KK=K1,NV MEJ=MEK+KK 470 PS=PS-R(MEJ)*Y(KK) 480 Y(K)=PS/RR(K) C C TEST FOR POSITIVITY C 500 DMU=-DBLE(ZERO*EPS) DO 530 K=1,NV IF(JC(K).EQ.1) GO TO 520 IF(Y(K).LE.DZERO) GO TO 550 GO TO 530 520 IF(Y(K).LE.DMU) GO TO 550 530 CONTINUE DO 540 K=1,NV 540 XPR(K)=Y(K) GO TO 200 C INTERPOLATING BETWEEN X AND Y 550 TETA=0.D0 K0=K DO 560 K=1,NV IF(Y(K).GE.0.D0) GO TO 560 PS=Y(K)/(Y(K)-XPR(K)) IF(TETA.GE.PS) GO TO 560 TETA=PS K0=K 560 CONTINUE DO 570 K=1,NV PS=TETA*XPR(K)+(1.D0-TETA)*Y(K) IF(PS.LE.DZERO) PS=0.D0 570 XPR(K)=PS IF(IMP.LE.6) GO TO 600 PS1=0.D0 PS2=0.D0 DO 580 K=1,NV DO 580 KK=1,NV J1=MAX0(JC(K),JC(KK)) J2=MIN0(JC(K),JC(KK)) JJ=(J1-1)*MM1+J2 PS1=PS1+XPR(K)*XPR(KK)*R(JJ) PS2=PS2+Y(K)*Y(KK)*R(JJ) 580 CONTINUE C C COMPRESSING THE CORRAL C 600 NV=NV-1 INCR=0 K00=K0 W1S=0.D0 W2S=0.D0 W12S=0.D0 L=JC(K0) IF(L.NE.1) NC=NC-1 IF (IMP.GT.6) WRITE(IO,1007) K0,L,Y(K0),PS1,PS2 IF(K0.GT.NV) GO TO 400 K1=K0-1 DO 620 K=K0,NV XPR(K)=XPR(K+1) IF(K0.EQ.1) GO TO 620 DO 610 KK=1,K1 MEK=(KK-1)*MM1+K 610 R(MEK)=R(MEK+1) 620 JC(K)=JC(K+1) XPR(NV+1)=0.D0 630 MEK=(K0-1)*MM1+K0+1 PS=R(MEK) PS12=RR(K0+1) PS0=DSQRT(PS*PS+PS12*PS12) PS=PS/PS0 PS12=PS12/PS0 RR(K0)=PS0 IF(K0.EQ.NV) GO TO 400 K1=K0+1 MEK01=(K0-1)*MM1 MEK=K0*MM1 MEKK=MEK-MM1 DO 640 K=K1,NV J1=MEKK+K J2=MEK+K R(J1)=PS*R(J1+1)+PS12*R(J2+1) IF(K.GT.K1) R(J2)=PS2 640 PS2=-PS12*R(J1+1)+PS*R(J2+1) R(J2+1)=PS2 K0=K0+1 GO TO 630 C C *** FINISHED *** C 900 IFLAG=0 DO 930 J=1,NTOT 930 AL(J)=0. DO 940 K=1,NV J=JC(K)-1 IF(J.NE.0) AL(J)=SNGL(XPR(K)) 940 CONTINUE U=SNGL(U1) IF (IMP.LE.5) RETURN WRITE(IO,1014) NC,NV,S2,SP,U1 WRITE(IO,1015) (JC(K),K=1,NV) WRITE(IO,1018) RETURN END SUBROUTINE EFINF1 (N,NV,JC,XPR,P,S) DIMENSION JC(NV),P(*),S(N) DOUBLE PRECISION XPR(NV) C C CETTE SUBROUTINE CALCULE S = SIGMA XPR(.)*P(.) C SACHANT QUE LES XPR ONT ETE CALCULES PAR EPRF2 C DO 920 I=1,N PS=0. DO 910 K=1,NV J=JC(K)-1 IF(J.EQ.0) GO TO 910 NIJ=(J-1)*N+I PS=PS+SNGL(XPR(K))*P(NIJ) 910 CONTINUE 920 S(I)=PS RETURN END SUBROUTINE EREMF1 (PROSCA,IFLAG,N,NTOT,NTA,MM1,P,ALFA, 1 E,A,R,IZS,RZS,DZS) EXTERNAL PROSCA DIMENSION P(*),ALFA(NTOT),IZS(*),RZS(*) DOUBLE PRECISION E(MM1),A(MM1),R(*),DZS(*),PS C C CETTE SUBROUTINE REMPLIT LES DONNEES POUR EPRF2 C (PRODUITS SCALAIRES ET 2 CONTRAINTES LINEAIRES) C C DE 1 A NTOT +1 SI IFLAG=0 C DE NTA+1 +1 A NTOT +1 SINON C C (LE +1 EST DU A L'ECART, PLACE EN PREMIER) C C P CONTIENT LES OPPOSES DES GRADIENTS A LA QUEUE LEU LEU C -G(1), -G(2),..., -G(NTOT) SOIT NTOT*N COORDONNEES C NT1=NTOT+1 NTA1=NTA+1 IF(IFLAG.GT.0) GO TO 50 C C REMPLISSAGE DES ANCIENNES DONNEES C (PRODUITS SCALAIRES, ECART ET CONTRAINTE D'EGALITE) C DO 10 J=1,NTOT JJ=(J-1)*MM1+1 10 R(JJ)=0.D0 A(1)=1.D0 E(1)=0.D0 IF (NTA1.EQ.1) GO TO 50 DO 30 J=2,NTA1 E(J)=1.D0 NJ=(J-2)*N MEJ=(J-1)*MM1 DO 30 I=2,J NI=(I-2)*N C C PRODUIT SCALAIRE DE G(I-1) AVEC G(J-1) C POUR J-1=1,NTA ET I-1=1,J-1 C CALL PROSCA (N,P(NI+1),P(NJ+1),PS,IZS,RZS,DZS) NIJ=MEJ+I C LE PRODUIT SCALAIRE CI-DESSUS VA DANS R((J-1)*MM1+I) R(NIJ)=PS 30 CONTINUE C C 50 NTA2=NTA+2 C C REMPLISSAGE DES NOUVELLES DONNEES C IF (NTA2.GT.NT1) GO TO 100 DO 70 KK=NTA2,NT1 MEKK=(KK-1)*MM1 E(KK)=1.D0 R(MEKK+1)=0.D0 NJ=(KK-2)*N DO 70 I=2,KK NI=(I-2)*N C C PRODUIT SCALAIRE DE G(KK-1) AVEC G(I-1) C POUR KK-1=NTA+1,NTOT ET I-1=1,KK-1 C CALL PROSCA (N,P(NI+1),P(NJ+1),PS,IZS,RZS,DZS) NIJ=MEKK+I C LE PRODUIT SCALAIRE CI-DESSUS VA DANS R((KK-1)*MM1+I) 70 R(NIJ)=PS C C REMPLISSAGE DE LA CONTRAINTE D'INEGALITE C (TOUT ENTIERE SAUF L'ECART) C DO 80 I=2,NT1 80 A(I)=DBLE(ALFA(I-1)) 100 RETURN END SUBROUTINE M1GC2 (SIMUL,PROSCA,N,X,F,G,DXMIN,DF1,EPSREL,IMP,IO, / MODE,NITER,NSIM,RZ,NRZ,IZS,RZS,DZS) C C LE SOUS-PROGRAMME M1GC2 (GRADIENT CONJUGUE A ENCOMBREMENT VARIABLE) C EST UNE INTERFACE ENTRE LE PROGRAMME APPELANT ET LE SOUS-PROGRAMME C M1GC2A, MINIMISEUR PROPREMENT DIT. C NRZ EST LA DIMENSION DECLAREE POUR LE TABLEAU DE TRAVAIL RZ. C RZ EST SUBDIVISE EN 4 VECTEURS DE DIMENSION N C ET UN TABLEAU DE DIMENSION MEMH. C MEMH EST LA DIMENSION ALLOUEE A LA MATRICE DE QUASI NEWTON H. C REAL ZERO, UN PARAMETER ( ZERO=0. , UN=1. ) C DECLARATION DES TABLEAUX REAL X(N), G(N), RZ(NRZ), RZS(*) DOUBLE PRECISION DZS(*) INTEGER IZS(*) C DECLARATION DES SCALAIRES REAL F, EPSREL, DXMIN, DF1 INTEGER N, NRZ, IMP, NSIM, MODE INTEGER ID, IX, IG, IAUX, IH, MEMH C EXTERNAL SIMUL, PROSCA C IF (IMP .GT. 0) WRITE(IO,1)N,NRZ,NITER,NSIM,IMP, / EPSREL,DF1,DXMIN 1 FORMAT(19H ENTREE DANS M1GC2:,6X,22HDIMENSION DU PROBLEME , /I6/2X,4HNRZ=,I7,4X,6HNITER=,I6,4X,5HNSIM=,I6,4X,4HIMP=,I3/2X, /7HEPSREL=,E8.2,4X,4HDF1=,E8.2,4X,6HDXMIN=,E8.2) IF ( N.LE.0 .OR. NITER.LE.0 .OR. NSIM.LE.0 .OR. / DXMIN.LE.ZERO .OR. DF1.LE.ZERO / .OR. EPSREL.LE.ZERO .OR. EPSREL.GT.UN ) THEN MODE=2 IF (IMP .GT. 0) WRITE(IO,3) RETURN ENDIF C C CALCULS DES POINTEURS DESTINES A SUBDIVISER LE TABLEAU RZ ID=1 IX=ID +N IG=IX +N IAUX=IG +N IH=IAUX + N C C CALCUL DU NOMBRE DE PLACES MEMOIRE AFFECTEES A H MEMH=NRZ - 4*N C IF (MEMH .LE. 0) THEN MODE=3 GOTO 100 ELSE CONTINUE ENDIF C C APPEL DU SOUS-PROGRAMME M1GC2A QUI EFFECTUE LA REELLE OPTIMISATION CALL M1GC2A(SIMUL,PROSCA,N,X,F,G,DXMIN,DF1,EPSREL,IMP,IO, / NITER,NSIM,MODE,MEMH,RZ(ID),RZ(IX),RZ(IG), / RZ(IAUX),RZ(IH),IZS,RZS,DZS) C 100 IF (IMP .GT. 0) THEN IF (MODE .EQ. 3) THEN WRITE(IO,2) ELSE IF (MODE .EQ. 6) THEN WRITE(IO,4) ELSE WRITE(IO,5)EPSREL,NITER,NSIM ENDIF ENDIF RETURN 2 FORMAT(/,38H M1GC2 RZ INSUFFISAMMENT DIMENSIONNE) 3 FORMAT(/,25H M1GC2 APPEL INCOHERENT) 4 FORMAT(/,22H M1GC2 FIN SUR DXMIN) 5 FORMAT(/,16H SORTIE DE M1GC2,7X,12HNORME DE G =,E15.9/9X, / 6HNITER=,I4,4X,5HNSIM=,I5) END SUBROUTINE M1GC2A(SIMUL,PROSCA,N,X,F,G,DX,DF1,EPS,IMP,IO, / NITER,NSIM,INFO,MEMH,D,XX,GG,TABAUX,H, / IZS,RZS,DZS) C C PARAMETRES REAL ZERO , UN , DEUX , RO PARAMETER ( ZERO=0. , UN=1. , DEUX=2. , RO=0.2 ) C DECLARATION DES TABLEAUX REAL X(N),G(N),D(N),XX(N),GG(N),TABAUX(N),RZS(*),H(*) DOUBLE PRECISION DZS(*) INTEGER IZS(*) C DECLARATION DES SCALAIRES REAL F, DX, EPS, DF1 REAL DG1, DG, ALPHA, NORMG0, AUX1, AUX2, MU, ETA, / OMEGA, NORMG, GCARRE, GGCARR, NU, SIGMA, SSCALG, USCALG /,SK, SSCAEK INTEGER N, MEMH, IMP, IO, NSIM, NITER, INFO INTEGER MEMUTI, NRZUTI, MEMSUP, M, INDIC, RETOUR, ITER, / NAPPEL, NTOTAP, NMISAJ, I, IU, IS, IETA, INU, J, KJ, K, KP1 DOUBLE PRECISION PS LOGICAL GC, ITERQN, INTFOR, REDFOR, REDEM, TERMI C EXTERNAL SIMUL, PROSCA C C ************************************************************* C PHASE I:DETERMINATION DE LA METHODE ( ET DE M LE CAS ECHEANT) C ************************************************************* C MEMUTI=N*(N+1)/2 C C MEMSUP EST AUSSI LA DIMENSION MINIMALE DE LA MATRICE H MEMSUP=2*N+2 C IF (MEMH .GE. MEMUTI) THEN GC=.FALSE. NRZUTI=MEMUTI+4*N IF (IMP .GT. 1) WRITE(IO,1) NRZUTI ELSE IF (MEMH .LT. MEMSUP) THEN INFO=3 RETURN ELSE GC=.TRUE. C M EST LE NOMBRE DE MISES A JOUR ADMISSIBLE M=MEMH / MEMSUP C MEMUTI EST ICI LE NOMBRE DE PLACES MEMOIRE UTILISEES POUR STOCKER H MEMUTI=M * MEMSUP NRZUTI=MEMUTI+4*N IF (IMP .GT. 1) WRITE(IO,2) M,NRZUTI ENDIF 1 FORMAT(40H METHODE DE QUASI-NEWTON. NRZ UTILE=,I7) 2 FORMAT(38H METHODE DU GRADIENT CONJUGUE AVEC,I3, / 14H MISES A JOUR.,11H NRZ UTILE=,I7) C C *********************************************** C PHASE II:INITIALISATIONS PROPRES A L'OPTIMISEUR C *********************************************** C C INITIALISATION DES COMPTEURS ITER=1 NTOTAP=1 C C ****************************************************************** C PHASE III:DEMARRAGE A PARTIR DE X(0) AVEC DESCENTE SUIVANT -(GRAD) C ****************************************************************** C 3000 I=0 NMISAJ=0 C C CALCUL DE LA DIRECTION DE DESCENTE DO 3100 J=1,N D(J)=-G(J) 3100 CONTINUE C CALL PROSCA(N,G,D,PS,IZS,RZS,DZS) DG1=SNGL(PS) NORMG0=SQRT(ABS(DG1)) IF (ITER .EQ. 1) THEN OMEGA=EPS * NORMG0 ENDIF C C ************************************************************ C PHASE IV:DEBUT DE L'ITERATION X(I-1) DONNE X(I) LE LONG DE D C ************************************************************ C 4000 IF (ITER .EQ. NITER) THEN INFO=4 GOTO 99999 ENDIF ITER=ITER + 1 I=I+1 C C DETERMINATION DU TYPE DE PAS IF (GC) THEN ITERQN=(I .LE. M) .AND. (2 .LE. I) ENDIF C C ******************************* C PHASE V:INITIALISATION DE ALPHA C ******************************* C IF (ITER .EQ. 2) THEN ALPHA=DEUX * DF1 /(-DG1) ELSE IF (GC) THEN IF (I.EQ.1) THEN ALPHA=UN / NORMG0 ELSE IF (ITERQN) THEN ALPHA=UN ELSE ALPHA=ALPHA * DG / DG1 ENDIF ENDIF ELSE ALPHA=UN ENDIF C C *************************** C PHASE VI:RECHERCHE LINEAIRE C *************************** C DG=DG1 INTFOR=( GC .AND. (.NOT.ITERQN)).OR. ((.NOT.GC) .AND.(I.EQ.1)) DO 6000 J=1,N XX(J)=X(J) GG(J)=G(J) 6000 CONTINUE CALL M1GC2B(N,SIMUL,PROSCA,XX,F,DG,ALPHA,D,X,G,IMP,IO,RETOUR / ,NTOTAP,NSIM,INTFOR,DX,EPS,IZS,RZS,DZS) C IF (IMP .GT. 3) WRITE(IO,6003) IF ((RETOUR .EQ. 4).OR.((RETOUR .EQ. 1).AND.(I .EQ. 1))) THEN INFO=6 RETURN ELSE IF (RETOUR .EQ. 1) THEN IF (IMP .GT. 1) WRITE(IO,6002) ITER,NTOTAP GOTO 3000 ELSE C CALCUL DE (G,G) IF((I .GT. 1) .AND. GC) GGCARR=GCARRE CALL PROSCA(N,G,G,PS,IZS,RZS,DZS) GCARRE=SNGL(PS) NORMG=SQRT(GCARRE) IF (IMP .GT. 2) WRITE(IO,6001)ITER,NTOTAP,F,NMISAJ IF (RETOUR .EQ. 2) THEN INFO=0 GOTO 99999 ELSE IF (RETOUR .EQ. 3)THEN INFO=5 GOTO 99999 ENDIF ENDIF 6001 FORMAT(4X,6H M1GC2,3X,I4,6H ITERS,3X,I4,7H SIMULS,3X,2HF=,E15.8, 1 I6,4H MAJ) 6002 FORMAT(4X,6H M1GC2,3X,I4,6H ITERS,3X,I4,7H SIMULS, / 33H NECESSITE D'UN REDEMARRAGE TOTAL) 6003 FORMAT() C C ****************************************************** C PHASE VII:TEST D'ARRET PAR OBTENTION DE LA CONVERGENCE C ****************************************************** C TERMI=NORMG .LT. OMEGA IF (TERMI) THEN INFO=1 GOTO 99999 ELSE CONTINUE ENDIF C C ******************************************* C PHASE VIII:TEST X(I) POINT DE REDEMARRAGE? C ******************************************* C C DOIT ON FORCER UN REDEMARRAGE? REDFOR=GC .AND. ((I .EQ. 1) .OR. (I .EQ. M+N)) IF (REDFOR) THEN REDEM=.TRUE. ELSE IF (GC .AND. .NOT. ITERQN) THEN CALL PROSCA(N,G,GG,PS,IZS,RZS,DZS) AUX1=SNGL(PS) REDEM=ABS(AUX1) .GT. ABS(RO * GGCARR) ELSE REDEM=.FALSE. ENDIF C C ******************** C PHASE IX:MISE A JOUR C ******************** C C CALCULS DE S STOCKE DANS D ET DE Y STOCKE DANS XX DO 9000 J=1,N D(J)=ALPHA * D(J) XX(J)=G(J)-GG(J) 9000 CONTINUE IF (REDEM) THEN C CAS OU X(I) EST UN POINT DE REDEMARRAGE I=1 NMISAJ=1 C SAUVEGARDE DE S QUI EST ACTUELLEMENT DANS D C U=H*Y=Y C NU=(Y,HY)=(Y,Y) C ETA=(S,Y) C CALCUL DES INDICES INU=1 IETA=INU + 1 IU=IETA IS=IU + N C DO 9100 J=1,N H(IU +J)=XX(J) H(IS +J)=D(J) 9100 CONTINUE CALL PROSCA(N,XX,XX,PS,IZS,RZS,DZS) NU=SNGL(PS) H(INU)=NU CALL PROSCA(N,D,XX,PS,IZS,RZS,DZS) ETA=SNGL(PS) H(IETA)=ETA C H1 EST MAINTENANT DEFINIE C CALCUL DE H1*G QUE L'ON RANGE DANS XX CALL EMULB1(N,H,G,XX,TABAUX,NMISAJ,PROSCA,IZS,RZS,DZS) C ELSE IF (GC) THEN C CAS DE GC SANS REDAMARRAGE C CALCUL DE H*Y RANGE DANS GG CALL EMULB1(N,H,XX,GG,TABAUX,NMISAJ,PROSCA,IZS,RZS,DZS) C CALCULS DE NU, ETA, SSCALG, USCALG CALL PROSCA(N,XX,GG,PS,IZS,RZS,DZS) NU=SNGL(PS) CALL PROSCA(N,D,XX,PS,IZS,RZS,DZS) ETA=SNGL(PS) CALL PROSCA(N,D,G,PS,IZS,RZS,DZS) SSCALG=SNGL(PS) CALL PROSCA(N,GG,G,PS,IZS,RZS,DZS) USCALG=SNGL(PS) C CALCUL DE SIGMA ET DE MU SIGMA=(USCALG -(UN + NU / ETA)* SSCALG) / ETA MU=SSCALG /ETA C CALCUL DE H*G QUE L'ON RANGE DANS XX CALL EMULB1(N,H,G,XX,TABAUX,NMISAJ,PROSCA,IZS,RZS,DZS) C CALCUL DE LA NOUVELLE DIRECTION DE RECHERCHE: C H*G - MU * U - SIGMA * S DO 9200 J=1,N XX(J)= XX(J) - MU * GG(J) - SIGMA * D(J) 9200 CONTINUE C C CAS D'UNE ITERATION DE TYPE QUASI NEWTON IF (ITERQN) THEN NMISAJ=NMISAJ + 1 C SAUVEGARDE DES TERMES UTILES POUR STOCKER LA MATRICE MISE A JOUR INU=INU + MEMSUP IETA=INU + 1 IU=IETA IS=IU + N DO 9300 J=1,N H(IU +J)=GG(J) H(IS +J)=D(J) 9300 CONTINUE H(INU)=NU H(IETA)=ETA ENDIF C CAS DE LA METHODE QUASI NEWTON ELSE C CALCUL DE ETA=(S,Y) CALL PROSCA(N,D,XX,PS,IZS,RZS,DZS) ETA=SNGL(PS) IF (I .EQ. 1) THEN C ETAPE INITIALE CALCUL DE L'APPROXIMATION INITIALE DE L'INVERSE DE LA C MATRICE HESSIENNE C CALCUL DE NU=(Y,H0*Y)=(Y,Y) CALL PROSCA(N,XX,XX,PS,IZS,RZS,DZS) NU=SNGL(PS) C STOCKAGE DE CETTE MATRICE H=(ETA / NU) * I KJ=1 AUX1=ETA / NU DO 9500 K=1,N H(KJ)=AUX1 KJ=KJ +1 KP1=K+1 IF (N .GE. KP1) THEN DO 9400 J=KP1,N H(KJ)=ZERO KJ=KJ +1 9400 CONTINUE ENDIF GG(K)=AUX1 * XX(K) 9500 CONTINUE NU=ETA ELSE CALL EMULS1(N,H,XX,GG) CALL PROSCA(N,XX,GG,PS,IZS,RZS,DZS) NU=SNGL(PS) ENDIF C CALCUL DE LA MATRICE MISE A JOUR (UTILISATION DE LA FORMULE BFGS ) C NU, ETA ET H*Y (STOCKE DANS GG) SONT CONNUS AUX1=UN + NU / ETA KJ=1 DO 9800 K=1,N C CALCUL DU VECTEUR CONTENANT LA KEME COLONNE DE H LK=K KM1=K-1 IF (K .GE. 2) THEN DO 9610 L=1,KM1 TABAUX(L)=H(LK) LK=LK + (N-L) 9610 CONTINUE ENDIF DO 9620 L=K,N TABAUX(L)=H(LK) LK=LK+1 9620 CONTINUE C CALL PROSCA(N,XX,TABAUX,PS,IZS,RZS,DZS) AUX2=SNGL(PS) DO 9630 L=1,N TABAUX(L)=ZERO 9630 CONTINUE TABAUX(K)=UN CALL PROSCA(N,TABAUX,D,PS,IZS,RZS,DZS) SSCAEK=SNGL(PS) KJ=K-N DO 9700 J=1,K KJ=KJ+N-J+1 H(KJ)=H(KJ) - ( (AUX2 - AUX1*SSCAEK)*D(J) + SSCAEK*GG(J) )/ETA 9700 CONTINUE 9800 CONTINUE ENDIF C C ***************************************************** C PHASE X :CALCUL DE LA NOUVELLE DIRECTION DE RECHERCHE C ***************************************************** C IF (GC) THEN C XX CONTIENT -D DO 10000 J=1,N D(J)=-XX(J) 10000 CONTINUE C ELSE C CAS DE LA METHODE DE QUASI NEWTON C LA NOUVELLE DIRECTION D EGALE -(H * G) CALL EMULS1(N,H,G,D) DO 10100 J=1,N D(J)=-D(J) 10100 CONTINUE ENDIF C C TEST:LA DIRECTION DE RECHERCHE EST ELLE BIEN DE DESCENTE CALL PROSCA(N,D,G,PS,IZS,RZS,DZS) DG1=SNGL(PS) IF (DG1 .GE. ZERO) THEN INFO=7 IF (IMP .GT. 1) WRITE(IO,10101) DG1 GOTO 99999 ELSE GOTO 4000 ENDIF C C RETOUR AU PROGRAMME APPELANT 99999 NITER=ITER NSIM=NTOTAP IF (I .EQ. 0) THEN EPS=NORMG0 ELSE EPS=NORMG ENDIF 10101 FORMAT(40H M1GC2A ERREUR DANS LA HESSIENNE DG=,E9.2) END SUBROUTINE EMULS1(N,H,X,HX) C C CE SOUS-PROGRAMME EFFECTUE LE PRODUIT H * X AVEC: C N (E) DIMENSION DU PROBLEME C H (E) DIMENSION N(N+1)/2. TIANGLE SUPERIEUR, COEFFICIENTS PAR LIGNE C X (E) VECTEUR DE DIMENSION N C HX (S) DIMENSION N. RESULTAT DU PRODUIT C C C PARAMETRE REAL ZERO PARAMETER ( ZERO=0. ) C DECLARATIONS REAL H(*), X(N), HX(N), AUX1 INTEGER N, K, KM1, KJ, J C DO 3000 K=1,N C CALCUL DE LA KEME COMPOSANTE DU PRODUIT H* X AUX1=ZERO C H(KJ) EST LE COEFFICIENT (K,J) DE LA MATRICE SYMETRIQUE COMPLETE KJ=K KM1=K-1 C CONTRIBUTION DU TRIANGLE INFERIEUR IF (KM1.GE.1) THEN DO 1000 J=1,KM1 AUX1=AUX1 + H(KJ) * X(J) KJ=KJ+(N-J) 1000 CONTINUE ENDIF C CONTRIBUTION DU TRIANGLE SUPERIEUR DO 2000 J=K,N AUX1=AUX1 + H(KJ) * X(J) KJ=KJ+1 2000 CONTINUE C HX(K)=AUX1 3000 CONTINUE C RETURN END SUBROUTINE EMULB1(N,H,X,HX,TABAUX,NMISAJ,PROSCA,IZS,RZS,DZS) C C PARAMETRES REAL UN , DEUX PARAMETER ( UN=1., DEUX=2. ) C DECLARATIONS REAL H(*), X(N), HX(N), TABAUX(N), RZS(*) DOUBLE PRECISION DZS(*), PS INTEGER IZS(*) REAL USCALX, SSCALX, NU, ETA, GAMMA, MU, SIGMA INTEGER N, NMISAJ, MEMSUP, PTNU, COMPT, IU, IS, K C C CALCUL DE LA LONGUEUR D'UN BLOC MEMSUP=2*N+2 C CALCUL DE H0*X=X=X DO 1000 K=1,N HX(K)=X(K) 1000 CONTINUE C IF (NMISAJ.EQ.0) THEN RETURN ELSE PTNU=1 COMPT=1 ENDIF C 2000 IU=PTNU+1 IS=IU+N DO 3000 K=1,N TABAUX(K)=H(IU+K) 3000 CONTINUE CALL PROSCA(N,TABAUX,X,PS,IZS,RZS,DZS) USCALX=SNGL(PS) DO 4000 K=1,N TABAUX(K)=H(IS+K) 4000 CONTINUE CALL PROSCA(N,TABAUX,X,PS,IZS,RZS,DZS) SSCALX=SNGL(PS) NU=H(PTNU) ETA=H(PTNU+1) C CALCUL DU NOUVEAU TERME ET ADDITION DANS HX IF (COMPT.EQ.1) THEN GAMMA=ETA / NU DO 5000 K=1,N HX(K)=GAMMA * HX(K) 5000 CONTINUE MU=SSCALX / NU SIGMA=-(DEUX * SSCALX / ETA)+(USCALX / NU) ELSE MU=SSCALX / ETA SIGMA=-(UN + NU / ETA)* MU + USCALX / ETA ENDIF C DO 6000 K=1,N HX(K)=HX(K) - MU * H(IU+K) - SIGMA * H(IS+K) 6000 CONTINUE C COMPT=COMPT+1 IF (COMPT.LE.NMISAJ) THEN PTNU=PTNU+MEMSUP GOTO 2000 ELSE RETURN ENDIF END SUBROUTINE M1GC2B(N,SIMUL,PROSCA,XINIT,F,DG,ALPHA,D, / XFINAL,GFINAL,IMP,IO,RETOUR,NTOTAP,NSIM, / INTFOR,DX,EPS,IZS,RZS,DZS) C PARAMETRES REAL ZERO ,DEUX ,TROIS PARAMETER ( ZERO=0. ,DEUX=2. ,TROIS=3. ) REAL DIXIEM , PETIT PARAMETER ( DIXIEM=.1 , PETIT=.0001 ) REAL UNPLUS ,UNMOIN ,ENVIR1 PARAMETER ( UNPLUS=1.01 ,UNMOIN=.99 ,ENVIR1=.9 ) C DECLARATIONS DES TABLEAUX REAL XINIT(N),D(N),XFINAL(N),GFINAL(N),RZS(*) DOUBLE PRECISION DZS(*) INTEGER IZS(*) C DECLARATIONS DES SCALAIRES REAL F, FINIT, DG, ALPHA, EPS, DX, AP, DP, FP, / AUX1, AUX2, PAS, AT, DAL, BSUP, DELTA INTEGER N, IMP, IO, RETOUR, NSIM, NTOTAP, NAPPEL, INDIC, J LOGICAL INTFOR, MAXPAS, RFINIE, ACCEPT, ENCADR, DEPAS DOUBLE PRECISION PS C EXTERNAL SIMUL, PROSCA C C INITIALISATIONS DEPAS=.FALSE. BSUP=ZERO FINIT=F NAPPEL=0 AP=ZERO FP=FINIT DP=DG IF (IMP .GT. 3) WRITE(IO,1) ALPHA, DG C CALCUL DE LA LONGUEUR DU PAS CALL PROSCA(N,D,D,PS,IZS,RZS,DZS) PAS=SNGL(PS) PAS=SQRT(PAS) C TEST D'ERREUR DANS LA RECHERCHE LINEAIRE 1000 CONTINUE IF (ALPHA * PAS .LE. DX) THEN IF (IMP .GT. 3) WRITE(IO,1001) RETOUR=1 RETURN ELSE IF (NTOTAP .EQ. NSIM) THEN RETOUR=3 RETURN ELSE CONTINUE ENDIF C CALCUL DU NOUVEAU POINT SUSCEPTIBLE D'ETRE XFINAL DO 2000 J=1,N XFINAL(J)=XINIT(J) + ALPHA * D(J) 2000 CONTINUE C CALCULS DE F ET G EN CE POINT INDIC=4 CALL SIMUL(INDIC,N,XFINAL,F,GFINAL,IZS,RZS,DZS) NAPPEL=NAPPEL + 1 NTOTAP=NTOTAP + 1 IF (INDIC .LT. 0) THEN DEPAS=.TRUE. IF (IMP . GT. 3) THEN WRITE(IO,2001) ALPHA,INDIC ENDIF DELTA=ALPHA - AP IF (DELTA .LE. DX) THEN RETOUR=4 RETURN ELSE BSUP=ALPHA ALPHA=DELTA * DIXIEM + AP GOTO 1000 ENDIF ENDIF C CALCUL DE LA DERIVEE SUIVANT D AU POINT XFINAL CALL PROSCA(N,D,GFINAL,PS,IZS,RZS,DZS) DAL=SNGL(PS) C IF (IMP .GT. 3) THEN AUX2=F - FINIT WRITE(IO,2002) ALPHA, AUX2, DAL ENDIF IF (INDIC .EQ. 0) THEN RETOUR=2 RETURN ENDIF MAXPAS=(F .GT. FINIT) .AND. (DAL .LT. ZERO) IF (MAXPAS) THEN ALPHA=ALPHA / TROIS AP=ZERO FP=FINIT DP=DG RFINIE=.FALSE. C ELSE C TEST:LE NOUVEAU POINT EST IL ACCEPTABLE AUX1=FINIT + PETIT * ALPHA * DG AUX2=ABS(DAL/DG) ACCEPT=(F .LE. AUX1) .AND. (AUX2 .LE. ENVIR1) IF (ACCEPT) THEN C DOIT ON FAIRE UNE INTERPOLATION RFINIE=(NAPPEL .GT. 1) .OR. (.NOT. INTFOR) .OR. (AUX2 .LE. EPS) ELSE RFINIE=.FALSE. ENDIF C IF (.NOT. RFINIE) THEN C LA RECHERCHE LINEAIRE N'EST PAS FINIE. INTERPOLATION AUX1=DP + DAL- TROIS*(FP-F)/(AP-ALPHA) AUX2=AUX1 * AUX1 - DP * DAL IF (AUX2 .LE. ZERO) THEN AUX2=ZERO ELSE AUX2=SQRT(AUX2) ENDIF IF (DAL-DP+ DEUX * AUX2 .EQ. ZERO) THEN RETOUR=4 RETURN ENDIF AT=ALPHA - (ALPHA-AP)*(DAL+AUX2-AUX1)/(DAL-DP+ DEUX * AUX2) C TEST:LE MINIMUM A T-IL ETE ENCADRE ENCADR=(DAL/DP) .LE. ZERO IF (ENCADR) THEN C LE MINIMUM A ETE ENCADRE IF (ABS(ALPHA - AP) .LE. DX) THEN RETOUR=4 RETURN ENDIF AUX1=UNPLUS * AMIN1(ALPHA,AP) AUX2=UNMOIN * AMAX1(ALPHA,AP) IF ((AT .LT. AUX1) .OR. (AT .GT. AUX2)) AT=(ALPHA + AP)/DEUX ELSE C LE MINIMUM N'A PAS ETE ENCADRE AUX1=UNMOIN * AMIN1(AP,ALPHA) IF ((DAL .LE. ZERO) .OR. (AT .LE. ZERO) .OR. (AT .GE. AUX1)) THEN AUX1=UNPLUS * AMAX1(AP,ALPHA) IF ((DAL .GT. ZERO) .OR. (AT .LE. AUX1)) THEN IF (DAL .LE. ZERO) THEN AT=DEUX * MAX(AP,ALPHA) ELSE AT=AMIN1(AP,ALPHA) / DEUX ENDIF ENDIF ENDIF ENDIF IF ( (DEPAS) .AND. (AT .GE. BSUP)) THEN DELTA=BSUP - ALPHA IF (DELTA .LE. DX) THEN RETOUR=4 RETURN ELSE AT=ALPHA + DELTA * DIXIEM ENDIF ENDIF AP=ALPHA FP=F DP=DAL ALPHA=AT ENDIF ENDIF IF (RFINIE) THEN RETOUR=0 RETURN ELSE GOTO 1000 ENDIF 1 FORMAT(7H M1GC2B,6X,5H PAS,E10.3,5H DG=,E9.2) 1001 FORMAT(21H M1GC2B FIN SUR DX) 2001 FORMAT(7H M1GC2B,20X,E10.3,8H INDIC=,I3) 2002 FORMAT(7H M1GC2B,20X,E10.3,2E11.3) END