Ignore:
Timestamp:
Oct 2, 2003 1:21:02 PM (17 years ago)
Author:
forrest
Message:

lots of stuff

File:
1 edited

Legend:

Unmodified
Added
Removed
  • branches/pre/ClpSimplexPrimal.cpp

    r210 r212  
    185185  // save data
    186186  ClpDataSave data = saveData();
    187 
    188 
    189     // initialize - maybe values pass and algorithm_ is +1
     187 
     188  // Save so can see if doing after dual
     189  int initialStatus=problemStatus_;
     190  // initialize - maybe values pass and algorithm_ is +1
    190191    if (!startup(ifValuesPass)) {
    191192
     
    242243      // If getting nowhere - why not give it a kick
    243244#if 1
    244       if (perturbation_<101&&numberIterations_>2*(numberRows_+numberColumns_))
     245      if (perturbation_<101&&numberIterations_>2*(numberRows_+numberColumns_)
     246          &&initialStatus!=10)
    245247        perturb(1);
    246248#endif
     
    262264     
    263265      // Iterate
    264       whileIterating();
     266      whileIterating(ifValuesPass);
    265267    }
    266268  }
     
    292294*/
    293295int
    294 ClpSimplexPrimal::whileIterating()
     296ClpSimplexPrimal::whileIterating(int valuesOption)
    295297{
    296298  // Say if values pass
    297299  int ifValuesPass=(firstFree_>=0) ? 1 : 0;
    298300  int returnCode=-1;
     301  int superBasicType=1;
     302  if (valuesOption>1)
     303    superBasicType=3;
    299304  // status stays at -1 while iterating, >=0 finished, -2 to invert
    300305  // status -3 to go to top without an invert
     
    352357    } else {
    353358      // in values pass
    354       int sequenceIn=nextSuperBasic();
     359      int sequenceIn=nextSuperBasic(superBasicType,columnArray_[0]);
     360      if (valuesOption>1)
     361        superBasicType=2;
    355362      if (sequenceIn<0) {
    356363        // end of values pass - initialize weights etc
     
    418425    }
    419426  }
     427  if (valuesOption>1)
     428    columnArray_[0]->setNumElements(0);
    420429  return returnCode;
    421430}
     
    604613          gutsOfSolution(NULL,NULL);
    605614          problemStatus_=-1; //continue
     615          goToDual=false;
    606616        } else {
    607617          // say infeasible
     
    13691379  if (perturbation_>100)
    13701380    return; //perturbed already
     1381  if (perturbation_==100)
     1382    perturbation_=50; // treat as normal
    13711383  int i;
     1384  if (!numberIterations_)
     1385    cleanStatus(); // make sure status okay
     1386  if (!numberIterations_&&perturbation_==50) {
     1387    // See if we need to perturb
     1388    double * sort = new double[numberRows_];
     1389    for (i=0;i<numberRows_;i++) {
     1390      double lo = fabs(lower_[i]);
     1391      double up = fabs(upper_[i]);
     1392      double value=0.0;
     1393      if (lo&&lo<1.0e20) {
     1394        if (up&&up<1.0e20)
     1395          value = 0.5*(lo+up);
     1396        else
     1397          value=lo;
     1398      } else {
     1399        if (up&&up<1.0e20)
     1400          value = up;
     1401      }
     1402      sort[i]=value;
     1403    }
     1404    std::sort(sort,sort+numberRows_);
     1405    int number=1;
     1406    double last = sort[0];
     1407    for (i=1;i<numberRows_;i++) {
     1408      if (last!=sort[i])
     1409        number++;
     1410      last=sort[i];
     1411    }
     1412    delete [] sort;
     1413    if (number*4>numberRows_)
     1414      return; // good enough
     1415  }
    13721416  // primal perturbation
    13731417  double perturbation=1.0e-20;
     1418  int numberNonZero=0;
    13741419  // maximum fraction of rhs/bounds to perturb
    1375   double maximumFraction = 1.0e-8;
     1420  double maximumFraction = 1.0e-7;
    13761421  if (perturbation_>=50) {
    13771422    perturbation = 1.0e-4;
     
    13871432        else
    13881433          upperValue=0.0;
    1389         double value = max(lowerValue,upperValue);
     1434        double value = max(fabs(lowerValue),fabs(upperValue));
    13901435        value = min(value,upper_[i]-lower_[i]);
     1436#if 1
     1437        if (value) {
     1438          perturbation += value;
     1439          numberNonZero++;
     1440        }
     1441#else
    13911442        perturbation = max(perturbation,value);
    1392       }
    1393     }
     1443#endif
     1444      }
     1445    }
     1446    if (numberNonZero)
     1447      perturbation /= (double) numberNonZero;
     1448    else
     1449      perturbation = 1.0e-1;
    13941450  } else if (perturbation_<100) {
    13951451    perturbation = pow(10.0,perturbation_);
     
    14011457  double largestPerCent=0.0;
    14021458  bool printOut=(handler_->logLevel()==63);
     1459  printOut=false; //off
    14031460  // Check if all slack
    14041461  int number=0;
     
    14801537                 i,lower_[i],lowerValue,upper_[i],upperValue);
    14811538      }
    1482       if (solution_[i]==lower_[i])
    1483         solution_[i]=lowerValue;
    1484       else if (solution_[i]==upper_[i])
    1485         solution_[i]=upperValue;
    14861539      lower_[i]=lowerValue;
    14871540      upper_[i]=upperValue;
     
    15181571        printf("row %d lower from %g to %g, upper from %g to %g\n",
    15191572               i-numberColumns_,lower_[i],lowerValue,upper_[i],upperValue);
    1520       if (solution_[i]==lower_[i])
    1521         solution_[i]=lowerValue;
    1522       else if (solution_[i]==upper_[i])
    1523         solution_[i]=upperValue;
    15241573      lower_[i]=lowerValue;
    15251574      upper_[i]=upperValue;
     1575    }
     1576  }
     1577  // Clean up
     1578  for (i=0;i<numberColumns_+numberRows_;i++) {
     1579    switch(getStatus(i)) {
     1580     
     1581    case basic:
     1582      break;
     1583    case atUpperBound:
     1584      solution_[i]=upper_[i];
     1585      break;
     1586    case isFixed:
     1587    case atLowerBound:
     1588      solution_[i]=lower_[i];
     1589      break;
     1590    case isFree:
     1591      break;
     1592    case superBasic:
     1593      break;
    15261594    }
    15271595  }
     
    18101878      }
    18111879      // here do part of steepest - ready for next iteration
    1812       primalColumnPivot_->updateWeights(rowArray_[1]);
     1880      if (!ifValuesPass)
     1881        primalColumnPivot_->updateWeights(rowArray_[1]);
    18131882    } else {
    18141883      if (pivotRow_==-1) {
     
    19552024  }
    19562025}
     2026/* Get next superbasic -1 if none,
     2027   Normal type is 1
     2028   If type is 3 then initializes sorted list
     2029   if 2 uses list.
     2030*/
     2031int
     2032ClpSimplexPrimal::nextSuperBasic(int superBasicType,CoinIndexedVector * columnArray)
     2033{
     2034  if (firstFree_>=0&&superBasicType) {
     2035    int returnValue=firstFree_;
     2036    int iColumn=firstFree_+1;
     2037    if (superBasicType>1) {
     2038      if (superBasicType>2) {
     2039        // Initialize list
     2040        // Wild guess that lower bound more natural than upper
     2041        int number=0;
     2042        double * work=columnArray->denseVector();
     2043        int * which=columnArray->getIndices();
     2044        for (iColumn=0;iColumn<numberRows_+numberColumns_;iColumn++) {
     2045          if (getStatus(iColumn)==superBasic) {
     2046            if (fabs(solution_[iColumn]-lower_[iColumn])<=primalTolerance_) {
     2047              solution_[iColumn]=lower_[iColumn];
     2048              setStatus(iColumn,atLowerBound);
     2049            } else if (fabs(solution_[iColumn]-upper_[iColumn])
     2050                       <=primalTolerance_) {
     2051              solution_[iColumn]=upper_[iColumn];
     2052              setStatus(iColumn,atUpperBound);
     2053            } else if (lower_[iColumn]<-1.0e20&&upper_[iColumn]>1.0e20) {
     2054              setStatus(iColumn,isFree);
     2055              break;
     2056            } else if (!flagged(iColumn)) {
     2057              // put ones near bounds at end after sorting
     2058              work[number]= - min(0.1*(solution_[iColumn]-lower_[iColumn]),
     2059                                  upper_[iColumn]-solution_[iColumn]);
     2060              which[number++] = iColumn;
     2061            }
     2062          }
     2063        }
     2064        CoinSort_2(work,work+number,which);
     2065        columnArray->setNumElements(number);
     2066        memset(work,0,number*sizeof(double));
     2067      }
     2068      int * which=columnArray->getIndices();
     2069      int number = columnArray->getNumElements();
     2070      if (!number) {
     2071        // finished
     2072        iColumn = numberRows_+numberColumns_;
     2073        returnValue=-1;
     2074      } else {
     2075        number--;
     2076        returnValue=which[number];
     2077        iColumn=returnValue;
     2078        columnArray->setNumElements(number);
     2079      }     
     2080    } else {
     2081      for (;iColumn<numberRows_+numberColumns_;iColumn++) {
     2082        if (getStatus(iColumn)==superBasic) {
     2083          if (fabs(solution_[iColumn]-lower_[iColumn])<=primalTolerance_) {
     2084            solution_[iColumn]=lower_[iColumn];
     2085            setStatus(iColumn,atLowerBound);
     2086          } else if (fabs(solution_[iColumn]-upper_[iColumn])
     2087                     <=primalTolerance_) {
     2088            solution_[iColumn]=upper_[iColumn];
     2089            setStatus(iColumn,atUpperBound);
     2090          } else if (lower_[iColumn]<-1.0e20&&upper_[iColumn]>1.0e20) {
     2091            setStatus(iColumn,isFree);
     2092            break;
     2093          } else {
     2094            break;
     2095          }
     2096        }
     2097      }
     2098    }
     2099    firstFree_ = iColumn;
     2100    if (firstFree_==numberRows_+numberColumns_)
     2101      firstFree_=-1;
     2102    return returnValue;
     2103  } else {
     2104    return -1;
     2105  }
     2106}
     2107 
    19572108 
    19582109
Note: See TracChangeset for help on using the changeset viewer.