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