Ignore:
Timestamp:
Jan 6, 2019 2:43:06 PM (3 months ago)
Author:
unxusr
Message:

formatting

File:
1 edited

Legend:

Unmodified
Added
Removed
  • trunk/Clp/src/IdiSolve.cpp

    r2271 r2385  
    1717#define HISTORY 7
    1818#endif
    19 #define NSOLVE HISTORY-1
    20 static void solveSmall(int nsolve, double **aIn, double **a, double * b)
     19#define NSOLVE HISTORY - 1
     20static void solveSmall(int nsolve, double **aIn, double **a, double *b)
    2121{
    22      int i, j;
    23      /* copy */
    24      for (i = 0; i < nsolve; i++) {
    25           for (j = 0; j < nsolve; j++) {
    26                a[i][j] = aIn[i][j];
    27           }
    28      }
    29      for (i = 0; i < nsolve; i++) {
    30           /* update using all previous */
    31           double diagonal;
    32           int j;
    33           for (j = i; j < nsolve; j++) {
    34                int k;
    35                double value = a[i][j];
    36                for (k = 0; k < i; k++) {
    37                     value -= a[k][i] * a[k][j];
    38                }
    39                a[i][j] = value;
    40           }
    41           diagonal = a[i][i];
    42           if (diagonal < 1.0e-20) {
    43                diagonal = 0.0;
    44           } else {
    45                diagonal = 1.0 / sqrt(diagonal);
    46           }
    47           a[i][i] = diagonal;
    48           for (j = i + 1; j < nsolve; j++) {
    49                a[i][j] *= diagonal;
    50           }
    51      }
    52      /* update */
    53      for (i = 0; i < nsolve; i++) {
    54           int j;
    55           double value = b[i];
    56           for (j = 0; j < i; j++) {
    57                value -= b[j] * a[j][i];
    58           }
    59           value *= a[i][i];
    60           b[i] = value;
    61      }
    62      for (i = nsolve - 1; i >= 0; i--) {
    63           int j;
    64           double value = b[i];
    65           for (j = i + 1; j < nsolve; j++) {
    66                value -= b[j] * a[i][j];
    67           }
    68           value *= a[i][i];
    69           b[i] = value;
    70      }
     22  int i, j;
     23  /* copy */
     24  for (i = 0; i < nsolve; i++) {
     25    for (j = 0; j < nsolve; j++) {
     26      a[i][j] = aIn[i][j];
     27    }
     28  }
     29  for (i = 0; i < nsolve; i++) {
     30    /* update using all previous */
     31    double diagonal;
     32    int j;
     33    for (j = i; j < nsolve; j++) {
     34      int k;
     35      double value = a[i][j];
     36      for (k = 0; k < i; k++) {
     37        value -= a[k][i] * a[k][j];
     38      }
     39      a[i][j] = value;
     40    }
     41    diagonal = a[i][i];
     42    if (diagonal < 1.0e-20) {
     43      diagonal = 0.0;
     44    } else {
     45      diagonal = 1.0 / sqrt(diagonal);
     46    }
     47    a[i][i] = diagonal;
     48    for (j = i + 1; j < nsolve; j++) {
     49      a[i][j] *= diagonal;
     50    }
     51  }
     52  /* update */
     53  for (i = 0; i < nsolve; i++) {
     54    int j;
     55    double value = b[i];
     56    for (j = 0; j < i; j++) {
     57      value -= b[j] * a[j][i];
     58    }
     59    value *= a[i][i];
     60    b[i] = value;
     61  }
     62  for (i = nsolve - 1; i >= 0; i--) {
     63    int j;
     64    double value = b[i];
     65    for (j = i + 1; j < nsolve; j++) {
     66      value -= b[j] * a[i][j];
     67    }
     68    value *= a[i][i];
     69    b[i] = value;
     70  }
    7171}
    7272IdiotResult
    73 Idiot::objval(int nrows, int ncols, double * rowsol , double * colsol,
    74               double * pi, double * /*djs*/, const double * cost ,
    75               const double * /*rowlower*/,
    76               const double * rowupper, const double * /*lower*/,
    77               const double * /*upper*/, const double * elemnt,
    78               const int * row, const CoinBigIndex * columnStart,
    79               const int * length, int extraBlock, int * rowExtra,
    80               double * solExtra, double * elemExtra, double * /*upperExtra*/,
    81               double * costExtra, double weight)
     73Idiot::objval(int nrows, int ncols, double *rowsol, double *colsol,
     74  double *pi, double * /*djs*/, const double *cost,
     75  const double * /*rowlower*/,
     76  const double *rowupper, const double * /*lower*/,
     77  const double * /*upper*/, const double *elemnt,
     78  const int *row, const CoinBigIndex *columnStart,
     79  const int *length, int extraBlock, int *rowExtra,
     80  double *solExtra, double *elemExtra, double * /*upperExtra*/,
     81  double *costExtra, double weight)
    8282{
    83      IdiotResult result;
    84      double objvalue = 0.0;
    85      double sum1 = 0.0, sum2 = 0.0;
    86      int i;
    87      for (i = 0; i < nrows; i++) {
    88           rowsol[i] = -rowupper[i];
    89      }
    90      for (i = 0; i < ncols; i++) {
    91           CoinBigIndex j;
    92           double value = colsol[i];
    93           if (value) {
    94                objvalue += value * cost[i];
    95                if (elemnt) {
    96                     for (j = columnStart[i]; j < columnStart[i] + length[i]; j++) {
    97                          int irow = row[j];
    98                          rowsol[irow] += elemnt[j] * value;
    99                     }
    100                } else {
    101                     for (j = columnStart[i]; j < columnStart[i] + length[i]; j++) {
    102                          int irow = row[j];
    103                          rowsol[irow] += value;
    104                     }
    105                }
    106           }
    107      }
    108      /* adjust to make as feasible as possible */
    109      /* no */
    110      if (extraBlock) {
    111           for (i = 0; i < extraBlock; i++) {
    112                double element = elemExtra[i];
    113                int irow = rowExtra[i];
    114                objvalue += solExtra[i] * costExtra[i];
    115                rowsol[irow] += solExtra[i] * element;
    116           }
    117      }
    118      for (i = 0; i < nrows; i++) {
    119           double value = rowsol[i];
    120           sum1 += fabs(value);
    121           sum2 += value * value;
    122           pi[i] = -2.0 * weight * value;
    123      }
    124      result.infeas = sum1;
    125      result.objval = objvalue;
    126      result.weighted = objvalue + weight * sum2;
    127      result.dropThis = 0.0;
    128      result.sumSquared = sum2;
    129      return result;
     83  IdiotResult result;
     84  double objvalue = 0.0;
     85  double sum1 = 0.0, sum2 = 0.0;
     86  int i;
     87  for (i = 0; i < nrows; i++) {
     88    rowsol[i] = -rowupper[i];
     89  }
     90  for (i = 0; i < ncols; i++) {
     91    CoinBigIndex j;
     92    double value = colsol[i];
     93    if (value) {
     94      objvalue += value * cost[i];
     95      if (elemnt) {
     96        for (j = columnStart[i]; j < columnStart[i] + length[i]; j++) {
     97          int irow = row[j];
     98          rowsol[irow] += elemnt[j] * value;
     99        }
     100      } else {
     101        for (j = columnStart[i]; j < columnStart[i] + length[i]; j++) {
     102          int irow = row[j];
     103          rowsol[irow] += value;
     104        }
     105      }
     106    }
     107  }
     108  /* adjust to make as feasible as possible */
     109  /* no */
     110  if (extraBlock) {
     111    for (i = 0; i < extraBlock; i++) {
     112      double element = elemExtra[i];
     113      int irow = rowExtra[i];
     114      objvalue += solExtra[i] * costExtra[i];
     115      rowsol[irow] += solExtra[i] * element;
     116    }
     117  }
     118  for (i = 0; i < nrows; i++) {
     119    double value = rowsol[i];
     120    sum1 += fabs(value);
     121    sum2 += value * value;
     122    pi[i] = -2.0 * weight * value;
     123  }
     124  result.infeas = sum1;
     125  result.objval = objvalue;
     126  result.weighted = objvalue + weight * sum2;
     127  result.dropThis = 0.0;
     128  result.sumSquared = sum2;
     129  return result;
    130130}
    131131IdiotResult
    132132Idiot::IdiSolve(
    133      int nrows, int ncols, double * COIN_RESTRICT rowsol , double * COIN_RESTRICT colsol,
    134      double * COIN_RESTRICT pi, double * COIN_RESTRICT djs, const double * COIN_RESTRICT origcost , double * COIN_RESTRICT rowlower,
    135      double * COIN_RESTRICT rowupper, const double * COIN_RESTRICT lower,
    136      const double * COIN_RESTRICT upper, const double * COIN_RESTRICT elemnt,
    137      const int * row, const CoinBigIndex * columnStart,
    138      const int * length, double * COIN_RESTRICT lambda,
    139      int maxIts, double mu, double drop,
    140      double maxmin, double offset,
    141      int strategy, double djTol, double djExit, double djFlag,
    142      CoinThreadRandom * randomNumberGenerator)
     133  int nrows, int ncols, double *COIN_RESTRICT rowsol, double *COIN_RESTRICT colsol,
     134  double *COIN_RESTRICT pi, double *COIN_RESTRICT djs, const double *COIN_RESTRICT origcost, double *COIN_RESTRICT rowlower,
     135  double *COIN_RESTRICT rowupper, const double *COIN_RESTRICT lower,
     136  const double *COIN_RESTRICT upper, const double *COIN_RESTRICT elemnt,
     137  const int *row, const CoinBigIndex *columnStart,
     138  const int *length, double *COIN_RESTRICT lambda,
     139  int maxIts, double mu, double drop,
     140  double maxmin, double offset,
     141  int strategy, double djTol, double djExit, double djFlag,
     142  CoinThreadRandom *randomNumberGenerator)
    143143{
    144      IdiotResult result;
    145      int i, k, iter;
    146      CoinBigIndex j;
    147      double value = 0.0, objvalue = 0.0, weightedObj = 0.0;
    148      double tolerance = 1.0e-8;
    149      double * history[HISTORY+1];
    150      int ncolx;
    151      int nChange;
    152      int extraBlock = 0;
    153      int * rowExtra = NULL;
    154      double * COIN_RESTRICT solExtra = NULL;
    155      double * COIN_RESTRICT elemExtra = NULL;
    156      double * COIN_RESTRICT upperExtra = NULL;
    157      double * COIN_RESTRICT costExtra = NULL;
    158      double * COIN_RESTRICT useCostExtra = NULL;
    159      double * COIN_RESTRICT saveExtra = NULL;
    160      double * COIN_RESTRICT cost = NULL;
    161      double saveValue = 1.0e30;
    162      double saveOffset = offset;
    163      double useOffset = offset;
    164      /*#define NULLVECTOR*/
     144  IdiotResult result;
     145  int i, k, iter;
     146  CoinBigIndex j;
     147  double value = 0.0, objvalue = 0.0, weightedObj = 0.0;
     148  double tolerance = 1.0e-8;
     149  double *history[HISTORY + 1];
     150  int ncolx;
     151  int nChange;
     152  int extraBlock = 0;
     153  int *rowExtra = NULL;
     154  double *COIN_RESTRICT solExtra = NULL;
     155  double *COIN_RESTRICT elemExtra = NULL;
     156  double *COIN_RESTRICT upperExtra = NULL;
     157  double *COIN_RESTRICT costExtra = NULL;
     158  double *COIN_RESTRICT useCostExtra = NULL;
     159  double *COIN_RESTRICT saveExtra = NULL;
     160  double *COIN_RESTRICT cost = NULL;
     161  double saveValue = 1.0e30;
     162  double saveOffset = offset;
     163  double useOffset = offset;
     164  /*#define NULLVECTOR*/
    165165#ifndef NULLVECTOR
    166      int nsolve = NSOLVE;
     166  int nsolve = NSOLVE;
    167167#else
    168      int nsolve = NSOLVE + 1; /* allow for null vector */
    169 #endif
    170      int nflagged;
    171      double * COIN_RESTRICT thetaX;
    172      double * COIN_RESTRICT djX;
    173      double * COIN_RESTRICT bX;
    174      double * COIN_RESTRICT vX;
    175      double ** aX;
    176      double **aworkX;
    177      double ** allsum;
    178      double * COIN_RESTRICT saveSol = 0;
    179      const double * COIN_RESTRICT useCost = cost;
    180      double bestSol = 1.0e60;
    181      double weight = 0.5 / mu;
    182      char * statusSave = new char[2*ncols];
    183      char * statusWork = statusSave + ncols;
     168  int nsolve = NSOLVE + 1; /* allow for null vector */
     169#endif
     170  int nflagged;
     171  double *COIN_RESTRICT thetaX;
     172  double *COIN_RESTRICT djX;
     173  double *COIN_RESTRICT bX;
     174  double *COIN_RESTRICT vX;
     175  double **aX;
     176  double **aworkX;
     177  double **allsum;
     178  double *COIN_RESTRICT saveSol = 0;
     179  const double *COIN_RESTRICT useCost = cost;
     180  double bestSol = 1.0e60;
     181  double weight = 0.5 / mu;
     182  char *statusSave = new char[2 * ncols];
     183  char *statusWork = statusSave + ncols;
    184184#define DJTEST 5
    185      double djSave[DJTEST];
    186      double largestDj = 0.0;
    187      double smallestDj = 1.0e60;
    188      double maxDj = 0.0;
    189      int doFull = 0;
     185  double djSave[DJTEST];
     186  double largestDj = 0.0;
     187  double smallestDj = 1.0e60;
     188  double maxDj = 0.0;
     189  int doFull = 0;
    190190#define SAVEHISTORY 10
    191 #define EVERY (2*SAVEHISTORY)
    192 #define AFTER SAVEHISTORY*(HISTORY+1)
     191#define EVERY (2 * SAVEHISTORY)
     192#define AFTER SAVEHISTORY *(HISTORY + 1)
    193193#define DROP 5
    194      double after = AFTER;
    195      double obj[DROP];
    196      double kbad = 0, kgood = 0;
    197      if (strategy & 128) after = 999999; /* no acceleration at all */
    198      for (i = 0; i < DROP; i++) {
    199           obj[i] = 1.0e70;
    200      }
    201      //#define FOUR_GOES 2
     194  double after = AFTER;
     195  double obj[DROP];
     196  double kbad = 0, kgood = 0;
     197  if (strategy & 128)
     198    after = 999999; /* no acceleration at all */
     199  for (i = 0; i < DROP; i++) {
     200    obj[i] = 1.0e70;
     201  }
     202  //#define FOUR_GOES 2
    202203#ifdef FOUR_GOES
    203      double * COIN_RESTRICT pi2 = new double [3*nrows];
    204      double * COIN_RESTRICT rowsol2 = new double [3*nrows];
    205      double * COIN_RESTRICT piX[4];
    206      double * COIN_RESTRICT rowsolX[4];
    207      int startsX[2][5];
    208      int nChangeX[4];
    209      double maxDjX[4];
    210      double objvalueX[4];
    211      int nflaggedX[4];
    212      piX[0]=pi;
    213      piX[1]=pi2;
    214      piX[2]=pi2+nrows;
    215      piX[3]=piX[2]+nrows;
    216      rowsolX[0]=rowsol;
    217      rowsolX[1]=rowsol2;
    218      rowsolX[2]=rowsol2+nrows;
    219      rowsolX[3]=rowsolX[2]+nrows;
    220 #endif
    221      allsum = new double * [nsolve];
    222      aX = new double * [nsolve];
    223      aworkX = new double * [nsolve];
    224      thetaX = new double[nsolve];
    225      vX = new double[nsolve];
    226      bX = new double[nsolve];
    227      djX = new double[nsolve];
    228      allsum[0] = pi;
    229      for (i = 0; i < nsolve; i++) {
    230           if (i) allsum[i] = new double[nrows];
    231           aX[i] = new double[nsolve];
    232           aworkX[i] = new double[nsolve];
    233      }
    234      /* check = rows */
    235      for (i = 0; i < nrows; i++) {
    236           if (rowupper[i] - rowlower[i] > tolerance) {
    237                extraBlock++;
    238           }
    239      }
    240      cost = new double[ncols];
    241      memset(rowsol, 0, nrows * sizeof(double));
    242      for (i = 0; i < ncols; i++) {
    243           CoinBigIndex j;
    244           double value = origcost[i] * maxmin;
    245           double value2 = colsol[i];
    246           if (elemnt) {
    247                for (j = columnStart[i]; j < columnStart[i] + length[i]; j++) {
    248                     int irow = row[j];
    249                     value += elemnt[j] * lambda[irow];
    250                     rowsol[irow] += elemnt[j] * value2;
    251                }
     204  double *COIN_RESTRICT pi2 = new double[3 * nrows];
     205  double *COIN_RESTRICT rowsol2 = new double[3 * nrows];
     206  double *COIN_RESTRICT piX[4];
     207  double *COIN_RESTRICT rowsolX[4];
     208  int startsX[2][5];
     209  int nChangeX[4];
     210  double maxDjX[4];
     211  double objvalueX[4];
     212  int nflaggedX[4];
     213  piX[0] = pi;
     214  piX[1] = pi2;
     215  piX[2] = pi2 + nrows;
     216  piX[3] = piX[2] + nrows;
     217  rowsolX[0] = rowsol;
     218  rowsolX[1] = rowsol2;
     219  rowsolX[2] = rowsol2 + nrows;
     220  rowsolX[3] = rowsolX[2] + nrows;
     221#endif
     222  allsum = new double *[nsolve];
     223  aX = new double *[nsolve];
     224  aworkX = new double *[nsolve];
     225  thetaX = new double[nsolve];
     226  vX = new double[nsolve];
     227  bX = new double[nsolve];
     228  djX = new double[nsolve];
     229  allsum[0] = pi;
     230  for (i = 0; i < nsolve; i++) {
     231    if (i)
     232      allsum[i] = new double[nrows];
     233    aX[i] = new double[nsolve];
     234    aworkX[i] = new double[nsolve];
     235  }
     236  /* check = rows */
     237  for (i = 0; i < nrows; i++) {
     238    if (rowupper[i] - rowlower[i] > tolerance) {
     239      extraBlock++;
     240    }
     241  }
     242  cost = new double[ncols];
     243  memset(rowsol, 0, nrows * sizeof(double));
     244  for (i = 0; i < ncols; i++) {
     245    CoinBigIndex j;
     246    double value = origcost[i] * maxmin;
     247    double value2 = colsol[i];
     248    if (elemnt) {
     249      for (j = columnStart[i]; j < columnStart[i] + length[i]; j++) {
     250        int irow = row[j];
     251        value += elemnt[j] * lambda[irow];
     252        rowsol[irow] += elemnt[j] * value2;
     253      }
     254    } else {
     255      for (j = columnStart[i]; j < columnStart[i] + length[i]; j++) {
     256        int irow = row[j];
     257        value += lambda[irow];
     258        rowsol[irow] += value2;
     259      }
     260    }
     261    cost[i] = value;
     262  }
     263  if (extraBlock) {
     264    rowExtra = new int[extraBlock];
     265    solExtra = new double[extraBlock];
     266    elemExtra = new double[extraBlock];
     267    upperExtra = new double[extraBlock];
     268    costExtra = new double[extraBlock];
     269    saveExtra = new double[extraBlock];
     270    extraBlock = 0;
     271    int nbad = 0;
     272    for (i = 0; i < nrows; i++) {
     273      if (rowupper[i] - rowlower[i] > tolerance) {
     274        double smaller, difference;
     275        double value;
     276        saveExtra[extraBlock] = rowupper[i];
     277        if (fabs(rowupper[i]) > fabs(rowlower[i])) {
     278          smaller = rowlower[i];
     279          value = -1.0;
     280        } else {
     281          smaller = rowupper[i];
     282          value = 1.0;
     283        }
     284        if (fabs(smaller) > 1.0e10) {
     285          if (!nbad)
     286            COIN_DETAIL_PRINT(printf("Can't handle rows where both bounds >1.0e10 %d %g\n",
     287              i, smaller));
     288          nbad++;
     289          if (rowupper[i] < 0.0 || rowlower[i] > 0.0)
     290            abort();
     291          if (fabs(rowupper[i]) > fabs(rowlower[i])) {
     292            rowlower[i] = -0.9e10;
     293            smaller = rowlower[i];
     294            value = -1.0;
    252295          } else {
    253                for (j = columnStart[i]; j < columnStart[i] + length[i]; j++) {
    254                     int irow = row[j];
    255                     value += lambda[irow];
    256                     rowsol[irow] += value2;
    257                }
    258           }
    259           cost[i] = value;
    260      }
    261      if (extraBlock) {
    262           rowExtra = new int[extraBlock];
    263           solExtra = new double[extraBlock];
    264           elemExtra = new double[extraBlock];
    265           upperExtra = new double[extraBlock];
    266           costExtra = new double[extraBlock];
    267           saveExtra = new double[extraBlock];
    268           extraBlock = 0;
    269           int nbad = 0;
    270           for (i = 0; i < nrows; i++) {
    271                if (rowupper[i] - rowlower[i] > tolerance) {
    272                     double smaller, difference;
    273                     double value;
    274                     saveExtra[extraBlock] = rowupper[i];
    275                     if (fabs(rowupper[i]) > fabs(rowlower[i])) {
    276                          smaller = rowlower[i];
    277                          value = -1.0;
    278                     } else {
    279                          smaller = rowupper[i];
    280                          value = 1.0;
    281                     }
    282                     if (fabs(smaller) > 1.0e10) {
    283                          if (!nbad)
    284                               COIN_DETAIL_PRINT(printf("Can't handle rows where both bounds >1.0e10 %d %g\n",
    285                                                        i, smaller));
    286                          nbad++;
    287                          if (rowupper[i] < 0.0 || rowlower[i] > 0.0)
    288                               abort();
    289                          if (fabs(rowupper[i]) > fabs(rowlower[i])) {
    290                               rowlower[i] = -0.9e10;
    291                               smaller = rowlower[i];
    292                               value = -1.0;
    293                          } else {
    294                               rowupper[i] = 0.9e10;
    295                               saveExtra[extraBlock] = rowupper[i];
    296                               smaller = rowupper[i];
    297                               value = 1.0;
    298                          }
    299                     }
    300                     difference = rowupper[i] - rowlower[i];
    301                     difference = CoinMin(difference, 1.0e31);
    302                     rowupper[i] = smaller;
    303                     elemExtra[extraBlock] = value;
    304                     solExtra[extraBlock] = (rowupper[i] - rowsol[i]) / value;
    305                     if (solExtra[extraBlock] < 0.0) solExtra[extraBlock] = 0.0;
    306                     if (solExtra[extraBlock] > difference) solExtra[extraBlock] = difference;
    307                     costExtra[extraBlock] = lambda[i] * value;
    308                     upperExtra[extraBlock] = difference;
    309                     rowsol[i] += value * solExtra[extraBlock];
    310                     rowExtra[extraBlock++] = i;
    311                }
    312           }
    313           if (nbad)
    314             COIN_DETAIL_PRINT(printf("%d bad values - results may be wrong\n", nbad));
    315      }
    316      for (i = 0; i < nrows; i++) {
    317           offset += lambda[i] * rowsol[i];
    318      }
    319      if ((strategy & 256) != 0) {
    320           /* save best solution */
    321           saveSol = new double[ncols];
    322           CoinMemcpyN(colsol, ncols, saveSol);
    323           if (extraBlock) {
    324                useCostExtra = new double[extraBlock];
    325                memset(useCostExtra, 0, extraBlock * sizeof(double));
    326           }
    327           useCost = origcost;
    328           useOffset = saveOffset;
    329      } else {
    330           useCostExtra = costExtra;
    331           useCost = cost;
    332           useOffset = offset;
    333      }
    334      ncolx = ncols + extraBlock;
    335      for (i = 0; i < HISTORY + 1; i++) {
    336           history[i] = new double[ncolx];
    337      }
    338      for (i = 0; i < DJTEST; i++) {
    339           djSave[i] = 1.0e30;
    340      }
     296            rowupper[i] = 0.9e10;
     297            saveExtra[extraBlock] = rowupper[i];
     298            smaller = rowupper[i];
     299            value = 1.0;
     300          }
     301        }
     302        difference = rowupper[i] - rowlower[i];
     303        difference = CoinMin(difference, 1.0e31);
     304        rowupper[i] = smaller;
     305        elemExtra[extraBlock] = value;
     306        solExtra[extraBlock] = (rowupper[i] - rowsol[i]) / value;
     307        if (solExtra[extraBlock] < 0.0)
     308          solExtra[extraBlock] = 0.0;
     309        if (solExtra[extraBlock] > difference)
     310          solExtra[extraBlock] = difference;
     311        costExtra[extraBlock] = lambda[i] * value;
     312        upperExtra[extraBlock] = difference;
     313        rowsol[i] += value * solExtra[extraBlock];
     314        rowExtra[extraBlock++] = i;
     315      }
     316    }
     317    if (nbad)
     318      COIN_DETAIL_PRINT(printf("%d bad values - results may be wrong\n", nbad));
     319  }
     320  for (i = 0; i < nrows; i++) {
     321    offset += lambda[i] * rowsol[i];
     322  }
     323  if ((strategy & 256) != 0) {
     324    /* save best solution */
     325    saveSol = new double[ncols];
     326    CoinMemcpyN(colsol, ncols, saveSol);
     327    if (extraBlock) {
     328      useCostExtra = new double[extraBlock];
     329      memset(useCostExtra, 0, extraBlock * sizeof(double));
     330    }
     331    useCost = origcost;
     332    useOffset = saveOffset;
     333  } else {
     334    useCostExtra = costExtra;
     335    useCost = cost;
     336    useOffset = offset;
     337  }
     338  ncolx = ncols + extraBlock;
     339  for (i = 0; i < HISTORY + 1; i++) {
     340    history[i] = new double[ncolx];
     341  }
     342  for (i = 0; i < DJTEST; i++) {
     343    djSave[i] = 1.0e30;
     344  }
    341345#ifndef OSI_IDIOT
    342      int numberColumns = model_->numberColumns();
    343      for (int i=0;i<numberColumns;i++) {
    344        if (model_->getColumnStatus(i)!=ClpSimplex::isFixed)
    345          statusSave[i] = 0;
    346        else
    347          statusSave[i] = 2;
    348      }
    349      memset(statusSave+numberColumns,0,ncols-numberColumns);
    350      if ((strategy_&131072)==0) {
    351        for (int i=0;i<numberColumns;i++) {
    352          if (model_->getColumnStatus(i)==ClpSimplex::isFixed) {
    353            assert (colsol[i]<lower[i]+tolerance||
    354                    colsol[i]>upper[i]-tolerance);
    355          }
    356        }
    357      }
     346  int numberColumns = model_->numberColumns();
     347  for (int i = 0; i < numberColumns; i++) {
     348    if (model_->getColumnStatus(i) != ClpSimplex::isFixed)
     349      statusSave[i] = 0;
     350    else
     351      statusSave[i] = 2;
     352  }
     353  memset(statusSave + numberColumns, 0, ncols - numberColumns);
     354  if ((strategy_ & 131072) == 0) {
     355    for (int i = 0; i < numberColumns; i++) {
     356      if (model_->getColumnStatus(i) == ClpSimplex::isFixed) {
     357        assert(colsol[i] < lower[i] + tolerance || colsol[i] > upper[i] - tolerance);
     358      }
     359    }
     360  }
    358361#else
    359      for (i = 0; i < ncols; i++) {
    360           if (upper[i] - lower[i]) {
    361                statusSave[i] = 0;
    362           } else {
    363                statusSave[i] = 1;
    364           }
    365      }
    366 #endif
    367      // for two pass method
    368      int start[2];
    369      int stop[2];
    370      int direction = -1;
    371      start[0] = 0;
    372      stop[0] = ncols;
    373      start[1] = 0;
    374      stop[1] = 0;
    375      iter = 0;
    376      for (; iter < maxIts; iter++) {
    377           double sum1 = 0.0, sum2 = 0.0;
    378           double lastObj = 1.0e70;
    379           int good = 0, doScale = 0;
    380           if (strategy & 16) {
    381                int ii = iter / EVERY + 1;
    382                ii = ii * EVERY;
    383                if (iter > ii - HISTORY * 2 && (iter & 1) == 0) {
    384                     double * COIN_RESTRICT x = history[HISTORY-1];
    385                     for (i = HISTORY - 1; i > 0; i--) {
    386                          history[i] = history[i-1];
    387                     }
    388                     history[0] = x;
    389                     CoinMemcpyN(colsol, ncols, history[0]);
    390                     CoinMemcpyN(solExtra, extraBlock, history[0] + ncols);
    391                }
    392           }
    393           if ((iter % SAVEHISTORY) == 0 || doFull) {
    394                if ((strategy & 16) == 0) {
    395                     double * COIN_RESTRICT x = history[HISTORY-1];
    396                     for (i = HISTORY - 1; i > 0; i--) {
    397                          history[i] = history[i-1];
    398                     }
    399                     history[0] = x;
    400                     CoinMemcpyN(colsol, ncols, history[0]);
    401                     CoinMemcpyN(solExtra, extraBlock, history[0] + ncols);
    402                }
    403           }
    404           /* start full try */
    405           if ((iter % EVERY) == 0 || doFull) {
    406                // for next pass
    407                direction = - direction;
    408                // randomize.
    409                // The cast is to avoid gcc compiler warning
    410                int kcol = static_cast<int>(ncols * randomNumberGenerator->randomDouble());
    411                if (kcol == ncols)
    412                     kcol = ncols - 1;
    413                if (direction > 0) {
    414                     start[0] = kcol;
    415                     stop[0] = ncols;
    416                     start[1] = 0;
    417                     stop[1] = kcol;
     362  for (i = 0; i < ncols; i++) {
     363    if (upper[i] - lower[i]) {
     364      statusSave[i] = 0;
     365    } else {
     366      statusSave[i] = 1;
     367    }
     368  }
     369#endif
     370  // for two pass method
     371  int start[2];
     372  int stop[2];
     373  int direction = -1;
     374  start[0] = 0;
     375  stop[0] = ncols;
     376  start[1] = 0;
     377  stop[1] = 0;
     378  iter = 0;
     379  for (; iter < maxIts; iter++) {
     380    double sum1 = 0.0, sum2 = 0.0;
     381    double lastObj = 1.0e70;
     382    int good = 0, doScale = 0;
     383    if (strategy & 16) {
     384      int ii = iter / EVERY + 1;
     385      ii = ii * EVERY;
     386      if (iter > ii - HISTORY * 2 && (iter & 1) == 0) {
     387        double *COIN_RESTRICT x = history[HISTORY - 1];
     388        for (i = HISTORY - 1; i > 0; i--) {
     389          history[i] = history[i - 1];
     390        }
     391        history[0] = x;
     392        CoinMemcpyN(colsol, ncols, history[0]);
     393        CoinMemcpyN(solExtra, extraBlock, history[0] + ncols);
     394      }
     395    }
     396    if ((iter % SAVEHISTORY) == 0 || doFull) {
     397      if ((strategy & 16) == 0) {
     398        double *COIN_RESTRICT x = history[HISTORY - 1];
     399        for (i = HISTORY - 1; i > 0; i--) {
     400          history[i] = history[i - 1];
     401        }
     402        history[0] = x;
     403        CoinMemcpyN(colsol, ncols, history[0]);
     404        CoinMemcpyN(solExtra, extraBlock, history[0] + ncols);
     405      }
     406    }
     407    /* start full try */
     408    if ((iter % EVERY) == 0 || doFull) {
     409      // for next pass
     410      direction = -direction;
     411      // randomize.
     412      // The cast is to avoid gcc compiler warning
     413      int kcol = static_cast< int >(ncols * randomNumberGenerator->randomDouble());
     414      if (kcol == ncols)
     415        kcol = ncols - 1;
     416      if (direction > 0) {
     417        start[0] = kcol;
     418        stop[0] = ncols;
     419        start[1] = 0;
     420        stop[1] = kcol;
    418421#ifdef FOUR_GOES
    419                     for (int itry=0;itry<2;itry++) {
    420                       int chunk=(stop[itry]-start[itry]+FOUR_GOES-1)/FOUR_GOES;
    421                       startsX[itry][0]=start[itry];
    422                       for (int i=1;i<5;i++)
    423                         startsX[itry][i]=CoinMin(stop[itry],startsX[itry][i-1]+chunk);
    424                     }
    425 #endif
    426                } else {
    427                     start[0] = kcol;
    428                     stop[0] = -1;
    429                     start[1] = ncols - 1;
    430                     stop[1] = kcol;
     422        for (int itry = 0; itry < 2; itry++) {
     423          int chunk = (stop[itry] - start[itry] + FOUR_GOES - 1) / FOUR_GOES;
     424          startsX[itry][0] = start[itry];
     425          for (int i = 1; i < 5; i++)
     426            startsX[itry][i] = CoinMin(stop[itry], startsX[itry][i - 1] + chunk);
     427        }
     428#endif
     429      } else {
     430        start[0] = kcol;
     431        stop[0] = -1;
     432        start[1] = ncols - 1;
     433        stop[1] = kcol;
    431434#ifdef FOUR_GOES
    432                     for (int itry=0;itry<2;itry++) {
    433                       int chunk=(start[itry]-stop[itry]+FOUR_GOES-1)/FOUR_GOES;
    434                       startsX[itry][0]=start[itry];
    435                       for (int i=1;i<5;i++)
    436                         startsX[itry][i]=CoinMax(stop[itry],startsX[itry][i-1]-chunk);
    437                     }
    438 #endif
    439                }
    440                int itry = 0;
    441                /*if ((strategy&16)==0) {
     435        for (int itry = 0; itry < 2; itry++) {
     436          int chunk = (start[itry] - stop[itry] + FOUR_GOES - 1) / FOUR_GOES;
     437          startsX[itry][0] = start[itry];
     438          for (int i = 1; i < 5; i++)
     439            startsX[itry][i] = CoinMax(stop[itry], startsX[itry][i - 1] - chunk);
     440        }
     441#endif
     442      }
     443      int itry = 0;
     444      /*if ((strategy&16)==0) {
    442445                double * COIN_RESTRICT x=history[HISTORY-1];
    443446                for (i=HISTORY-1;i>0;i--) {
     
    448451                 CoinMemcpyN(solExtra,extraBlock,history[0]+ncols);
    449452                }*/
    450                while (!good) {
    451                     itry++;
     453      while (!good) {
     454        itry++;
    452455#define MAXTRY 5
    453                     if (iter > after && doScale < 2 && itry < MAXTRY) {
    454                          /* now full one */
    455                          for (i = 0; i < nrows; i++) {
    456                               rowsol[i] = -rowupper[i];
    457                          }
    458                          sum2 = 0.0;
    459                          objvalue = 0.0;
    460                          memset(pi, 0, nrows * sizeof(double));
    461                          {
    462                               double * COIN_RESTRICT theta = thetaX;
    463                               double * COIN_RESTRICT dj = djX;
    464                               double * COIN_RESTRICT b = bX;
    465                               double ** a = aX;
    466                               double ** awork = aworkX;
    467                               double * COIN_RESTRICT v = vX;
    468                               double c;
     456        if (iter > after && doScale < 2 && itry < MAXTRY) {
     457          /* now full one */
     458          for (i = 0; i < nrows; i++) {
     459            rowsol[i] = -rowupper[i];
     460          }
     461          sum2 = 0.0;
     462          objvalue = 0.0;
     463          memset(pi, 0, nrows * sizeof(double));
     464          {
     465            double *COIN_RESTRICT theta = thetaX;
     466            double *COIN_RESTRICT dj = djX;
     467            double *COIN_RESTRICT b = bX;
     468            double **a = aX;
     469            double **awork = aworkX;
     470            double *COIN_RESTRICT v = vX;
     471            double c;
    469472#ifdef FIT
    470                               int ntot = 0, nsign = 0, ngood = 0, mgood[4] = {0, 0, 0, 0};
    471                               double diff1, diff2, val0, val1, val2, newValue;
    472                               CoinMemcpyN(colsol, ncols, history[HISTORY-1]);
    473                               CoinMemcpyN(solExtra, extraBlock, history[HISTORY-1] + ncols);
    474 #endif
    475                               dj[0] = 0.0;
    476                               for (i = 1; i < nsolve; i++) {
    477                                    dj[i] = 0.0;
    478                                    memset(allsum[i], 0, nrows * sizeof(double));
    479                               }
    480                               for (i = 0; i < ncols; i++) {
    481                                    double value2 = colsol[i];
    482                                    if (value2 > lower[i] + tolerance) {
    483                                         if(value2 < (upper[i] - tolerance)) {
    484                                              int k;
    485                                              objvalue += value2 * cost[i];
     473            int ntot = 0, nsign = 0, ngood = 0, mgood[4] = { 0, 0, 0, 0 };
     474            double diff1, diff2, val0, val1, val2, newValue;
     475            CoinMemcpyN(colsol, ncols, history[HISTORY - 1]);
     476            CoinMemcpyN(solExtra, extraBlock, history[HISTORY - 1] + ncols);
     477#endif
     478            dj[0] = 0.0;
     479            for (i = 1; i < nsolve; i++) {
     480              dj[i] = 0.0;
     481              memset(allsum[i], 0, nrows * sizeof(double));
     482            }
     483            for (i = 0; i < ncols; i++) {
     484              double value2 = colsol[i];
     485              if (value2 > lower[i] + tolerance) {
     486                if (value2 < (upper[i] - tolerance)) {
     487                  int k;
     488                  objvalue += value2 * cost[i];
    486489#ifdef FIT
    487                                              ntot++;
    488                                              val0 = history[0][i];
    489                                              val1 = history[1][i];
    490                                              val2 = history[2][i];
    491                                              diff1 = val0 - val1;
    492                                              diff2 = val1 - val2;
    493                                              if (diff1*diff2 >= 0.0) {
    494                                                   nsign++;
    495                                                   if (fabs(diff1) < fabs(diff2)) {
    496                                                        int ii = static_cast<int>(fabs(4.0 * diff1 / diff2));
    497                                                        if (ii == 4) ii = 3;
    498                                                        mgood[ii]++;
    499                                                        ngood++;
    500                                                   }
    501                                                   if (fabs(diff1) < 0.75 * fabs(diff2)) {
    502                                                        newValue = val1 + (diff1 * diff2) / (diff2 - diff1);
    503                                                   } else {
    504                                                        newValue = val1 + 4.0 * diff1;
    505                                                   }
    506                                              } else {
    507                                                   newValue = 0.333333333 * (val0 + val1 + val2);
    508                                              }
    509                                              if (newValue > upper[i] - tolerance) {
    510                                                   newValue = upper[i];
    511                                              } else if (newValue < lower[i] + tolerance) {
    512                                                   newValue = lower[i];
    513                                              }
    514                                              history[HISTORY-1][i] = newValue;
    515 #endif
    516                                              for (k = 0; k < HISTORY - 1; k++) {
    517                                                   value = history[k][i] - history[k+1][i];
    518                                                   dj[k] += value * cost[i];
    519                                                   v[k] = value;
    520                                              }
    521                                              if (elemnt) {
    522                                                   for (j = columnStart[i]; j < columnStart[i] + length[i]; j++) {
    523                                                        int irow = row[j];
    524                                                        for (k = 0; k < HISTORY - 1; k++) {
    525                                                             allsum[k][irow] += elemnt[j] * v[k];
    526                                                        }
    527                                                        rowsol[irow] += elemnt[j] * value2;
    528                                                   }
    529                                              } else {
    530                                                   for (j = columnStart[i]; j < columnStart[i] + length[i]; j++) {
    531                                                        int irow = row[j];
    532                                                        for (k = 0; k < HISTORY - 1; k++) {
    533                                                             allsum[k][irow] += v[k];
    534                                                        }
    535                                                        rowsol[irow] += value2;
    536                                                   }
    537                                              }
    538                                         } else {
    539                                              /* at ub */
    540                                              colsol[i] = upper[i];
    541                                              value2 = colsol[i];
    542                                              objvalue += value2 * cost[i];
    543                                              if (elemnt) {
    544                                                   for (j = columnStart[i]; j < columnStart[i] + length[i]; j++) {
    545                                                        int irow = row[j];
    546                                                        rowsol[irow] += elemnt[j] * value2;
    547                                                   }
    548                                              } else {
    549                                                   for (j = columnStart[i]; j < columnStart[i] + length[i]; j++) {
    550                                                        int irow = row[j];
    551                                                        rowsol[irow] += value2;
    552                                                   }
    553                                              }
    554                                         }
    555                                    } else {
    556                                         /* at lb */
    557                                         if (value2) {
    558                                              objvalue += value2 * cost[i];
    559                                              if (elemnt) {
    560                                                   for (j = columnStart[i]; j < columnStart[i] + length[i]; j++) {
    561                                                        int irow = row[j];
    562                                                        rowsol[irow] += elemnt[j] * value2;
    563                                                   }
    564                                              } else {
    565                                                   for (j = columnStart[i]; j < columnStart[i] + length[i]; j++) {
    566                                                        int irow = row[j];
    567                                                        rowsol[irow] += value2;
    568                                                   }
    569                                              }
    570                                         }
    571                                    }
    572                               }
     490                  ntot++;
     491                  val0 = history[0][i];
     492                  val1 = history[1][i];
     493                  val2 = history[2][i];
     494                  diff1 = val0 - val1;
     495                  diff2 = val1 - val2;
     496                  if (diff1 * diff2 >= 0.0) {
     497                    nsign++;
     498                    if (fabs(diff1) < fabs(diff2)) {
     499                      int ii = static_cast< int >(fabs(4.0 * diff1 / diff2));
     500                      if (ii == 4)
     501                        ii = 3;
     502                      mgood[ii]++;
     503                      ngood++;
     504                    }
     505                    if (fabs(diff1) < 0.75 * fabs(diff2)) {
     506                      newValue = val1 + (diff1 * diff2) / (diff2 - diff1);
     507                    } else {
     508                      newValue = val1 + 4.0 * diff1;
     509                    }
     510                  } else {
     511                    newValue = 0.333333333 * (val0 + val1 + val2);
     512                  }
     513                  if (newValue > upper[i] - tolerance) {
     514                    newValue = upper[i];
     515                  } else if (newValue < lower[i] + tolerance) {
     516                    newValue = lower[i];
     517                  }
     518                  history[HISTORY - 1][i] = newValue;
     519#endif
     520                  for (k = 0; k < HISTORY - 1; k++) {
     521                    value = history[k][i] - history[k + 1][i];
     522                    dj[k] += value * cost[i];
     523                    v[k] = value;
     524                  }
     525                  if (elemnt) {
     526                    for (j = columnStart[i]; j < columnStart[i] + length[i]; j++) {
     527                      int irow = row[j];
     528                      for (k = 0; k < HISTORY - 1; k++) {
     529                        allsum[k][irow] += elemnt[j] * v[k];
     530                      }
     531                      rowsol[irow] += elemnt[j] * value2;
     532                    }
     533                  } else {
     534                    for (j = columnStart[i]; j < columnStart[i] + length[i]; j++) {
     535                      int irow = row[j];
     536                      for (k = 0; k < HISTORY - 1; k++) {
     537                        allsum[k][irow] += v[k];
     538                      }
     539                      rowsol[irow] += value2;
     540                    }
     541                  }
     542                } else {
     543                  /* at ub */
     544                  colsol[i] = upper[i];
     545                  value2 = colsol[i];
     546                  objvalue += value2 * cost[i];
     547                  if (elemnt) {
     548                    for (j = columnStart[i]; j < columnStart[i] + length[i]; j++) {
     549                      int irow = row[j];
     550                      rowsol[irow] += elemnt[j] * value2;
     551                    }
     552                  } else {
     553                    for (j = columnStart[i]; j < columnStart[i] + length[i]; j++) {
     554                      int irow = row[j];
     555                      rowsol[irow] += value2;
     556                    }
     557                  }
     558                }
     559              } else {
     560                /* at lb */
     561                if (value2) {
     562                  objvalue += value2 * cost[i];
     563                  if (elemnt) {
     564                    for (j = columnStart[i]; j < columnStart[i] + length[i]; j++) {
     565                      int irow = row[j];
     566                      rowsol[irow] += elemnt[j] * value2;
     567                    }
     568                  } else {
     569                    for (j = columnStart[i]; j < columnStart[i] + length[i]; j++) {
     570                      int irow = row[j];
     571                      rowsol[irow] += value2;
     572                    }
     573                  }
     574                }
     575              }
     576            }
    573577#ifdef FIT
    574                               /*printf("total %d, same sign %d, good %d %d %d %d %d\n",
     578            /*printf("total %d, same sign %d, good %d %d %d %d %d\n",
    575579                                ntot,nsign,ngood,mgood[0],mgood[1],mgood[2],mgood[3]);*/
    576580#endif
    577                               if (extraBlock) {
    578                                    for (i = 0; i < extraBlock; i++) {
    579                                         double element = elemExtra[i];
    580                                         int irow = rowExtra[i];
    581                                         objvalue += solExtra[i] * costExtra[i];
    582                                         if (solExtra[i] > tolerance
    583                                                   && solExtra[i] < (upperExtra[i] - tolerance)) {
    584                                              double value2 = solExtra[i];
    585                                              int k;
    586                                              for (k = 0; k < HISTORY - 1; k++) {
    587                                                   value = history[k][i+ncols] - history[k+1][i+ncols];
    588                                                   dj[k] += value * costExtra[i];
    589                                                   allsum[k][irow] += element * value;
    590                                              }
    591                                              rowsol[irow] += element * value2;
    592                                         } else {
    593                                              double value2 = solExtra[i];
    594                                              double element = elemExtra[i];
    595                                              int irow = rowExtra[i];
    596                                              rowsol[irow] += element * value2;
    597                                         }
    598                                    }
    599                               }
     581            if (extraBlock) {
     582              for (i = 0; i < extraBlock; i++) {
     583                double element = elemExtra[i];
     584                int irow = rowExtra[i];
     585                objvalue += solExtra[i] * costExtra[i];
     586                if (solExtra[i] > tolerance
     587                  && solExtra[i] < (upperExtra[i] - tolerance)) {
     588                  double value2 = solExtra[i];
     589                  int k;
     590                  for (k = 0; k < HISTORY - 1; k++) {
     591                    value = history[k][i + ncols] - history[k + 1][i + ncols];
     592                    dj[k] += value * costExtra[i];
     593                    allsum[k][irow] += element * value;
     594                  }
     595                  rowsol[irow] += element * value2;
     596                } else {
     597                  double value2 = solExtra[i];
     598                  double element = elemExtra[i];
     599                  int irow = rowExtra[i];
     600                  rowsol[irow] += element * value2;
     601                }
     602              }
     603            }
    600604#ifdef NULLVECTOR
    601                               if ((strategy & 64)) {
    602                                    double djVal = dj[0];
    603                                    for (i = 0; i < ncols - nrows; i++) {
    604                                         double value2 = colsol[i];
    605                                         if (value2 > lower[i] + tolerance && value2 < upper[i] - tolerance) {
    606                                              value = history[0][i] - history[1][i];
    607                                         } else {
    608                                              value = 0.0;
    609                                         }
    610                                         history[HISTORY][i] = value;
    611                                    }
    612                                    for (; i < ncols; i++) {
    613                                         double value2 = colsol[i];
    614                                         double delta;
    615                                         int irow = i - (ncols - nrows);
    616                                         double oldSum = allsum[0][irow];;
    617                                         if (value2 > lower[i] + tolerance && value2 < upper[i] - tolerance) {
    618                                              delta = history[0][i] - history[1][i];
    619                                         } else {
    620                                              delta = 0.0;
    621                                         }
    622                                         djVal -= delta * cost[i];
    623                                         oldSum -= delta;
    624                                         delta = - oldSum;
    625                                         djVal += delta * cost[i];
    626                                         history[HISTORY][i] = delta;
    627                                    }
    628                                    dj[HISTORY-1] = djVal;
    629                                    djVal = 0.0;
    630                                    for (i = 0; i < ncols; i++) {
    631                                         double value2 = colsol[i];
    632                                         if (value2 > lower[i] + tolerance && value2 < upper[i] - tolerance ||
    633                                                   i >= ncols - nrows) {
    634                                              int k;
    635                                              value = history[HISTORY][i];
    636                                              djVal += value * cost[i];
    637                                              for (j = columnStart[i]; j < columnStart[i] + length[i]; j++) {
    638                                                   int irow = row[j];
    639                                                   allsum[nsolve-1][irow] += value;
    640                                              }
    641                                         }
    642                                    }
    643                                    printf("djs %g %g\n", dj[HISTORY-1], djVal);
    644                               }
    645 #endif
    646                               for (i = 0; i < nsolve; i++) {
    647                                    int j;
    648                                    b[i] = 0.0;
    649                                    for (j = 0; j < nsolve; j++) {
    650                                         a[i][j] = 0.0;
    651                                    }
    652                               }
    653                               c = 0.0;
    654                               for (i = 0; i < nrows; i++) {
    655                                    double value = rowsol[i];
    656                                    for (k = 0; k < nsolve; k++) {
    657                                         v[k] = allsum[k][i];
    658                                         b[k] += v[k] * value;
    659                                    }
    660                                    c += value * value;
    661                                    for (k = 0; k < nsolve; k++) {
    662                                         for (j = k; j < nsolve; j++) {
    663                                              a[k][j] += v[k] * v[j];
    664                                         }
    665                                    }
    666                               }
    667                               sum2 = c;
    668                               if (itry == 1) {
    669                                    lastObj = objvalue + weight * sum2;
    670                               }
    671                               for (k = 0; k < nsolve; k++) {
    672                                    b[k] = - (weight * b[k] + 0.5 * dj[k]);
    673                                    for (j = k; j < nsolve; j++) {
    674                                         a[k][j] *= weight;
    675                                         a[j][k] = a[k][j];
    676                                    }
    677                               }
    678                               c *= weight;
    679                               for (k = 0; k < nsolve; k++) {
    680                                    theta[k] = b[k];
    681                               }
    682                               solveSmall(nsolve, a, awork, theta);
    683                               if ((strategy & 64) != 0) {
    684                                    value = 10.0;
    685                                    for (k = 0; k < nsolve; k++) {
    686                                         value = CoinMax(value, fabs(theta[k]));
    687                                    }
    688                                    if (value > 10.0 && ((logLevel_ & 4) != 0)) {
    689                                         printf("theta %g %g %g\n", theta[0], theta[1], theta[2]);
    690                                    }
    691                                    value = 10.0 / value;
    692                                    for (k = 0; k < nsolve; k++) {
    693                                         theta[k] *= value;
    694                                    }
    695                               }
    696                               for (i = 0; i < ncolx; i++) {
    697                                    double valueh = 0.0;
    698                                    for (k = 0; k < HISTORY - 1; k++) {
    699                                         value = history[k][i] - history[k+1][i];
    700                                         valueh += value * theta[k];
    701                                    }
     605            if ((strategy & 64)) {
     606              double djVal = dj[0];
     607              for (i = 0; i < ncols - nrows; i++) {
     608                double value2 = colsol[i];
     609                if (value2 > lower[i] + tolerance && value2 < upper[i] - tolerance) {
     610                  value = history[0][i] - history[1][i];
     611                } else {
     612                  value = 0.0;
     613                }
     614                history[HISTORY][i] = value;
     615              }
     616              for (; i < ncols; i++) {
     617                double value2 = colsol[i];
     618                double delta;
     619                int irow = i - (ncols - nrows);
     620                double oldSum = allsum[0][irow];
     621                ;
     622                if (value2 > lower[i] + tolerance && value2 < upper[i] - tolerance) {
     623                  delta = history[0][i] - history[1][i];
     624                } else {
     625                  delta = 0.0;
     626                }
     627                djVal -= delta * cost[i];
     628                oldSum -= delta;
     629                delta = -oldSum;
     630                djVal += delta * cost[i];
     631                history[HISTORY][i] = delta;
     632              }
     633              dj[HISTORY - 1] = djVal;
     634              djVal = 0.0;
     635              for (i = 0; i < ncols; i++) {
     636                double value2 = colsol[i];
     637                if (value2 > lower[i] + tolerance && value2 < upper[i] - tolerance || i >= ncols - nrows) {
     638                  int k;
     639                  value = history[HISTORY][i];
     640                  djVal += value * cost[i];
     641                  for (j = columnStart[i]; j < columnStart[i] + length[i]; j++) {
     642                    int irow = row[j];
     643                    allsum[nsolve - 1][irow] += value;
     644                  }
     645                }
     646              }
     647              printf("djs %g %g\n", dj[HISTORY - 1], djVal);
     648            }
     649#endif
     650            for (i = 0; i < nsolve; i++) {
     651              int j;
     652              b[i] = 0.0;
     653              for (j = 0; j < nsolve; j++) {
     654                a[i][j] = 0.0;
     655              }
     656            }
     657            c = 0.0;
     658            for (i = 0; i < nrows; i++) {
     659              double value = rowsol[i];
     660              for (k = 0; k < nsolve; k++) {
     661                v[k] = allsum[k][i];
     662                b[k] += v[k] * value;
     663              }
     664              c += value * value;
     665              for (k = 0; k < nsolve; k++) {
     666                for (j = k; j < nsolve; j++) {
     667                  a[k][j] += v[k] * v[j];
     668                }
     669              }
     670            }
     671            sum2 = c;
     672            if (itry == 1) {
     673              lastObj = objvalue + weight * sum2;
     674            }
     675            for (k = 0; k < nsolve; k++) {
     676              b[k] = -(weight * b[k] + 0.5 * dj[k]);
     677              for (j = k; j < nsolve; j++) {
     678                a[k][j] *= weight;
     679                a[j][k] = a[k][j];
     680              }
     681            }
     682            c *= weight;
     683            for (k = 0; k < nsolve; k++) {
     684              theta[k] = b[k];
     685            }
     686            solveSmall(nsolve, a, awork, theta);
     687            if ((strategy & 64) != 0) {
     688              value = 10.0;
     689              for (k = 0; k < nsolve; k++) {
     690                value = CoinMax(value, fabs(theta[k]));
     691              }
     692              if (value > 10.0 && ((logLevel_ & 4) != 0)) {
     693                printf("theta %g %g %g\n", theta[0], theta[1], theta[2]);
     694              }
     695              value = 10.0 / value;
     696              for (k = 0; k < nsolve; k++) {
     697                theta[k] *= value;
     698              }
     699            }
     700            for (i = 0; i < ncolx; i++) {
     701              double valueh = 0.0;
     702              for (k = 0; k < HISTORY - 1; k++) {
     703                value = history[k][i] - history[k + 1][i];
     704                valueh += value * theta[k];
     705              }
    702706#ifdef NULLVECTOR
    703                                    value = history[HISTORY][i];
    704                                    valueh += value * theta[HISTORY-1];
    705 #endif
    706                                    history[HISTORY][i] = valueh;
    707                               }
    708                          }
     707              value = history[HISTORY][i];
     708              valueh += value * theta[HISTORY - 1];
     709#endif
     710              history[HISTORY][i] = valueh;
     711            }
     712          }
    709713#ifdef NULLVECTOR
    710                          if ((strategy & 64)) {
    711                               for (i = 0; i < ncols - nrows; i++) {
    712                                    if (colsol[i] <= lower[i] + tolerance
    713                                              || colsol[i] >= (upper[i] - tolerance)) {
    714                                         history[HISTORY][i] = 0.0;;
    715                                    }
    716                               }
    717                               tolerance = -tolerance; /* switch off test */
    718                          }
    719 #endif
    720                          if (!doScale) {
    721                               for (i = 0; i < ncols; i++) {
    722                                    if (colsol[i] > lower[i] + tolerance
    723                                              && colsol[i] < (upper[i] - tolerance)) {
    724                                         value = history[HISTORY][i];
    725                                         colsol[i] += value;
    726                                         if (colsol[i] < lower[i] + tolerance) {
    727                                              colsol[i] = lower[i];
    728                                         } else if (colsol[i] > upper[i] - tolerance) {
    729                                              colsol[i] = upper[i];
    730                                         }
    731                                    }
    732                               }
    733                               if (extraBlock) {
    734                                    for (i = 0; i < extraBlock; i++) {
    735                                         if (solExtra[i] > tolerance
    736                                                   && solExtra[i] < (upperExtra[i] - tolerance)) {
    737                                              value = history[HISTORY][i+ncols];
    738                                              solExtra[i] += value;
    739                                              if (solExtra[i] < 0.0) {
    740                                                   solExtra[i] = 0.0;
    741                                              } else if (solExtra[i] > upperExtra[i]) {
    742                                                   solExtra[i] = upperExtra[i];
    743                                              }
    744                                         }
    745                                    }
    746                               }
    747                          } else {
    748                               double theta = 1.0;
    749                               double saveTheta = theta;
    750                               for (i = 0; i < ncols; i++) {
    751                                    if (colsol[i] > lower[i] + tolerance
    752                                              && colsol[i] < (upper[i] - tolerance)) {
    753                                         value = history[HISTORY][i];
    754                                         if (value > 0) {
    755                                              if (theta * value + colsol[i] > upper[i]) {
    756                                                   theta = (upper[i] - colsol[i]) / value;
    757                                              }
    758                                         } else if (value < 0) {
    759                                              if (colsol[i] + theta * value < lower[i]) {
    760                                                   theta = (lower[i] - colsol[i]) / value;
    761                                              }
    762                                         }
    763                                    }
    764                               }
    765                               if (extraBlock) {
    766                                    for (i = 0; i < extraBlock; i++) {
    767                                         if (solExtra[i] > tolerance
    768                                                   && solExtra[i] < (upperExtra[i] - tolerance)) {
    769                                              value = history[HISTORY][i+ncols];
    770                                              if (value > 0) {
    771                                                   if (theta * value + solExtra[i] > upperExtra[i]) {
    772                                                        theta = (upperExtra[i] - solExtra[i]) / value;
    773                                                   }
    774                                              } else if (value < 0) {
    775                                                   if (solExtra[i] + theta * value < 0.0) {
    776                                                        theta = -solExtra[i] / value;
    777                                                   }
    778                                              }
    779                                         }
    780                                    }
    781                               }
    782                               if ((iter % 100 == 0) && (logLevel_ & 8) != 0) {
    783                                    if (theta < saveTheta) {
    784                                         printf(" - modified theta %g\n", theta);
    785                                    }
    786                               }
    787                               for (i = 0; i < ncols; i++) {
    788                                    if (colsol[i] > lower[i] + tolerance
    789                                              && colsol[i] < (upper[i] - tolerance)) {
    790                                         value = history[HISTORY][i];
    791                                         colsol[i] += value * theta;
    792                                    }
    793                               }
    794                               if (extraBlock) {
    795                                    for (i = 0; i < extraBlock; i++) {
    796                                         if (solExtra[i] > tolerance
    797                                                   && solExtra[i] < (upperExtra[i] - tolerance)) {
    798                                              value = history[HISTORY][i+ncols];
    799                                              solExtra[i] += value * theta;
    800                                         }
    801                                    }
    802                               }
    803                          }
     714          if ((strategy & 64)) {
     715            for (i = 0; i < ncols - nrows; i++) {
     716              if (colsol[i] <= lower[i] + tolerance
     717                || colsol[i] >= (upper[i] - tolerance)) {
     718                history[HISTORY][i] = 0.0;
     719                ;
     720              }
     721            }
     722            tolerance = -tolerance; /* switch off test */
     723          }
     724#endif
     725          if (!doScale) {
     726            for (i = 0; i < ncols; i++) {
     727              if (colsol[i] > lower[i] + tolerance
     728                && colsol[i] < (upper[i] - tolerance)) {
     729                value = history[HISTORY][i];
     730                colsol[i] += value;
     731                if (colsol[i] < lower[i] + tolerance) {
     732                  colsol[i] = lower[i];
     733                } else if (colsol[i] > upper[i] - tolerance) {
     734                  colsol[i] = upper[i];
     735                }
     736              }
     737            }
     738            if (extraBlock) {
     739              for (i = 0; i < extraBlock; i++) {
     740                if (solExtra[i] > tolerance
     741                  && solExtra[i] < (upperExtra[i] - tolerance)) {
     742                  value = history[HISTORY][i + ncols];
     743                  solExtra[i] += value;
     744                  if (solExtra[i] < 0.0) {
     745                    solExtra[i] = 0.0;
     746                  } else if (solExtra[i] > upperExtra[i]) {
     747                    solExtra[i] = upperExtra[i];
     748                  }
     749                }
     750              }
     751            }
     752          } else {
     753            double theta = 1.0;
     754            double saveTheta = theta;
     755            for (i = 0; i < ncols; i++) {
     756              if (colsol[i] > lower[i] + tolerance
     757                && colsol[i] < (upper[i] - tolerance)) {
     758                value = history[HISTORY][i];
     759                if (value > 0) {
     760                  if (theta * value + colsol[i] > upper[i]) {
     761                    theta = (upper[i] - colsol[i]) / value;
     762                  }
     763                } else if (value < 0) {
     764                  if (colsol[i] + theta * value < lower[i]) {
     765                    theta = (lower[i] - colsol[i]) / value;
     766                  }
     767                }
     768              }
     769            }
     770            if (extraBlock) {
     771              for (i = 0; i < extraBlock; i++) {
     772                if (solExtra[i] > tolerance
     773                  && solExtra[i] < (upperExtra[i] - tolerance)) {
     774                  value = history[HISTORY][i + ncols];
     775                  if (value > 0) {
     776                    if (theta * value + solExtra[i] > upperExtra[i]) {
     777                      theta = (upperExtra[i] - solExtra[i]) / value;
     778                    }
     779                  } else if (value < 0) {
     780                    if (solExtra[i] + theta * value < 0.0) {
     781                      theta = -solExtra[i] / value;
     782                    }
     783                  }
     784                }
     785              }
     786            }
     787            if ((iter % 100 == 0) && (logLevel_ & 8) != 0) {
     788              if (theta < saveTheta) {
     789                printf(" - modified theta %g\n", theta);
     790              }
     791            }
     792            for (i = 0; i < ncols; i++) {
     793              if (colsol[i] > lower[i] + tolerance
     794                && colsol[i] < (upper[i] - tolerance)) {
     795                value = history[HISTORY][i];
     796                colsol[i] += value * theta;
     797              }
     798            }
     799            if (extraBlock) {
     800              for (i = 0; i < extraBlock; i++) {
     801                if (solExtra[i] > tolerance
     802                  && solExtra[i] < (upperExtra[i] - tolerance)) {
     803                  value = history[HISTORY][i + ncols];
     804                  solExtra[i] += value * theta;
     805                }
     806              }
     807            }
     808          }
    804809#ifdef NULLVECTOR
    805                          tolerance = fabs(tolerance); /* switch back on */
    806 #endif
    807                          if ((iter % 100) == 0 && (logLevel_ & 8) != 0) {
    808                               printf("\n");
    809                          }
    810                     }
    811                     good = 1;
    812                     result = objval(nrows, ncols, rowsol, colsol, pi, djs, useCost,
    813                                     rowlower, rowupper, lower, upper,
    814                                     elemnt, row, columnStart, length, extraBlock, rowExtra,
    815                                     solExtra, elemExtra, upperExtra, useCostExtra,
    816                                     weight);
    817                     weightedObj = result.weighted;
    818                     if (!iter) saveValue = weightedObj;
    819                     objvalue = result.objval;
    820                     sum1 = result.infeas;
    821                     if (saveSol) {
    822                          if (result.weighted < bestSol) {
    823                               COIN_DETAIL_PRINT(printf("%d %g better than %g\n", iter,
    824                                                        result.weighted * maxmin - useOffset, bestSol * maxmin - useOffset));
    825                               bestSol = result.weighted;
    826                               CoinMemcpyN(colsol, ncols, saveSol);
    827                          }
    828                     }
     810          tolerance = fabs(tolerance); /* switch back on */
     811#endif
     812          if ((iter % 100) == 0 && (logLevel_ & 8) != 0) {
     813            printf("\n");
     814          }
     815        }
     816        good = 1;
     817        result = objval(nrows, ncols, rowsol, colsol, pi, djs, useCost,
     818          rowlower, rowupper, lower, upper,
     819          elemnt, row, columnStart, length, extraBlock, rowExtra,
     820          solExtra, elemExtra, upperExtra, useCostExtra,
     821          weight);
     822        weightedObj = result.weighted;
     823        if (!iter)
     824          saveValue = weightedObj;
     825        objvalue = result.objval;
     826        sum1 = result.infeas;
     827        if (saveSol) {
     828          if (result.weighted < bestSol) {
     829            COIN_DETAIL_PRINT(printf("%d %g better than %g\n", iter,
     830              result.weighted * maxmin - useOffset, bestSol * maxmin - useOffset));
     831            bestSol = result.weighted;
     832            CoinMemcpyN(colsol, ncols, saveSol);
     833          }
     834        }
    829835#ifdef FITz
    830                     if (iter > after) {
    831                          IdiotResult result2;
    832                          double ww, oo, ss;
    833                          if (extraBlock) abort();
    834                          result2 = objval(nrows, ncols, row2, sol2, pi2, djs, cost,
    835                                           rowlower, rowupper, lower, upper,
    836                                           elemnt, row, columnStart, extraBlock, rowExtra,
    837                                           solExtra, elemExtra, upperExtra, costExtra,
    838                                           weight);
    839                          ww = result2.weighted;
    840                          oo = result2.objval;
    841                          ss = result2.infeas;
    842                          printf("wobj %g obj %g inf %g last %g\n", ww, oo, ss, lastObj);
    843                          if (ww < weightedObj && ww < lastObj) {
    844                               printf(" taken");
    845                               ntaken++;
    846                               saving += weightedObj - ww;
    847                               weightedObj = ww;
    848                               objvalue = oo;
    849                               sum1 = ss;
    850                               CoinMemcpyN(row2, nrows, rowsol);
    851                               CoinMemcpyN(pi2, nrows, pi);
    852                               CoinMemcpyN(sol2, ncols, colsol);
    853                               result = objval(nrows, ncols, rowsol, colsol, pi, djs, cost,
    854                                               rowlower, rowupper, lower, upper,
    855                                               elemnt, row, columnStart, extraBlock, rowExtra,
    856                                               solExtra, elemExtra, upperExtra, costExtra,
    857                                               weight);
    858                               weightedObj = result.weighted;
    859                               objvalue = result.objval;
    860                               sum1 = result.infeas;
    861                               if (ww < weightedObj) abort();
    862                          } else {
    863                               printf(" not taken");
    864                               nottaken++;
    865                          }
    866                     }
    867 #endif
    868                     /*printf("%d %g %g %g %g\n",itry,lastObj,weightedObj,objvalue,sum1);*/
    869                     if (weightedObj > lastObj + 1.0e-4 && itry < MAXTRY) {
    870                          if((logLevel_ & 16) != 0 && doScale) {
    871                               printf("Weighted objective from %g to %g **** bad move\n",
    872                                      lastObj, weightedObj);
    873                          }
    874                          if (doScale) {
    875                               good = 1;
    876                          }
    877                          if ((strategy & 3) == 1) {
    878                               good = 0;
    879                               if (weightedObj > lastObj + djExit) {
    880                                    if ((logLevel_ & 16) != 0) {
    881                                         printf("Weighted objective from %g to %g ?\n", lastObj, weightedObj);
    882                                    }
    883                                    CoinMemcpyN(history[0], ncols, colsol);
    884                                    CoinMemcpyN(history[0] + ncols, extraBlock, solExtra);
    885                                    good = 1;
    886                               }
    887                          } else if ((strategy & 3) == 2) {
    888                               if (weightedObj > lastObj + 0.1 * maxDj) {
    889                                    CoinMemcpyN(history[0], ncols, colsol);
    890                                    CoinMemcpyN(history[0] + ncols, extraBlock, solExtra);
    891                                    doScale++;
    892                                    good = 0;
    893                               }
    894                          } else if ((strategy & 3) == 3) {
    895                               if (weightedObj > lastObj + 0.001 * maxDj) {
    896                                    /*doScale++;*/
    897                                    good = 0;
    898                               }
    899                          }
    900                     }
    901                }
    902                if ((iter % checkFrequency_) == 0) {
    903                     double best = weightedObj;
    904                     double test = obj[0];
    905                     for (i = 1; i < DROP; i++) {
    906                          obj[i-1] = obj[i];
    907                          if (best > obj[i]) best = obj[i];
    908                     }
    909                     obj[DROP-1] = best;
    910                     if (test - best < drop && (strategy & 8) == 0) {
    911                          if ((logLevel_ & 8) != 0) {
    912                               printf("Exiting as drop in %d its is %g after %d iterations\n",
    913                                      DROP * checkFrequency_, test - best, iter);
    914                          }
    915                          goto RETURN;
    916                     }
    917                }
    918                if ((iter % logFreq_) == 0) {
    919                     double piSum = 0.0;
    920                     for (i = 0; i < nrows; i++) {
    921                          piSum += (rowsol[i] + rowupper[i]) * pi[i];
    922                     }
    923                     if ((logLevel_ & 2) != 0) {
    924                          printf("%d Infeas %g, obj %g - wtObj %g dual %g maxDj %g\n",
    925                                 iter, sum1, objvalue * maxmin - useOffset, weightedObj - useOffset,
    926                                 piSum * maxmin - useOffset, maxDj);
    927                     }
    928                }
    929                CoinMemcpyN(statusSave, ncols, statusWork);
    930                nflagged = 0;
    931           }
    932           nChange = 0;
    933           doFull = 0;
    934           maxDj = 0.0;
    935           // go through forwards or backwards and starting at odd places
     836        if (iter > after) {
     837          IdiotResult result2;
     838          double ww, oo, ss;
     839          if (extraBlock)
     840            abort();
     841          result2 = objval(nrows, ncols, row2, sol2, pi2, djs, cost,
     842            rowlower, rowupper, lower, upper,
     843            elemnt, row, columnStart, extraBlock, rowExtra,
     844            solExtra, elemExtra, upperExtra, costExtra,
     845            weight);
     846          ww = result2.weighted;
     847          oo = result2.objval;
     848          ss = result2.infeas;
     849          printf("wobj %g obj %g inf %g last %g\n", ww, oo, ss, lastObj);
     850          if (ww < weightedObj && ww < lastObj) {
     851            printf(" taken");
     852            ntaken++;
     853            saving += weightedObj - ww;
     854            weightedObj = ww;
     855            objvalue = oo;
     856            sum1 = ss;
     857            CoinMemcpyN(row2, nrows, rowsol);
     858            CoinMemcpyN(pi2, nrows, pi);
     859            CoinMemcpyN(sol2, ncols, colsol);
     860            result = objval(nrows, ncols, rowsol, colsol, pi, djs, cost,
     861              rowlower, rowupper, lower, upper,
     862              elemnt, row, columnStart, extraBlock, rowExtra,
     863              solExtra, elemExtra, upperExtra, costExtra,
     864              weight);
     865            weightedObj = result.weighted;
     866            objvalue = result.objval;
     867            sum1 = result.infeas;
     868            if (ww < weightedObj)
     869              abort();
     870          } else {
     871            printf(" not taken");
     872            nottaken++;
     873          }
     874        }
     875#endif
     876        /*printf("%d %g %g %g %g\n",itry,lastObj,weightedObj,objvalue,sum1);*/
     877        if (weightedObj > lastObj + 1.0e-4 && itry < MAXTRY) {
     878          if ((logLevel_ & 16) != 0 && doScale) {
     879            printf("Weighted objective from %g to %g **** bad move\n",
     880              lastObj, weightedObj);
     881          }
     882          if (doScale) {
     883            good = 1;
     884          }
     885          if ((strategy & 3) == 1) {
     886            good = 0;
     887            if (weightedObj > lastObj + djExit) {
     888              if ((logLevel_ & 16) != 0) {
     889                printf("Weighted objective from %g to %g ?\n", lastObj, weightedObj);
     890              }
     891              CoinMemcpyN(history[0], ncols, colsol);
     892              CoinMemcpyN(history[0] + ncols, extraBlock, solExtra);
     893              good = 1;
     894            }
     895          } else if ((strategy & 3) == 2) {
     896            if (weightedObj > lastObj + 0.1 * maxDj) {
     897              CoinMemcpyN(history[0], ncols, colsol);
     898              CoinMemcpyN(history[0] + ncols, extraBlock, solExtra);
     899              doScale++;
     900              good = 0;
     901            }
     902          } else if ((strategy & 3) == 3) {
     903            if (weightedObj > lastObj + 0.001 * maxDj) {
     904              /*doScale++;*/
     905              good = 0;
     906            }
     907          }
     908        }
     909      }
     910      if ((iter % checkFrequency_) == 0) {
     911        double best = weightedObj;
     912        double test = obj[0];
     913        for (i = 1; i < DROP; i++) {
     914          obj[i - 1] = obj[i];
     915          if (best > obj[i])
     916            best = obj[i];
     917        }
     918        obj[DROP - 1] = best;
     919        if (test - best < drop && (strategy & 8) == 0) {
     920          if ((logLevel_ & 8) != 0) {
     921            printf("Exiting as drop in %d its is %g after %d iterations\n",
     922              DROP * checkFrequency_, test - best, iter);
     923          }
     924          goto RETURN;
     925        }
     926      }
     927      if ((iter % logFreq_) == 0) {
     928        double piSum = 0.0;
     929        for (i = 0; i < nrows; i++) {
     930          piSum += (rowsol[i] + rowupper[i]) * pi[i];
     931        }
     932        if ((logLevel_ & 2) != 0) {
     933          printf("%d Infeas %g, obj %g - wtObj %g dual %g maxDj %g\n",
     934            iter, sum1, objvalue * maxmin - useOffset, weightedObj - useOffset,
     935            piSum * maxmin - useOffset, maxDj);
     936        }
     937      }
     938      CoinMemcpyN(statusSave, ncols, statusWork);
     939      nflagged = 0;
     940    }
     941    nChange = 0;
     942    doFull = 0;
     943    maxDj = 0.0;
     944    // go through forwards or backwards and starting at odd places
    936945#ifdef FOUR_GOES
    937           for (int i=1;i<FOUR_GOES;i++) {
    938             cilk_spawn memcpy(piX[i], pi, nrows * sizeof(double));
    939             cilk_spawn memcpy(rowsolX[i], rowsol, nrows * sizeof(double));
    940           }
    941           for (int i=0;i<FOUR_GOES;i++) {
    942             nChangeX[i]=0;
    943             maxDjX[i]=0.0;
    944             objvalueX[i]=0.0;
    945             nflaggedX[i]=0;
    946           }
    947           cilk_sync;
    948 #endif
    949           //printf("PASS\n");
     946    for (int i = 1; i < FOUR_GOES; i++) {
     947      cilk_spawn memcpy(piX[i], pi, nrows * sizeof(double));
     948      cilk_spawn memcpy(rowsolX[i], rowsol, nrows * sizeof(double));
     949    }
     950    for (int i = 0; i < FOUR_GOES; i++) {
     951      nChangeX[i] = 0;
     952      maxDjX[i] = 0.0;
     953      objvalueX[i] = 0.0;
     954      nflaggedX[i] = 0;
     955    }
     956    cilk_sync;
     957#endif
     958    //printf("PASS\n");
    950959#ifdef FOUR_GOES
    951                cilk_for (int iPar = 0; iPar < FOUR_GOES; iPar++) {
    952                     double * COIN_RESTRICT pi=piX[iPar];
    953                     double * COIN_RESTRICT rowsol = rowsolX[iPar];
    954           for (int itry = 0; itry < 2; itry++) {
    955                     int istop;
    956                     int istart;
     960    cilk_for(int iPar = 0; iPar < FOUR_GOES; iPar++)
     961    {
     962      double *COIN_RESTRICT pi = piX[iPar];
     963      double *COIN_RESTRICT rowsol = rowsolX[iPar];
     964      for (int itry = 0; itry < 2; itry++) {
     965        int istop;
     966        int istart;
    957967#if 0
    958968                    int chunk = (start[itry]+stop[itry]+FOUR_GOES-1)/FOUR_GOES;
     
    970980                           startsX[itry][iPar],startsX[itry][iPar+1]);
    971981#endif
    972                     istart=startsX[itry][iPar];
    973                     istop=startsX[itry][iPar+1];
     982        istart = startsX[itry][iPar];
     983        istop = startsX[itry][iPar + 1];
    974984#else
    975           for (int itry = 0; itry < 2; itry++) {
    976                int istart = start[itry];
    977                int istop = stop[itry];
    978 #endif
    979                     for (int icol=istart; icol != istop; icol += direction) {
    980                          if (!statusWork[icol]) {
    981                               CoinBigIndex j;
    982                               double value = colsol[icol];
    983                               double djval = cost[icol];
    984                               double djval2, value2;
    985                               double theta, a, b, c;
    986                               if (elemnt) {
    987                                    for (j = columnStart[icol]; j < columnStart[icol] + length[icol]; j++) {
    988                                         int irow = row[j];
    989                                         djval -= elemnt[j] * pi[irow];
    990                                    }
    991                               } else {
    992                                    for (j = columnStart[icol]; j < columnStart[icol] + length[icol]; j++) {
    993                                         int irow = row[j];
    994                                         djval -= pi[irow];
    995                                    }
    996                               }
    997                               /*printf("xx iter %d seq %d djval %g value %g\n",
     985    for (int itry = 0; itry < 2; itry++) {
     986      int istart = start[itry];
     987      int istop = stop[itry];
     988#endif
     989        for (int icol = istart; icol != istop; icol += direction) {
     990          if (!statusWork[icol]) {
     991            CoinBigIndex j;
     992            double value = colsol[icol];
     993            double djval = cost[icol];
     994            double djval2, value2;
     995            double theta, a, b, c;
     996            if (elemnt) {
     997              for (j = columnStart[icol]; j < columnStart[icol] + length[icol]; j++) {
     998                int irow = row[j];
     999                djval -= elemnt[j] * pi[irow];
     1000              }
     1001            } else {
     1002              for (j = columnStart[icol]; j < columnStart[icol] + length[icol]; j++) {
     1003                int irow = row[j];
     1004                djval -= pi[irow];
     1005              }
     1006            }
     1007            /*printf("xx iter %d seq %d djval %g value %g\n",
    9981008                                iter,i,djval,value);*/
    999                               if (djval > 1.0e-5) {
    1000                                    value2 = (lower[icol] - value);
    1001                               } else {
    1002                                    value2 = (upper[icol] - value);
    1003                               }
    1004                               djval2 = djval * value2;
    1005                               djval = fabs(djval);
    1006                               if (djval > djTol) {
    1007                                    if (djval2 < -1.0e-4) {
     1009            if (djval > 1.0e-5) {
     1010              value2 = (lower[icol] - value);
     1011            } else {
     1012              value2 = (upper[icol] - value);
     1013            }
     1014            djval2 = djval * value2;
     1015            djval = fabs(djval);
     1016            if (djval > djTol) {
     1017              if (djval2 < -1.0e-4) {
    10081018#ifndef FOUR_GOES
    1009                                         nChange++;
    1010                                         if (djval > maxDj) maxDj = djval;
     1019                nChange++;
     1020                if (djval > maxDj)
     1021                  maxDj = djval;
    10111022#else
    1012                                         nChangeX[iPar]++;
    1013                                         if (djval > maxDjX[iPar]) maxDjX[iPar] = djval;
    1014 #endif
    1015                                         /*if (djval>3.55e6) {
     1023              nChangeX[iPar]++;
     1024              if (djval > maxDjX[iPar])
     1025                maxDjX[iPar] = djval;
     1026#endif
     1027                /*if (djval>3.55e6) {
    10161028                                                        printf("big\n");
    10171029                                                        }*/
    1018                                         a = 0.0;
    1019                                         b = 0.0;
    1020                                         c = 0.0;
    1021                                         djval2 = cost[icol];
    1022                                         if (elemnt) {
    1023                                              for (j = columnStart[icol]; j < columnStart[icol] + length[icol]; j++) {
    1024                                                   int irow = row[j];
    1025                                                   double value = rowsol[irow];
    1026                                                   c += value * value;
    1027                                                   a += elemnt[j] * elemnt[j];
    1028                                                   b += value * elemnt[j];
    1029                                              }
    1030                                         } else {
    1031                                              for (j = columnStart[icol]; j < columnStart[icol] + length[icol]; j++) {
    1032                                                   int irow = row[j];
    1033                                                   double value = rowsol[irow];
    1034                                                   c += value * value;
    1035                                                   a += 1.0;
    1036                                                   b += value;
    1037                                              }
    1038                                         }
    1039                                         a *= weight;
    1040                                         b = b * weight + 0.5 * djval2;
    1041                                         c *= weight;
    1042                                         /* solve */
    1043                                         theta = -b / a;
     1030                a = 0.0;
     1031                b = 0.0;
     1032                c = 0.0;
     1033                djval2 = cost[icol];
     1034                if (elemnt) {
     1035                  for (j = columnStart[icol]; j < columnStart[icol] + length[icol]; j++) {
     1036                    int irow = row[j];
     1037                    double value = rowsol[irow];
     1038                    c += value * value;
     1039                    a += elemnt[j] * elemnt[j];
     1040                    b += value * elemnt[j];
     1041                  }
     1042                } else {
     1043                  for (j = columnStart[icol]; j < columnStart[icol] + length[icol]; j++) {
     1044                    int irow = row[j];
     1045                    double value = rowsol[irow];
     1046                    c += value * value;
     1047                    a += 1.0;
     1048                    b += value;
     1049                  }
     1050                }
     1051                a *= weight;
     1052                b = b * weight + 0.5 * djval2;
     1053                c *= weight;
     1054                /* solve */
     1055                theta = -b / a;
    10441056#ifndef FOUR_GOES
    1045                                         if ((strategy & 4) != 0) {
    1046                                           double valuep, thetap;
    1047                                              value2 = a * theta * theta + 2.0 * b * theta;
    1048                                              thetap = 2.0 * theta;
    1049                                              valuep = a * thetap * thetap + 2.0 * b * thetap;
    1050                                              if (valuep < value2 + djTol) {
    1051                                                   theta = thetap;
    1052                                                   kgood++;
    1053                                              } else {
    1054                                                   kbad++;
    1055                                              }
    1056                                         }
    1057 #endif
    1058                                         if (theta > 0.0) {
    1059                                              if (theta < upper[icol] - colsol[icol]) {
    1060                                                   value2 = theta;
    1061                                              } else {
    1062                                                   value2 = upper[icol] - colsol[icol];
    1063                                              }
    1064                                         } else {
    1065                                              if (theta > lower[icol] - colsol[icol]) {
    1066                                                   value2 = theta;
    1067                                              } else {
    1068                                                   value2 = lower[icol] - colsol[icol];
    1069                                              }
    1070                                         }
    1071                                         colsol[icol] += value2;
     1057                if ((strategy & 4) != 0) {
     1058                  double valuep, thetap;
     1059                  value2 = a * theta * theta + 2.0 * b * theta;
     1060                  thetap = 2.0 * theta;
     1061                  valuep = a * thetap * thetap + 2.0 * b * thetap;
     1062                  if (valuep < value2 + djTol) {
     1063                    theta = thetap;
     1064                    kgood++;
     1065                  } else {
     1066                    kbad++;
     1067                  }
     1068                }
     1069#endif
     1070                if (theta > 0.0) {
     1071                  if (theta < upper[icol] - colsol[icol]) {
     1072                    value2 = theta;
     1073                  } else {
     1074                    value2 = upper[icol] - colsol[icol];
     1075                  }
     1076                } else {
     1077                  if (theta > lower[icol] - colsol[icol]) {
     1078                    value2 = theta;
     1079                  } else {
     1080                    value2 = lower[icol] - colsol[icol];
     1081                  }
     1082                }
     1083                colsol[icol] += value2;
    10721084#ifndef FOUR_GOES
    1073                                         objvalue += cost[icol] * value2;
     1085                objvalue += cost[icol] * value2;
    10741086#else
    1075                                         objvalueX[iPar] += cost[icol] * value2;
    1076 #endif
    1077                                         if (elemnt) {
    1078                                              for (j = columnStart[icol]; j < columnStart[icol] + length[icol]; j++) {
    1079                                                   int irow = row[j];
    1080                                                   double value;
    1081                                                   rowsol[irow] += elemnt[j] * value2;
    1082                                                   value = rowsol[irow];
    1083                                                   pi[irow] = -2.0 * weight * value;
    1084                                              }
    1085                                         } else {
    1086                                              for (j = columnStart[icol]; j < columnStart[icol] + length[icol]; j++) {
    1087                                                   int irow = row[j];
    1088                                                   double value;
    1089                                                   rowsol[irow] += value2;
    1090                                                   value = rowsol[irow];
    1091                                                   pi[irow] = -2.0 * weight * value;
    1092                                              }
    1093                                         }
    1094                                    } else {
    1095                                         /* dj but at bound */
    1096                                         if (djval > djFlag) {
    1097                                              statusWork[icol] = 1;
     1087              objvalueX[iPar] += cost[icol] * value2;
     1088#endif
     1089                if (elemnt) {
     1090                  for (j = columnStart[icol]; j < columnStart[icol] + length[icol]; j++) {
     1091                    int irow = row[j];
     1092                    double value;
     1093                    rowsol[irow] += elemnt[j] * value2;
     1094                    value = rowsol[irow];
     1095                    pi[irow] = -2.0 * weight * value;
     1096                  }
     1097                } else {
     1098                  for (j = columnStart[icol]; j < columnStart[icol] + length[icol]; j++) {
     1099                    int irow = row[j];
     1100                    double value;
     1101                    rowsol[irow] += value2;
     1102                    value = rowsol[irow];
     1103                    pi[irow] = -2.0 * weight * value;
     1104                  }
     1105                }
     1106              } else {
     1107                /* dj but at bound */
     1108                if (djval > djFlag) {
     1109                  statusWork[icol] = 1;
    10981110#ifndef FOUR_GOES
    1099                                              nflagged++;
     1111                  nflagged++;
    11001112#else
    1101                                              nflaggedX[iPar]++;
    1102 #endif
    1103                                         }
    1104                                    }
    1105                               }
    1106                          }
    1107                     }
     1113                nflaggedX[iPar]++;
     1114#endif
     1115                }
     1116              }
     1117            }
     1118          }
     1119        }
    11081120#ifdef FOUR_GOES
    1109                }
    1110 #endif
    1111           }
     1121      }
     1122#endif
     1123    }
    11121124#ifdef FOUR_GOES
    1113           for (int i=0;i<FOUR_GOES;i++) {
    1114             nChange += nChangeX[i];
    1115             maxDj = CoinMax(maxDj,maxDjX[i]);
    1116             objvalue += objvalueX[i];
    1117             nflagged += nflaggedX[i];
    1118           }
    1119           cilk_for (int i = 0; i < nrows; i++) {
    1120 #if FOUR_GOES==2
    1121             rowsol[i] = 0.5 * (rowsolX[0][i] + rowsolX[1][i]);
    1122             pi[i] = 0.5 * (piX[0][i] + piX[1][i]);
    1123 #elif FOUR_GOES==3
    1124             pi[i] = 0.33333333333333 * (piX[0][i] + piX[1][i]+piX[2][i]);
    1125             rowsol[i] = 0.3333333333333 * (rowsolX[0][i] + rowsolX[1][i]+rowsolX[2][i]);
     1125    for (int i = 0; i < FOUR_GOES; i++) {
     1126      nChange += nChangeX[i];
     1127      maxDj = CoinMax(maxDj, maxDjX[i]);
     1128      objvalue += objvalueX[i];
     1129      nflagged += nflaggedX[i];
     1130    }
     1131    cilk_for(int i = 0; i < nrows; i++)
     1132    {
     1133#if FOUR_GOES == 2
     1134      rowsol[i] = 0.5 * (rowsolX[0][i] + rowsolX[1][i]);
     1135      pi[i] = 0.5 * (piX[0][i] + piX[1][i]);
     1136#elif FOUR_GOES == 3
     1137      pi[i] = 0.33333333333333 * (piX[0][i] + piX[1][i] + piX[2][i]);
     1138      rowsol[i] = 0.3333333333333 * (rowsolX[0][i] + rowsolX[1][i] + rowsolX[2][i]);
    11261139#else
    1127             pi[i] = 0.25 * (piX[0][i] + piX[1][i]+piX[2][i]+piX[3][i]);
    1128             rowsol[i] = 0.25 * (rowsolX[0][i] + rowsolX[1][i]+rowsolX[2][i]+rowsolX[3][i]);
    1129 #endif
    1130           }
    1131 #endif
    1132           if (extraBlock) {
    1133                for (int i = 0; i < extraBlock; i++) {
    1134                     double value = solExtra[i];
    1135                     double djval = costExtra[i];
    1136                     double djval2, value2;
    1137                     double element = elemExtra[i];
    1138                     double theta, a, b, c;
    1139                     int irow = rowExtra[i];
    1140                     djval -= element * pi[irow];
    1141                     /*printf("xxx iter %d extra %d djval %g value %g\n",
     1140      pi[i] = 0.25 * (piX[0][i] + piX[1][i] + piX[2][i] + piX[3][i]);
     1141      rowsol[i] = 0.25 * (rowsolX[0][i] + rowsolX[1][i] + rowsolX[2][i] + rowsolX[3][i]);
     1142#endif
     1143    }
     1144#endif
     1145    if (extraBlock) {
     1146      for (int i = 0; i < extraBlock; i++) {
     1147        double value = solExtra[i];
     1148        double djval = costExtra[i];
     1149        double djval2, value2;
     1150        double element = elemExtra[i];
     1151        double theta, a, b, c;
     1152        int irow = rowExtra[i];
     1153        djval -= element * pi[irow];
     1154        /*printf("xxx iter %d extra %d djval %g value %g\n",
    11421155                          iter,irow,djval,value);*/
    1143                     if (djval > 1.0e-5) {
    1144                          value2 = -value;
    1145                     } else {
    1146                          value2 = (upperExtra[i] - value);
    1147                     }
    1148                     djval2 = djval * value2;
    1149                     if (djval2 < -1.0e-4 && fabs(djval) > djTol) {
    1150                          nChange++;
    1151                          a = 0.0;
    1152                          b = 0.0;
    1153                          c = 0.0;
    1154                          djval2 = costExtra[i];
    1155                          value = rowsol[irow];
    1156                          c += value * value;
    1157                          a += element * element;
    1158                          b += element * value;
    1159                          a *= weight;
    1160                          b = b * weight + 0.5 * djval2;
    1161                          c *= weight;
    1162                          /* solve */
    1163                          theta = -b / a;
    1164                          if (theta > 0.0) {
    1165                               value2 = CoinMin(theta, upperExtra[i] - solExtra[i]);
    1166                          } else {
    1167                               value2 = CoinMax(theta, -solExtra[i]);
    1168                          }
    1169                          solExtra[i] += value2;
    1170                          rowsol[irow] += element * value2;
    1171                          value = rowsol[irow];
    1172                          pi[irow] = -2.0 * weight * value;
    1173                     }
    1174                }
    1175           }
    1176           if ((iter % 10) == 2) {
    1177                for (int i = DJTEST - 1; i > 0; i--) {
    1178                     djSave[i] = djSave[i-1];
    1179                }
    1180                djSave[0] = maxDj;
    1181                largestDj = CoinMax(largestDj, maxDj);
    1182                smallestDj = CoinMin(smallestDj, maxDj);
    1183                for (int i = DJTEST - 1; i > 0; i--) {
    1184                     maxDj += djSave[i];
    1185                }
    1186                maxDj = maxDj / static_cast<double> (DJTEST);
    1187                if (maxDj < djExit && iter > 50) {
    1188                     //printf("Exiting on low dj %g after %d iterations\n",maxDj,iter);
    1189                     break;
    1190                }
    1191                if (nChange < 100) {
    1192                     djTol *= 0.5;
    1193                }
    1194           }
    1195      }
     1156        if (djval > 1.0e-5) {
     1157          value2 = -value;
     1158        } else {
     1159          value2 = (upperExtra[i] - value);
     1160        }
     1161        djval2 = djval * value2;
     1162        if (djval2 < -1.0e-4 && fabs(djval) > djTol) {
     1163          nChange++;
     1164          a = 0.0;
     1165          b = 0.0;
     1166          c = 0.0;
     1167          djval2 = costExtra[i];
     1168          value = rowsol[irow];
     1169          c += value * value;
     1170          a += element * element;
     1171          b += element * value;
     1172          a *= weight;
     1173          b = b * weight + 0.5 * djval2;
     1174          c *= weight;
     1175          /* solve */
     1176          theta = -b / a;
     1177          if (theta > 0.0) {
     1178            value2 = CoinMin(theta, upperExtra[i] - solExtra[i]);
     1179          } else {
     1180            value2 = CoinMax(theta, -solExtra[i]);
     1181          }
     1182          solExtra[i] += value2;
     1183          rowsol[irow] += element * value2;
     1184          value = rowsol[irow];
     1185          pi[irow] = -2.0 * weight * value;
     1186        }
     1187      }
     1188    }
     1189    if ((iter % 10) == 2) {
     1190      for (int i = DJTEST - 1; i > 0; i--) {
     1191        djSave[i] = djSave[i - 1];
     1192      }
     1193      djSave[0] = maxDj;
     1194      largestDj = CoinMax(largestDj, maxDj);
     1195      smallestDj = CoinMin(smallestDj, maxDj);
     1196      for (int i = DJTEST - 1; i > 0; i--) {
     1197        maxDj += djSave[i];
     1198      }
     1199      maxDj = maxDj / static_cast< double >(DJTEST);
     1200      if (maxDj < djExit && iter > 50) {
     1201        //printf("Exiting on low dj %g after %d iterations\n",maxDj,iter);
     1202        break;
     1203      }
     1204      if (nChange < 100) {
     1205        djTol *= 0.5;
     1206      }
     1207    }
     1208  }
    11961209RETURN:
    1197      if (kgood || kbad) {
    1198        COIN_DETAIL_PRINT(printf("%g good %g bad\n", kgood, kbad));
    1199      }
    1200      result = objval(nrows, ncols, rowsol, colsol, pi, djs, useCost,
    1201                      rowlower, rowupper, lower, upper,
    1202                      elemnt, row, columnStart, length, extraBlock, rowExtra,
    1203                      solExtra, elemExtra, upperExtra, useCostExtra,
    1204                      weight);
    1205      result.djAtBeginning = largestDj;
    1206      result.djAtEnd = smallestDj;
    1207      result.dropThis = saveValue - result.weighted;
    1208      if (saveSol) {
    1209           if (result.weighted < bestSol) {
    1210                bestSol = result.weighted;
    1211                CoinMemcpyN(colsol, ncols, saveSol);
    1212           } else {
    1213                COIN_DETAIL_PRINT(printf("restoring previous - now %g best %g\n",
    1214                                         result.weighted * maxmin - useOffset, bestSol * maxmin - useOffset));
    1215           }
    1216      }
    1217      if (saveSol) {
    1218           if (extraBlock) {
    1219                delete [] useCostExtra;
    1220           }
    1221           CoinMemcpyN(saveSol, ncols, colsol);
    1222           delete [] saveSol;
    1223      }
    1224      for (i = 0; i < nsolve; i++) {
    1225           if (i) delete [] allsum[i];
    1226           delete [] aX[i];
    1227           delete [] aworkX[i];
    1228      }
    1229      delete [] thetaX;
    1230      delete [] djX;
    1231      delete [] bX;
    1232      delete [] vX;
    1233      delete [] aX;
    1234      delete [] aworkX;
    1235      delete [] allsum;
    1236      delete [] cost;
     1210  if (kgood || kbad) {
     1211    COIN_DETAIL_PRINT(printf("%g good %g bad\n", kgood, kbad));
     1212  }
     1213  result = objval(nrows, ncols, rowsol, colsol, pi, djs, useCost,
     1214    rowlower, rowupper, lower, upper,
     1215    elemnt, row, columnStart, length, extraBlock, rowExtra,
     1216    solExtra, elemExtra, upperExtra, useCostExtra,
     1217    weight);
     1218  result.djAtBeginning = largestDj;
     1219  result.djAtEnd = smallestDj;
     1220  result.dropThis = saveValue - result.weighted;
     1221  if (saveSol) {
     1222    if (result.weighted < bestSol) {
     1223      bestSol = result.weighted;
     1224      CoinMemcpyN(colsol, ncols, saveSol);
     1225    } else {
     1226      COIN_DETAIL_PRINT(printf("restoring previous - now %g best %g\n",
     1227        result.weighted * maxmin - useOffset, bestSol * maxmin - useOffset));
     1228    }
     1229  }
     1230  if (saveSol) {
     1231    if (extraBlock) {
     1232      delete[] useCostExtra;
     1233    }
     1234    CoinMemcpyN(saveSol, ncols, colsol);
     1235    delete[] saveSol;
     1236  }
     1237  for (i = 0; i < nsolve; i++) {
     1238    if (i)
     1239      delete[] allsum[i];
     1240    delete[] aX[i];
     1241    delete[] aworkX[i];
     1242  }
     1243  delete[] thetaX;
     1244  delete[] djX;
     1245  delete[] bX;
     1246  delete[] vX;
     1247  delete[] aX;
     1248  delete[] aworkX;
     1249  delete[] allsum;
     1250  delete[] cost;
    12371251#ifdef FOUR_GOES
    1238      delete [] pi2 ;
    1239      delete [] rowsol2 ;
    1240 #endif
    1241      for (i = 0; i < HISTORY + 1; i++) {
    1242           delete [] history[i];
    1243      }
    1244      delete [] statusSave;
    1245      /* do original costs objvalue*/
    1246      result.objval = 0.0;
    1247      for (i = 0; i < ncols; i++) {
    1248           result.objval += colsol[i] * origcost[i];
    1249      }
    1250      if (extraBlock) {
    1251           for (i = 0; i < extraBlock; i++) {
    1252                int irow = rowExtra[i];
    1253                rowupper[irow] = saveExtra[i];
    1254           }
    1255           delete [] rowExtra;
    1256           delete [] solExtra;
    1257           delete [] elemExtra;
    1258           delete [] upperExtra;
    1259           delete [] costExtra;
    1260           delete [] saveExtra;
    1261      }
    1262      result.iteration = iter;
    1263      result.objval -= saveOffset;
    1264      result.weighted = result.objval + weight * result.sumSquared;
    1265      return result;
     1252  delete[] pi2;
     1253  delete[] rowsol2;
     1254#endif
     1255  for (i = 0; i < HISTORY + 1; i++) {
     1256    delete[] history[i];
     1257  }
     1258  delete[] statusSave;
     1259  /* do original costs objvalue*/
     1260  result.objval = 0.0;
     1261  for (i = 0; i < ncols; i++) {
     1262    result.objval += colsol[i] * origcost[i];
     1263  }
     1264  if (extraBlock) {
     1265    for (i = 0; i < extraBlock; i++) {
     1266      int irow = rowExtra[i];
     1267      rowupper[irow] = saveExtra[i];
     1268    }
     1269    delete[] rowExtra;
     1270    delete[] solExtra;
     1271    delete[] elemExtra;
     1272    delete[] upperExtra;
     1273    delete[] costExtra;
     1274    delete[] saveExtra;
     1275  }
     1276  result.iteration = iter;
     1277  result.objval -= saveOffset;
     1278  result.weighted = result.objval + weight * result.sumSquared;
     1279  return result;
    12661280}
     1281
     1282/* vi: softtabstop=2 shiftwidth=2 expandtab tabstop=2
     1283*/
Note: See TracChangeset for help on using the changeset viewer.