1 : C SIMPLE SOURCE CHAT 05/01/13 03:18:48 5004 2 : SUBROUTINE SIMPLE 3 : ********************************************************************** 4 : * * 5 : * IMPLEMENTATION D'UNE METHODE DU SIMPLEX DANS CASTEM 2000 * 6 : * * 7 : * SUBROUTINE (ESOPE) UTILISEE : TACVEC * 8 : * * 9 : * SUBROUTINE (FORTRAN) UTILISEES : SIMPLX,SIMP1,SIMP2,SIMP3 * 10 : * * 11 : * REFERENCE : NUMERICAL RECIPES pp 312-325 * 12 : * * 13 : * PROGRAMEUR : P.PEGON 31/8/92 * 14 : ********************************************************************** 15 : IMPLICIT INTEGER(I-N) 16 : IMPLICIT REAL*8(A-H,O-Z) 17 : -INC CCOPTIO 18 : -INC CCREEL 19 : -INC SMLREEL 20 : -INC SMLENTI 21 : -INC SMTABLE 22 : C 23 : SEGMENT WORK 24 : REAL*8 A(MP,NP) 25 : INTEGER IZROV(NN),IPOSV(MM) 26 : INTEGER L1(MMAX),L2(MMAX),L3(MMAX) 27 : ENDSEGMENT 28 : C 29 : LOGICAL LOGIN,LOGRE 30 : CHARACTER*8 TYPOBJ 31 : CHARACTER*1 CHARIN,CHARRE 32 : * 33 : * VX VINEG= SIMP VF CINEG CEGAL; 34 : * 35 : * lecture des tables 36 : * 36 : CALL LIRTAB('VECTEUR',KVF,1,IRETOU) 37 : IF(IERR.NE.0) RETURN 37 : CALL LIROBJ('TABLE',KNEG,1,IRETOU) 38 : IF(IERR.NE.0) RETURN 38 : CALL LIROBJ('TABLE',KEGA,1,IRETOU) 39 : IF(IERR.NE.0) RETURN 40 : * 41 : * et du reel 42 : * 42 : CALL LIRREE(XXCONV,0,IRETOU) 43 : IF(IRETOU.EQ.0)XXCONV=1.D-10 44 : * 45 : * On cherche a determiner la dimension du tableau 46 : * 47 : * nb de ligne MM 48 : * 49 : JG=1 50 : SEGINI MLENTI 51 : LECT(1)=KVF 52 : * 53 : * nb de contrainte d'inegalite 54 : * 55 : MTABLE=KNEG 56 : SEGACT MTABLE 57 : IF(MLOTAB.EQ.0) GOTO 11 58 : DO 10 I=1,MLOTAB 59 : TYPOBJ=' ' 59 : CALL ACCTAB(MTABLE,'ENTIER ',I,XVALIN,CHARIN,LOGIN,IOBIN, 60 : $ TYPOBJ,IVALRE,XVALRE,CHARRE,LOGRE,IOBRE) 61 : SEGACT MTABLE 62 : IF(TYPOBJ.NE.' ') GOTO 10 63 : M1M2=I-1 64 : GOTO 12 65 : 10 CONTINUE 66 : 11 M1M2=MLOTAB 67 : 12 CONTINUE 68 : IF(M1M2.EQ.0)THEN 68 : CALL ERREUR(-284) 69 : ELSE 70 : JG=JG+M1M2 71 : SEGADJ MLENTI 72 : DO 15 I=1,M1M2 73 : TYPOBJ='TABLE ' 73 : CALL ACCTAB(MTABLE,'ENTIER ',I,XVALIN,CHARIN,LOGIN,IOBIN, 74 : $ TYPOBJ,IVALRE,XVALRE,CHARRE,LOGRE,IOBRE) 75 : IF(IERR.NE.0)GOTO 9999 76 : SEGACT MTABLE 77 : LECT(I+1)=IOBRE 78 : 15 CONTINUE 79 : ENDIF 80 : * 81 : * nb de contrainte d'egalite 82 : * 83 : MTABLE=KEGA 84 : SEGACT MTABLE 85 : IF(MLOTAB.EQ.0) GOTO 21 86 : DO 20 I=1,MLOTAB 87 : TYPOBJ=' ' 87 : CALL ACCTAB(MTABLE,'ENTIER ',I,XVALIN,CHARIN,LOGIN,IOBIN, 88 : $ TYPOBJ,IVALRE,XVALRE,CHARRE,LOGRE,IOBRE) 89 : SEGACT MTABLE 90 : IF(TYPOBJ.NE.' ') GOTO 20 91 : M3=I-1 92 : GOTO 22 93 : 20 CONTINUE 94 : 21 M3=MLOTAB 95 : 22 CONTINUE 96 : IF(M3.EQ.0)THEN 96 : CALL ERREUR(-285) 97 : ELSE 98 : JG=JG+M3 99 : SEGADJ MLENTI 100 : DO 25 I=1,M3 101 : TYPOBJ='TABLE ' 101 : CALL ACCTAB(MTABLE,'ENTIER ',I,XVALIN,CHARIN,LOGIN,IOBIN, 102 : $ TYPOBJ,IVALRE,XVALRE,CHARRE,LOGRE,IOBRE) 103 : IF(IERR.NE.0)GOTO 9998 104 : SEGACT MTABLE 105 : LECT(I+M1M2+1)=IOBRE 106 : 25 CONTINUE 107 : ENDIF 108 : MM=JG-1 109 : * 110 : IF (MM.EQ.0)THEN 110 : CALL ERREUR(619) 111 : GOTO 9997 112 : ENDIF 113 : * 114 : * nb de colonne NN 115 : * 116 : JG=0 117 : DO 30 I=1,MM+1 118 : MTABLE=LECT(I) 119 : SEGACT MTABLE 120 : JJG=0 121 : IXXXXX=0 121 : CALL TACVEC (MTABLE,IXXXXX,JJG) 122 : IF(JJG.EQ.0)THEN 122 : CALL ERREUR(620) 123 : GOTO 9996 124 : ENDIF 125 : JG=MAX(JG,JJG) 126 : 30 CONTINUE 127 : NN=JG 128 : * 129 : * allocation du vecteur de lecture 130 : * 131 : JG=NN+1 132 : SEGINI MLREEL 133 : * 134 : * allocation de la zone de travail 135 : * 136 : NP=NN+1 137 : MP=MM+2 138 : MMAX=MP+NP 139 : SEGINI WORK 140 : * 141 : * lecture du tableau 142 : * 143 : * 1ere ligne = la fonction 144 : * 145 : JLECT=LECT(1) 145 : CALL TACVEC (JLECT,MLREEL,NN+1) 146 : DO 35 J=1,NN+1 147 : A(1,J)=PROG(J) 148 : 35 CONTINUE 149 : * 150 : * M1M2 lignes = les inegalites 151 : * 152 : M1=0 153 : M2=0 154 : IF(M1M2.GT.0)THEN 155 : JG=M1M2 156 : SEGINI MLENT1 157 : DO 45 I=1,M1M2 158 : JLECT=LECT(I+1) 158 : CALL TACVEC (JLECT,MLREEL,NN+1) 159 : IF(PROG(1).GE.0.D0)THEN 160 : M1=M1+1 161 : MLENT1.LECT(I)=M1 162 : A(M1+1,1)=PROG(1) 163 : DO 40 J=2,NN+1 164 : A(M1+1,J)=-PROG(J) 165 : 40 CONTINUE 166 : ELSE 167 : M2=M2+1 168 : MLENT1.LECT(I)=M1M2-M2+1 169 : A(M1M2-M2+1+1,1)=-PROG(1) 170 : DO 42 J=2,NN+1 171 : A(M1M2-M2+1+1,J)=PROG(J) 172 : 42 CONTINUE 173 : ENDIF 174 : 45 CONTINUE 175 : * 176 : * on "inverse MLENT1" 177 : * 178 : JG=M1M2 179 : SEGINI MLENT2 180 : DO 46 I=1,M1M2 181 : MLENT2.LECT(MLENT1.LECT(I))=I 182 : 46 CONTINUE 183 : SEGSUP MLENT1 184 : MLENT1=MLENT2 185 : ENDIF 186 : * 187 : * M3 lignes = les egalites 188 : * 189 : IF(M3.GT.0)THEN 190 : DO 50 I=1,M3 191 : JLECT=LECT(I+M1M2+1) 191 : CALL TACVEC (JLECT,MLREEL,NN+1) 192 : IF(PROG(1).GE.0.D0)THEN 193 : IISIGN=-1 194 : ELSE 195 : IISIGN=+1 196 : ENDIF 197 : A(I+M1M2+1,1)=-IISIGN*PROG(1) 198 : DO 47 J=2,NN+1 199 : A(I+M1M2+1,J)=IISIGN*PROG(J) 200 : 47 CONTINUE 201 : 50 CONTINUE 202 : ENDIF 203 : * 204 : * on libere les tables d'entrees 205 : * 206 : DO 60 I=1,MM+1 207 : MTABLE=LECT(I) 208 : SEGDES MTABLE 209 : 60 CONTINUE 210 : MTABLE=KEGA 211 : SEGDES MTABLE 212 : MTABLE=KNEG 213 : SEGDES MTABLE 214 : * 215 : * on veut voir les entres 216 : * 217 : IF(IIMPI.EQ.1789)THEN 218 : WRITE(6,*)'MM,NN,MP,NP',MM,NN,MP,NP 219 : WRITE(6,*)'M1,M2,M3',M1,M2,M3 220 : WRITE(6,*)'la fonction' 221 : WRITE(6,*)(A(1,J),J=1,NN+1) 222 : IF(M1M2.GT.0)THEN 223 : WRITE(6,*)'les inegalites' 224 : DO 200 I=1,M1M2 225 : WRITE(6,*)'numero ',I,' index',MLENT1.LECT(I) 226 : WRITE(6,*)(A(I+1,J),J=1,NN+1) 227 : 200 CONTINUE 228 : ENDIF 229 : IF(M3.GT.0)THEN 230 : WRITE(6,*)'les egalites' 231 : DO 201 I=1,M3 232 : WRITE(6,*)'numero ',I 233 : WRITE(6,*)(A(I+1+M1M2,J),J=1,NN+1) 234 : 201 CONTINUE 235 : ENDIF 236 : ENDIF 237 : * 238 : * on fait le simplex 239 : * 239 : CALL SIMPLX(A,MM,NN,MP,NP,M1,M2,M3,ICASE,IZROV,IPOSV, 240 : > XXCONV,L1,L2,L3,MMAX) 241 : * 242 : * on sort en cas d'erreur 243 : * 244 : IF (ICASE.NE.0)THEN 245 : IF(ICASE.EQ.-1)THEN 245 : CALL ERREUR(621) 246 : ELSE 246 : CALL ERREUR(622) 247 : ENDIF 248 : M=0 249 : SEGINI MTABLE 250 : SEGDES MTABLE 251 : KX=MTABLE 252 : KSEG=MTABLE 253 : GOTO 100 254 : ENDIF 255 : * 256 : * on elabore la solution 257 : * 258 : * nb de variable principale et secondaire non nulle 259 : * 260 : NPR=0 261 : NSEG=0 262 : DO 65 I=1,MM 263 : IF(IPOSV(I).LE.NN)THEN 264 : NPR=NPR+1 265 : ELSEIF(IPOSV(I).LE.NN+M1M2)THEN 266 : NSEG=NSEG+1 267 : ENDIF 268 : 65 CONTINUE 269 : * 270 : * table des variables principale et de la valeur de la fonction 271 : * 272 : M=NPR+1 273 : SEGINI MTABLE 274 : KX=MTABLE 275 : MLOTAB=M 276 : DO 70 I=1,MLOTAB 277 : MTABTI(I)='ENTIER ' 278 : MTABTV(I)='FLOTTANT' 279 : 70 CONTINUE 280 : NPR=1 281 : MTABII(NPR)=0 282 : RMTABV(NPR)=A(1,1) 283 : DO 71 I=1,MM 284 : IF(IPOSV(I).LE.NN)THEN 285 : NPR=NPR+1 286 : MTABII(NPR)=IPOSV(I) 287 : RMTABV(NPR)=A(I+1,1) 288 : ENDIF 289 : 71 CONTINUE 290 : SEGDES MTABLE 291 : * 292 : * table des variables secondaires 293 : * 294 : M=NSEG 295 : SEGINI MTABLE 296 : KSEG=MTABLE 297 : MLOTAB=M 298 : IF(NSEG.GT.0)THEN 299 : DO 80 I=1,MLOTAB 300 : MTABTI(I)='ENTIER ' 301 : MTABTV(I)='FLOTTANT' 302 : 80 CONTINUE 303 : NSEG=0 304 : DO 81 I=1,MM 305 : C ERR? IF(IPOSV(I).GT.NN)THEN 306 : IF(IPOSV(I).GT.NN.AND.IPOSV(I).LE.NN+M1M2)THEN 307 : NSEG=NSEG+1 308 : MTABII(NSEG)=MLENT1.LECT(IPOSV(I)-NN) 309 : RMTABV(NSEG)=A(I+1,1) 310 : ENDIF 311 : 81 CONTINUE 312 : ENDIF 313 : SEGDES MTABLE 314 : * 315 : * bye bye 316 : * 316 : 100 CALL ECROBJ('TABLE ',KSEG) 316 : CALL ECROBJ('TABLE ',KX) 316 : CALL ECRENT(ICASE) 317 : IF(M1M2.GT.0)SEGSUP MLENT1 318 : SEGSUP WORK 319 : SEGSUP MLREEL,MLENTI 320 : RETURN 321 : * 322 : * erreur 323 : * 324 : 9996 CONTINUE 325 : DO 9896 II=1,I 326 : MTABLE=LECT(II) 327 : SEGDES MTABLE 328 : 9896 CONTINUE 329 : 9997 CONTINUE 330 : 9998 CONTINUE 331 : MTABLE=KEGA 332 : SEGDES MTABLE 333 : 9999 CONTINUE 334 : MTABLE=KNEG 335 : SEGDES MTABLE 336 : SEGSUP MLENTI 337 : RETURN 338 : * 339 : END 340 : 341 : 342 : 343 :
© Cast3M 2003 - All rights reserved.
Disclaimer