Changeset 190 for branches


Ignore:
Timestamp:
Aug 5, 2003 4:57:40 AM (16 years ago)
Author:
forrest
Message:

Stuff

Location:
branches/pre
Files:
7 edited

Legend:

Unmodified
Added
Removed
  • branches/pre/ClpMessage.cpp

    r180 r190  
    6161  {CLP_CRASH,28,1,"Crash put %d variables in basis, %d dual infeasibilities"},
    6262  {CLP_END_VALUES_PASS,29,1,"End of values pass after %d iterations"},
     63  {CLP_QUADRATIC_BOTH,108,32,"%s %d (%g) and %d (%g) both basic"},
     64  {CLP_QUADRATIC_PRIMAL_DETAILS,109,32,"coeff %g, %g, %g - dj %g - deriv zero at %g, sj at %g"},
    6365  {CLP_DUMMY_END,999999,0,""}
    6466};
  • branches/pre/ClpNonLinearCost.cpp

    r185 r190  
    661661  int end = start_[iPivot+1]-1;
    662662  if (!bothWays_) {
    663     for (iRange=start; iRange<end;iRange++) {
    664       if (value<lower_[iRange+1]+primalTolerance) {
    665         // put in better range
    666         if (value>=lower_[iRange+1]-primalTolerance&&infeasible(iRange)&&iRange==start)
    667           iRange++;
    668         break;
     663    // If fixed try and get feasible
     664    if (lower_[start+1]==lower_[start+2]&&fabs(value-lower_[start+1])<1.001*primalTolerance) {
     665      iRange =start+1;
     666    } else {
     667      for (iRange=start; iRange<end;iRange++) {
     668        if (value<=lower_[iRange+1]+primalTolerance) {
     669          // put in better range
     670          if (value>=lower_[iRange+1]-primalTolerance&&infeasible(iRange)&&iRange==start)
     671            iRange++;
     672          break;
     673        }
    669674      }
    670675    }
  • branches/pre/ClpPrimalQuadraticDantzig.cpp

    r183 r190  
    107107  const double * objective = quadraticInfo_->originalObjective();
    108108  for (iSequence=0;iSequence<originalNumberColumns;iSequence++) {
     109    if (model_->flagged(iSequence))
     110      continue;
    109111    int jSequence = which[iSequence];
    110112    double value;
     
    163165  for (jSequence=0;jSequence<originalNumberRows;jSequence++) {
    164166    int iSequence  = jSequence + firstSlack;
     167    if (model_->flagged(iSequence))
     168      continue;
    165169    int iPi = jSequence+originalNumberColumns;
    166170    // for slacks either pi zero or row at bound
  • branches/pre/ClpSimplex.cpp

    r182 r190  
    10811081    }
    10821082  }
    1083 #if 0
    1084   int status=-99;
    1085   int * rowIsBasic = new int[numberRows_];
    1086   int * columnIsBasic = new int[numberColumns_];
    1087   //returns 0 -okay, -1 singular, -2 too many in basis, -99 memory */
    1088   while (status<-98) {
    1089  
     1083  if (0)  {
     1084    static int k=0;
     1085    printf("start basis\n");
    10901086    int i;
    1091     int numberBasic=0;
    1092     for (i=0;i<numberRows_;i++) {
    1093       if (getRowStatus(i) == basic) {
    1094         rowIsBasic[i]=1;
    1095         numberBasic++;
    1096       } else {
    1097         rowIsBasic[i]=-1;
    1098       }
    1099     }
    1100     for (i=0;i<numberColumns_;i++) {
    1101       if (getColumnStatus(i) == basic) {
    1102         columnIsBasic[i]=1;
    1103         numberBasic++;
    1104       } else {
    1105         columnIsBasic[i]=-1;
    1106       }
    1107     }
    1108     assert (numberBasic<=numberRows_);
    1109     while (status==-99) {
    1110       status =  factorization_->factorizeOld(this,matrix_,
    1111                                           numberRows_,numberColumns_,
    1112                                           rowIsBasic, columnIsBasic,
    1113                                           0.0);
    1114       if (status==-99) {
    1115         // get more memory
    1116         factorization_->areaFactor(2.0*factorization_->areaFactor());
    1117       }
    1118     }
    1119     if (!status) {
    1120       // See whether to redo pivot order
    1121       if (factorization_->needToReorder()||
    1122           progress_->lastIterationNumber(0)<0) {
    1123         // do pivot information
    1124         for (i=0;i<numberRows_;i++) {
    1125           if (getRowStatus(i) == basic) {
    1126             pivotVariable_[rowIsBasic[i]]=i+numberColumns_;
    1127           }
    1128         }
    1129         for (i=0;i<numberColumns_;i++) {
    1130           if (getColumnStatus(i) == basic) {
    1131             pivotVariable_[columnIsBasic[i]]=i;
    1132           }
    1133         }
    1134       }
    1135     } else {
    1136       // leave pivotVariable_ in useful form for cleaning basis
    1137       for (i=0;i<numberRows_;i++) {
    1138         pivotVariable_[i]=-1;
    1139       }
    1140       for (i=0;i<numberRows_;i++) {
    1141         if (getRowStatus(i) == basic) {
    1142           int iPivot = rowIsBasic[i];
    1143           if (iPivot>=0)
    1144             pivotVariable_[iPivot]=i+numberColumns_;
    1145         }
    1146       }
    1147       for (i=0;i<numberColumns_;i++) {
    1148         if (getColumnStatus(i) == basic) {
    1149           int iPivot = columnIsBasic[i];
    1150           if (iPivot>=0)
    1151             pivotVariable_[iPivot]=i;
    1152         }
    1153       }
    1154     }
    1155     if (status==-1) {
    1156       if (!solveType) {
    1157         //redo basis - first take ALL columns out
    1158         for (iColumn=0;iColumn<numberColumns_;iColumn++) {
    1159           if (getColumnStatus(iColumn)==
    1160               basic) {
    1161             // take out
    1162             if (!valuesPass) {
    1163               double lower=columnLowerWork_[iColumn];
    1164               double upper=columnUpperWork_[iColumn];
    1165               double value=columnActivityWork_[iColumn];
    1166               if (lower>-largeValue_||upper<largeValue_) {
    1167                 if (fabs(value-lower)<fabs(value-upper)) {
    1168                   setColumnStatus(iColumn,atLowerBound);
    1169                   columnActivityWork_[iColumn]=lower;
    1170                 } else {
    1171                   setColumnStatus(iColumn,atUpperBound);
    1172                   columnActivityWork_[iColumn]=upper;
    1173                 }
    1174               } else {
    1175                 setColumnStatus(iColumn,isFree);
    1176               }
    1177             } else {
    1178               setColumnStatus(iColumn,superBasic);
    1179             }
    1180           }
    1181         }
    1182         for (iRow=0;iRow<numberRows_;iRow++) {
    1183           int iSequence=pivotVariable_[iRow];
    1184           if (iSequence>=0) {
    1185             // basic
    1186             if (iSequence>=numberColumns_) {
    1187               // slack in - leave
    1188               //assert (iSequence-numberColumns_==iRow);
    1189             } else {
    1190               // put back structural
    1191               setColumnStatus(iSequence,basic);
    1192             }
    1193           } else {
    1194             // put in slack
    1195             setRowStatus(iRow,basic);
    1196           }
    1197         }
    1198         // signal repeat
    1199         status=-99;
    1200         // set fixed if they are
    1201         for (iRow=0;iRow<numberRows_;iRow++) {
    1202           if (getRowStatus(iRow)!=basic ) {
    1203             if (rowLowerWork_[iRow]==rowUpperWork_[iRow]) {
    1204               rowActivityWork_[iRow]=rowLowerWork_[iRow];
    1205               setRowStatus(iRow,isFixed);
    1206             }
    1207           }
    1208         }
    1209         for (iColumn=0;iColumn<numberColumns_;iColumn++) {
    1210           if (getColumnStatus(iColumn)!=basic ) {
    1211             if (columnLowerWork_[iColumn]==columnUpperWork_[iColumn]) {
    1212               columnActivityWork_[iColumn]=columnLowerWork_[iColumn];
    1213               setColumnStatus(iColumn,isFixed);
    1214             }
    1215           }
    1216         }
    1217       }
    1218     }
    1219   }
    1220   delete [] rowIsBasic;
    1221   delete [] columnIsBasic;
    1222 #else
     1087    for (i=0;i<numberRows_;i++)
     1088      printf ("xx %d %d\n",i,pivotVariable_[i]);
     1089    for (i=0;i<numberRows_+numberColumns_;i++)
     1090      if (getColumnStatus(i)==basic)
     1091        printf ("yy %d basic\n",i);
     1092    if (k>20)
     1093      exit(0);
     1094    k++;
     1095  }
    12231096  int status = factorization_->factorize(this, solveType,valuesPass);
    1224 #endif
    12251097  if (status) {
    12261098    handler_->message(CLP_SIMPLEX_BADFACTOR,messages_)
  • branches/pre/ClpSimplexPrimal.cpp

    r180 r190  
    12511251    double change = theta*alpha;
    12521252    value -= change;
    1253 
     1253    // But make sure one going out is feasible
    12541254    if (change>0.0) {
    12551255      // going down
    1256       if (value<=lower(iPivot)+primalTolerance_) {
     1256      if (value<lower(iPivot)+primalTolerance_) {
     1257        if (iPivot==sequenceOut_&&value>lower(iPivot)-1.001*primalTolerance_)
     1258          value=lower(iPivot);
    12571259        double difference = nonLinearCost_->setOne(iPivot,value);
    12581260        work[iRow] = difference;
     
    12671269    } else {
    12681270      // going up
    1269       if (value>=upper(iPivot)-primalTolerance_) {
     1271      if (value>upper(iPivot)-primalTolerance_) {
     1272        if (iPivot==sequenceOut_&&value<upper(iPivot)+1.001*primalTolerance_)
     1273          value=upper(iPivot);
    12701274        double difference = nonLinearCost_->setOne(iPivot,value);
    12711275        work[iRow] = difference;
     
    16851689             valueOut_<=upperValue+primalTolerance_);
    16861690      // may not be exactly at bound and bounds may have changed
     1691      // Make sure outgoing looks feasible
     1692      double useValue=valueOut_;
    16871693      if (valueOut_<=lowerValue+primalTolerance_) {
    16881694        directionOut_=1;
     1695        useValue= lower_[sequenceOut_];
    16891696      } else if (valueOut_>=upperValue-primalTolerance_) {
    16901697        directionOut_=-1;
     1698        useValue= upper_[sequenceOut_];
    16911699      } else {
    16921700        printf("*** variable wandered off bound %g %g %g!\n",
    16931701               lowerValue,valueOut_,upperValue);
    16941702        if (valueOut_-lowerValue<=upperValue-valueOut_) {
     1703          useValue= lower_[sequenceOut_];
    16951704          valueOut_=lowerValue;
    16961705          directionOut_=1;
    16971706        } else {
     1707          useValue= upper_[sequenceOut_];
    16981708          valueOut_=upperValue;
    16991709          directionOut_=-1;
     
    17011711      }
    17021712      solution_[sequenceOut_]=valueOut_;
    1703       nonLinearCost_->setOne(sequenceOut_,valueOut_);
     1713      nonLinearCost_->setOne(sequenceOut_,useValue);
    17041714    }
    17051715    // change cost and bounds on incoming if primal
  • branches/pre/ClpSimplexPrimalQuadratic.cpp

    r187 r190  
    563563      if (lastGoodIteration_==numberIterations_)
    564564        factorType=3;
    565       if (info->currentPhase())
    566         info->setCurrentSolution(solution_);
    567565      // may factorize, checks if problem finished
    568566      statusOfProblemInPrimal(lastCleaned,factorType,progress_,info);
    569       if (info->currentPhase())
    570         memcpy(solution_,info->currentSolution(),
    571                (numberRows_+numberColumns_)*sizeof(double));
    572567     
    573568      // Compute objective function from scratch
     
    610605            if ((value>lower+primalTolerance_&&value<upper-primalTolerance_)
    611606                ||getColumnStatus(iColumn)==superBasic) {
    612               if (getColumnStatus(iColumn)!=basic) {
     607              if (getColumnStatus(iColumn)!=basic&&!flagged(iColumn)) {
    613608                if (phase!=2) {
    614609                  phase=2;
     
    636631          if ((value>lower+primalTolerance_&&value<upper-primalTolerance_)
    637632            ||getColumnStatus(iColumn)==superBasic) {
    638             if (getColumnStatus(iColumn)!=basic) {
     633            if (getColumnStatus(iColumn)!=basic&&!flagged(iColumn)) {
    639634              if (phase!=2) {
    640635                phase=2;
     
    659654
    660655      // exit if victory declared
    661       if (!info->currentPhase()&&primalColumnPivot_->pivotColumn(rowArray_[1],
     656      if (!info->currentPhase()&&info->sequenceIn()<0&&primalColumnPivot_->pivotColumn(rowArray_[1],
    662657                                          rowArray_[2],rowArray_[3],
    663658                                          columnArray_[0],
     
    693688  checkComplementarity (info,rowArray_[0],rowArray_[1]);
    694689  int crucialSj = info->crucialSj();
     690  if (info->crucialSj()>=0) {
     691    printf("after inv %d Sj value %g\n",crucialSj,solution_[info->crucialSj()]);
     692  }
    695693  int returnCode=-1;
    696694  int phase = info->currentPhase();
    697695  double saveObjective = objectiveValue_;
    698696  int numberXColumns = info->numberXColumns();
     697  int numberXRows = info->numberXRows();
    699698  const int * which = info->quadraticSequence();
    700699  const double * pi = solution_+numberXColumns;
     
    706705  int nextSequenceIn=info->sequenceIn();
    707706  int oldSequenceIn=nextSequenceIn;
     707  int saveSequenceIn = sequenceIn_;
    708708  // status stays at -1 while iterating, >=0 finished, -2 to invert
    709709  // status -3 to go to top without an invert
     
    730730    // Only do if one not already chosen
    731731    bool cleanupIteration;
    732     if (phase==2) {
    733       // values pass
    734       if (nextSequenceIn<0) {
     732    if (nextSequenceIn<0) {
     733      cleanupIteration=false;
     734      if (phase==2) {
     735        // values pass
    735736        // get next
    736737        int iColumn;
     
    741742          double upper = upper_[iColumn];
    742743          if (value>lower+primalTolerance_&&value<upper-primalTolerance_) {
    743             if (getColumnStatus(iColumn)!=basic) {
     744            if (getColumnStatus(iColumn)!=basic&&!flagged(iColumn)) {
    744745              nextSequenceIn=iColumn;
    745746              break;
     
    756757            double upper = upper_[iColumn];
    757758            if (value>lower+primalTolerance_&&value<upper-primalTolerance_) {
    758               if (getColumnStatus(iColumn)!=basic) {
     759              if (getColumnStatus(iColumn)!=basic&&!flagged(iColumn)) {
    759760                nextSequenceIn=iColumn;
    760761                break;
     
    763764          }
    764765        }
    765       }
    766       oldSequenceIn=nextSequenceIn;
    767       sequenceIn_ = nextSequenceIn;
    768       cleanupIteration=false;
    769       lowerIn_=lower_[sequenceIn_];
    770       upperIn_=upper_[sequenceIn_];
    771       if (sequenceIn_<numberXColumns) {
    772         int jSequence = which[sequenceIn_];
    773         if (jSequence>=0) {
    774           dualIn_ = solution_[jSequence+numberRows_];
    775         } else {
    776           // need coding
    777           dualIn_=cost_[sequenceIn_];
    778           int j;
    779           for (j=columnStart[sequenceIn_];j<columnStart[sequenceIn_]+columnLength[sequenceIn_]; j++) {
    780             int iRow = row[j];
    781             dualIn_ -= element[j]*pi[iRow];
    782           }
    783         }
     766        oldSequenceIn=nextSequenceIn;
     767        saveSequenceIn=sequenceIn_;
     768        sequenceIn_ = nextSequenceIn;
    784769      } else {
    785         dualIn_ = solution_[sequenceIn_-numberColumns_+numberXColumns];
    786       }
    787       valueIn_=solution_[sequenceIn_];
    788       if (dualIn_>0.0)
    789         directionIn_ = -1;
    790       else
    791         directionIn_ = 1;
    792     } else {
    793       if (nextSequenceIn<0) {
     770        saveSequenceIn=sequenceIn_;
    794771        primalColumn(rowArray_[1],rowArray_[2],rowArray_[3],
    795772                     columnArray_[0],columnArray_[1]);
     773      }
     774    } else {
     775      saveSequenceIn=sequenceIn_;
     776      sequenceIn_ = nextSequenceIn;
     777      if (sequenceIn_<numberXColumns||sequenceIn_>=numberColumns_)
    796778        cleanupIteration=false;
    797       } else {
    798         sequenceIn_ = nextSequenceIn;
     779      else
    799780        cleanupIteration=true;
    800       }
    801781    }
    802782    pivotRow_=-1;
     
    811791      // do second half of iteration
    812792      while (chosen>=0) {
     793        {
     794          // Compute objective function from scratch
     795          const CoinPackedMatrix * quadratic =
     796            info->originalModel()->quadraticObjective();
     797          const int * columnQuadratic = quadratic->getIndices();
     798          const int * columnQuadraticStart = quadratic->getVectorStarts();
     799          const int * columnQuadraticLength = quadratic->getVectorLengths();
     800          const double * quadraticElement = quadratic->getElements();
     801          const double * originalCost = info->originalObjective();
     802          int iColumn;
     803          objectiveValue_=0.0;
     804          double infeasCost=0.0;
     805          double linearCost=0.0;
     806          for (iColumn=0;iColumn<numberXColumns;iColumn++) {
     807            double valueI = solution_[iColumn];
     808            linearCost += valueI*originalCost[iColumn];
     809            double diff =cost_[iColumn]-originalCost[iColumn];
     810            // WE are always feasible so look at low not up!
     811            if (diff>0.0) {
     812              double above = valueI-lower_[iColumn];
     813              assert(above>=0.0);
     814              infeasCost += diff*above;
     815            } else if (diff<0.0) {
     816              double below = upper_[iColumn]-valueI;
     817              assert(below>=0.0);
     818              infeasCost -= diff*below;
     819            }
     820            int j;
     821            for (j=columnQuadraticStart[iColumn];
     822                 j<columnQuadraticStart[iColumn]+columnQuadraticLength[iColumn];j++) {
     823              int jColumn = columnQuadratic[j];
     824              double valueJ = solution_[jColumn];
     825              double elementValue = 0.5*quadraticElement[j];
     826              objectiveValue_ += valueI*valueJ*elementValue;
     827            }
     828          }
     829          int iRow;
     830          for (iRow=0;iRow<numberXRows;iRow++) {
     831            int iColumn = iRow + numberColumns_;
     832            double diff =cost_[iColumn];
     833            double valueI = solution_[iColumn];
     834            if (diff>0.0) {
     835              double above = valueI-lower_[iColumn];
     836              assert(above>=0.0);
     837              infeasCost += diff*above;
     838            } else if (diff<0.0) {
     839              double below = upper_[iColumn]-valueI;
     840              assert(below>=0.0);
     841              infeasCost -= diff*below;
     842            }
     843          }
     844          objectiveValue_ += linearCost;
     845          assert (infeasCost>=0.0);
     846          printf("True objective is %g, infeas cost %g, sum %g\n",
     847                 objectiveValue_,infeasCost,objectiveValue_+infeasCost);
     848        }
    813849        objectiveValue_=saveObjective;
    814850        returnCode=-1;
     
    841877                iSjRow = iRow;
    842878                double d2 = solution_[crucialSj]/tj;
     879                if (fabs(solution_[crucialSj])<1.0e-8)
     880                  printf("zero crucialSj - pivot out - get right way\n");
     881                //if (sequenceIn_>=numberColumns_)
     882                //d2=-d2; // Try both ways? - makes no difference!!!
     883                //if (sequenceIn_>=numberXColumns&&sequenceIn_<numberColumns_)
     884                //d2=-d2; // Try both ways? - makes no difference!!!
    843885                // see which way to go
    844886                if (d2>0)
     
    865907            directionIn_ = 1;
    866908        } else {
     909          // not clean up
     910          lowerIn_=lower_[sequenceIn_];
     911          upperIn_=upper_[sequenceIn_];
     912          if (sequenceIn_<numberXColumns) {
     913            int jSequence = which[sequenceIn_];
     914            if (jSequence>=0) {
     915              dualIn_ = solution_[jSequence+numberRows_];
     916            } else {
     917              // need coding
     918              dualIn_=cost_[sequenceIn_];
     919              int j;
     920              for (j=columnStart[sequenceIn_];j<columnStart[sequenceIn_]+columnLength[sequenceIn_]; j++) {
     921                int iRow = row[j];
     922                dualIn_ -= element[j]*pi[iRow];
     923              }
     924            }
     925          } else {
     926            dualIn_ = solution_[sequenceIn_-numberColumns_+numberXColumns];
     927          }
     928          valueIn_=solution_[sequenceIn_];
     929          if (dualIn_>0.0)
     930            directionIn_ = -1;
     931          else
     932            directionIn_ = 1;
    867933          if (sequenceIn_<numberColumns_) {
    868934            // Set dj as value of slack
     
    878944          info->setCrucialSj(crucialSj);
    879945        }
     946        double oldSj=1.0e30;
     947        if (info->crucialSj()>=0&&cleanupIteration)
     948          oldSj= solution_[info->crucialSj()];
    880949        // save reduced cost
    881950        //double saveDj = dualIn_;
     
    886955                             info,
    887956                             cleanupIteration);
     957        if (pivotRow_==-1&&phase==2&&fabs(dualIn_)<1.0e-3) {
     958          // try other way
     959          dualIn_=-dualIn_;
     960          directionIn_=-directionIn_;
     961          if (info->crucialSj()>=0)
     962            setColumnStatus(info->crucialSj(),basic);
     963          rowArray_[3]->checkClear();
     964          rowArray_[2]->checkClear();
     965          rowArray_[0]->checkClear();
     966          result=primalRow(rowArray_[1],rowArray_[3],
     967                           rowArray_[2],rowArray_[0],
     968                           info,
     969                           cleanupIteration);
     970        }
    888971        saveObjective = objectiveValue_;
    889972        if (pivotRow_>=0) {
     
    913996          } else if (updateStatus==2) {
    914997            // major error
     998            // Reset sequenceIn_
     999            // sequenceIn_=saveSequenceIn;
     1000            pivotRow_=-1;
    9151001            // better to have small tolerance even if slower
    9161002            factorization_->zeroTolerance(1.0e-15);
     
    10201106                 valueOut_<=upperValue+primalTolerance_);
    10211107          // may not be exactly at bound and bounds may have changed
     1108          // Make sure outgoing looks feasible
     1109          double useValue=valueOut_;
    10221110          if (valueOut_<=lowerValue+primalTolerance_) {
    10231111            directionOut_=1;
     1112            useValue= lower_[sequenceOut_];
    10241113          } else if (valueOut_>=upperValue-primalTolerance_) {
    10251114            directionOut_=-1;
     1115            useValue= upper_[sequenceOut_];
    10261116          } else {
    10271117            printf("*** variable wandered off bound %g %g %g!\n",
    10281118                   lowerValue,valueOut_,upperValue);
    10291119            if (valueOut_-lowerValue<=upperValue-valueOut_) {
     1120              useValue= lower_[sequenceOut_];
    10301121              valueOut_=lowerValue;
    10311122              directionOut_=1;
    10321123            } else {
     1124              useValue= upper_[sequenceOut_];
    10331125              valueOut_=upperValue;
    10341126              directionOut_=-1;
     
    10361128          }
    10371129          solution_[sequenceOut_]=valueOut_;
    1038           nonLinearCost_->setOne(sequenceOut_,valueOut_);
     1130          nonLinearCost_->setOne(sequenceOut_,useValue);
    10391131        }
    10401132        // change cost and bounds on incoming if primal
    10411133        nonLinearCost_->setOne(sequenceIn_,valueIn_);
    10421134        int whatNext=housekeeping(objectiveChange);
     1135        if (info->crucialSj()>=0) {
     1136          assert(fabs(oldSj)>= fabs(solution_[info->crucialSj()]));
     1137          printf("%d Sj value went from %g to %g\n",crucialSj,oldSj,solution_[info->crucialSj()]);
     1138        }
     1139        {
     1140          // Compute objective function from scratch
     1141          const CoinPackedMatrix * quadratic =
     1142            info->originalModel()->quadraticObjective();
     1143          const int * columnQuadratic = quadratic->getIndices();
     1144          const int * columnQuadraticStart = quadratic->getVectorStarts();
     1145          const int * columnQuadraticLength = quadratic->getVectorLengths();
     1146          const double * quadraticElement = quadratic->getElements();
     1147          const double * originalCost = info->originalObjective();
     1148          int iColumn;
     1149          objectiveValue_=0.0;
     1150          double infeasCost=0.0;
     1151          double linearCost=0.0;
     1152          for (iColumn=0;iColumn<numberXColumns;iColumn++) {
     1153            double valueI = solution_[iColumn];
     1154            linearCost += valueI*originalCost[iColumn];
     1155            double diff =cost_[iColumn]-originalCost[iColumn];
     1156            // WE are always feasible so look at low not up!
     1157            if (diff>0.0) {
     1158              double above = valueI-lower_[iColumn];
     1159              assert(above>=0.0);
     1160              infeasCost += diff*above;
     1161            } else if (diff<0.0) {
     1162              double below = upper_[iColumn]-valueI;
     1163              assert(below>=0.0);
     1164              infeasCost -= diff*below;
     1165            }
     1166            int j;
     1167            for (j=columnQuadraticStart[iColumn];
     1168                 j<columnQuadraticStart[iColumn]+columnQuadraticLength[iColumn];j++) {
     1169              int jColumn = columnQuadratic[j];
     1170              double valueJ = solution_[jColumn];
     1171              double elementValue = 0.5*quadraticElement[j];
     1172              objectiveValue_ += valueI*valueJ*elementValue;
     1173            }
     1174          }
     1175          int iRow;
     1176          for (iRow=0;iRow<numberXRows;iRow++) {
     1177            int iColumn = iRow + numberColumns_;
     1178            double diff =cost_[iColumn];
     1179            double valueI = solution_[iColumn];
     1180            if (diff>0.0) {
     1181              double above = valueI-lower_[iColumn];
     1182              assert(above>=0.0);
     1183              infeasCost += diff*above;
     1184            } else if (diff<0.0) {
     1185              double below = upper_[iColumn]-valueI;
     1186              assert(below>=0.0);
     1187              infeasCost -= diff*below;
     1188            }
     1189          }
     1190          objectiveValue_ += linearCost;
     1191          assert (infeasCost>=0.0);
     1192          printf("After Housekeeping True objective is %g, infeas cost %g, sum %g\n",
     1193                 objectiveValue_,infeasCost,objectiveValue_+infeasCost);
     1194        }
    10431195        if (sequenceOut_>=numberXColumns&&sequenceOut_<numberColumns_) {
    10441196          // really free but going to zero
    10451197          setStatus(sequenceOut_,isFree);
     1198          if (sequenceOut_==info->crucialSj())
     1199            info->setCrucialSj(-1);
    10461200        }
    10471201        checkComplementarity (info,rowArray_[2],rowArray_[3]);
     
    10601214        cleanupIteration = true;
    10611215        // may not be correct on second time
    1062         if (result) {
     1216        if (result&&(returnCode==-1||returnCode==-2)) {
    10631217          assert(sequenceOut_<numberXColumns||
    10641218                 sequenceOut_>=numberColumns_);
     1219          int crucialSj=info->crucialSj();
    10651220          if (sequenceOut_<numberXColumns) {
    10661221            chosen =sequenceOut_ + numberRows_; // sj variable
     
    10701225            if (iRow<numberRows_-numberXColumns) {
    10711226              int iPi = iRow+numberXColumns;
    1072               printf("pi for row %d is %g\n",
    1073                      iRow,solution_[iPi]);
     1227              //printf("pi for row %d is %g\n",
     1228              //     iRow,solution_[iPi]);
    10741229              chosen=iPi;
    10751230            } else {
     
    10771232              abort();
    10781233            }
     1234          }
     1235          // But double check as original incoming might have gone out
     1236          if (chosen==crucialSj) {
     1237            chosen=-1;
     1238            break; // means original has gone out
     1239          }
     1240          if (returnCode==-2) {
     1241            // factorization
     1242            nextSequenceIn = chosen;
     1243            break;
    10791244          }
    10801245        } else {
     
    12121377  // sj
    12131378  int iSjRow=-1;
     1379  // sj for incoming
     1380  int iSjRow2=-1,crucialSj2=-1;
     1381  if (sequenceIn_<numberColumns_) {
     1382    const int * which = info->quadraticSequence();
     1383    crucialSj2 = which[sequenceIn_];
     1384    if (crucialSj2>=0)
     1385      crucialSj2 += numberRows_; // sj2 which should go to 0.0
     1386  } else {
     1387    crucialSj2 = sequenceIn_-numberColumns_+numberXColumns; // pi which should go to 0.0
     1388  }
    12141389  double tj = 0.0;
    12151390  // Change in objective will be theta*coeff1 + theta*theta*coeff2
     
    12341409      }
    12351410    }
     1411    if (iPivot==crucialSj2)
     1412      iSjRow2=iRow;
    12361413  }
    12371414  // Incoming
     
    12391416    index2[number2++]=sequenceIn_;
    12401417    rhs[sequenceIn_]=way;
    1241     printf("incoming col %d alpha %g solution %g\n",sequenceIn_,way,valueIn_);
     1418    //printf("incoming col %d alpha %g solution %g\n",sequenceIn_,way,valueIn_);
    12421419    coeff1 += way*cost_[sequenceIn_];
    12431420  } else {
    1244     printf("incoming new col %d alpha %g solution %g\n",sequenceIn_,way,valueIn_);
     1421    //printf("incoming new col %d alpha %g solution %g\n",sequenceIn_,way,valueIn_);
    12451422    coeffSlack += way*cost_[sequenceIn_];
    12461423  }
    12471424  rhsArray->setNumElements(number2);
    1248   if (numberIterations_>=0||cleanupIteration) {
     1425  if (numberIterations_!=-701) {  //if (numberIterations_>=0||cleanupIteration) {
    12491426    for (iIndex=0;iIndex<number2;iIndex++) {
    12501427      int iColumn=index2[iIndex];
     
    13071484  }
    13081485  coeff2 *= 0.5;
    1309   printf("coefficients %g %g %g - dualIn %g\n",coeff1,coeff2,coeffSlack,dualIn_);
    13101486  if (coeffSlack)
    13111487    printf("trouble\n");
     
    13141490  int accuracyFlag=0;
    13151491  if (!cleanupIteration) {
    1316     assert (fabs(way*coeff1-dualIn_)<1.0e-4*(1.0+fabs(dualIn_)));
    1317     assert (way*coeff1*dualIn_>=0.0);
     1492    //assert (fabs(way*coeff1-dualIn_)<1.0e-4*(1.0+fabs(dualIn_)));
     1493    //assert (way*coeff1*dualIn_>=0.0);
    13181494    if (way*coeff1*dualIn_<0.0) {
    1319       // bad
    1320       accuracyFlag=2;
    1321     } else if (fabs(way*coeff1-dualIn_)>1.0e-2*(1.0e-3+fabs(dualIn_))) {
    13221495      // bad
    13231496      accuracyFlag=2;
     
    13261499      accuracyFlag=1;
    13271500    }
    1328     dualIn_ = way*coeff1;
    13291501    if (crucialSj>=0) {
    1330       solution_[crucialSj]=dualIn_;
     1502      solution_[crucialSj]=way*coeff1;
    13311503    }
    13321504  } else if (dualIn_<0.0&&way*coeff1>1.0e-2) {
     
    13451517  else
    13461518    d1 = 0.0;
    1347   if (fabs(tj)<1.0e-7) {
    1348     if (sequenceIn_<numberXColumns)
    1349       printf("column %d is basically linear\n",sequenceIn_);
     1519  if (fabs(tj)<1.0e-7||!cleanupIteration) {
     1520    //if (sequenceIn_<numberXColumns)
     1521      //printf("column %d is basically linear\n",sequenceIn_);
    13501522    //assert(!columnQuadraticLength[sequenceIn_]);
    13511523  } else {
     
    13561528    }
    13571529  }
    1358   printf("derivative zero at %g, sj zero at %g\n",d1,d2);
    13591530  if (d1>1.0e10&&d2>1.0e10) {
    13601531    // linear variable entering
    13611532    // maybe we should have done dual iteration to force sj to 0.0
    1362     printf("linear variable\n");
    1363   }
    1364   maximumMovement = min(maximumMovement,d2);
    1365   if (d1>=0.0) {
    1366     maximumMovement = min(maximumMovement,d1);
    1367     d2 = min(d1,d2);
    1368   }
    1369    
     1533    //printf("linear variable\n");
     1534  }
     1535  handler_->message(CLP_QUADRATIC_PRIMAL_DETAILS,messages_)
     1536    <<coeff1<<coeff2<<coeffSlack<<dualIn_<<d1<<d2
     1537    <<CoinMessageEol;
     1538  if (!cleanupIteration)
     1539    dualIn_ = way*coeff1;
     1540  double mind1d2=1.0e50;
     1541  if (cleanupIteration) {
     1542    // we are only interested in d1 if smaller than d2
     1543    mind1d2 = min(maximumMovement,d2);
     1544    //if (d1>1.0e-8&&d1<0.999*d2)
     1545    //mind1d2=d1;
     1546  } else {
     1547    // There is no d2 - d1 refers to crucialSj
     1548    if (d1>1.0e-5) {
     1549      mind1d2 = min(maximumMovement,d1);
     1550      //assert (d1>=0.0);
     1551    }
     1552  }
     1553  maximumMovement = min(maximumMovement,mind1d2);
    13701554  rhsArray->clear();
    13711555  double tentativeTheta = maximumMovement;
    13721556  double upperTheta = maximumMovement;
    1373 
     1557  bool throwAway=false;
     1558  if (numberIterations_==700||numberIterations_==694) {
     1559    printf("Bad iteration coming up after iteration %d\n",numberIterations_);
     1560  }
    13741561
    13751562  for (iIndex=0;iIndex<number;iIndex++) {
     
    14581645      maximumSwapped = max(maximumSwapped,numberSwapped);
    14591646
    1460       double dualCheck = - 2.0*coeff2*tentativeTheta - coeff1 - 99999999;
     1647      double dualCheck = - 2.0*coeff2*tentativeTheta - coeff1 - 0.1*infeasibilityCost_;
    14611648      // but make a bit more pessimistic
    14621649      dualCheck=max(dualCheck-100.0*dualTolerance_,0.99*dualCheck);
     
    14751662        memset(spare+numberRemaining,0,
    14761663               (saveNumber-numberRemaining)*sizeof(double));
     1664        double lastThru = totalThru;
    14771665        int iTry;
    14781666#define MAXTRY 100
     
    15391727            break;
    15401728          } else {
    1541             dualCheck = - 2.0*coeff2*theta_ - coeff1-9999999;
     1729            dualCheck = - 2.0*coeff2*theta_ - coeff1;
    15421730            //dualCheck =0.0;
    1543             if (totalThru>=dualCheck) {
     1731            if (totalThru>=dualCheck&&(sequenceIn_<numberXColumns||sequenceIn_>=numberColumns_)) {
     1732              if (!cleanupIteration) {
     1733                // We can pivot out sj
     1734                if (upperTheta==maximumMovement) {
     1735                  printf("figures\n");
     1736                } else {
     1737                  if (lastThru>=dualCheck)
     1738                    printf("totalThru %g, lastThru %g - taking sj to zero?\n",totalThru,lastThru);
     1739                }
     1740              } else {
     1741                if (pivotRow_<0&&lastPivotRow<0)
     1742                  throwAway=true;
     1743              }
     1744              break; // no point trying another loop
     1745            } else if (totalThru>=dualCheck||upperTheta==maximumMovement) {
    15441746              break; // no point trying another loop
    15451747            } else {
     
    15541756              if (bestPivot>bestEverPivot)
    15551757                bestEverPivot=bestPivot;
     1758              lastThru = totalThru;
    15561759            }
    15571760          }
     
    16371840    dualOut_ = reducedCost(sequenceOut_);
    16381841  } else {
     1842    // If no pivot but bad move then throw away
     1843    if (throwAway&&(sequenceIn_<numberXColumns||sequenceIn_>=numberColumns_)) {
     1844      assert (pivotRow_<0);
     1845      accuracyFlag=2;
     1846    }
    16391847    double trueMaximumMovement;
    16401848    if (way>0.0)
     
    16631871      }
    16641872      // we may still have sj to get rid of
    1665     } else if (fabs(maximumMovement-d2)<dualTolerance_) {
    1666       // sj going to zero
    1667       result=0;
     1873    } else if (fabs(maximumMovement-mind1d2)<dualTolerance_) {
     1874      // crucialSj local copy i.e. dj going to zero
     1875      if (maximumMovement==d2) {
     1876        // sj going to zero
     1877        result=0;
     1878      } else {
     1879        // incoming dj going to zero
     1880        assert(!cleanupIteration);
     1881        if (!cleanupIteration)
     1882          result=0;
     1883        iSjRow=iSjRow2;
     1884        crucialSj = pivotVariable_[iSjRow];
     1885      }
    16681886      assert (pivotRow_<0);
    16691887      setColumnStatus(crucialSj,isFree);
     
    16721890      // translate to sequence
    16731891      sequenceOut_ = pivotVariable_[pivotRow_];
     1892      assert (sequenceOut_==crucialSj);
    16741893      valueOut_ = solution(sequenceOut_);
    1675       theta_ = d2;
     1894      theta_ = maximumMovement;
    16761895      if (way<0.0)
    16771896        theta_ = - theta_;
     
    16971916  // Change in objective will be theta*coeff1 + theta*theta*coeff2
    16981917  objectiveValue_ += theta_*coeff1+theta_*theta_*coeff2;
    1699   printf("New objective value %g\n",objectiveValue_);
    17001918  {
    17011919    int iColumn;
     
    17221940      }
    17231941    }
    1724     printf("Objective value %g\n",objectiveValue());
    17251942  }
    17261943  if (accuracyFlag<2) {
     
    17371954                                                   ClpQuadraticInfo * info)
    17381955{
     1956  if (info->currentPhase())
     1957    info->setCurrentSolution(solution_);
     1958  //info->saveStatus();
     1959 
    17391960  if (type==2) {
    17401961    // trouble - restore solution
     
    17611982    if (type) {
    17621983      // is factorization okay?
     1984      int numberPivots = factorization_->pivots();
    17631985      if (internalFactorize(1)) {
    17641986        if (solveType_==2) {
     
    17671989          return;
    17681990        }
    1769 
     1991        numberIterations_ = progress_->lastIterationNumber(0);
    17701992        // no - restore previous basis
    17711993        assert (type==1);
     
    17792001        assert (internalFactorize(1)==0);
    17802002        info->restoreStatus();
     2003        // flag incoming
     2004        if (numberPivots==1) {
     2005          int crucialSj = info->crucialSj();
     2006          if (crucialSj<0) {
     2007            setFlagged(sequenceIn_);
     2008          } else {
     2009            int iSequence;
     2010            int numberXColumns = info->numberXColumns();
     2011            if (crucialSj<numberRows_) {
     2012              // row
     2013              iSequence = crucialSj-numberXColumns+numberColumns_;
     2014            } else {
     2015              // column
     2016              iSequence = crucialSj-numberRows_;
     2017            }
     2018            if (getColumnStatus(iSequence)!=basic) {
     2019              setFlagged(iSequence);
     2020              info->setSequenceIn(-1);
     2021            } else {
     2022              printf("**** %d is basic ?\n",iSequence);
     2023            }
     2024          }
     2025        }
    17812026        changeMade_++; // say change made
    17822027      }
     
    17892034  // put back original costs and then check
    17902035  createRim(4);
     2036  if (info->currentPhase()) {
     2037    memcpy(solution_,info->currentSolution(),
     2038           (numberRows_+numberColumns_)*sizeof(double));
     2039  }
    17912040  gutsOfSolution(NULL,NULL);
     2041  if (info->currentPhase()) {
     2042    memcpy(solution_,info->currentSolution(),
     2043           (numberRows_+numberColumns_)*sizeof(double));
     2044    nonLinearCost_->checkInfeasibilities(primalTolerance_);
     2045  }
    17922046  checkComplementarity(info,rowArray_[2],rowArray_[3]);
    17932047  // Double check reduced costs if no action
     
    17992053  }
    18002054  // Check if looping
     2055  // Take out for now
    18012056  int loop;
    1802   if (type!=2)
     2057  if (type!=2&&0)
    18032058    loop = progress->looping();
    18042059  else
     
    18522107      // put back original costs
    18532108      createRim(4);
     2109      // put all non-basic variables to bounds
     2110      int iSequence;
     2111      for (iSequence=0;iSequence<numberRows_+numberColumns_;iSequence++) {
     2112        Status status = getStatus(iSequence);
     2113        if (status==atLowerBound||status==isFixed)
     2114          solution_[iSequence]=lower_[iSequence];
     2115        else if (status==atUpperBound)
     2116          solution_[iSequence]=upper_[iSequence];
     2117      }
    18542118      nonLinearCost_->checkInfeasibilities(primalTolerance_);
    18552119      gutsOfSolution(NULL,NULL);
     2120      nonLinearCost_->checkInfeasibilities(primalTolerance_);
    18562121      checkComplementarity(info,rowArray_[2],rowArray_[3]);
    18572122
     
    18652130        problemStatus_ = -1;
    18662131        infeasibilityCost_=saveWeight;
     2132        nonLinearCost_->checkInfeasibilities(primalTolerance_);
     2133        checkComplementarity(info,rowArray_[2],rowArray_[3]);
    18672134      } else {
    18682135        nonLinearCost_=NULL;
     
    18882155          infeasibilityCost_=0.0;
    18892156          createRim(4);
     2157          // put all non-basic variables to bounds
     2158          int iSequence;
     2159          for (iSequence=0;iSequence<numberRows_+numberColumns_;iSequence++) {
     2160            Status status = getStatus(iSequence);
     2161            if (status==atLowerBound||status==isFixed)
     2162              solution_[iSequence]=lower_[iSequence];
     2163            else if (status==atUpperBound)
     2164              solution_[iSequence]=upper_[iSequence];
     2165          }
    18902166          nonLinearCost_->checkInfeasibilities(primalTolerance_);
    18912167          gutsOfSolution(NULL,NULL);
     
    19402216          // put back original costs and then check
    19412217          createRim(4);
     2218          // put all non-basic variables to bounds
     2219          int iSequence;
     2220          for (iSequence=0;iSequence<numberRows_+numberColumns_;iSequence++) {
     2221            Status status = getStatus(iSequence);
     2222            if (status==atLowerBound||status==isFixed)
     2223              solution_[iSequence]=lower_[iSequence];
     2224            else if (status==atUpperBound)
     2225              solution_[iSequence]=upper_[iSequence];
     2226          }
    19422227          nonLinearCost_->checkInfeasibilities(oldTolerance);
    19432228          gutsOfSolution(NULL,NULL);
     
    20542339    // zero out basic values
    20552340    int iRow;
    2056     double * save = new double [numberRows_];
    20572341    for (iRow=0;iRow<numberRows_;iRow++) {
    20582342      int iPivot=pivotVariable_[iRow];
    2059       save[iRow]=solution_[iPivot];
     2343      // save
     2344      dj_[iPivot]=solution_[iPivot];
    20602345      solution_[iPivot] = 0.0;
    20612346    }
     
    21032388      int iSequence  = iRow + numberColumns_;
    21042389      double value = cost_[iSequence];
    2105       if (value) {
    2106         assert(getRowStatus(iRow)==basic);
     2390      // Its possible a fixed vector may have had this set
     2391      if (value&&getRowStatus(iRow)==basic) {
    21072392        // add in pi column
    21082393        matrix_->add(this,array1,iRow+numberXColumns,value);
     
    21162401      int iRow = index[k];
    21172402      double value=modifiedCost[iRow];
    2118       int iSequence = pivotVariable_[iRow];
    2119       if (iSequence>=numberXColumns&&iSequence<numberColumns_)
    2120         solution_[iSequence] = value;
     2403      int iPivot = pivotVariable_[iRow];
     2404      if (iPivot>=numberXColumns&&iPivot<numberColumns_)
     2405        solution_[iPivot] = value;
    21212406      else
    2122         solution_[iSequence] = save[iRow];
    2123       //solution_[iSequence] = value;
     2407        solution_[iPivot] = dj_[iPivot];
     2408      //solution_[iPivot] = value;
    21242409     
    21252410    }
    2126     delete [] save;
    21272411    array1->clear();
    21282412    // Now modify pi values by slack costs to make djs
     
    21442428    for (k=0;k<numberRows_;k++) {
    21452429      if (k<numberXRows) {
    2146         if (fabs(modifiedCost[k]-rowSolution[k])>1.0e-3*(1.0+fabs(rowSolution[k])))
    2147           printf("bad row %d %g %g\n",k,modifiedCost[k],rowSolution[k]);
     2430        //if (fabs(modifiedCost[k]-rowSolution[k])>1.0e-3*(1.0+fabs(rowSolution[k])))
     2431        //printf("bad row %d %g %g\n",k,modifiedCost[k],rowSolution[k]);
    21482432      } else {
    21492433        rowSolution[k]=modifiedCost[k];
     
    22042488        jSequence += start;
    22052489        if (getColumnStatus(jSequence)==basic)
    2206           printf("Struct %d (%g) and %d (%g) both basic\n",
    2207                  iSequence,solution_[iSequence],jSequence,solution_[jSequence]);
     2490          handler_->message(CLP_QUADRATIC_BOTH,messages_)
     2491            <<"Structural"<<iSequence
     2492            <<solution_[iSequence]<<jSequence<<solution_[jSequence]
     2493            <<CoinMessageEol;
    22082494      }
    22092495      break;
     
    28503136  int numberRows = numberXRows_+numberQuadraticColumns_;
    28513137  int size = numberRows+numberColumns;
    2852   if (!validSolution_)
    2853     validSolution_ = new double [size];
    2854   memcpy(validSolution_,currentSolution_,size*sizeof(double));
     3138  if (currentSolution_) {
     3139    if (!validSolution_)
     3140      validSolution_ = new double [size];
     3141    memcpy(validSolution_,currentSolution_,size*sizeof(double));
     3142  }
    28553143}
    28563144// Restore previous
  • branches/pre/include/ClpMessage.hpp

    r180 r190  
    6161  CLP_CRASH,
    6262  CLP_END_VALUES_PASS,
     63  CLP_QUADRATIC_BOTH,
     64  CLP_QUADRATIC_PRIMAL_DETAILS,
    6365  CLP_DUMMY_END
    6466};
Note: See TracChangeset for help on using the changeset viewer.