home *** CD-ROM | disk | FTP | other *** search
/ Gold Fish 1 / GoldFishApril1994_CD1.img / d1xx / d199 / asimplex / phase2.c < prev    next >
Encoding:
Text File  |  1989-03-31  |  17.5 KB  |  530 lines
  1. /*****************************************************************************
  2.  * Modul                : phase2.c                                           *
  3.  * Zweck                : Phase II des Simplexalgorithmus                    *
  4.  * Autor                : Stefan Förster                                     *
  5.  *                                                                           *
  6.  * Datum      | Version | Bemerkung                                          *
  7.  * -----------|---------|--------------------------------------------------- *
  8.  * 06.02.1989 | 0.0     |                                                    *
  9.  * 07.02.1989 | 0.1     | Anpassung an PHASE_I                               *
  10.  * 25.02.1989 | 0.2     | CHUZR(), WRETA(), 1. Testlauf erfolgreich          *
  11.  * 26.02.1989 | 0.3     | PRICE() verbessert                                 *
  12.  * 26.02.1989 | 0.4     | VERBOSE-Option, Reinvertierung von AB1             *
  13.  * 02.03.1989 | 0.5     | Bug in CHUZR(): dqi<-EPS_NULL statt dqi<EPS_NULL   *
  14.  *            |         | Bug in PRICE(): sum-=pi[qs-nm] statt sum-=pi[qs]   *
  15.  *            |         | Bug in PRICE(): Bei PHASE_I hat A phys. nm Spalten *
  16.  *            |         | Bug in FTRAN(): Bei PHASE_I hat A phys. nm Spalten *
  17.  *            |         | Bug in WRETA(): M(A,phase==PHASE_I?nm:n,i,j,phase) *
  18.  * 03.03.1989 |         | Bug in CHUZR(): Abfrage if(lambdaX==INFINITE||...) *
  19.  *            | 1.0     | Bug in FTRAN(): ptr1 += m  statt  *ptr1 += m  !!!! *
  20.  *            |         | Abfrage auf Invertierbarkeit von AB                *
  21.  * 04.03.1989 | 1.1     | Bug in WRETA(): c0-Änderung, wenn Nq[s-1]<0 ist    *
  22.  * 06.03.1989 | 1.2     | Bug in WRETA(): else ptr1+=m; bei AB1-Update       *
  23.  *            |         | Neustrukturierung von WRETA()                      *
  24.  * 10.03.1989 | 1.3     | c0start, lower zusätzlich übergeben                *
  25.  *****************************************************************************/
  26.  
  27.  
  28. IMPORT VOID   Mult();
  29. IMPORT USHORT Invers();
  30. IMPORT DOUBLE M();
  31. IMPORT VOID   SetM();
  32.  
  33. GLOBAL DOUBLE INFINITE;
  34. GLOBAL BOOL   minimize;
  35. GLOBAL FILE   *file[2];
  36.  
  37.  
  38. /*****************************************************************************
  39.  * USHORT PhaseII()                                                          *
  40.  * - OPTIMAL, falls optimale Lösung gefunden                                 *
  41.  * - UNBOUNDED, falls das Problem unbeschränkt ist                           *
  42.  * - NOT_INVERTABLE, falls (unzulässigerweise!) AB nicht invertierbar ist    *
  43.  *                                                                           *
  44.  * Input:     B, Nq,                                                         *
  45.  *            A, AB1,                                                        *
  46.  *            b, b2q,                                                        *
  47.  *            c, c0, c0start,                                                *
  48.  *            flag,                                                          *
  49.  *            upper, lower bzw. NULL (für Phase1())                          *
  50.  * Output:    x, c0, B, Nq                                                   *
  51.  *****************************************************************************/
  52.  
  53. USHORT PhaseII(m,n,B,Nq,A,AB1,b,b2q,c,c0,c0start,flags,upper,lower,x,cq,pi,dq,
  54.                Nminus,help,iter)
  55.  
  56. SHORT   m,n;          /* m = Zeilenzahl, n = Spaltenzahl des Problems        */
  57. SHORT   *B;           /* (m)-Vektor     - Basis                              */
  58. SHORT   *Nq;          /* (n-m)-Vektor   - E-Nichtbasis                       */
  59. DOUBLE  *A;           /* (m,n)-Matrix   - Matrix A                           */
  60. DOUBLE  *AB1;         /* (m,m)-Matrix   - Inverse von AB                     */
  61. DOUBLE  *b;           /* (m)-Vektor     - RHS (rechte Seite)                 */
  62. DOUBLE  *b2q;         /* (m)-Vektor     - bq-AB1*AN*uN_                      */
  63. DOUBLE  *c;           /* (n)-Vektor     - Vektor c                           */
  64. DOUBLE  *c0;          /* Zahl (Zeiger!) - Zielfunktionswert                  */
  65. DOUBLE  c0start;      /* Korrekturwert für Zielfunktionswert                 *
  66. USHORT  flags;        /* Flags          - PHASE_I/PHASE_II,VERBOSE           */
  67. DOUBLE  *upper;       /* (n)-Vektor     - Obergrenzen für x                  */
  68. DOUBLE  *lower;       /* (n)-Vektor     - Untergrenzen für x                 */
  69. DOUBLE  *x;           /* (n)-Vektor     - Lösung (falls existiert)           */
  70. DOUBLE  *cq;          /* (n-m)-Vektor   - reduzierte Kosten                  */
  71. DOUBLE  *pi;          /* (m)-Vektor     - pi = cB*AB1                        */
  72. DOUBLE  *dq;          /* (m)-Vektor     - dq = AB1*A.qs                      */
  73. SHORT   *Nminus;      /* (n-m)-Vektor   - Negativer Anteil von Nq            */
  74. DOUBLE  *help;        /* (m)-Vektor     - Hilfsvektor                        */
  75. ULONG   *iter;        /* Zeiger auf "Anzahl Iterationen"                     */
  76.  
  77. {
  78.   USHORT  method = STEEPEST_ASCENT, rowflag = LAMBDA_0;
  79.   USHORT  phase = flags & (PHASE_I | PHASE_II), verbose = flags & VERBOSE;
  80.   SHORT   r ,s ,nm = n-m, i, qs, lastpos = nm-1; /* lastpos: 0,...,(n-m)-1 */
  81.   DOUBLE  c0_old;
  82.   VOID    BTRAN(), FTRAN();
  83.   USHORT  PRICE(), CHUZR(), WRETA();
  84.   ULONG   count = 1L;
  85.  
  86.  
  87.   for( ; ; count++) {
  88.  
  89.     c0_old = *c0;
  90.  
  91.     if(rowflag!=LAMBDA_2) BTRAN(m,AB1,B,c,pi,x);
  92.  
  93.     if(PRICE(m,n,A,Nq,c,cq,pi,&s,method|phase,&lastpos)==OPTIMAL) {
  94.       for(i=0; i<nm; ++i) {
  95.         qs = ABS(Nq[i])-1;
  96.         if(lower) x[qs] = Nq[i]>0 ? lower[qs] : upper[qs]+lower[qs];
  97.         else      x[qs] = Nq[i]>0 ? 0.0 : upper[qs];
  98.       }
  99.       for(i=0; i<m; ++i) {
  100.         if(lower) x[B[i]-1] = b2q[i]+lower[qs];
  101.         else      x[B[i]-1] = b2q[i];
  102.       }
  103.       *iter = count;
  104.       return(OPTIMAL);
  105.     }
  106.  
  107.     FTRAN(m,n,AB1,A,Nq[s-1],dq,phase);
  108.  
  109.     if(CHUZR(m,dq,b2q,upper,B,Nq[s-1],s,&r,&rowflag)==UNBOUNDED) {
  110.       *iter = count;
  111.       return(UNBOUNDED);
  112.     }
  113.  
  114.     /* x als eta-Vektor mißbrauchen */
  115.     if(WRETA(m,n,B,Nq,A,AB1,b,b2q,c,cq,c0,c0start,dq,upper,r,s,x,help,
  116.        Nminus,rowflag|phase|verbose,count) == NOT_INVERTABLE)
  117.       return(NOT_INVERTABLE);
  118.  
  119.     if(verbose) {
  120.       printf("@@ ");
  121.       if(phase == PHASE_I) printf("phase I : ");
  122.       else                 printf("phase II: ");
  123.       printf("iter# %5ld : s =%4d, r =",count,s);
  124.       if(rowflag == LAMBDA_2) printf(" np , ");
  125.       else                    printf("%4d, ",r);
  126.       printf("c0 = %16.10lg\n",minimize ? -(*c0) : *c0);
  127.       if(method==STEEPEST_ASCENT) puts("   steepest-ascent-rule");
  128.       else                        puts("   cyclic smallest-index-rule");
  129.       if(file[1]) {
  130.         fprintf(file[1],"@@ ");
  131.         if(phase == PHASE_I) fprintf(file[1],"phase I : ");
  132.         else                 fprintf(file[1],"phase II: ");
  133.         fprintf(file[1],"iter# %5ld : s =%4d, r =",count,s);
  134.         if(rowflag == LAMBDA_2) fprintf(file[1]," np , ");
  135.         else                    fprintf(file[1],"%4d, ",r);
  136.         fprintf(file[1],"c0 = %16.10lg\n",minimize ? -(*c0) : *c0);
  137.         if(method==STEEPEST_ASCENT) fputs("   steepest-ascent-rule\n",file[1]);
  138.         else fputs("   cyclic smallest-index-rule\n",file[1]);
  139.       }
  140.     }
  141.  
  142.     if(fabs(c0_old)>EPS_NULL)
  143.       method = (*c0-c0_old)/fabs(c0_old) <= PERCENT ? STEEPEST_ASCENT :
  144.                                                       SMALLEST_INDEX;
  145.   }
  146. }
  147.  
  148.  
  149.  
  150. /*****************************************************************************
  151.  * VOID BTRAN()                                                              *
  152.  * Backward Transformation                                                   *
  153.  *****************************************************************************/
  154.  
  155. VOID BTRAN(m,AB1,B,c,pi,dummy)
  156.  
  157. SHORT   m;
  158. DOUBLE  *AB1;
  159. SHORT   *B;
  160. DOUBLE  *c, *pi, *dummy; /* dummy : (m)-Hilfsvektor */
  161.  
  162. {
  163.   DOUBLE *ptr=dummy;
  164.   SHORT  i;
  165.  
  166.   for(i=0; i<m; ++i) *ptr++ = c[B[i]-1];
  167.   Mult(dummy,AB1,pi,1,m,m);
  168.  
  169. }
  170.  
  171.  
  172.  
  173. /*****************************************************************************
  174.  * USHORT PRICE() -> OPTIMAL/NOT_OPTIMAL                                     *
  175.  * method & STEEPEST_ASCENT => Steilster-Anstieg-Regel                       *
  176.  * method & SMALLEST_INDEX  => (zyklische) Kleinster-Index-Regel             *
  177.  *****************************************************************************/
  178.  
  179. USHORT PRICE(m,n,A,Nq,c,cq,pi,s,flags,lastpos)
  180.  
  181. SHORT   m, n;
  182. DOUBLE  *A;
  183. SHORT   *Nq;
  184. DOUBLE  *c, *cq, *pi;
  185. SHORT   *s;
  186. USHORT  flags;       /* PHASE_I/PHASE_II, SMALLEST_INDEX/STEEPEST_ASCENT */
  187. SHORT   *lastpos;
  188.  
  189. {
  190.   REGISTER DOUBLE *ptr1, *ptr2, sum;
  191.   REGISTER SHORT  j, i;
  192.   SHORT           qs, nm = n-m, pos;
  193.   USHORT          flag = OPTIMAL;
  194.   USHORT          phase = flags & (PHASE_I | PHASE_II);
  195.   USHORT          method = flags & (SMALLEST_INDEX | STEEPEST_ASCENT);
  196.   SHORT           columns = phase==PHASE_I ? nm : n;
  197.   DOUBLE          dummy;
  198.  
  199.   if(method == STEEPEST_ASCENT) {
  200.     for(i=0; i<nm; ++i) {
  201.       ptr1 = pi;
  202.       qs = ABS(Nq[i])-1;
  203.       sum = c[qs];
  204.       if(phase == PHASE_I && qs>=nm) sum -= pi[qs-nm];
  205.       else {
  206.         ptr2 = A+qs;
  207.         for(j=0; j<m; ++j) {
  208.           sum -= (*ptr1++)*(*ptr2);
  209.           ptr2 += columns;
  210.         }
  211.       }
  212.       cq[i] = sum;
  213.     }
  214.     for(i=0; i<nm; ++i) {
  215.       if((dummy=fabs(cq[i])) > EPS_NULL && cq[i]*Nq[i] > 0.0) {
  216.         if(flag==OPTIMAL) {
  217.           flag = NOT_OPTIMAL;
  218.           *s = i+1;           /* Index s : 1,...,n-m */
  219.           sum = dummy;
  220.         }
  221.         else {
  222.           if(dummy>sum) {
  223.             *s = i+1;
  224.             sum = dummy;
  225.           }
  226.         }
  227.       }
  228.     }
  229.   }
  230.  
  231.   else {                              /* method == SMALLEST_INDEX */
  232.     pos = ((*lastpos)+1) % nm;        /* lastpos+1 modulo nm */
  233.     for(i=0; i<nm; ++i) {
  234.       ptr1 = pi;
  235.       qs = ABS(Nq[pos])-1;
  236.       sum = c[qs];
  237.       if(phase == PHASE_I && qs>=nm) sum -= pi[qs-nm];
  238.       else {
  239.         ptr2 = A+qs;
  240.         for(j=0; j<m; ++j) {
  241.           sum -= (*ptr1++)*(*ptr2);
  242.           ptr2 += columns;
  243.         }
  244.       }
  245.       if(fabs(sum)>EPS_NULL && sum*Nq[pos]>0.0) {
  246.         cq[pos] = sum;
  247.         *s = pos+1;                   /* s: Index 1,...,nm */
  248.         *lastpos = pos;
  249.         return(NOT_OPTIMAL);
  250.       }
  251.       ++pos;
  252.       pos %= nm;
  253.     }
  254.   }
  255.  
  256. return(flag);
  257. }
  258.  
  259.  
  260.  
  261. /*****************************************************************************
  262.  * VOID FTRAN()                                                              *
  263.  * Forward Transformation                                                    *
  264.  *****************************************************************************/
  265.  
  266. VOID FTRAN(m,n,AB1,A,qs,dq,phase)
  267.  
  268. SHORT   m, n;
  269. DOUBLE  *AB1, *A;
  270. SHORT   qs;
  271. DOUBLE  *dq;
  272. USHORT  phase;
  273.  
  274. {
  275.   REGISTER DOUBLE *ptr1, sum;
  276.   REGISTER SHORT  j, i;
  277.   SHORT           nm = n-m, columns = phase==PHASE_I ? nm : n;
  278.  
  279.   qs = ABS(qs)-1;
  280.   if(phase == PHASE_I && qs>=nm) {
  281.     ptr1 = AB1 + (qs-nm);
  282.     for(i=0; i<m; ++i) {
  283.       dq[i] = *ptr1;
  284.       ptr1 += m;
  285.     }
  286.   }
  287.   else {
  288.     for(i=0; i<m; ++i) {
  289.       ptr1 = A+qs;    
  290.       sum = 0.0;
  291.       for(j=0; j<m; ++j) {
  292.         sum += (*AB1++)*(*ptr1);
  293.         ptr1 += columns;
  294.       }
  295.       dq[i] = sum;
  296.     }
  297.   }
  298.  
  299. }
  300.  
  301.  
  302.  
  303. /*****************************************************************************
  304.  * USHORT CHUZR() -> UNBOUNDED/NOT_UNBOUNDED                                 *
  305.  * Ergebnis: Zeile r, rowflag (Art des Pivotschritts)                        *
  306.  *****************************************************************************/
  307.  
  308. USHORT  CHUZR(m,dq,b2q,upper,B,qs,s,r,rowflag)
  309.  
  310. SHORT   m;
  311. DOUBLE  *dq, *b2q, *upper;
  312. SHORT   *B, qs, s, *r;      /* qs: 0,...,n-1   s: 1,...,n-m   r: 1,...,m  */
  313. USHORT  *rowflag;           /* LAMBDA_0, LAMBDA_1 oder LAMBDA_2 */
  314.  
  315. {
  316.   REGISTER DOUBLE lambda0 = INFINITE, lambda1 = INFINITE;
  317.   REGISTER SHORT  i;
  318.   DOUBLE          lambda2 = upper[ABS(qs)-1];
  319.   REGISTER DOUBLE dqi, dummy, min = lambda2;
  320.   SHORT           sigma = SGN(qs), min_i;
  321.  
  322.   *rowflag = LAMBDA_2;
  323.  
  324.   for(i=0; i<m; ++i) {
  325.     dqi = sigma*dq[i];
  326.     if(dqi > EPS_NULL) {
  327.       dummy = b2q[i]/dqi;
  328.       if(lambda0 == INFINITE || dummy<lambda0) lambda0 = dummy;
  329.       if(min==INFINITE || dummy<min || dummy==min &&
  330.          fabs(dqi)>fabs(dq[min_i-1])) {
  331.         lambda0 = min = dummy;
  332.         min_i = i+1;
  333.         *rowflag = LAMBDA_0;
  334.       }
  335.     }
  336.     else if(dqi < -EPS_NULL && (dummy = upper[B[i]-1])!=INFINITE) {
  337.       dummy = (b2q[i]-dummy)/dqi;
  338.       if(lambda1 == INFINITE || dummy<lambda1) lambda1 = dummy;
  339.       if(min==INFINITE || dummy<min || dummy==min &&
  340.          fabs(dqi)>fabs(dq[min_i-1])) {
  341.         lambda1 = min = dummy;
  342.         min_i = i+1;
  343.         *rowflag = LAMBDA_1;
  344.       }
  345.     }
  346.   }
  347.  
  348.   *r = min_i;
  349.  
  350.   if(lambda0==INFINITE && lambda1==INFINITE && lambda2==INFINITE)
  351.     return(UNBOUNDED);
  352.  
  353.   return(NOT_UNBOUNDED);
  354. }
  355.  
  356.  
  357.  
  358. /*****************************************************************************
  359.  * USHORT WRETA() -> INVERTABLE/NOT_INVERTABLE                               *
  360.  * je nachdem, ob AB invertierbar ist (sollte es zumindest) oder nicht       *
  361.  *****************************************************************************/
  362.  
  363. USHORT  WRETA(m,n,B,Nq,A,AB1,b,b2q,c,cq,c0,c0start,dq,upper,r,s,eta,
  364.               help,Nminus,flags,iter)
  365.  
  366. SHORT   m, n;
  367. SHORT   *B, *Nq;
  368. DOUBLE  *A, *AB1, *b, *b2q, *c, *cq, *c0, c0start, *dq, *upper;
  369. SHORT   r, s;
  370. DOUBLE  *eta, *help;         /* eta : (n)-Vektor (=x), als (m)-Vektor verw. */
  371. SHORT   *Nminus;             /* Nminus (n-m)-Vektor, negativer Anteil von Nq */
  372. USHORT  flags;               /* PHASE_I/PHASE_II, rowflag, verbose */
  373. ULONG   iter;                /* Zahl der bewältigten Iterationen */
  374.  
  375. {
  376.   REGISTER DOUBLE *ptr1, *ptr2, eta_r = 1.0/dq[r-1], eta_dummy;
  377.   REGISTER SHORT  j, i;
  378.   LONG            offset = (LONG)(r-1)*(LONG)m;
  379.   SHORT           anzneg = 0; /* Anzahl negativer Elemente in Nq */
  380.   SHORT           nm = n-m, swap;
  381.   USHORT          phase = flags & (PHASE_I | PHASE_II);
  382.   USHORT          rowflag = flags & (LAMBDA_0 | LAMBDA_1 | LAMBDA_2);
  383.   SHORT           columns = phase==PHASE_I ? nm : n;
  384.  
  385.  
  386.   iter %= INVERT_FREQUENCY;
  387.  
  388.   if(rowflag != LAMBDA_2) {
  389.  
  390.     /* Falls Nq[s-1]<0 ist, ändert sich c0 zusätzlich !!!! */
  391.     /* Falls AB1 reinvertiert wird, wird c0 jedoch völlig neu berechnet */
  392.     if(Nq[s-1]<0) *c0 -= cq[s-1]*upper[ABS(Nq[s-1])-1];
  393.  
  394.     swap = B[r-1];
  395.     B[r-1] = ABS(Nq[s-1]);
  396.     Nq[s-1] = rowflag == LAMBDA_0 ? swap : -swap;
  397.  
  398.     for(j=0; j<nm; ++j) {
  399.       if(Nq[j]<0) Nminus[anzneg++] = ABS(Nq[j]);
  400.     }
  401.  
  402.     if(!iter) {  /* AB1 reinvertieren */
  403.  
  404.       if(flags & VERBOSE) {
  405.         printf("@@ ");
  406.         if(phase == PHASE_I) printf("phase I : ");
  407.         else                 printf("phase II: ");
  408.         puts("reinversion");
  409.         if(file[1]) {
  410.           fprintf(file[1],"@@ ");
  411.           if(phase == PHASE_I) fprintf(file[1],"phase I : ");
  412.           else                 fprintf(file[1],"phase II: ");
  413.           fputs("reinversion\n",file[1]);
  414.         }
  415.       }
  416.  
  417.       for(j=1; j<=m; ++j) {
  418.         swap = B[j-1];        /* B[j].te Spalte von A */
  419.         for(i=1; i<=m; ++i) SetM(AB1,m,i,j,M(A,columns,i,swap,phase));
  420.       }
  421.       /* x==eta als (SHORT *)-Permutationsvektor !! */
  422.       if(Invers(AB1,m,(SHORT *)eta) == NOT_INVERTABLE) return(NOT_INVERTABLE);
  423.     }
  424.  
  425.     else {  /* AB1 nicht reinvertieren, sondern nur mit eta mult. */
  426.  
  427.       ptr1 = eta;
  428.       ptr2 = dq;
  429.  
  430.       for(j=1; j<=m; ++j) {
  431.         *ptr1++ = j==r ? eta_r : -(*ptr2)*eta_r;
  432.         ++ptr2;
  433.       }
  434.  
  435.       ptr1 = AB1;
  436.  
  437.       for(i=0; i<m; ++i) {
  438.         ptr2 = AB1+offset;
  439.         if(i+1==r) ptr1 += m; /* r.te Zeile überspringen */
  440.         else {
  441.           if((eta_dummy = eta[i]) != 0.0) {
  442.             for(j=0; j<m; ++j) *ptr1++ += eta_dummy*(*ptr2++);
  443.           }
  444.           else ptr1 += m;
  445.         }
  446.       }
  447.  
  448.       ptr2 = AB1+offset;
  449.  
  450.       for(j=0; j<m; ++j) *ptr2++ *= eta_r;  /* r.te Spalte mit eta_r mult. */
  451.     }
  452.  
  453.  
  454.     for(i=0; i<m; ++i) {
  455.       eta_dummy = b[i];
  456.       for(j=0; j<anzneg; ++j)
  457.         eta_dummy -= M(A,columns,i+1,Nminus[j],phase)*upper[Nminus[j]-1];
  458.       eta[i] = eta_dummy;
  459.     }
  460.  
  461.     Mult(AB1,eta,b2q,m,m,1); /* eta-Vektor als Zwischenspeicher */
  462.  
  463.     if(iter) *c0 += cq[s-1]*b2q[r-1];
  464.     else {
  465.       eta_dummy = 0.0;
  466.       ptr1 = b2q;
  467.       for(i=0; i<m; ++i) eta_dummy += c[B[i]-1]*(*ptr1++);
  468.       for(i=0; i<anzneg; ++i) {
  469.         j = Nminus[i]-1;
  470.         eta_dummy += c[j]*upper[j];
  471.       }
  472.       *c0 = eta_dummy+c0start;
  473.     }
  474.  
  475.   }
  476.  
  477.  
  478.  
  479.   else {   /* rowflag == LAMBDA_2 */
  480.  
  481.     if(!iter) {  /* AB1 reinvertieren */
  482.  
  483.       if(flags & VERBOSE) {
  484.         printf("@@ ");
  485.         if(phase == PHASE_I) printf("phase I : ");
  486.         else                 printf("phase II: ");
  487.         puts("reinversion");
  488.         if(file[1]) {
  489.           fprintf(file[1],"@@ ");
  490.           if(phase == PHASE_I) fprintf(file[1],"phase I : ");
  491.           else                 fprintf(file[1],"phase II: ");
  492.           fputs("reinversion\n",file[1]);
  493.         }
  494.       }
  495.  
  496.       for(j=1; j<=m; ++j) {
  497.         swap = B[j-1];        /* B[j].te Spalte von A */
  498.         for(i=1; i<=m; ++i) SetM(AB1,m,i,j,M(A,columns,i,swap,phase));
  499.       }
  500.       /* x==eta als (SHORT *)-Permutationsvektor !! */
  501.       if(Invers(AB1,m,(SHORT *)eta) == NOT_INVERTABLE) return(NOT_INVERTABLE);
  502.     }
  503.  
  504.     swap = Nq[s-1];               /* swap      : qs_quer */
  505.     anzneg = ABS(swap);           /* anzneg    : qs      */
  506.     eta_dummy = upper[anzneg-1];  /* eta_dummy : u(qs)   */
  507.  
  508.     ptr1 = eta;
  509.     for(i=1; i<=m; ++i) *ptr1++ = M(A,columns,i,anzneg,phase)*eta_dummy;
  510.  
  511.     Mult(AB1,eta,help,m,m,1);
  512.  
  513.     ptr1 = b2q;
  514.     ptr2 = help;
  515.  
  516.     if(swap>0) {    /* wenn also das alte qs_quer positiv war */
  517.       for(i=0; i<m; ++i) *ptr1++ -= *ptr2++;
  518.       *c0 += cq[s-1]*upper[anzneg-1];
  519.     }
  520.     else {
  521.       for(i=0; i<m; ++i) *ptr1++ += *ptr2++;
  522.       *c0 -= cq[s-1]*upper[anzneg-1];
  523.     }
  524.  
  525.     Nq[s-1] *= -1;  /* qs_quer auf neuen Stand bringen */
  526.   }
  527.  
  528.   return(INVERTABLE);
  529. }
  530.