PIPS
r1.c
Go to the documentation of this file.
1 /* ========================================================================= */
2 /* SIMPLEXE for integer variables */
3 /* ALL-INTEGER METHODS */
4 /* Jean Claude SOGNO */
5 /* Projet CHLOE -- INRIA ROCQUENCOURT */
6 /* Juin 1994 */
7 /* ========================================================================= */
8 
9 /* ========================================================================= */
10 /* Duong NGUYEN QUE */
11 /* Adaption to abstract computation: janusvalue */
12 /* CRI-ENSMP */
13 /* ========================================================================= */
14 #ifdef HAVE_CONFIG_H
15  #include "config.h"
16 #endif
17 
18 #include <stdio.h>
19 #include "rproblem.h"
20 
21 #define FT RR->ftrace
22 #define TRACE RR->ntrace
23 #define A(i,j) RR->a[i][j]
24 #define D RR->d
25 #define E RR->e
26 #define M1 RR->m
27 #define M2 RR->m2
28 #define MB RR->mb
29 #define N1 RR->n
30 #define NB RR->nb
31 #define NP RR->np
32 #define J0 RR->j0
33 #define TESTP RR->testp
34 #define DUAL RR->meth
35 #define COPIE RR->copie
36 #define BASE RR->base
37 #define COST RR->cost
38 #define RHS RR->rhs
39 #define VRESULT RR->vresult
40 #define ITER RR->iter
41 #define ITEMAX RR->tmax
42 #define X RR->x
43 #define PI RR->pi
44 #define B RR->b
45 #define G RR->g
46 #define INF RR->inf
47 #define CONORM(A) *(RR->pconor+A)
48 #define INFCOL(A) *(RR->pinfcolonn+A)
49 #define VRNUL VREPS
50 #define VARIABLELIBRE(v) v>NP && v<=N1
51 static float absf(float vf) /* ************************************************* */
52 { if (vf<0) return(-vf) ;
53  return(vf) ;
54 }
55 static int fraction(struct rproblem *RR,float vf) /* ***************************** */
56 { int vi ;
57  float vfp ;
58  if (vf<0) vfp= -vf ; else vfp= vf ;
59  vi=vfp+RR->epss ;
60  if (vfp-vi<RR->epss) return(0) ; return(1) ;
61 }
62 static int voir4(struct rproblem *RR, int v, int i11, int j11) /* ********************** */
63 { int i,j ;
64  if (TRACE<2) return(0) ; if (TRACE<3&&v!=30) return(0) ; /*99*/
65  if (v==0) { NB=N1 ; M2=M1 ; } ;/* appel avant les initialisations */
66  fprintf(FT,"m=%2d mb=%2d nb=%2d t=%2d ",M1,MB,NB,ITER) ;
67  if (v==0) fprintf(FT,"tout debut") ;
68  if (v==1)fprintf(FT,"situation initiale, contraintes <= , = , >=");
69  if (v==2) fprintf(FT,"forme standard, contraintes <= , =") ;
70  if (v==3) fprintf(FT,"lignes ont ete normalisees") ;
71  if (v==4) fprintf(FT,"colonnes ont ete normalisees") ;
72  if (v==5) fprintf(FT,"apres norm") ;
73  if (v==6) fprintf(FT,"forme standard normalisee") ;
74  if (v==7) fprintf(FT,"egalites eliminees, forme canonique <=") ;
75 
76  if (v==8) { fprintf(FT,"methode primale\n") ; return(0) ; } ;
77  if (v==9) fprintf(FT,"phase 1 primal") ;
78  if (v==10) fprintf(FT,"variable artificielle primal introduite") ;
79  if (v==11) {fprintf(FT,"algorithme simplexe primal ph. 1\n"); return(0);};
80  if (v==12) fprintf(FT,"phase 2 primal") ;
81  if (v==13) { fprintf(FT,"methode duale\n") ; return(0) ; } ;
82  if (v==14) fprintf(FT,"phase 1 dual") ;
83  if (v==15) fprintf(FT,"variable artificielle phase 1 dual introduite");
84  if (v==16) {fprintf(FT,"algorithme simplexe dual ph. 1\n"); return(0);};
85  if (v==17) fprintf(FT,"ligne artificielle inutile deplacee") ;
86  if (v==18) fprintf(FT,"retour necessaire de dual a primal") ;
87  if (v==19) fprintf(FT,"phase 2 dual") ;
88 /*if(v==30){fprintf(FT,"pppivot\‍(%2d,%2d\‍):%6.3f",i11,j11,1/ A(i11,j11));*/
89 if (v==30) {fprintf(FT,"pivot(%2d,%2d):%6.3f",i11,j11,1/ A(i11,j11));
90  fprintf(FT," variables %2d <-> %2d",B[j11],G[i11]) ; };
91  /*fprintf(FT,"d=%6.3f a=%15.9f\\\\\n",D[i11],A(0,j11)) ;*provisoire*/
92  fprintf(FT,"\\\\\n") ; /*99*/
93 if (v==1)
94 { for ( i=0 ; i <= M1 ; i++)
95  { if (i==0)
96  { if (E[i] ==1) fprintf(FT,"min ") ;
97  else if (E[i] == -1) fprintf(FT,"max ") ;
98  else fprintf(FT,"\?\? ") ;
99  } else
100  if (i+N1 < 10) fprintf(FT," x%1d : ",i+N1) ;
101  else fprintf(FT,"x%2d : ",i+N1) ;
102  for ( j=1 ; j <= N1 ; j++)
103  { if (A(i,j) == 0) continue ;
104  if (A(i,j) > 0) fprintf(FT,"+") ;
105  if (j < 10) fprintf(FT,"%6.3fx%1d ",A(i,j),j) ;
106  else fprintf(FT,"%6.3fx%2d ",A(i,j),j) ;
107  } ;
108  if (i>0)
109  { if (E[i] ==1) fprintf(FT,"<= ") ;
110  else if (E[i] == -1) fprintf(FT,">= ") ;
111  else if (E[i] ==0) fprintf(FT,"= ") ;
112  fprintf(FT,"%7.3f",D[i]) ;
113  }
114  else
115  if (D[i]!=0)
116  { if (D[i]>0) fprintf(FT,"+") ;
117  fprintf(FT,"%7.3f",D[i]) ;
118  }
119  fprintf(FT,"\n") ;
120  } ;
121  if (NP==N1) fprintf(FT,"variables hors-base toutes >=0\n") ;
122  else if (NP==0)
123  fprintf(FT,"variables hors-base toutes de signe quelconque\n") ;
124  else fprintf
125 (FT,"%3d variables hors-base >=0 ,%3d variables libres \n",NP,N1-NP) ;
126 };
127  if (TRACE<4) return(0) ;
128  if (TRACE<5 && (v!=1 ) ) return(0) ;
129  if (TRACE<6 && (v!=1 && v!=6 ) ) return(0) ;
130  if (TRACE<7 && (v!=1 && v!=6 && v!=30) ) return(0) ;
131  fprintf(FT,"b ") ;
132  for ( j=0 ; j <= NB ; j++)
133  { fprintf(FT,"[%2d]=%2d ",j,B[j]) ;
134  } ;
135  fprintf(FT,"\n") ;
136  for ( i=0 ; i <= M2 ; i++)
137  { fprintf(FT,"g[%2d]=%2d ",i,G[i]) ;
138  for ( j=0 ; j <= NB ; j++)
139  { fprintf(FT,"%6.3f ",A(i,j)) ;
140  } ;
141  fprintf(FT,"e=%4d d=%7.3f",E[i],D[i]) ;
142  fprintf(FT,"\n") ;
143  } ;
144  return(0);
145 } /* fin voir */
146 static int voir(struct rproblem *RR, int v) /* ********************** */
147 {// int i,j ;
148  if (TRACE<3) return(0) ; voir4(RR,v, 0,0); return(0);
149 }
150 /***************************************************************/
151 /***************************************************************/
152 int realsimplex(Prproblem RR)
153 { int mee,nvt;
154 /*********************** initialisations **********************/
155  M1=RR->mcontr ; N1=RR->nvar ; NP=RR->nvpos ;
156  mee=M1; nvt=N1;
157  if (M1>RMAXLIGNES) return(VRDEB) ;
158  if (N1>RMAXCOLONNES) return(VRDEB) ;
159  if (NP>nvt) return(VRDEB) ;
160  VRESULT=0 ;
161  RR->pconor= &RR->tconor[1] ; /* pconor utilise par macro CONORM */
162  RR->pinfcolonn= &RR->tinfcolonn[1] ; /* utilise par macro INFCOL */
163  RR->eps=0.0000001 ; /* peut-etre a modifier 0.0000001 0.0000003 */
164  NB=N1 ; M2=M1 ;
165  if (TRACE>0)
166  { fprintf(FT,"%3d contraintes%3d variables%3d variables non-negatives\n",
167  M1,N1,NP) ;
168  if (DUAL)fprintf(FT," dual"); else fprintf(FT," primal");
169  if (COPIE) fprintf(FT,", doublage");
170  /*else fprintf(FT,",memes colonnes" );*/
171  /*if (TESTP==0) fprintf(FT,", tests precision pousses\n" );
172  if (TESTP==1) fprintf(FT,", tests precision moyens\n" );
173  if (TESTP==2) fprintf(FT,", tests precision sommaires\n" );
174  if (TESTP>2) fprintf(FT,",sommaire,entiers inconnus\n");*/
175  if (TESTP==0) fprintf(FT,", tests precision pousses" );
176  if (TESTP==1) fprintf(FT,", tests precision moyens" );
177  if (TESTP==2) fprintf(FT,", tests precision sommaires" );
178  if (TESTP>2) fprintf(FT,",sommaire,entiers inconnus");
179  fprintf(FT,", niveau trace: %d\\\\\n",TRACE); /*99*/
180  }
181  if (COPIE)
182  { int ii,jj,jk ;
183  if(TRACE)fprintf(FT,"systeme original contraintes\n"); voir(RR,1) ;
184  if (TRACE) fprintf(FT,"systeme modifie contraintes double:\n");
185  if ((N1=(2*nvt)-NP) >RMAXCOLONNES) return(VRDEB) ;
186  for (ii=0 ; ii<= M1 ; ii++)
187  { jk=N1 ;
188  for (jj=nvt ; jj>= 1 ; jj--)
189  { float s; s=A(ii,jj);
190  if (jj>NP) A(ii,jk--)= -s;
191  A(ii,jk--)= s;
192  };
193  } ;
194  NP=N1 ; /* variables rendues toutes >=0 */
195  }
196  /*if (simbad(RR)) return(VRESULT);*/
197  VRESULT=simbad(RR);
198  if (TRACE)
199  { fprintf(FT,"resultat simplexe apres %d iterations: ",ITER) ;
200  if (VRESULT==VRFIN) fprintf(FT,"solution finie\n") ;
201  if (VRESULT==VRVID) fprintf(FT,"polyedre vide\n") ;
202  if (VRESULT==VRINF) fprintf(FT,"solution infinie\n") ;
203  if (VRESULT==VRINS) fprintf(FT,"nombre insuffisant\n") ;
204  if (VRESULT==VRDEB) fprintf(FT,"debordement tableaux\n") ;
205  if (VRESULT==VRCAL) fprintf(FT,"appel errone\n") ;
206  if (VRESULT==VREPS) fprintf(FT,"pivot anormalement petit\n") ;
207  if (VRESULT==VRBUG) fprintf(FT,"bug de programmation\n") ;
208  fprintf(FT,"\\\\\n") ; /*99*/
209  }
210  if (VRESULT) return(VRESULT);
211 /******************** cas ou une solution finie est obtenue ******************/
212  if (TRACE>1)
213  { int jj,jk ; float xf ;
214  for (jj=0 ; jj<= N1+M1 ; jj++)
215  { if (X[jj] !=0)
216  { if (jj<= N1) fprintf(FT," " );
217  else fprintf(FT,"%2d: ",jj-N1 );
218  /*fprintf(FT,"x[%2d] =%9.3f\n",jj,X[jj] );*/
219  fprintf(FT,"x[%2d] =%13.6f\n",jj,X[jj] );
220  } ;
221  } ;
222  if (TESTP>2 && MB!=M1)
223  { fprintf(FT," sans signification pour variables libres:\n");
224  for (jj=MB+1 ; jj<= M1 ; jj++) fprintf(FT," %3d",G[jj]);
225  fprintf(FT,"\n");
226  }
227  if (COPIE && (nvt-NP)) /* if (NP!=nvp) */
228  { jk=0 ;
229  for (jj=1 ; jj<= nvt+mee ; jj++)
230  { jk++ ; xf=X[jk] ;
231  if (jj>NP && jj<=nvt)
232  { jk++ ; if (xf==0) xf= -X[jk] ;
233  };
234  X[jj] = xf ;
235  };
236  fprintf(FT,"resultats probleme initial:\n") ;
237  for (jj=0 ; jj<= nvt+mee ; jj++)
238  if (X[jj] !=0)
239  { if (jj <= nvt) fprintf(FT," " );
240  else fprintf(FT,"%2d: ",jj-nvt );
241  fprintf(FT,"x[%2d] =%9.3f\n",jj,X[jj] );
242  }
243  };
244  fprintf(FT,"\\\\\n"); /*99*/
245  }
246  VRESULT=VRFIN ;
247  if (RR->tfrac==0) return(VRESULT) ;
248  /* examen si variables fractionnaires utilisation marginale */
249  if (TESTP<3)
250  /* if (mx==M1)*/
251  { int bfract, jj ; bfract=0 ;/*FAUX*/
252  for (jj=0 ; jj<= N1 ; jj++)
253  if (fraction(RR,X[jj]))
254  { if (TRACE>0)
255  { if (!bfract) fprintf(FT," variables non entieres: ");
256  fprintf(FT," %2d",jj);
257  } ;
258  bfract=1 ;/*VRAI*/
259  } ;
260  if (TRACE>0){
261  if (bfract) fprintf(FT,"\n");
262  else fprintf(FT," certitude variables toutes entieres\n");
263  }
264  RR->rfrac=0 ; if (bfract) RR->rfrac=1 ;
265  } ;
266  /*if (TRACE) fprintf(FT,"\\end{document}\n");*/
267 return(VRESULT) ;
268 }
269 /***************************************************************/
270 /***************************************************************/
271 static void norm(struct rproblem *RR)
272 /* cette procedure normalise la fonction cout, calcule les valeurs des
273  seconds membres resultant d'une normalisation du tableau a, ou
274  effectue ces deux operations */
275 {
276  float s; int i,j;
277  if (COST || ! COST && ! RHS)
278  {
279  s=0 ;
280  for (j=1 ; j<=N1 ; j++)
281  {
282  A(0,j)= A(0,j)* CONORM(j) ;
283  if (s < absf(A(0,j))) s=absf(A(0,j)) ;
284  } ;
285  if (s==0) s=1 ;
286  for (j=1 ; j<=N1 ; j++) A(0,j)= A(0,j)/s ;
287  CONORM(0)=s ;
288  } ;
289  if (RHS || ! COST && ! RHS)
290  {
291  s=0 ;
292  for (i=1 ; i<=M1 ; i++)
293  {
294  D[i]= D[i]/ CONORM(N1+i) ;
295  if (s < absf(D[i])) s=absf(D[i]) ;
296  } ;
297  CONORM(-1)=s ;
298  } ;
299  D[0]= D[0]/ CONORM(0) ;
300 } /* fin norm() */ ;
301 /*static*/ void elimination(Prproblem RR,int j1)
302 /* cette procedure permute les colonnes j1 (qui sera eliminee) et NB */
303 {
304  if (j1 != NB)
305  {
306  int i ;
307  float s ;
308  s= INFCOL(j1) ; INFCOL(j1)= INFCOL(NB) ;
309  INFCOL(NB)=s ;
310  i=B[j1] ; B[j1]=B[NB] ; B[NB]=i ;
311  for (i=0 ; i<=M1 ; i++)
312  {
313  s= A(i,j1) ; A(i,j1)= A(i,NB) ; A(i,NB)=s ;
314  } ;
315  } ;
316  NB=NB-1 ;
317 } /* fin elimination */ ;
318 
319 static void retrait(struct rproblem *RR, int i1)
320 /* cette procedure permute les lignes i1 (qui sera eliminee) et mb */
321 {
322  if (i1 != MB)
323  {
324  int j ; float s ;
325  s=INF[i1] ; INF[i1]=INF[MB] ; INF[MB]=s ;
326  s= D[i1] ; D[i1]= D[MB] ; D[MB]=s ;
327  j=G[i1] ; G[i1]=G[MB] ; G[MB]=j ;
328  for (j=0 ; j<=N1 ; j++)
329  {
330  s= A(i1,j) ; A(i1,j)= A(MB,j) ; A(MB,j)=s ;
331  } ;
332  } ;
333  MB=MB-1 ;
334 } /* fin retrait */ ;
335 
336 static void reduction(struct rproblem *RR)
337 /* cette procedure groupe les vecteurs provenant de l'elimination des
338  egalites dans les colonnes NB+1 a n */
339 { int j;
340  for (j=N1 ; j>=1 ; j--)
341  if (B[j] > N1)
342  {
343  if (E[B[j]-N1] == 0) elimination(RR,j) ;
344  } ;
345  for (j=M1 ; j>=1 ; j--)
346  if (VARIABLELIBRE(G[j])) retrait(RR,j) ;
347 } /* fin reduction() */ ;
348 
349 static int validite (struct rproblem *RR)
350 /* cette procedure verifie, lorsque la base est donnee ou lorsqu'on
351  utilise une matrice des contraintes ayant deja pivote, que le
352  tableau g ou que les tableaux g et b fournis sont compatibles */
353 { int i,j;
354  for (i=1 ; i<=M1 ; i++)
355  {
356  if (G[i]<1 || G[i]>N1+M1)
357  return(VRESULT=VRCAL) ;
358  for (j=i+1 ; j<=M1 ; j++)
359  if (G[j]==G[i]) return(VRESULT=VRCAL) ;
360  if (COST || RHS)
361  {
362  for (j=1 ; j<=N1 ; j++)
363  if (B[j]==G[i]) return(VRESULT=VRCAL);
364  }
365  else
366  if (G[i]>N1 && G[i]!=N1+i)
367  {
368  j=G[G[i]-N1] ; G[G[i]-N1]=G[i] ; G[i]=j ;
369  i=i-1 ;
370  } ;
371  } ;
372  if (COST || RHS)
373  for (j=1 ; j<=N1 ; j++)
374  {
375  if (B[j]<1 || B[j]>N1+M1)
376  return(VRESULT=VRCAL) ;
377  for (i=j+1 ; i<=N1 ; i++)
378  if (B[i]==B[j]) return(VRESULT=VRCAL) ;
379  }
380  return(0) ;
381 } /* fin validite */ ;
382 
383 static void precision(struct rproblem *RR,int i,int nx)
384 /*
385  cette procedure evalue l'incertitude moyenne dont sont affectes
386  les termes de la ligne i du tableau a */
387 {
388  int j;
389  float s ;
390  if (i==0) s=1
391  ; else
392  s=(i>M2 || G[i]>N1) ? 1 : 0 ;
393  for (j=J0 ; j<=nx ; j++)
394  if (B[j]>N1) s=s+absf(A(i,j)) ;
395  s=s*RR->eps1 ;
396  if (s>INF[i]) INF[i]=s ;
397 } /* fin precision */ ;
398 
399 static void precisioncolonne(struct rproblem *RR, int j, int mx)
400 /*
401  cette procedure evalue l'incertitude moyenne dont sont affectes
402  les termes de la colonne j du tableau a, ou ceux de la colonne
403  des seconds membres (j=-1) */
404 {
405  int i ;
406  float s ;
407  if (j != -1)
408  {
409  s=(B[j]>N1) ? 0 : 1 ;
410  for (i=1 ; i<=mx ; i++)
411  if (G[i] <= N1)
412  {
413  if (A(i,j) != 0) s=s+absf(A(i,j)) ;
414  }
415  }
416  else
417  {
418  s= CONORM(-1) ;
419  for (i=1 ; i<=mx ; i++)
420  if (G[i] <= N1)
421  {
422  if (D[i] != 0) s=s+absf(D[i]) ;
423  } ;
424  } ;
425  s=s*RR->eps1 ;
426  if (s > INFCOL(j)) INFCOL(j)=s ;
427 } /* fin precisioncolonne */ ;
428 static int pivotage(struct rproblem *RR, int i1, int j1)
429 /* le pivotage de la matrice a s'effectue a partir de la ligne i1 et de
430  la colonne j1. Les numeros de la variable d'ecart (basique) de la
431  ligne i1, et de la variable independante (hors-base) de la colonne j1,
432  contenus respectivement dans les tableaux g et b, sont permutes. Dans
433  le cas de l'elimination d'une egalite, l'operation equivaut a chasser
434  de la base une variable d'ecart identiquement nulle */
435 {
436  float p1,p2,s ;
437  int i,j,mx,nx ;
438  if (ITER==ITEMAX) return(VRESULT=VRINS);
439  s= A(i1,j1); if (absf(s)<RR->eps2) return(VRESULT=VREPS);
440  nx=N1 ; mx=M1 ; if (TESTP>0) nx=NB ; if (TESTP>2) mx=MB ;
441  ITER=ITER+1 ; p1=1/s ;
442  j=B[j1] ; B[j1]=G[i1] ; G[i1]=j ;
443  for (i=0 ; i<=mx ; i++)
444  if (i != i1 && A(i,j1) !=0)
445  {
446  p2= A(i,j1)*p1 ; D[i]= D[i]-D[i1]*p2 ;
447  if (RR->sommaire) if (absf(D[i])<RR->epss) D[i]=0 ;
448  for (j=J0 ; j<=nx ; j++)
449  if (j==j1) A(i,j)= -p2
450  ; else
451  if (A(i1,j) != 0)
452  {
453  A(i,j)= A(i,j)-p2* A(i1,j) ;
454  if (RR->sommaire)
455  if (absf(A(i,j))<RR->epss) A(i,j)=0 ;
456  } ;
457  if (! RR->sommaire) precision(RR,i,nx) ;
458  } ;
459  D[i1]= D[i1]*p1 ;
460  for (j=J0 ; j<=nx ; j++)
461  A(i1,j)=(j==j1) ? p1 : A(i1,j)*p1 ;
462  if (RR->sommaire) goto finpivotage ;
463  p1=absf(p1) ;
464  INF[i1]=INF[i1]*p1 ;
465  /* les quantites inferieures en valeurs absolue a la tolerance
466  correspondante sont annulees */
467  if (D[i1] != 0)
468  { precisioncolonne(RR,-1,mx) ;
469  s= INFCOL(-1)* CONORM(-1) ;
470  for (i=0 ; i<=mx ; i++)
471  if (D[i] != 0)
472  {
473  if (D[i]* D[i] <= INF[i]*s) D[i]=0 ;
474  } ;
475  } ;
476  for (j=J0 ; j<=nx ; j++)
477  if (j==j1) INFCOL(j)= INFCOL(j)*p1
478  ; else
479  if (A(i1,j) != 0)
480  { precisioncolonne(RR,j,mx) ;
481  s= INFCOL(j) ;
482  for (i=0 ; i<=mx ; i++)
483  if (A(i,j) != 0)
484  {
485  if (A(i,j)* A(i,j) <= INF[i]*s) A(i,j)=0 ;
486  } ;
487  } ;
488 finpivotage:
489  voir4(RR,30,i1,j1) ;
490  return(0);
491 } /* fin pivotage */ ;
492 static int simplexe(struct rproblem *RR,int i2,int bphase2,int *pj1,int *pbool1)
493 /* Cette procedure minimise le cout defini dans la ligne i2
494  (s'il s'agit du cout artificiel, il est represente en valeur opposee)*/
495 {
496  int i,j;
497  int i1= 0,i3 =0;
498  float gamma,teta,c1,c2,zeta,tet,pivot ;
499  float alpha =0;
500  /* en phase 1, bphase2 est faux */
501 choix:
502  /*
503  Le pivot est choisi de facon a assurer la diminution la plus grande du
504  cout. Si plusieurs pivots realisent cette condition (ce qui se produit
505  essentiellement dans le cas d'un sommet multiple, la plus grande
506  diminution possible pouvant etre nulle), on prendra celui situe dans
507  la colonne dont le facteur de cout relatif est le plus negatif, et
508  dans une meme colonne, le plus grand pivot. Si au cours de
509  l'exploration d'une colonne aucun coefficient positif n'est rencontre,
510  on se renvoie a l'etiquette solution infinie */
511  *pbool1=1 ; c1=0 ; tet=0 ;
512  for (j=J0 ; j<=NB ; j++)
513  {
514  teta=(bphase2) ? A(i2,j) : -A(i2,j) ;
515  if (teta < 0)
516  {
517  *pbool1=1 ; pivot=0 ;
518  for (i=1 ; i<=MB ; i++)
519  if (A(i,j) > 0)
520  {
521  zeta= D[i]/ A(i,j) ;
522  if (*pbool1
523  || zeta<alpha
524  || zeta==alpha && pivot<A(i,j))
525  {
526  *pbool1=0;
527  alpha=zeta ; i3=i ; pivot= A(i,j) ;
528  }
529  } /* fin i */ ;
530  if (*pbool1) return(VRESULT=VRINF);
531  c2=alpha*teta ;
532  if (c2<c1 || c2==c1 && teta<tet)
533  {
534  c1=c2 ; *pj1=j ; i1=i3 ; tet=teta ;
535  }
536  }
537  } /* fin j */ ;
538  if (*pbool1) goto optimum ;
539  if (pivotage(RR, i1,*pj1)) return(VRESULT);
540  /* si la variable artificielle quitte la base, on sort de la
541  procedure simplexe */
542  if (i1==i2) goto sortie ;
543  if (bphase2 || D[i2]>0) goto choix ;
544 optimum:
545  /* le cout minimum est atteint */
546  if (bphase2) goto sortie ;
547  /* on chasse la variable artificielle de la base */
548  gamma=0 ;
549  for (j=J0 ; j<=NB ; j++)
550  if (absf(A(i2,j))>gamma)
551  {
552  gamma=absf(A(i2,j)) ; *pj1=j ;
553  } ;
554  if (pivotage(RR, i2,*pj1)) return(VRESULT);
555 sortie:
556  /* en phase 1 le booleen bool1 indique si le polyedre est vide */
557 } /* fin simplexe */
558 static int simplexedual(struct rproblem *RR, int j2, int bphase2, int *pi1, int *pbool1)
559 /* cette procedure maximise dans le systeme dual soit le cout defini
560  par la colonne des seconds membres (lorsque j2=0), soit le cout
561  defini par la colonne j2 */
562 {
563  int i,j;
564  int j1 = 0 ,j3 = 0 ;
565  float gamma,teta,c1,c2,zeta,tet,pivot ;
566  float alpha=0;
567  J0=1;
568 choix:
569  *pbool1=1 ; c1=0 ; tet=0 ;
570  for (i=1 ; i<=MB ; i++)
571  {
572  teta=(bphase2) ? D[i] : A(i,j2) ;
573  if (teta < 0)
574  {
575  *pbool1=1; pivot=0 ;
576  for (j=1 ; j<=NB ; j++)
577  if (A(i,j)<0) {
578  zeta= A(0,j)/ A(i,j) ;
579  /*if (A(0,j)<0.0001) zeta=0; provisoire*/
580  /*fprintf(FT,"zetafoir d=%6.3f a=%15.9f\\\\\n",D[i],A(0,j)) ;*/
581  if (*pbool1
582  || zeta>alpha
583  || zeta==alpha && pivot>A(i,j))
584  {
585  *pbool1=0 ;
586  alpha=zeta ; j3=j ; pivot= A(i,j) ;
587  }
588  } /* fin j */ ;
589  if (*pbool1) return(VRESULT=VRVID);
590  c2=alpha*teta ;
591  if (c2>c1 || c2==c1 && teta<tet)
592  {
593  c1=c2 ; *pi1=i ; j1=j3 ; tet=teta ;
594  } ;
595  }
596  } /* fin i */ ;
597  if (*pbool1) goto optimum ;
598  if (pivotage(RR,*pi1,j1)) return(VRESULT);
599  if (j1==j2) goto sortie ;
600  if (bphase2 || A(0,j2)>0) goto choix ;
601 optimum:
602  if (bphase2) goto sortie ;
603  /* nulle ou non, la variable artificielle est chassee de la base
604  duale */
605  gamma=0 ;
606  for (i=1 ; i<=MB ; i++)
607  if (gamma < absf(A(i,j2)))
608  {
609  gamma=absf(A(i,j2)) ; *pi1=i ;
610  } ;
611  if (pivotage(RR,*pi1,j2)) return(VRESULT);
612 sortie:
613 } /* fin simplexedual */ ;
614 static int simbad(struct rproblem *RR)
615 {
616 int i,j,i1,j1,bphase2,bool1,risqueinfini; float s,gamma,pivot;
617  RR->sommaire=(TESTP>1) ? 1 : 0 ;
618  if (NP<0)
619  { if (TRACE) fprintf(FT,"Negative value for nvpos:%d\n",NP);
620  return(VRCAL);
621  }
622  if (N1<NP)
623  { if (TRACE)fprintf(FT,"nvar (%d) < nvpos (%d)\n",N1,NP);
624  return(VRCAL);
625  }
626  if (M1<0)
627  { if (TRACE) fprintf(FT,"Negative value for mcontr:%d\n",M1);
628  return(VRCAL);
629  }
630  ITER=0 ; ITEMAX=2*M1+N1 ; RR->eps1=RR->eps*10 ;
631  RR->eps2=RR->eps*100 ;RR->epss=RR->eps*10 ;
632  J0=1 ; bphase2=1 ; NB=N1 ; M2=M1 ;
633  MB=M1 ; VRESULT=0;
634  risqueinfini=0 ;
635  if (E[0]>0) E[0]=1;
636  if (E[0]<0) E[0]= -1;
637  if (COST<0) /* apres modification directe de la matrice */
638  {
639  }
640  else
641  if (! COST && ! RHS)
642  {
643  for (j=1 ; j<=N1 ; j++) B[j]=j ;
644  voir(RR,1) ;
645  for (i=0 ; i<=M1 ; i++)
646  if (E[i]== -1)
647  {
648  /* on donne aux inequations le meme sens et on
649  se ramene eventuellement a un probleme de
650  minimisation */
651  D[i]= -D[i] ;
652  for (j=1 ; j<=N1 ; j++) A(i,j)= -A(i,j) ;
653  } ;
654  voir(RR,2) ;
655  if (! RR->sommaire)
656  {
657  for (i=1 ; i<=M1 ; i++)
658  {
659  /* normalisation des lignes */
660  s=0 ;
661  for (j=1 ; j<=N1 ; j++)
662  if (s < absf(A(i,j))) s=absf(A(i,j)) ;
663  if (s==0) s=1 ;
664  for (j=1 ; j<=N1 ; j++)
665  A(i,j)= A(i,j)/s ;
666  CONORM(N1+i)=s ;
667  } ;
668  voir(RR,3) ;
669  for (j=1 ; j<=N1 ; j++)
670  {
671  /* normalisation des colonnes */
672  s=0 ;
673  for (i=1 ; i<=M1 ; i++)
674  if (s < absf(A(i,j))) s=absf(A(i,j)) ;
675  if (s==0) s=1 ;
676  if (s<1)
677  {
678  for (i=1 ; i<=M1 ; i++)
679  A(i,j)= A(i,j)/s ;
680  CONORM(j)=1/s ;
681  }
682  else CONORM(j)=1 ;
683  } ;
684  voir(RR,4) ;
685  norm(RR) ;
686  voir(RR,5) ;
687  for (i=0 ; i<=M1+1 ; i++) INF[i]=0 ;
688  for (j= -1 ; j<=N1 ; j++) INFCOL(j)=0 ;
689  } ;
690 /*#include "base.c"*/
691  if (BASE)
692  {
693  if (validite(RR)) return(VRESULT) ;
694  for (i=1 ; i<=M1 ; i++)
695  if (G[i] <= N1)
696  {
697  s=0 ;
698  for (i1=i ; i1<=M1 ; i1++)
699  if (G[i1] <= N1)
700  {
701  if (s <= absf(A(i,G[i1])))
702  {
703  s=absf(A(i,G[i1])) ; j=i1 ;
704  } ;
705  } ;
706  j1=G[j] ; G[j]=G[i] ; G[i]=N1+i ;
707  if (s==0)
708  {
709  for (i1=i+1 ; i1<=M1 ; i1++)
710  G[i1]=N1+i1 ;
711  return(VRESULT=VRNUL) ;
712  } ;
713  M2=i ;
714  if (pivotage(RR,i,j1)) return(VRESULT);
715  } ;
716  reduction(RR) ;
717  }
718  else
719 /* reprendre ................................................ */
720  for (i=1 ; i<=M1 ; i++) G[i]=N1+i ;
721 /*#include "corh.c"*/
722  }
723  else
724  {
725  if (validite(RR)) return(VRESULT) ;
726  reduction(RR) ;
727  norm(RR) ;
728  if (E[0]== -1) D[0]= -D[0] ;
729  if (COST)
730  {
731  /* calcul de la ligne de la nouvelle fonction
732  cout a l'aide de l'inverse de la base du systeme dual.
733  pi est utilise comme tableau auxiliaire */
734  for (j=1 ; j<=N1 ; j++)
735  PI[j]=(E[0]== -1) ? -A(0,j) : A(0,j) ;
736  for (j=1 ; j<=N1 ; j++)
737  {
738  A(0,j)=(B[j]>N1) ? 0 : PI[B[j]] ;
739  for (i=1 ; i<=M1 ; i++)
740  if (G[i] <= N1)
741  A(0,j)= A(0,j)-A(i,j)*PI[G[i]];
742  } ;
743  INF[0]=0 ; precision(RR,0,N1) ;
744  s=INF[0] ;
745  for (j=1 ; j<=N1 ; j++)
746  if (A(0,j)* A(0,j) <= s* INFCOL(j)) A(0,j)=0 ;
747  } ;
748  if (RHS)
749  {
750  /* calcul de la colonne des nouveaux seconds
751  membres a l'aide de l'inverse de la base (systeme
752  primal). x est utilise comme tableau auxiliaire */
753  for (i=1 ; i<=M1 ; i++)
754  X[i]=(E[i]== -1) ? -D[i] : D[i] ;
755  for (i=0 ; i<=M1 ; i++)
756  {
757  if (i != 0)
758  D[i]=(G[i]>N1) ? X[G[i]-N1] : 0 ;
759  for (j=1 ; j<=N1 ; j++)
760  if (B[j]>N1)
761  D[i]= D[i]+A(i,j)*X[B[j]-N1] ;
762  } ;
763  INFCOL(-1)=0 ; precisioncolonne(RR,-1,M1) ;
764  s= INFCOL(-1)* CONORM(-1) ;
765  for (i=1 ; i<=M1 ; i++)
766  if (D[i]* D[i] <= s*INF[i]) D[i]=0 ;
767  }
768  else
769  for (i=1 ; i<=M1 ; i++)
770  if (G[i] <= N1)
771  D[0]= D[0]-D[i]*PI[G[i]] ;
772 /* reprendre ................................................ */
773  } ;
774  voir(RR,6) ;
775  /* elimination des egalites */
776  for (i=M1 ; i>=1 ; i--)
777  if (E[i]==0 && G[i]==N1+i)
778  {
779  pivot=0 ;
780  for (j=NB ; j>=1 ; j--)
781  if (pivot>=0)
782  {
783  s=absf(A(i,j)) ;
784  if (s>pivot)
785  {
786  pivot=s ; j1=j ;
787  /* la premiere variable libre possible est
788  choisie d'autorite,signalee par pivot*/
789  if (VARIABLELIBRE(B[j])) pivot= -1 ;
790  }
791  } ;
792  if (pivot==0)
793  {
794  /* l'equation est surabondante ou incompatible */
795  if (D[i] != 0) return(VRESULT=VRVID);
796  /* l'equation surabondante est ignoree */
797  }
798  else
799  {
800  /* l'equation est changee en inequation dont la
801  variable d'ecart provient du vecteur j1, qui
802  entre dans la base */
803  if (pivotage(RR,i,j1)) return(VRESULT);
804  elimination(RR,j1) ;
805  if (pivot<0) retrait(RR,i) ;
806  } ;
807  } ;
808  voir(RR,7) ;
809  /* le systeme auquel on s'est ramene ne comporte que des
810  inequations, dont le nombre de variables est egal a NB. dans
811  les colonnes NB+1 a n ont ete groupes les vecteurs provenant
812  des pivotages effectues au cours des eliminations des egalites
813  ces vecteurs ne sont pas detruits, car ils appartiennent a
814  l'inverse de la base, utilisee dans les tests de precision, et
815  eventuellement dans le calcul de la colonne des nouveaux
816  seconds membres */
817  /* traitement variables libres */
818  for (j=NB ; j>=1 ; j--)
819  if (VARIABLELIBRE(B[j]))
820  {
821  pivot=0 ;
822  for (i=MB ; i>=1 ; i--)
823  {
824  s=absf(A(i,j)) ;
825  if (s>pivot)
826  {
827  pivot=s ; i1=i ;
828  } ;
829  } ;
830  if (pivot==0)
831  {
832  /* la variable n'intervient pas dans le domaine
833  mais si le cout relatif n'est pas nul, il ne
834  pourra exister de solution finie */
835  if (A(0,j) != 0) risqueinfini=1 ;
836  elimination(RR,j) ;
837  }
838  else
839  {
840  /* on fait entrer dans la base la variable libre
841  la variable d'ecart resultante sera desormais
842  ignoree */
843  if (pivotage(RR,i1,j)) return(VRESULT);
844  retrait(RR,i1) ;
845  } ;
846  } ;
847  if (DUAL) goto methodeduale ;
848 methodeprimale:
849  voir(RR,8) ;
850  gamma=0 ;
851  for (i=1 ; i<=MB ; i++)
852  if (D[i]<gamma)
853  {
854  bphase2=0 ; gamma= D[i] ; i1=i ;
855  } ;
856  if (bphase2) goto phase2 ;
857  /* une premiere phase est necessaire. une (seule) variable
858  artificielle est introduite dans le systeme au moyen de la
859  colonne 0 du tableau a */
860  voir(RR,9) ;
861  J0=0 ;
862  B[0]=0 ; A(0,0)=0 ; INFCOL(0)=0 ;
863  for (i=1 ; i<=M1 ; i++)
864  A(i,0)=(D[i]<0 && i<=MB) ? -1 : 0 ;
865  voir(RR,10) ;
866  if (pivotage(RR,i1,0)) return(VRESULT);
867  voir(RR,11) ;
868  (simplexe(RR,i1,bphase2,&j1,&bool1)) ; if (VRESULT) return(VRESULT);
869  B[j1]=B[0] ; INFCOL(j1)= INFCOL(0) ;
870  for (i=0 ; i<=M1 ; i++) A(i,j1)= A(i,0) ;
871  if (bool1) return(VRESULT=VRVID);
872  J0=1 ; bphase2=1 ;
873  /* fin de la premiere phase avec obtention d'une solution
874  realisable */
875 phase2:
876  voir(RR,12) ;
877  (simplexe(RR,0,bphase2,&j1,&bool1)) ; if (VRESULT) return(VRESULT);
878  goto resultats ;
879 methodeduale:
880  voir(RR,13) ;
881  gamma=0 ;
882  for (j=1 ; j<=NB ; j++)
883  if (A(0,j)<gamma)
884  {
885  bphase2=0 ; gamma= A(0,j) ; j1=j ;
886  } ;
887  if (bphase2) goto phase2dual ;
888  /* une (seule) variable artificielle est ajoutee au systeme dual
889  au moyen de la ligne mb+1 du tableau a */
890  voir(RR,14) ;
891  M2=M1=M1+1 ;
892  MB=MB+1 ;
893  /* ligne conservee uniquement pour tests de precision */
894  if (MB!=M1)
895  {
896  G[M1]=G[MB] ; D[M1]= D[MB] ; INF[M1]=INF[MB] ;
897  for (j=1 ; j<=N1 ; j++) A(M1,j)= A(MB,j) ;
898  } ;
899  G[MB]=M1+N1 ; D[MB]=0 ; INF[MB]=0 ;
900  for (j=1 ; j<=NB ; j++)
901  A(MB,j)=(A(0,j)<0) ? 1 : 0 ;
902  for (j=NB+1 ; j<=N1 ; j++) A(MB,j)=0 ;
903  voir(RR,15) ;
904  if (pivotage(RR,MB,j1)) return(VRESULT);
905  voir(RR,16) ;
906  simplexedual(RR,j1,bphase2,&i1,&bool1); if (VRESULT) return(VRESULT);
907  if (i1 != MB)
908  {
909  D[i1]= D[MB] ; G[i1]=G[MB] ; INF[i1]=INF[MB] ;
910  for (j=1 ; j<=N1 ; j++) A(i1,j)= A(MB,j) ;
911  } ;
912  if (MB != M1)
913  {
914  D[MB]= D[M1] ; G[MB]=G[M1] ; INF[MB]=INF[M1] ;
915  for (j=1 ; j<=N1 ; j++) A(MB,j)= A(M1,j) ;
916  } ;
917  M2=M1=M1-1 ;
918  MB=MB-1 ;
919  voir(RR,17) ;
920  bphase2=1 ;
921  if (bool1)
922  {
923  /* il n'existe pas de solution realisable dans le domaine
924  dual, mais on ignore si cela correspond a un polyedre
925  vide ou a une solution infinie pour le probleme
926  primal*/
927  voir(RR,18) ;
928  ITEMAX=ITEMAX+ITER ; goto methodeprimale ;
929  } ;
930 phase2dual:
931  voir(RR,19) ;
932  simplexedual(RR,0,bphase2,&i1,&bool1); if (VRESULT) return(VRESULT);
933 resultats:
934  if (risqueinfini) return(VRESULT=VRINF);
935  for (i=1 ; i<=M1+N1 ; i++) X[i]=0 ;
936  if (RR->sommaire)
937  {
938  for (i=1 ; i<=M1 ; i++) X[G[i]]= D[i] ;
939  X[0]= D[0] ;
940  }
941  else
942  {
943  for (i=1 ; i<=M1 ; i++)
944  X[G[i]]= D[i]* CONORM(G[i]) ;
945  X[0]= D[0]* CONORM(0) ;
946  } ;
947  if (E[0]== -1) X[0]= -X[0] ;
948  for (j=1 ; j<=N1+M1 ; j++) PI[j]=0 ;
949  if (RR->sommaire)
950  for (j=1 ; j<=N1 ; j++)
951  PI[B[j]]= A(0,j) ;
952  else
953  for (j=1 ; j<=N1 ; j++)
954  PI[B[j]]= A(0,j)* CONORM(0)/ CONORM(B[j]) ;
955  return(VRESULT=VRFIN) ;
956 } /* fin simbad */ ;
RR
#define RR
Definition: genspec.h:95
VRVID
#define VRVID
Definition: iproblem.h:114
VRBUG
#define VRBUG
Definition: iproblem.h:119
VRCAL
#define VRCAL
Definition: iproblem.h:118
VRDEB
#define VRDEB
Definition: iproblem.h:117
VRINS
#define VRINS
Definition: iproblem.h:116
VRINF
#define VRINF
Definition: iproblem.h:115
VRFIN
#define VRFIN
.....................
Definition: iproblem.h:113
phase2
static int phase2(Pproblem XX, struct rproblem *RR, int rvar, int mini)
Definition: isolve.c:1614
validite
static int validite(struct rproblem *RR)
fin reduction()
Definition: r1.c:349
COPIE
#define COPIE
Definition: r1.c:35
VRESULT
#define VRESULT
Definition: r1.c:39
E
#define E
Definition: r1.c:25
realsimplex
int realsimplex(Prproblem RR)
Definition: r1.c:152
COST
#define COST
Definition: r1.c:37
J0
#define J0
Definition: r1.c:32
DUAL
#define DUAL
Definition: r1.c:34
B
#define B
Definition: r1.c:44
INF
#define INF
Definition: r1.c:46
X
#define X
Definition: r1.c:42
CONORM
#define CONORM(A)
Definition: r1.c:47
VRNUL
#define VRNUL
Definition: r1.c:49
M2
#define M2
Definition: r1.c:27
precisioncolonne
static void precisioncolonne(struct rproblem *RR, int j, int mx)
fin precision
Definition: r1.c:399
simplexedual
static int simplexedual(struct rproblem *RR, int j2, int bphase2, int *pi1, int *pbool1)
fin simplexe
Definition: r1.c:558
VARIABLELIBRE
#define VARIABLELIBRE(v)
Definition: r1.c:50
norm
static void norm(struct rproblem *RR)
cette procedure normalise la fonction cout, calcule les valeurs des seconds membres resultant d'une n...
Definition: r1.c:271
absf
static float absf(float vf)
Definition: r1.c:51
precision
static void precision(struct rproblem *RR, int i, int nx)
fin validite
Definition: r1.c:383
TESTP
#define TESTP
Definition: r1.c:33
NB
#define NB
Definition: r1.c:30
PI
#define PI
Definition: r1.c:43
elimination
void elimination(Prproblem RR, int j1)
fin norm()
Definition: r1.c:301
BASE
#define BASE
Definition: r1.c:36
voir4
static int voir4(struct rproblem *RR, int v, int i11, int j11)
Definition: r1.c:62
MB
#define MB
Definition: r1.c:28
pivotage
static int pivotage(struct rproblem *RR, int i1, int j1)
fin precisioncolonne
Definition: r1.c:428
fraction
static int fraction(struct rproblem *RR, float vf)
Definition: r1.c:55
TRACE
#define TRACE
Definition: r1.c:22
NP
#define NP
Definition: r1.c:31
retrait
static void retrait(struct rproblem *RR, int i1)
fin elimination
Definition: r1.c:319
ITER
#define ITER
Definition: r1.c:40
INFCOL
#define INFCOL(A)
Definition: r1.c:48
voir
static int voir(struct rproblem *RR, int v)
fin voir
Definition: r1.c:146
N1
#define N1
Definition: r1.c:29
ITEMAX
#define ITEMAX
Definition: r1.c:41
simplexe
static int simplexe(struct rproblem *RR, int i2, int bphase2, int *pj1, int *pbool1)
fin pivotage
Definition: r1.c:492
M1
#define M1
Definition: r1.c:26
reduction
static void reduction(struct rproblem *RR)
fin retrait
Definition: r1.c:336
A
#define A(i, j)
Definition: r1.c:23
simbad
static int simbad(struct rproblem *RR)
fin simplexedual
Definition: r1.c:614
G
#define G
Definition: r1.c:45
D
#define D
Definition: r1.c:24
FT
#define FT
=========================================================================
Definition: r1.c:21
RHS
#define RHS
Definition: r1.c:38
RMAXCOLONNES
#define RMAXCOLONNES
=========================================================================
Definition: rproblem.h:16
VREPS
#define VREPS
Definition: rproblem.h:76
RMAXLIGNES
#define RMAXLIGNES
Definition: rproblem.h:17
fprintf
int fprintf()
test sc_min : ce test s'appelle par : programme fichier1.data fichier2.data ...
pivot
void pivot(frac *X, frac A, frac B, frac C, frac D, bool ofl_ctrl)
Definition: sc_simplex_feasibility_fixprec.c:770
rproblem
define RMAXCOLONNES 800
Definition: rproblem.h:21