C SIMPLE SOURCE CB215821 17/07/21 21:15:30 9513 SUBROUTINE SIMPLE ********************************************************************** * * * IMPLEMENTATION D'UNE METHODE DU SIMPLEX DANS CASTEM 2000 * * * * SUBROUTINE (ESOPE) UTILISEE : TACVEC * * * * SUBROUTINE (FORTRAN) UTILISEES : SIMPLX,SIMP1,SIMP2,SIMP3 * * * * REFERENCE : NUMERICAL RECIPES pp 312-325 * * * * PROGRAMEUR : P.PEGON 31/8/92 * ********************************************************************** IMPLICIT INTEGER(I-N) IMPLICIT REAL*8(A-H,O-Z) -INC PPARAM -INC CCOPTIO -INC CCREEL -INC SMLREEL -INC SMLENTI -INC SMTABLE C SEGMENT WORK REAL*8 A(MP,NP) INTEGER IZROV(NN),IPOSV(MM) INTEGER L1(MMAX),L2(MMAX),L3(MMAX) ENDSEGMENT C LOGICAL LOGIN,LOGRE CHARACTER*8 TYPOBJ CHARACTER*1 CHARIN,CHARRE * * VX VINEG= SIMP VF CINEG CEGAL; * * lecture des tables * CALL LIRTAB('VECTEUR',KVF,1,IRETOU) IF(IERR.NE.0) RETURN CALL LIROBJ('TABLE',KNEG,1,IRETOU) IF(IERR.NE.0) RETURN CALL LIROBJ('TABLE',KEGA,1,IRETOU) IF(IERR.NE.0) RETURN * * et du reel * CALL LIRREE(XXCONV,0,IRETOU) IF(IRETOU.EQ.0)XXCONV=1.D-10 * * On cherche a determiner la dimension du tableau * * nb de ligne MM * JG=1 SEGINI MLENTI LECT(1)=KVF * * nb de contrainte d'inegalite * MTABLE=KNEG SEGACT MTABLE IF(MLOTAB.EQ.0) GOTO 11 DO 10 I=1,MLOTAB TYPOBJ=' ' CALL ACCTAB(MTABLE,'ENTIER ',I,XVALIN,CHARIN,LOGIN,IOBIN, $ TYPOBJ,IVALRE,XVALRE,CHARRE,LOGRE,IOBRE) SEGACT MTABLE IF(TYPOBJ.NE.' ') GOTO 10 M1M2=I-1 GOTO 12 10 CONTINUE 11 M1M2=MLOTAB 12 CONTINUE IF(M1M2.EQ.0)THEN CALL ERREUR(-284) ELSE JG=JG+M1M2 SEGADJ MLENTI DO 15 I=1,M1M2 TYPOBJ='TABLE ' CALL ACCTAB(MTABLE,'ENTIER ',I,XVALIN,CHARIN,LOGIN,IOBIN, $ TYPOBJ,IVALRE,XVALRE,CHARRE,LOGRE,IOBRE) IF(IERR.NE.0)GOTO 9999 SEGACT MTABLE LECT(I+1)=IOBRE 15 CONTINUE ENDIF * * nb de contrainte d'egalite * MTABLE=KEGA SEGACT MTABLE IF(MLOTAB.EQ.0) GOTO 21 DO 20 I=1,MLOTAB TYPOBJ=' ' CALL ACCTAB(MTABLE,'ENTIER ',I,XVALIN,CHARIN,LOGIN,IOBIN, $ TYPOBJ,IVALRE,XVALRE,CHARRE,LOGRE,IOBRE) SEGACT MTABLE IF(TYPOBJ.NE.' ') GOTO 20 M3=I-1 GOTO 22 20 CONTINUE 21 M3=MLOTAB 22 CONTINUE IF(M3.EQ.0)THEN CALL ERREUR(-285) ELSE JG=JG+M3 SEGADJ MLENTI DO 25 I=1,M3 TYPOBJ='TABLE ' CALL ACCTAB(MTABLE,'ENTIER ',I,XVALIN,CHARIN,LOGIN,IOBIN, $ TYPOBJ,IVALRE,XVALRE,CHARRE,LOGRE,IOBRE) IF(IERR.NE.0)GOTO 9998 SEGACT MTABLE LECT(I+M1M2+1)=IOBRE 25 CONTINUE ENDIF MM=JG-1 * IF (MM.EQ.0)THEN CALL ERREUR(619) GOTO 9997 ENDIF * * nb de colonne NN * JG=0 DO 30 I=1,MM+1 MTABLE=LECT(I) SEGACT MTABLE JJG=0 IXXXXX=0 CALL TACVEC (MTABLE,IXXXXX,JJG) IF(JJG.EQ.0)THEN CALL ERREUR(620) GOTO 9996 ENDIF JG=MAX(JG,JJG) 30 CONTINUE NN=JG * * allocation du vecteur de lecture * JG=NN+1 SEGINI MLREEL * * allocation de la zone de travail * NP=NN+1 MP=MM+2 MMAX=MP+NP SEGINI WORK * * lecture du tableau * * 1ere ligne = la fonction * JLECT=LECT(1) NNP1 = NN+1 CALL TACVEC (JLECT,MLREEL,NNP1) DO 35 J=1,NN+1 A(1,J)=PROG(J) 35 CONTINUE * * M1M2 lignes = les inegalites * M1=0 M2=0 IF(M1M2.GT.0)THEN JG=M1M2 SEGINI MLENT1 DO 45 I=1,M1M2 JLECT=LECT(I+1) CALL TACVEC (JLECT,MLREEL,NNP1) IF(PROG(1).GE.0.D0)THEN M1=M1+1 MLENT1.LECT(I)=M1 A(M1+1,1)=PROG(1) DO 40 J=2,NN+1 A(M1+1,J)=-PROG(J) 40 CONTINUE ELSE M2=M2+1 MLENT1.LECT(I)=M1M2-M2+1 A(M1M2-M2+1+1,1)=-PROG(1) DO 42 J=2,NN+1 A(M1M2-M2+1+1,J)=PROG(J) 42 CONTINUE ENDIF 45 CONTINUE * * on "inverse MLENT1" * JG=M1M2 SEGINI MLENT2 DO 46 I=1,M1M2 MLENT2.LECT(MLENT1.LECT(I))=I 46 CONTINUE SEGSUP MLENT1 MLENT1=MLENT2 ENDIF * * M3 lignes = les egalites * IF(M3.GT.0)THEN DO 50 I=1,M3 JLECT=LECT(I+M1M2+1) CALL TACVEC (JLECT,MLREEL,NNP1) IF(PROG(1).GE.0.D0)THEN IISIGN=-1 ELSE IISIGN=+1 ENDIF A(I+M1M2+1,1)=-IISIGN*PROG(1) DO 47 J=2,NN+1 A(I+M1M2+1,J)=IISIGN*PROG(J) 47 CONTINUE 50 CONTINUE ENDIF * * on libere les tables d'entrees * DO 60 I=1,MM+1 MTABLE=LECT(I) SEGDES MTABLE 60 CONTINUE MTABLE=KEGA SEGDES MTABLE MTABLE=KNEG SEGDES MTABLE * * on veut voir les entres * IF(IIMPI.EQ.1789)THEN WRITE(6,*)'MM,NN,MP,NP',MM,NN,MP,NP WRITE(6,*)'M1,M2,M3',M1,M2,M3 WRITE(6,*)'la fonction' WRITE(6,*)(A(1,J),J=1,NN+1) IF(M1M2.GT.0)THEN WRITE(6,*)'les inegalites' DO 200 I=1,M1M2 WRITE(6,*)'numero ',I,' index',MLENT1.LECT(I) WRITE(6,*)(A(I+1,J),J=1,NN+1) 200 CONTINUE ENDIF IF(M3.GT.0)THEN WRITE(6,*)'les egalites' DO 201 I=1,M3 WRITE(6,*)'numero ',I WRITE(6,*)(A(I+1+M1M2,J),J=1,NN+1) 201 CONTINUE ENDIF ENDIF * * on fait le simplex * CALL SIMPLX(A,MM,NN,MP,NP,M1,M2,M3,ICASE,IZROV,IPOSV, > XXCONV,L1,L2,L3,MMAX) * * on sort en cas d'erreur * IF (ICASE.NE.0)THEN IF(ICASE.EQ.-1)THEN CALL ERREUR(621) ELSE CALL ERREUR(622) ENDIF M=0 SEGINI MTABLE SEGDES MTABLE KX=MTABLE KSEG=MTABLE GOTO 100 ENDIF * * on elabore la solution * * nb de variable principale et secondaire non nulle * NPR=0 NSEG=0 DO 65 I=1,MM IF(IPOSV(I).LE.NN)THEN NPR=NPR+1 ELSEIF(IPOSV(I).LE.NN+M1M2)THEN NSEG=NSEG+1 ENDIF 65 CONTINUE * * table des variables principale et de la valeur de la fonction * M=NPR+1 SEGINI MTABLE KX=MTABLE MLOTAB=M DO 70 I=1,MLOTAB MTABTI(I)='ENTIER ' MTABTV(I)='FLOTTANT' 70 CONTINUE NPR=1 MTABII(NPR)=0 RMTABV(NPR)=A(1,1) DO 71 I=1,MM IF(IPOSV(I).LE.NN)THEN NPR=NPR+1 MTABII(NPR)=IPOSV(I) RMTABV(NPR)=A(I+1,1) ENDIF 71 CONTINUE SEGDES MTABLE * * table des variables secondaires * M=NSEG SEGINI MTABLE KSEG=MTABLE MLOTAB=M IF(NSEG.GT.0)THEN DO 80 I=1,MLOTAB MTABTI(I)='ENTIER ' MTABTV(I)='FLOTTANT' 80 CONTINUE NSEG=0 DO 81 I=1,MM C ERR? IF(IPOSV(I).GT.NN)THEN IF(IPOSV(I).GT.NN.AND.IPOSV(I).LE.NN+M1M2)THEN NSEG=NSEG+1 MTABII(NSEG)=MLENT1.LECT(IPOSV(I)-NN) RMTABV(NSEG)=A(I+1,1) ENDIF 81 CONTINUE ENDIF SEGDES MTABLE * * bye bye * 100 CALL ECROBJ('TABLE ',KSEG) CALL ECROBJ('TABLE ',KX) CALL ECRENT(ICASE) IF(M1M2.GT.0)SEGSUP MLENT1 SEGSUP WORK SEGSUP MLREEL,MLENTI RETURN * * erreur * 9996 CONTINUE DO 9896 II=1,I MTABLE=LECT(II) SEGDES MTABLE 9896 CONTINUE 9997 CONTINUE 9998 CONTINUE MTABLE=KEGA SEGDES MTABLE 9999 CONTINUE MTABLE=KNEG SEGDES MTABLE SEGSUP MLENTI RETURN * END