home *** CD-ROM | disk | FTP | other *** search
/ Dream 52 / Amiga_Dream_52.iso / Amiga / Applications / Mathematiques / ASimplex.lha / ASimplex / source / phase1.c < prev    next >
Text File  |  1989-09-16  |  13KB  |  430 lines
  1. /*****************************************************************************
  2.  * Modul                : phase1.c                                           *
  3.  * Zweck                : Phase I des Simplexalgorithmus                     *
  4.  * Autor                : Stefan F÷rster                                     *
  5.  *                                                                           *
  6.  * Datum      | Version | Bemerkung                                          *
  7.  * -----------|---------|--------------------------------------------------- *
  8.  * 01.03.1989 | 0.0     |                                                    *
  9.  * 02.03.1989 | 0.1     | TryPivot(), Calc1(), Calc2()                       *
  10.  * 03.03.1989 | 0.2     | Bug in TryPivot(): sum = eta[k] statt sum = eta[i] *
  11.  *            |         | Abfrage auf Invertierbarkeit von AB                *
  12.  * 04.03.1989 | 1.0     | Bug in Calc2(): zusΣtzl.: b[i] = b[j]              *
  13.  * 06.03.1989 | 1.1     | Bug in TryPivot(): else ptr1+=mm; bei AB1-Update   *
  14.  * 07.03.1989 |         | Ein Aufruf von Mult() gespart                      *
  15.  * 10.03.1989 | 1.2     | ─nd. in Calc1(): *c0 = dummy + c0start;            *
  16.  * 14.03.1989 | 1.3     | BUG in PhaseI(): c0 neuberechnen, falls CLEAR_CUT  *
  17.  * 06.05.1989 | 1.3a    | BUG in TryPivot(): Pivot auch, wenn Nq[i]<0 !!     *
  18.  *            |         | ─nd. in Calc1(): b2q neuberechnen                  *
  19.  * 20.05.1989 | 1.4     | ─nd. in PhaseI(): bsum globale Variable            *
  20.  * 06.08.1989 | 1.5     | ▄berflⁿssige Parameter gestrichen                  *
  21.  *****************************************************************************/
  22.  
  23.  
  24. #define PIVOT     (USHORT)1
  25. #define NO_PIVOT  (USHORT)0
  26.  
  27. IMPORT VOID   Mult();
  28. IMPORT USHORT Invers();
  29. IMPORT DOUBLE M();
  30. IMPORT VOID   SetM();
  31. IMPORT USHORT PhaseII();
  32. IMPORT VOID   CopyMemQuick();
  33.  
  34. GLOBAL DOUBLE INFINITE;
  35. GLOBAL BOOL   minimize;
  36. GLOBAL DOUBLE *lower, *upper;      /*         T  */
  37. GLOBAL DOUBLE bsum;                /* bsum = 1 b */
  38. GLOBAL DOUBLE *A, *AB1, *b, *b2q;
  39. GLOBAL SHORT  *B, *Nq;
  40.  
  41.  
  42.  
  43. /*****************************************************************************
  44.  * USHORT PhaseI()                                                           *
  45.  * - EMPTY, falls das Polyeder leer ist                                      *
  46.  * - CLEAR_CUT, falls das Polyeder einelementig ist                          *
  47.  * - OPTIMAL, falls eine zulΣssige Basis gefunden wurde                      *
  48.  * - UNBOUNDED, falls (unzulΣssigerweise !) das Hilfsprogramm unbeschr. ist  *
  49.  * - NOT_INVERTABLE, falls AB (unzulΣssigerweise!) nicht invertierbar ist    *
  50.  *                                                                           *
  51.  * Input:     m, n, A, b, c, c0start, upper  (b>=0)                          *
  52.  *                                                                           *
  53.  * Output:    EMPTY          : -                                             *
  54.  *            CLEAR_CUT      : x                                             *
  55.  *            OPTIMAL        : m, A, b, B, Nq, AB1, b2q                      *
  56.  *            UNBOUNDED      : -                                             *
  57.  *            NOT_INVERTABLE : -                                             *
  58.  *****************************************************************************/
  59.  
  60.  
  61.  
  62. USHORT PhaseI(m,n,c,c2,c0,c0start,flags,x,cq,pi,dq,Nminus,help,iter)
  63.  
  64. SHORT   *m;         /* Zeiger auf m           */
  65. SHORT   *n;         /* Zeiger auf n           */
  66. /* *B                  (m)-Vektor             */
  67. /* *Nq                 (n)-Vektor: (n-m)+m    */
  68. /* *b                  (m)-Vektor             */
  69. /* *b2q                (m)-Vektor             */
  70. DOUBLE  *c;         /* (n)-Vektor             */
  71. DOUBLE  *c2;        /* (n+m)-Vektor           */
  72. DOUBLE  *c0;        /* Zeiger auf c0          */
  73. DOUBLE  c0start;    /* Korr.wert f. Zielfktsw.*/
  74. USHORT  flags;      /* VERBOSE                */
  75. /* *upper              (n+m)-Vektor           */
  76. DOUBLE  *x;         /* (n+m)-Vektor           */
  77. DOUBLE  *cq;        /* (n)-Vektor: (n-m)+m    */
  78. DOUBLE  *pi;        /* (m)-Vektor             */
  79. DOUBLE  *dq;        /* (m)-Vektor             */
  80. SHORT   *Nminus;    /* (n)-Vektor: (n-m)+m    */
  81. DOUBLE  *help;      /* (m)-Vektor             */
  82. ULONG   *iter;      /* Anzahl Iterationen     */
  83.  
  84. {
  85.   DOUBLE  sum, *ptr1;
  86.   SHORT   mm = *m, nn = *n, i, j, pivots, no_of_bad;
  87.   USHORT  result;
  88.   VOID    PrepareData(), Calc1();
  89.   USHORT  TryPivot(), Calc2();
  90.  
  91.  
  92.   PrepareData(mm,nn,&bsum,c2,c0);
  93.  
  94.  
  95.   /* Das Hilfsprogramm hat eine optimale L÷sung (theoretisches Ergebnis) */
  96.  
  97.   result = PhaseII(mm,mm+nn,c2,c0,0.0,PHASE_I | flags,x,cq,pi,dq,Nminus,
  98.                    help,iter);
  99.  
  100.   if(result == UNBOUNDED) {
  101.     puts("\n\n!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!");
  102.     puts("!! Hilfsprogramm von Phase I unbeschrΣnkt !!");
  103.     puts("!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!\n\n");
  104.     return(UNBOUNDED);
  105.   }
  106.  
  107.   /*                        T      T      */
  108.   /* P(A,b) == EMPTY, wenn 1 Ax < 1 b ist */
  109.  
  110.   if(bsum-(*c0) > EPS_NULL) return(EMPTY);
  111.  
  112.   /* Versuch, kⁿnstliche Variablen aus der Basis zu werfen */
  113.  
  114.   do {
  115.     pivots = no_of_bad = 0;
  116.     for(i=1; i<=mm; ++i) {
  117.       if(B[i-1]>nn) {         /* help als eta-Vektor */
  118.         if(TryPivot(i,mm,nn,help,dq) == PIVOT) ++pivots;
  119.         else                                   ++no_of_bad;
  120.       }
  121.     }
  122.   } while(pivots>0);
  123.  
  124.  
  125.   if(no_of_bad == 0) {  /* keine kⁿnstlichen Variablen in der Basis (Yeah!) */
  126.     if(nn == mm) {      /* trivialer Fall */
  127.       Mult(b,help,mm,_FALSE);
  128.       for(i=0; i<mm; ++i) x[B[i]-1] = help[i]; /* umsortieren */
  129.       sum = c0start;                           /* c0 berechnen */
  130.       ptr1 = help;
  131.       for(i=0; i<mm; ++i) sum += c[B[i]-1]*(*ptr1++);
  132.       for(i=0; i<nn; ++i) {
  133.         if(Nq[i]<0) {
  134.           j = ABS(Nq[i])-1;
  135.           sum += c[j]*upper[j];
  136.         }
  137.       }
  138.       *c0 = sum;
  139.       return(CLEAR_CUT);
  140.     }
  141.     else {              /* nn > mm */
  142.       Calc1(mm,nn,c,c0,c0start,Nminus,help);
  143.       return(OPTIMAL);
  144.     }
  145.   }
  146.  
  147.   else {                /* kⁿnstl. Variablen in der Basis */
  148.     if(Calc2(m,nn,help) == NOT_INVERTABLE) return(NOT_INVERTABLE);
  149.     mm = *m;
  150.     if(nn == mm) {      /* trivialer Fall */
  151.       Mult(b,help,mm,_FALSE);
  152.       for(i=0; i<mm; ++i) x[B[i]-1] = help[i]; /* umsortieren */
  153.       sum = c0start;                           /* c0 berechnen */
  154.       ptr1 = help;
  155.       for(i=0; i<mm; ++i) sum += c[B[i]-1]*(*ptr1++);
  156.       for(i=0; i<nn; ++i) {
  157.         if(Nq[i]<0) {
  158.           j = ABS(Nq[i])-1;
  159.           sum += c[j]*upper[j];
  160.         }
  161.       }
  162.       *c0 = sum;
  163.       return(CLEAR_CUT);
  164.     }
  165.     else {
  166.       Calc1(mm,nn,c,c0,c0start,Nminus,help);
  167.       return(OPTIMAL);
  168.     }
  169.   }
  170. }
  171.  
  172.  
  173.  
  174. /*****************************************************************************
  175.  * VOID PrepareData()                                                        *
  176.  * bereitet die Daten fⁿr den Aufruf von PhaseII() mit dem Hilfsprogramm vor *
  177.  *****************************************************************************/
  178.  
  179. VOID    PrepareData(mm,nn,bsum,c2,c0)
  180.  
  181. SHORT   mm, nn;
  182. DOUBLE  *bsum, *c2, *c0;
  183.  
  184. {
  185.   REGISTER DOUBLE sum, *ptr1, *ptr2;
  186.   REGISTER SHORT  j, i, *ptr;
  187.  
  188.  
  189.   *c0 = 0.0;
  190.  
  191.   /*                     T          */
  192.   /* Neue Zielfunktion: 1 A  ->  c2 */
  193.  
  194.   ptr1 = c2;
  195.   for(i=0; i<nn; ++i) {
  196.     ptr2 = A+i;
  197.     for(sum=0.0,j=0; j<mm; ++j) {
  198.       sum += *ptr2;
  199.       ptr2 += nn;
  200.     }
  201.     *ptr1++ = sum;
  202.   }
  203.   for(i=0; i<mm; ++i) *ptr1++ = 0.0;
  204.  
  205.   /*         T  */
  206.   /* bsum = 1 b */
  207.  
  208.   ptr1 = b;
  209.   for(sum=0.0,i=0; i<mm; ++i) sum += *ptr1++;
  210.   *bsum = sum;
  211.  
  212.   /* Fⁿr das Hilfsprogramm gilt b[] == b2q[] */
  213.  
  214.   ptr1 = b;
  215.   ptr2 = b2q;
  216.   for(i=0; i<mm; ++i) *ptr2++ = *ptr1++;
  217.  
  218.   /* AB1 ist die Einheitsmatrix (nur Diagonale setzen, wegen MEMF_CLEAR) */
  219.  
  220.   ptr1 = AB1;
  221.   for(i=0; i<mm; ++i) {
  222.     *ptr1++ = 1.0;
  223.     ptr1 += mm;
  224.   }
  225.  
  226.   /* Basis B = { n+1, ... , n+m } */
  227.  
  228.   ptr = B;
  229.   for(i=0,j=nn; i<mm; ++i) *ptr++ = ++j;
  230.  
  231.   /* Nichtbasis Nq = { 1, ... , n } */
  232.  
  233.   ptr = Nq;
  234.   for(i=1; i<=nn; ++i) *ptr++ = i;
  235.  
  236. }
  237.  
  238.  
  239.  
  240. /*****************************************************************************
  241.  * USHORT TryPivot() -> PIVOT/NO_PIVOT                                       *
  242.  * Versuch, die kⁿnstliche Variable B[pos-1] aus der Basis zu entfernen      *
  243.  *****************************************************************************/
  244.  
  245. USHORT  TryPivot(pos,mm,nn,eta,dq)
  246.  
  247. SHORT   pos, mm, nn;
  248. DOUBLE  *eta, *dq;
  249.  
  250. {
  251.   REGISTER SHORT  j, k, i;
  252.   REGISTER DOUBLE *ptr1, *ptr2, sum;
  253.   SHORT           nandm = nn+mm, column;
  254.   LONG            offset = (LONG)(--pos)*(LONG)mm; /* pos: 0,...,m-1 */
  255.  
  256.   for(i=0; i<nandm; ++i) {
  257.  
  258.     if((column=ABS(Nq[i])-1) < nn) {  /* m÷gl. Kandidat gefunden */
  259.  
  260.       ptr1 = AB1+offset;      /* dq[pos] Pivotelement berechnen */
  261.       ptr2 = A+column;        /* column : 0,...,n-1 */
  262.       for(sum=0.0,j=0; j<mm; ++j,ptr2+=nn) sum += (*ptr1++)*(*ptr2);
  263.  
  264.       if(fabs(sum) > EPS_NULL) {   /* Pivotelement != 0.0 => Pivotschritt */
  265.  
  266.         dq[pos] = sum;            /* dq[] jetzt ganz berechnen */
  267.         ptr1 = AB1;
  268.         for(k=0; k<mm; ++k) {
  269.           if(k != pos) {
  270.             ptr2 = A+column;
  271.             for(sum=0.0,j=0; j<mm; ++j,ptr2+=nn) sum += (*ptr1++)*(*ptr2);
  272.             dq[k] = sum;
  273.           }
  274.           else ptr1 += mm;
  275.         }
  276.  
  277.         ptr1 = eta;                 /* eta-Vektor berechnen */
  278.         ptr2 = dq;
  279.         sum = 1.0/dq[pos];
  280.         for(j=0; j<mm; ++j,++ptr2) *ptr1++ = j == pos ? sum : -(*ptr2)*sum;
  281.  
  282.         ptr1 = AB1;                 /* AB1 update */
  283.         for(k=0; k<mm; ++k) {
  284.           ptr2 = AB1+offset;
  285.           if(k==pos) ptr1 += mm; /* pos.te Zeile ⁿberspringen */
  286.           else {
  287.             if((sum = eta[k]) != 0.0) {
  288.               for(j=0; j<mm; ++j) *ptr1++ += sum*(*ptr2++);
  289.             }
  290.             else ptr1 += mm;
  291.           }
  292.         }
  293.         ptr2 = AB1+offset;
  294.         sum = 1.0/dq[pos];
  295.         for(j=0; j<mm; ++j) *ptr2++ *= sum;  /* pos.te Spalte * 1/dq[pos] */
  296.  
  297.  
  298.         k = B[pos];
  299.         B[pos] = ABS(Nq[i]);
  300.         Nq[i] = k;
  301.  
  302.         return(PIVOT);
  303.  
  304.       }    /* if(Pivotelement > 0.0) */
  305.     }      /* if(m÷glicher Kandidat) */
  306.   }        /* for(i= ...             */
  307.  
  308.   return(NO_PIVOT);
  309. }
  310.  
  311.  
  312.  
  313. /*****************************************************************************
  314.  * VOID Calc1()                                                              *
  315.  * Entfernt die kⁿnstlichen Variablen aus Nq und berechnet b2q und c0 neu    *
  316.  *****************************************************************************/
  317.  
  318. VOID    Calc1(mm,nn,c,c0,c0start,Nminus,help)
  319.  
  320. SHORT   mm, nn;
  321. DOUBLE  *c, *c0, c0start;
  322. SHORT   *Nminus;
  323. DOUBLE  *help;
  324.  
  325. {
  326.   REGISTER DOUBLE dummy, *ptr1, *ptr3;
  327.   REGISTER SHORT  j, i, *ptr2;
  328.   SHORT           nm = nn-mm, anzneg = 0;
  329.  
  330.   /* kⁿnstliche Variable aus Nq entfernen und gleichzeitig Nminus erstellen */
  331.  
  332.   for(ptr2=Nq,i=0; i<nn && *ptr2<=nn; ++i,++ptr2) {
  333.     if(*ptr2 < 0) Nminus[anzneg++] = ABS(*ptr2);
  334.   }
  335.  
  336.   for(j=i+1; j<nn; ++j) {
  337.     if(Nq[j]<=nn) {
  338.       Nq[i++] = Nq[j];
  339.       if(Nq[j] < 0) Nminus[anzneg++] = ABS(Nq[j]);
  340.     }
  341.   }
  342.  
  343.   /* b2q neuberechnen */  
  344.  
  345.   ptr1 = help;
  346.   ptr3 = b;
  347.   for(j=0; j<mm; ++j) *ptr1++ = *ptr3++;
  348.  
  349.   for(j=0; j<anzneg; ++j) {
  350.     ptr1 = help;
  351.     ptr3 = A+(Nminus[j]-1);
  352.     dummy = upper[Nminus[j]-1];
  353.     for(i=0; i<mm; ++i) {
  354.       *ptr1++ -= (*ptr3)*dummy;
  355.       ptr3 += nn;
  356.     }
  357.   }
  358.  
  359.   Mult(help,b2q,mm,_FALSE);
  360.  
  361.  
  362.   /* fehlt noch c0 fⁿr die gefundene zulΣssige Basis */
  363.   /*        T       T              */
  364.   /* c0 = cB xB + cN_uN_ + c0start */
  365.  
  366.   ptr1 = b2q;
  367.   for(dummy=0.0,i=0; i<mm; ++i) dummy += c[B[i]-1]*(*ptr1++);
  368.   ptr2 = Nminus;
  369.   for(i=0; i<anzneg; ++i) {
  370.     j = (*ptr2++)-1;
  371.     dummy += c[j]*upper[j];
  372.   }
  373.   *c0 = dummy + c0start;
  374.  
  375. }
  376.  
  377.  
  378.  
  379. /*****************************************************************************
  380.  * USHORT Calc2() ->INVERTABLE/NOT_INVERTABLE                                *
  381.  * je nachdem, ob AB invertierbar ist (sollte es) oder nicht                 *
  382.  * Streichen der ⁿberflⁿssigen Zeilen von A, Entfernen der kⁿnstlichen       *
  383.  * Variablen aus B, Update von *m und Neuberechnung von AB1                  *
  384.  *****************************************************************************/
  385.  
  386. USHORT  Calc2(m,nn,help)
  387.  
  388. SHORT   *m, nn;
  389. DOUBLE  *help;
  390.  
  391. {
  392.   REGISTER SHORT  j, i, mm = *m, *ptr;
  393.   DOUBLE          *dest, *src;
  394.   SHORT           column;
  395.   LONG            length = (LONG)nn*(LONG)sizeof(DOUBLE);
  396.  
  397.   for(ptr=B,i=0; i<mm && (*ptr++)<=nn; ++i);
  398.  
  399.   if(i<mm) {
  400.  
  401.     dest = A+(LONG)i*(LONG)nn;
  402.     src = dest + nn;
  403.  
  404.     for(j=i+1; j<mm; ++j,src+=nn) {
  405.       if(B[j]<=nn) {
  406.         CopyMemQuick(src,dest,length);
  407.         dest += nn;
  408.         b2q[i] = b2q[j];    /* Auch Elemente aus b2q */
  409.         b[i] = b[j];        /* und b streichen       */
  410.         B[i++] = B[j];
  411.       }
  412.     }
  413.  
  414.     *m = mm = i;
  415.  
  416.     /* AB1 neu berechnen */
  417.  
  418.     for(ptr=B,j=0; j<mm; ++j) {
  419.       src = A+((*ptr++)-1);        /* B[j].te Spalte von A */
  420.       dest = AB1+j;
  421.       for(i=0; i<mm; ++i,dest+=mm,src+=nn) *dest = *src;
  422.     }
  423.     /* help als (SHORT *)-Permutationsvektor !! */
  424.     if(Invers(mm,(SHORT *)help) == NOT_INVERTABLE) return(NOT_INVERTABLE);
  425.   }
  426.  
  427.   return(INVERTABLE);
  428.  
  429. }
  430.