Re-Volt: /Users/frarees/Desktop/rv1/source/ACM

00001 #include "ReVolt.h"
00002 #include "Geom.h"
00003 #include "Gaussian.h"
00004 
00005 /*C     ALGORITHM 478 COLLECTED ALGORITHMS FROM ACM.
00006 C     ALGORITHM APPEARED IN COMM. ACM, VOL. 17, NO. 06,
00007 C     P. 319.
00008 C THIS SUBROUTINE USES A MODIFICATION OF THE SIMPLEX METHOD
00009 C OF LINEAR PROGRAMMING TO CALCULATE AN L1 SOLUTION TO AN
00010 C OVER-DETERMINED SYSTEM OF LINEAR EQUATIONS.
00011 C DESCRIPTION OF PARAMETERS.
00012 C M      NUMBER OF EQUATIONS.
00013 C N      NUMBER OF UNKNOWNS (M.GE.N).
00014 C M2     SET EQUAL TO M+2 FOR ADJUSTABLE DIMENSIONS.
00015 C N2     SET EQUAL TO N+2 FOR ADJUSTABLE DIMENSIONS.
00016 C A      TWO DIMENSIONAL REAL ARRAY OF SIZE (M2,N2).
00017 C        ON ENTRY, THE COEFFICIENTS OF THE MAT MUST BE
00018 C        STORED IN THE FIRST M ROWS AND N COLUMNS OF A.
00019 C        THESE VALUES ARE DESTROYED BY THE SUBROUTINE.
00020 C B      ONE DIMENSIONAL REAL ARRAY OF SIZE M. ON ENTRY, B
00021 C        MUST CONTAIN THE RIGHT HAND SIDE OF THE EQUATIONS.
00022 C        THESE VALUES ARE DESTROYED BY THE SUBROUTINE.
00023 C TOLER  A SMALL POSITIVE TOLERANCE. EMPIRICAL EVIDENCE
00024 C        SUGGESTS TOLER=10**(-D*2/3) WHERE D REPRESENTS
00025 C        THE NUMBER OF DECIMAL DIGITS OF ACCURACY AVALABLE
00026 C        (SEE DESCRIPTION).
00027 C X      ONE DIMENSIONAL REAL ARRAY OF SIZE N. ON EXIT, THIS
00028 C        ARRAY CONTAINS A SOLUTION TO THE L1 PROBLEM.
00029 C E      ONE DIMENSIONAL REAL ARRAY OF SIZE M. ON EXIT, THIS
00030 C        ARRAY CONTAINS THE RESIDUALS IN THE EQUATIONS.
00031 C S      INTEGER ARRAY OF SIZE M USED FOR WORKSPACE.
00032 C ON EXIT FROM THE SUBROUTINE, THE ARRAY A CONTAINS THE
00033 C FOLLOWING INFORMATION.
00034 C A(M+1,N+1)  THE MINIMUM SUM OF THE ABSOLUTE VALUES OF
00035 C             THE RESIDUALS.
00036 C A(M+1,N+2)  THE RANK OF THE MAT OF COEFFICIENTS.
00037 C A(M+2,N+1)  EXIT CODE WITH VALUES.
00038 C             0 - OPTIMAL SOLUTION WHICH IS PROBABLY NON-
00039 C                 UNIQUE (SEE DESCRIPTION).
00040 C             1 - UNIQUE OPTIMAL SOLUTION.
00041 C             2 - CALCULATIONS TERMINATED PREMATURELY DUE TO
00042 C                 ROUNDING ERRORS.
00043 C A(M+2,N+2)  NUMBER OF SIMPLEX ITERATIONS PERFORMED.
00044 */
00045 
00046 void L1(BIGMAT *a, BIGVEC *b, REAL toler, BIGVEC *x, BIGVEC *e, int *s)
00047 {
00048   int   m, m1, m2, n, n1, n2, out, i, j, k, l, kount, kr, kl, in;
00049   REAL  sum;
00050   REAL  min, max, d, pivot;
00051   bool  stage, test;
00052 
00053   // big must be set equal to any very large real constant.
00054     const REAL big = 1.0e7f;
00055   
00056   m = NRows(a);
00057   n = NCols(a);
00058 
00059 
00060   // initialization.
00061   m1 = m + 1;
00062   n1 = n + 1;
00063   m2 = m + 2;
00064   n2 = n + 2;
00065   
00066   for (j = 0; j < n; j++) {
00067     a->m[m1][j] = (REAL)j;
00068         x->v[j] = ZERO;
00069   }
00070   for (i = 0; i < m; i++) {
00071         a->m[i][n1] = (REAL)(n + i + 1);    // possible error (probably not)
00072     a->m[i][n] = b->v[i];       
00073     if (b->v[i] < ZERO) {
00074       for (j = 0; j < n2; j++) {
00075         a->m[i][j] = -a->m[i][j];
00076       }
00077     }
00078     e->v[i] = ZERO;
00079   }
00080 
00081   // compute the marginal costs.
00082   for (j = 0; j < n1; j++) {
00083     sum = ZERO;
00084     for (i = 0; i < m; i++) {
00085       sum += a->m[i][j];
00086     }
00087         a->m[m][j] = sum;
00088   }
00089 
00090   // stage i.
00091   // determine the vector to enter the basis.
00092   stage = TRUE;
00093   test = FALSE;
00094   kount = -1;           // possible error
00095   kr = 0;
00096   kl = 0;
00097 
00098 l70:
00099   max = -ONE;
00100   for (j = kr; j < n; j++) {
00101     if (abs(a->m[m1][j]) <= (REAL)n) {
00102       d = abs(a->m[m][j]);
00103       if (d > max) {
00104         max = d;
00105         in = j;
00106       }
00107     }
00108   }
00109   if (a->m[m][in] < ZERO) {
00110     for (i=0; i < m2; i++) {
00111       a->m[i][in] = -a->m[i][in];
00112     }
00113   }
00114 
00115   // determine the vector to leave the basis.
00116 l100:
00117   k = -1;
00118   for (i = kl; i < m; i++) {
00119     d = a->m[i][in];
00120     if (d > toler) {
00121       k++;
00122       b->v[k] = a->m[i][n1] / d;
00123       s[k] = i;
00124       test = TRUE;
00125     }
00126   }
00127 l120:
00128   if (k >= 0) {
00129     min = big;
00130     for (i = 0; i <= k; i++) {
00131       if (b->v[i] < min) {
00132         j = i;
00133         min = b->v[i];
00134         out = s[i];
00135       }
00136     }
00137     b->v[j] = b->v[k];
00138     s[j] = s[k];
00139     k = k - 1;      //possible problem
00140   } else {
00141     test = FALSE;
00142   }
00143 
00144   // check for linear dependence in stage i.
00145   if (!(test || !stage)) {
00146     for (i = 0; i < m2; i++) {
00147       d = a->m[i][kr];
00148       a->m[i][kr] = a->m[i][in];
00149       a->m[i][in] = d;
00150     }
00151     kr++;         // possible problem
00152     goto l260;
00153   }
00154   if (!test) {
00155     a->m[m1][n] = Real(2);
00156     goto l350;
00157   }
00158   pivot = a->m[out][in];
00159   if (a->m[m][in] - pivot - pivot > toler) {
00160     for (j = kr; j < n1; j++) {
00161       d = a->m[out][j];
00162       a->m[m][j] = a->m[m][j] - d - d;
00163       a->m[out][j] = -d;
00164     }
00165     a->m[out][n1] = -a->m[out][n1];
00166     goto l120;
00167   }
00168 
00169   // pivot on a(out,in).
00170   for (j = kr; j < n1; j++) {
00171     if (j != in) {
00172       a->m[out][j] = a->m[out][j] / pivot;
00173     }
00174   }
00175   for (i = 0; i < m1; i++) {
00176     if (i != out) {
00177       d = a->m[i][in];
00178       for (j = kr; j < n1; j++) {
00179         if (j != in) {
00180           a->m[i][j] -= d * a->m[out][j];
00181         }
00182       }
00183     }
00184   }
00185   for (i = 0; i < m1; i++) {
00186     if (i != out) {
00187       a->m[i][in] = -a->m[i][in] / pivot;
00188     }
00189   }
00190   a->m[out][in] = ONE / pivot;
00191   d = a->m[out][n1];
00192   a->m[out][n1] = a->m[m1][in];
00193   a->m[m1][in] = d;
00194   kount++;            // possible error
00195   if (stage) {
00196   // interchange rows in stage i.
00197     kl++;
00198     for (j = kr; j < n2; j++) {
00199       d = a->m[out][j];
00200       a->m[out][j] = a->m[kount][j];
00201       a->m[kount][j] = d;
00202     }
00203 l260:
00204     if (kount + kr != n1 - 2) goto l70;   // possible error
00205     // stage ii.
00206     stage = FALSE;
00207   }
00208   
00209   // determine the vector to enter the basis.
00210   max = -big;
00211   for (j = kr; j < n; j++) {
00212     d = a->m[m][j];
00213     if (d >= ZERO) goto l280;
00214     if (d > Real(-2)) goto l290;
00215         d = - (d + Real(2));
00216 l280:
00217     if (d <= max) goto l290;
00218     max = d;
00219     in = j;
00220   }
00221 l290:
00222   if (max > toler) {
00223     if (a->m[m][in] <= ZERO) {
00224       for (i = 0; i < m2; i++) {
00225         a->m[i][in] = -a->m[i][in];
00226       }
00227       a->m[m][in] = a->m[m][in] - Real(2);
00228     }
00229     goto l100;
00230   }
00231 
00232   //prepare output.
00233   l = kl - 1;       // possible error
00234   for (i = 0; i < l; i++) {
00235     if (a->m[i][n] < ZERO) {
00236       for (j = kr; j < n2; j++) {
00237         a->m[i][j] = -a->m[i][j];
00238       }
00239     }
00240   }
00241   a->m[m1][n] = ZERO;
00242   if (kr == 1) {
00243     for (j = 0; j < n; j++) {
00244       d = abs(a->m[m][j]);
00245       if ((d <= toler) || (Real(2 - d) <= toler)) goto l350;
00246     }
00247     a->m[m][n1] = ONE;
00248   }
00249 l350:
00250   for (i = 0; i < m; i++) {
00251     k = NearestInt(a->m[i][n1]);
00252     d = a->m[i][n];
00253     if (k <= 0) {
00254       k = -k;
00255       d = -d;
00256     }
00257     if (i < kl) {
00258       x->v[k] = d;
00259     } else {
00260       k = k - n;
00261       e->v[k] = d;
00262     }
00263   }
00264   a->m[m1][n1] = (REAL)kount;
00265   a->m[m][n1] = (REAL)n1 - kr;
00266   sum = ZERO;
00267   for (i = kl; i < m; i++) {
00268     sum += a->m[i][n];
00269   }
00270   a->m[m][n] = sum;
00271 
00272 }
00273 
00274 
00275 
00277 //
00278 // TestSolver: do a few tests on the linear problem solver
00279 //
00281 
00282 static BIGMAT Coef, OrigCoef;
00283 static BIGVEC Res, NewRes, OrigRes;
00284 static BIGVEC Soln;
00285 static int    Work[BIG_NMAX], OrigRow[BIG_NMAX], OrigCol[BIG_NMAX];
00286 
00287 
00288 #ifdef _PC
00289 void TestL1Solver()
00290 {
00291   int   ii, kk;
00292   REAL  xx;
00293 
00295   // Test 1: Simple solution
00296   //
00297   SetBigMatSize(&Coef, 4, 4);
00298   SetBigVecSize(&Res, 4);
00299 
00300   Coef.m[0][0] = 1.0f;
00301   Coef.m[0][1] = 0.0f;
00302   Coef.m[0][2] = 1.0f;
00303   Coef.m[0][3] = 0.0f;
00304 
00305   Coef.m[1][0] = 1.0f;
00306   Coef.m[1][1] = 0.0f;
00307   Coef.m[1][2] = 1.0f;
00308   Coef.m[1][3] = 0.0f;
00309 
00310   Coef.m[2][0] = 0.0f;
00311   Coef.m[2][1] = 1.0f;
00312   Coef.m[2][2] = 0.0f;
00313   Coef.m[2][3] = 1.0f;
00314 
00315   Coef.m[3][0] = 0.0f;
00316   Coef.m[3][1] = 1.0f;
00317   Coef.m[3][2] = 1.0f;
00318   Coef.m[3][3] = 1.0f;
00319 
00320   Res.v[0] = 2.0f;
00321   Res.v[1] = 2.0f;
00322   Res.v[2] = 2.0f;
00323   Res.v[3] = 2.0f;
00324 
00325   CopyBigMat(&Coef, &OrigCoef);
00326   CopyBigVec(&Res, &OrigRes);
00327 
00328   L1(&Coef, &Res, 0.0001f, &Soln, &NewRes, Work);
00329   SolveLinearEquations(&OrigCoef, &OrigRes, 0.0001f, 0.0001f, OrigRow, OrigCol, &NewRes, &Soln);
00330 
00332   // Test 2:
00333   //
00334   SetBigMatSize(&Coef, 2, 2);
00335   SetBigVecSize(&Res, 2);
00336 
00337   Coef.m[0][0] = 1.0f;
00338   Coef.m[0][1] = 0.0f;
00339   Coef.m[1][0] = 1.01f;
00340   Coef.m[1][1] = 0.0f;
00341 
00342   Res.v[0] = 2.0f;
00343   Res.v[1] = 1.0f;
00344 
00345   CopyBigMat(&Coef, &OrigCoef);
00346   CopyBigVec(&Res, &OrigRes);
00347 
00348   L1(&Coef, &Res, 0.0001f, &Soln, &NewRes, Work);
00349   SolveLinearEquations(&OrigCoef, &OrigRes, 0.0001f, 0.0001f, OrigRow, OrigCol, &NewRes, &Soln);
00350 
00351 
00353   // Test 3:
00354   //
00355   SetBigMatSize(&Coef, 5, 5);
00356   SetBigVecSize(&Res, 5);
00357   SetBigVecSize(&Soln, 5);
00358   for (ii = 0; ii < NRows(&Coef); ii++) {
00359     xx = (ii + 1) * Real(0.25);
00360     for (kk = 0; kk < NRows(&Coef); kk++) {
00361       Coef.m[ii][kk] = (REAL)pow(xx, kk + 1);
00362     }
00363     Res.v[ii] = xx * (REAL)exp(-xx);
00364   }
00365 
00366   CopyBigMat(&Coef, &OrigCoef);
00367   CopyBigVec(&Res, &OrigRes);
00368 
00369   L1(&Coef, &Res, 0.0001f, &Soln, &NewRes, Work);
00370   SolveLinearEquations(&OrigCoef, &OrigRes, 0.0001f, 0.0001f, OrigRow, OrigCol, &NewRes, &Soln);
00371 
00373   // Test 4: 
00374   //
00375   SetBigMatSize(&Coef, 9, 9);
00376   SetBigVecSize(&Res, 9);
00377   SetBigVecSize(&Soln, 9);
00378   ClearBigMat(&Coef);
00379   ClearBigVec(&Res);
00380 
00381   for (ii = 0; ii < NRows(&Coef); ii++) {
00382     xx = ii / Real(8.0);
00383     Coef.m[ii][0] = ONE;
00384     Coef.m[ii][1] = xx;
00385     Coef.m[ii][2] = xx - ONE;
00386     Coef.m[ii][3] = (REAL)pow(xx, 2);
00387     Coef.m[ii][4] = (REAL)pow(xx, 2) - xx;
00388     Coef.m[ii][5] = (REAL)pow(xx, 3);
00389     Coef.m[ii][6] = (REAL)pow(xx, 3) - (REAL)pow(xx, 2);
00390     Coef.m[ii][7] = (REAL)pow(xx, 4);
00391     Coef.m[ii][8] = (REAL)pow(xx, 4) - (REAL)pow(xx, 3);
00392     Res.v[ii] = (REAL)exp(xx);
00393   }
00394 
00395   CopyBigMat(&Coef, &OrigCoef);
00396   CopyBigVec(&Res, &OrigRes);
00397 
00398   L1(&Coef, &Res, 0.0001f, &Soln, &NewRes, Work);
00399   SolveLinearEquations(&OrigCoef, &OrigRes, 0.0001f, 0.0001f, OrigRow, OrigCol, &NewRes, &Soln);
00400 
00402   // Test 5: 
00403   //
00404   SetBigMatSize(&Coef, 2, 2);
00405   SetBigVecSize(&Res, 2);
00406   SetBigVecSize(&Soln, 2);
00407   ClearBigMat(&Coef);
00408   ClearBigVec(&Res);
00409 
00410   Coef.m[0][0] = ONE;
00411   Coef.m[0][1] = ZERO;
00412   Coef.m[1][0] = ZERO;
00413   Coef.m[1][1] = 0.00001f;
00414   Res.v[0] = ONE;
00415   Res.v[1] = 1.00001f;
00416 
00417   CopyBigMat(&Coef, &OrigCoef);
00418   CopyBigVec(&Res, &OrigRes);
00419 
00420   L1(&Coef, &Res, 0.0001f, &Soln, &NewRes, Work);
00421   SolveLinearEquations(&OrigCoef, &OrigRes, 0.0001f, 0.0001f, OrigRow, OrigCol, &NewRes, &Soln);
00422 
00423 }
00424 #endif
00425
/Users/frarees/Desktop/rv1/source/ACM_Simplex..cpp
