Télécharger la plaquette Cast3M


fr en

Télécharger simple.eso

Retour à la liste

Numérotation des lignes : 

simple

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

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