Changeset 201 for branches


Ignore:
Timestamp:
Aug 27, 2003 4:06:01 AM (16 years ago)
Author:
forrest
Message:

Still trying

Location:
branches/pre
Files:
3 edited

Legend:

Unmodified
Added
Removed
  • branches/pre/ClpSimplexPrimalQuadratic.cpp

    r200 r201  
    659659        phase=0;
    660660        int nextSequenceIn=-1;
     661        int numberQuadraticRows = info->numberQuadraticRows();
     662        const int * whichColumn = info->quadraticSequence();
    661663        for (iColumn=0;iColumn<numberXColumns;iColumn++) {
    662664          double value = solution_[iColumn];
     
    666668              ||getColumnStatus(iColumn)==superBasic) {
    667669            if (getColumnStatus(iColumn)!=basic&&!flagged(iColumn)) {
    668               if (phase!=2) {
    669                 phase=2;
    670                 nextSequenceIn=iColumn;
     670              int jSequence=whichColumn[iColumn];
     671              if (fabs(dj_[iColumn])>10.0*dualTolerance_||
     672                       (jSequence>=0&&
     673                        getColumnStatus(jSequence+numberXColumns+numberQuadraticRows)==basic)) {
     674                if (phase!=2) {
     675                  phase=2;
     676                  nextSequenceIn=iColumn;
     677                }
    671678              }
    672679            }
    673680          }
     681#if 0
    674682          if (which[iColumn]>=0) {
    675683            if (getColumnStatus(iColumn)==basic) {
     
    684692            }
    685693          }
     694#endif
    686695        }
    687696        const int * whichRow = info->quadraticRowSequence();
     
    692701          if ((value>lower+primalTolerance_&&value<upper-primalTolerance_)
    693702            ||getColumnStatus(iColumn)==superBasic) {
    694             if (getColumnStatus(iColumn)!=basic&&!flagged(iColumn)) {
     703            if (getColumnStatus(iColumn)!=basic&&!flagged(iColumn)&&fabs(dj_[iColumn])>10.0*dualTolerance_) {
    695704              if (phase!=2) {
    696705                phase=2;
     
    699708            }
    700709          }
     710#if 0
     711          // may work later when whichRow working??
    701712          int iRow = iColumn-numberColumns_;
    702713          if (whichRow[iRow]>=0) {
     
    712723            }
    713724          }
     725#endif
    714726        }
    715727        info->setSequenceIn(nextSequenceIn);
     
    779791  // status -3 to go to top without an invert
    780792  while (problemStatus_==-1) {
     793    if (info->crucialSj()<0&&factorization_->pivots()>=10) {
     794      returnCode =-2; // refactorize
     795      break;
     796    }
    781797#ifdef CLP_DEBUG
    782798    {
     
    799815    // Initially Dantzig and look at s variables
    800816    // Only do if one not already chosen
    801     bool cleanupIteration;
     817    int cleanupIteration;
    802818    if (nextSequenceIn<0) {
    803       cleanupIteration=false;
     819      cleanupIteration=0;
    804820      if (phase==2) {
    805821        // values pass
     
    813829          if (value>lower+primalTolerance_&&value<upper-primalTolerance_) {
    814830            if (getColumnStatus(iColumn)!=basic&&!flagged(iColumn)) {
    815               nextSequenceIn=iColumn;
    816               break;
     831              int jSequence=whichColumn[iColumn];
     832              if (fabs(dj_[iColumn])>10.0*dualTolerance_||
     833                       (jSequence>=0&&
     834                        getColumnStatus(jSequence+numberXColumns+numberQuadraticRows)==basic)) {
     835                nextSequenceIn=iColumn;
     836                break;
     837              }
    817838            }
    818839          }
     
    821842          iStart=max(iStart,numberColumns_);
    822843          int numberXRows = info->numberXRows();
    823           for (iColumn=numberColumns_;iColumn<numberColumns_+numberXRows;
     844          for (iColumn=iStart;iColumn<numberColumns_+numberXRows;
    824845               iColumn++) {
    825846            double value = solution_[iColumn];
     
    827848            double upper = upper_[iColumn];
    828849            if (value>lower+primalTolerance_&&value<upper-primalTolerance_) {
    829               if (getColumnStatus(iColumn)!=basic&&!flagged(iColumn)) {
     850              if (getColumnStatus(iColumn)!=basic&&!flagged(iColumn)&&fabs(dj_[iColumn])>10.0*dualTolerance_) {
    830851                nextSequenceIn=iColumn;
    831852                break;
     
    849870      sequenceIn_ = nextSequenceIn;
    850871      if (phase==2) {
    851         if (sequenceIn_<numberXColumns||sequenceIn_>=numberColumns_)
    852           cleanupIteration=false;
     872        if ((sequenceIn_<numberXColumns||sequenceIn_>=numberColumns_)&&info->crucialSj()<0)
     873          cleanupIteration=0;
    853874        else
    854           cleanupIteration=true;
     875          cleanupIteration=1;
    855876      } else {
    856         cleanupIteration=true;
    857       }
     877        cleanupIteration=1;
     878      }
     879      if ((sequenceIn_>=numberXColumns&&sequenceIn_<numberColumns_)&&info->crucialSj()<0)
     880        cleanupIteration=2;
    858881    }
    859882    pivotRow_=-1;
     
    866889      // do second half of iteration
    867890      while (chosen>=0) {
     891        int saveSequenceInInfo=chosen;
     892        int saveCrucialSjInfo=info->crucialSj();
    868893        rowArray_[1]->clear();
    869894        checkComplementarity (info,rowArray_[3],rowArray_[1]);
     
    911936          // should keep pivot row of crucialSj as well (speed)
    912937          int iSjRow=-1;
    913           {
     938          if (crucialSj>=0) {
    914939            double * work=rowArray_[1]->denseVector();
    915940            int number=rowArray_[1]->getNumElements();
     
    945970            //assert(tj);
    946971            }
     972          } else {
     973            assert (cleanupIteration==2);
     974            // get correct dj
     975            int iSequence = sequenceIn_-numberXColumns;
     976            if (iSequence<numberQuadraticRows) {
     977              int iRow=backRow[iSequence];
     978              dj_[sequenceIn_] = - dj_[numberColumns_+iRow];
     979            } else {
     980              iSequence = backColumn[iSequence-numberQuadraticRows];
     981              dj_[sequenceIn_]= - dj_[iSequence];
     982            }
    947983          }
    948984
     
    962998            if (jSequence>=0) {
    963999              dualIn_ = solution_[jSequence+numberXColumns+numberQuadraticRows];
     1000              dualIn_=dj_[sequenceIn_];
    9641001            } else {
    9651002              // need coding
     
    9741011            int iRow = sequenceIn_-numberColumns_;
    9751012            assert (iRow>=0);
    976             if (whichRow[iRow]>=0)
     1013            if (whichRow[iRow]>=0) {
    9771014              dualIn_ = solution_[numberXColumns+whichRow[iRow]];
     1015              dualIn_=dj_[sequenceIn_];
     1016              // temporary?
     1017              const int * columnLength = matrix_->getVectorLengths();
     1018              if (!columnLength[numberXColumns+whichRow[iRow]])
     1019                dualIn_ = dual_[iRow]+cost_[sequenceIn_];
     1020            }
    9781021          }
    9791022          valueIn_=solution_[sequenceIn_];
     
    9941037            else
    9951038              crucialSj=-1;
     1039            // temporary?
     1040            const int * columnLength = matrix_->getVectorLengths();
     1041            if (!columnLength[numberXColumns+iRow])
     1042              crucialSj=-1;
    9961043          }
     1044          if (crucialSj>=0&&getColumnStatus(crucialSj)!=basic)
     1045            crucialSj=-1;
    9971046          info->setCrucialSj(crucialSj);
    9981047        }
    9991048        double oldSj=1.0e30;
     1049        double oldObjectiveValue = objectiveValue_;
    10001050        if (info->crucialSj()>=0&&cleanupIteration)
    10011051          oldSj= solution_[info->crucialSj()];
     
    10481098            // Reset sequenceIn_
    10491099            // sequenceIn_=saveSequenceIn;
     1100            nextSequenceIn=saveSequenceInInfo;
     1101            info->setCrucialSj(saveCrucialSjInfo);
    10501102            pivotRow_=-1;
    10511103            // better to have small tolerance even if slower
     
    11811233        }
    11821234        // change cost and bounds on incoming if primal
     1235        if (sequenceIn_==sequenceOut_&&(lowerOut_!=lowerIn_||upperOut_!=upperIn_)) {
     1236          // linear variable superbasic
     1237          setStatus(sequenceOut_,superBasic);
     1238        }
    11831239        nonLinearCost_->setOne(sequenceIn_,valueIn_);
    11841240        int whatNext=housekeeping(objectiveChange);
     1241        // and again as housekeeping may have changed
     1242        if (sequenceIn_==sequenceOut_&&(lowerOut_!=lowerIn_||upperOut_!=upperIn_)) {
     1243          // linear variable superbasic
     1244          setStatus(sequenceOut_,superBasic);
     1245        }
    11851246        if (info->crucialSj()>=0) {
    11861247          assert(fabs(oldSj)>= fabs(solution_[info->crucialSj()]));
    11871248          printf("%d Sj value went from %g to %g\n",crucialSj,oldSj,solution_[info->crucialSj()]);
    11881249        }
    1189         double oldObjectiveValue = objectiveValue_;
    11901250        checkComplementarity (info,rowArray_[2],rowArray_[3]);
    11911251        printf("After Housekeeping True objective is %g, infeas cost %g, sum %g\n",
     
    12111271          if (sequenceOut_==info->crucialSj())
    12121272            info->setCrucialSj(-1);
    1213         }
    1214 
     1273        } else if (info->crucialSj()>=0) {
     1274          // Just possible crucialSj still in but original out again
     1275          int crucialSj=info->crucialSj();
     1276          if (crucialSj>=numberXColumns+numberQuadraticRows) {
     1277            if (sequenceOut_==backColumn[crucialSj-numberXColumns-numberQuadraticRows])
     1278              info->setCrucialSj(-1);
     1279          } else {
     1280            if (sequenceOut_-numberColumns_==whichRow[crucialSj-numberXColumns])
     1281              info->setCrucialSj(-1);
     1282          }
     1283        }
     1284        if (info->crucialSj()<0)
     1285          result=0;
    12151286        if (whatNext==1) {
    12161287          returnCode =-2; // refactorize
     
    12241295        }
    12251296        // may need to go round again
    1226         cleanupIteration = true;
     1297        cleanupIteration = 1;
    12271298        // may not be correct on second time
    12281299        if (result&&(returnCode==-1||returnCode==-2||returnCode==-3)) {
     
    12871358            int best=-1;
    12881359            double bestAlpha=1.0e-6;
    1289             if (numberIterations_==29)
     1360            if (numberIterations_==1180)
    12901361              printf("iteration %d coming up\n",numberIterations_);
    12911362            for (k=0;k<number;k++) {
     
    13011372              } else if (iSequence<numberXColumns) {
    13021373                double djValue = dj_[iSequence];
     1374                // take out if other one basic
     1375                int sj = whichColumn[iSequence];
     1376                if (sj>=0&&getColumnStatus(sj+numberXColumns+numberQuadraticRows)==basic)
     1377                  value=0.0;
    13031378                if (value>0.0)
    13041379                  printf("X%d going up alpha %g, sol %g, dj %g\n",iSequence,value,solution_[crucialSj],djValue);
    1305                 else
     1380                else if (value<0.0)
    13061381                  printf("X%d going down alpha %g, sol %g, dj %g\n",iSequence,value,solution_[crucialSj],djValue);
    13071382                ClpSimplex::Status status = getColumnStatus(iSequence);
     
    13211396                  break;
    13221397                case ClpSimplex::atUpperBound:
    1323                   value = -value;
     1398                  if (djValue>=-1.0e-8)
     1399                    value = -value;
     1400                  else
     1401                    value = -1.0e-4*value;
    13241402                  break;
    13251403                case ClpSimplex::atLowerBound:
     1404                  if (djValue>=1.0e-8)
     1405                    value = 1.0e-4*value;
    13261406                  break;
    13271407                }
     
    13501430              } else if (iSequence<numberXRows) {
    13511431                double djValue = dj_[iSequence+numberColumns_];
     1432                // take out if other one basic
     1433                int sj = whichRow[iSequence];
     1434                if (sj>=0&&getColumnStatus(sj+numberXColumns)==basic)
     1435                  value=0.0;
    13521436                if (value>0.0)
    13531437                  printf("R%d going up alpha %g, sol %g, dj %g\n",iSequence,value,solution_[crucialSj],djValue);
    1354                 else
     1438                else if (value<0.0)
    13551439                  printf("R%d going down alpha %g, sol %g, dj %g\n",iSequence,value,solution_[crucialSj],djValue);
    13561440                ClpSimplex::Status status = getRowStatus(iSequence);
     
    13701454                  break;
    13711455                case ClpSimplex::atUpperBound:
    1372                   value = -value;
     1456                  if (djValue>=-1.0e-8)
     1457                    value = -value;
     1458                  else
     1459                    value = -1.0e-4*value;
    13731460                  break;
    13741461                case ClpSimplex::atLowerBound:
     1462                  if (djValue>=1.0e-8)
     1463                    value = 1.0e-4*value;
    13751464                  break;
    13761465                }
     
    13971486          }
    13981487        } else {
     1488          // may need to bring sj in
     1489          if (0&&(returnCode==-1||returnCode==-2||returnCode==-3)) {
     1490            assert(info->crucialSj()<0);
     1491            chosen=-1;
     1492            if (sequenceOut_<numberXColumns&&whichColumn[sequenceOut_]>=0) {
     1493              chosen =whichColumn[sequenceOut_] + numberQuadraticRows + numberXColumns; // sj variable
     1494              if (getColumnStatus(chosen)==basic)
     1495                chosen=-1; // already in
     1496            } else if (sequenceOut_>=numberColumns_&&whichRow[sequenceOut_-numberColumns_]>=0) {
     1497              chosen = whichRow[sequenceOut_-numberColumns_]+numberXColumns;
     1498              if (getColumnStatus(chosen)==basic)
     1499                chosen=-1; // already in
     1500              // temporary?
     1501              const int * columnLength = matrix_->getVectorLengths();
     1502              if (chosen>=0&&!columnLength[chosen])
     1503                chosen=-1;
     1504            }
     1505            if (chosen>=0) {
     1506              printf("what now\n");
     1507              nextSequenceIn=chosen;
     1508            }
     1509          }
    13991510          break;
    14001511        }
     
    14401551                                     CoinIndexedVector * spareArray2,
    14411552                                     ClpQuadraticInfo * info,
    1442                                      bool cleanupIteration)
     1553                                     int cleanupIteration)
    14431554{
    14441555  int result=1;
     
    16131724  coeff2 *= 0.5;
    16141725  if (coeffSlack)
    1615     printf("trouble\n");
     1726    printf("slack has non-zero cost - trouble?\n");
    16161727  coeff1 += coeffSlack;
    16171728  //coeff1 = way*dualIn_;
     
    16311742      accuracyFlag=2;
    16321743    } else if (fabs(way*coeff1-dualIn_)>1.0e-3*(1.0+fabs(dualIn_))) {
    1633       // not wondeful
     1744      // not wonderful
    16341745      accuracyFlag=1;
    16351746    }
     
    16711782    <<coeff1<<coeff2<<coeffSlack<<dualIn_<<d1<<d2
    16721783    <<CoinMessageEol;
     1784  // reality check
     1785  {
     1786    if (crucialSj2>=0) {
     1787      // see if works out
     1788      if (getColumnStatus(crucialSj2)==basic) {
     1789        if (iSjRow2>=0) {
     1790          //corresponding alpha nonzero
     1791          double alpha=work[iSjRow2];
     1792          printf("Sj alpha %g sol %g ratio %g\n",alpha,solution_[crucialSj2],solution_[crucialSj2]/alpha);
     1793          if (fabs(fabs(d1)-fabs(solution_[crucialSj2]/alpha))>1.0e-3) {
     1794            printf("bad test\n");
     1795            if (factorization_->pivots())
     1796              accuracyFlag=2;
     1797          }
     1798          d1=fabs(solution_[crucialSj2]/alpha);
     1799        } else {
     1800          printf("Sj basic but no alpha\n");
     1801        }
     1802      } else {
     1803        printf("Sj not basic\n");
     1804      }
     1805    }
     1806  }
    16731807  //if (!cleanupIteration)
    16741808  //dualIn_ = way*coeff1;
     
    16961830  double upperTheta = maximumMovement;
    16971831  bool throwAway=false;
    1698   if (numberIterations_==1453) {
     1832  if (numberIterations_==931) {
    16991833    printf("Bad iteration coming up after iteration %d\n",numberIterations_);
    17001834  }
     
    20782212          result=0;
    20792213          iSjRow=iSjRow2;
    2080           crucialSj = pivotVariable_[iSjRow];
     2214          if (iSjRow>=0)
     2215            crucialSj = pivotVariable_[iSjRow];
     2216          else
     2217            crucialSj=-1;
    20812218          assert (pivotRow_<0);
     2219          // If crucialSj <0 then make superbasic?
     2220          if (crucialSj<0) {
     2221            printf("**ouch\n");
     2222            theta_ = maximumMovement;
     2223            pivotRow_ = -2; // so we can tell its a flip
     2224            result=0;
     2225            sequenceOut_ = sequenceIn_;
     2226            valueOut_ = valueIn_;
     2227            dualOut_ = dualIn_;
     2228            lowerOut_ = lowerIn_;
     2229            upperOut_ = upperIn_;
     2230            alpha_ = 0.0;
     2231            if (accuracyFlag<2)
     2232              accuracyFlag=0;
     2233            if (way<0.0) {
     2234              directionOut_=1;      // to lower bound
     2235              lowerIn_ = valueOut_-theta_;
     2236              theta_ = lowerIn_ - valueOut_;
     2237            } else {
     2238              directionOut_=-1;      // to upper bound
     2239              upperIn_ = valueOut_+theta_;
     2240              theta_ = upperIn_ - valueOut_;
     2241            }
     2242          }
    20822243        } else {
    20832244          assert (fabs(dualIn_)<1.0e-3);
     
    21012262        }
    21022263      }
    2103       if (!result) {
     2264      if (!result&&crucialSj>=0) {
    21042265        setColumnStatus(crucialSj,isFree);
    21052266        pivotRow_ = iSjRow;
     
    22122373        numberIterations_ = progress_->lastIterationNumber(0);
    22132374        // no - restore previous basis
     2375        assert (info->crucialSj()<0); // need to work out how to unwind
    22142376        assert (type==1);
    22152377        memcpy(status_ ,saveStatus_,(numberColumns_+numberRows_)*sizeof(char));
     
    22512413      problemStatus_=-3;
    22522414  }
    2253   // at this stage status is -3 or -5 if looks unbounded
    2254   // get primal and dual solutions
    2255   // put back original costs and then check
    2256   createRim(4);
    2257   if (info->currentPhase()) {
    2258     memcpy(solution_,info->currentSolution(),
    2259            (numberRows_+numberColumns_)*sizeof(double));
    2260   }
    2261   gutsOfSolution(NULL,NULL);
    2262   if (info->currentPhase()) {
    2263     memcpy(solution_,info->currentSolution(),
    2264            (numberRows_+numberColumns_)*sizeof(double));
    2265     nonLinearCost_->checkInfeasibilities(primalTolerance_);
    2266   }
    2267   checkComplementarity(info,rowArray_[2],rowArray_[3]);
    2268   // Double check reduced costs if no action
    2269   if (progress->lastIterationNumber(0)==numberIterations_) {
    2270     if (primalColumnPivot_->looksOptimal()) {
     2415  if (info->crucialSj()<0) {
     2416    // at this stage status is -3 or -5 if looks unbounded
     2417    // get primal and dual solutions
     2418    // put back original costs and then check
     2419    createRim(4);
     2420    if (info->currentPhase()) {
     2421      memcpy(solution_,info->currentSolution(),
     2422             (numberRows_+numberColumns_)*sizeof(double));
     2423    }
     2424    gutsOfSolution(NULL,NULL);
     2425    if (info->currentPhase()) {
     2426      memcpy(solution_,info->currentSolution(),
     2427             (numberRows_+numberColumns_)*sizeof(double));
     2428      nonLinearCost_->checkInfeasibilities(primalTolerance_);
     2429    }
     2430    checkComplementarity(info,rowArray_[2],rowArray_[3]);
     2431    // Double check reduced costs if no action
     2432    if (progress->lastIterationNumber(0)==numberIterations_) {
     2433      if (primalColumnPivot_->looksOptimal()) {
     2434        numberDualInfeasibilities_ = 0;
     2435        sumDualInfeasibilities_ = 0.0;
     2436      }
     2437    }
     2438    // Check if looping
     2439    int loop;
     2440    if (type!=2)
     2441      loop = progress->looping(); // saves iteration number
     2442    else
     2443      loop=-1;
     2444    //loop=-1;
     2445    if (loop>=0) {
     2446      problemStatus_ = loop; //exit if in loop
     2447      return ;
     2448    } else if (loop<-1) {
     2449      // something may have changed
     2450      gutsOfSolution(NULL,NULL);
     2451      checkComplementarity(info,rowArray_[2],rowArray_[3]);
     2452    }
     2453    // really for free variables in
     2454    //if((progressFlag_&2)!=0)
     2455    //problemStatus_=-1;;
     2456    progressFlag_ = 0; //reset progress flag
     2457   
     2458    handler_->message(CLP_SIMPLEX_STATUS,messages_)
     2459      <<numberIterations_<<objectiveValue();
     2460    handler_->printing(sumPrimalInfeasibilities_>0.0)
     2461      <<sumPrimalInfeasibilities_<<numberPrimalInfeasibilities_;
     2462    handler_->printing(sumDualInfeasibilities_>0.0)
     2463      <<sumDualInfeasibilities_<<numberDualInfeasibilities_;
     2464    handler_->printing(numberDualInfeasibilitiesWithoutFree_
     2465                       <numberDualInfeasibilities_)
     2466                         <<numberDualInfeasibilitiesWithoutFree_;
     2467    handler_->message()<<CoinMessageEol;
     2468    assert (primalFeasible());
     2469    // we may wish to say it is optimal even if infeasible
     2470    bool alwaysOptimal = (specialOptions_&1)!=0;
     2471    // give code benefit of doubt
     2472    if (sumOfRelaxedDualInfeasibilities_ == 0.0&&
     2473        sumOfRelaxedPrimalInfeasibilities_ == 0.0&&info->sequenceIn()<0) {
     2474      // say optimal (with these bounds etc)
    22712475      numberDualInfeasibilities_ = 0;
    22722476      sumDualInfeasibilities_ = 0.0;
    2273     }
    2274   }
    2275   // Check if looping
    2276   int loop;
    2277   if (type!=2)
    2278     loop = progress->looping(); // saves iteration number
    2279   else
    2280     loop=-1;
    2281   //loop=-1;
    2282   if (loop>=0) {
    2283     problemStatus_ = loop; //exit if in loop
    2284     return ;
    2285   } else if (loop<-1) {
    2286     // something may have changed
    2287     gutsOfSolution(NULL,NULL);
    2288     checkComplementarity(info,rowArray_[2],rowArray_[3]);
    2289   }
    2290   // really for free variables in
    2291   //if((progressFlag_&2)!=0)
    2292   //problemStatus_=-1;;
    2293   progressFlag_ = 0; //reset progress flag
    2294 
    2295   handler_->message(CLP_SIMPLEX_STATUS,messages_)
    2296     <<numberIterations_<<objectiveValue();
    2297   handler_->printing(sumPrimalInfeasibilities_>0.0)
    2298     <<sumPrimalInfeasibilities_<<numberPrimalInfeasibilities_;
    2299   handler_->printing(sumDualInfeasibilities_>0.0)
    2300     <<sumDualInfeasibilities_<<numberDualInfeasibilities_;
    2301   handler_->printing(numberDualInfeasibilitiesWithoutFree_
    2302                      <numberDualInfeasibilities_)
    2303                        <<numberDualInfeasibilitiesWithoutFree_;
    2304   handler_->message()<<CoinMessageEol;
    2305   assert (primalFeasible());
    2306   // we may wish to say it is optimal even if infeasible
    2307   bool alwaysOptimal = (specialOptions_&1)!=0;
    2308   // give code benefit of doubt
    2309   if (sumOfRelaxedDualInfeasibilities_ == 0.0&&
    2310       sumOfRelaxedPrimalInfeasibilities_ == 0.0&&info->sequenceIn()<0) {
    2311     // say optimal (with these bounds etc)
    2312     numberDualInfeasibilities_ = 0;
    2313     sumDualInfeasibilities_ = 0.0;
    2314     numberPrimalInfeasibilities_ = 0;
    2315     sumPrimalInfeasibilities_ = 0.0;
    2316   }
    2317   // had ||(type==3&&problemStatus_!=-5) -- ??? why ????
    2318   if (dualFeasible()||problemStatus_==-4) {
    2319     if (nonLinearCost_->numberInfeasibilities()&&!alwaysOptimal) {
    2320       //may need infeasiblity cost changed
    2321       // we can see if we can construct a ray
    2322       // make up a new objective
    2323       double saveWeight = infeasibilityCost_;
    2324       // save nonlinear cost as we are going to switch off costs
    2325       ClpNonLinearCost * nonLinear = nonLinearCost_;
    2326       // do twice to make sure Primal solution has settled
    2327       // put non-basics to bounds in case tolerance moved
    2328       // put back original costs
    2329       createRim(4);
    2330       // put all non-basic variables to bounds
    2331       int iSequence;
    2332       for (iSequence=0;iSequence<numberRows_+numberColumns_;iSequence++) {
    2333         Status status = getStatus(iSequence);
    2334         if (status==atLowerBound||status==isFixed)
    2335           solution_[iSequence]=lower_[iSequence];
    2336         else if (status==atUpperBound)
    2337           solution_[iSequence]=upper_[iSequence];
    2338       }
    2339       nonLinearCost_->checkInfeasibilities(primalTolerance_);
    2340       gutsOfSolution(NULL,NULL);
    2341       nonLinearCost_->checkInfeasibilities(primalTolerance_);
    2342       checkComplementarity(info,rowArray_[2],rowArray_[3]);
    2343 
    2344       infeasibilityCost_=1.0e100;
    2345       // put back original costs
    2346       createRim(4);
    2347       nonLinearCost_->checkInfeasibilities(primalTolerance_);
    2348       // may have fixed infeasibilities - double check
    2349       if (nonLinearCost_->numberInfeasibilities()==0) {
    2350         // carry on
    2351         problemStatus_ = -1;
    2352         infeasibilityCost_=saveWeight;
     2477      numberPrimalInfeasibilities_ = 0;
     2478      sumPrimalInfeasibilities_ = 0.0;
     2479    }
     2480    // had ||(type==3&&problemStatus_!=-5) -- ??? why ????
     2481    if (dualFeasible()||problemStatus_==-4) {
     2482      if (nonLinearCost_->numberInfeasibilities()&&!alwaysOptimal) {
     2483        //may need infeasiblity cost changed
     2484        // we can see if we can construct a ray
     2485        // make up a new objective
     2486        double saveWeight = infeasibilityCost_;
     2487        // save nonlinear cost as we are going to switch off costs
     2488        ClpNonLinearCost * nonLinear = nonLinearCost_;
     2489        // do twice to make sure Primal solution has settled
     2490        // put non-basics to bounds in case tolerance moved
     2491        // put back original costs
     2492        createRim(4);
     2493        // put all non-basic variables to bounds
     2494        int iSequence;
     2495        for (iSequence=0;iSequence<numberRows_+numberColumns_;iSequence++) {
     2496          Status status = getStatus(iSequence);
     2497          if (status==atLowerBound||status==isFixed)
     2498            solution_[iSequence]=lower_[iSequence];
     2499          else if (status==atUpperBound)
     2500            solution_[iSequence]=upper_[iSequence];
     2501        }
     2502        nonLinearCost_->checkInfeasibilities(primalTolerance_);
     2503        gutsOfSolution(NULL,NULL);
    23532504        nonLinearCost_->checkInfeasibilities(primalTolerance_);
    23542505        checkComplementarity(info,rowArray_[2],rowArray_[3]);
     2506       
     2507        infeasibilityCost_=1.0e100;
     2508        // put back original costs
     2509        createRim(4);
     2510        nonLinearCost_->checkInfeasibilities(primalTolerance_);
     2511        // may have fixed infeasibilities - double check
     2512        if (nonLinearCost_->numberInfeasibilities()==0) {
     2513          // carry on
     2514          problemStatus_ = -1;
     2515          infeasibilityCost_=saveWeight;
     2516          nonLinearCost_->checkInfeasibilities(primalTolerance_);
     2517          checkComplementarity(info,rowArray_[2],rowArray_[3]);
     2518        } else {
     2519          nonLinearCost_=NULL;
     2520          // scale
     2521          int i;
     2522          for (i=0;i<numberRows_+numberColumns_;i++)
     2523            cost_[i] *= 1.0e-100;
     2524          gutsOfSolution(NULL,NULL);
     2525          nonLinearCost_=nonLinear;
     2526          infeasibilityCost_=saveWeight;
     2527          if ((infeasibilityCost_>=1.0e18||
     2528               numberDualInfeasibilities_==0)&&perturbation_==101) {
     2529            unPerturb(); // stop any further perturbation
     2530            numberDualInfeasibilities_=1; // carry on
     2531          }
     2532          if (infeasibilityCost_>=1.0e20||
     2533              numberDualInfeasibilities_==0) {
     2534            // we are infeasible - use as ray
     2535            delete [] ray_;
     2536            ray_ = new double [numberRows_];
     2537            memcpy(ray_,dual_,numberRows_*sizeof(double));
     2538            // and get feasible duals
     2539            infeasibilityCost_=0.0;
     2540            createRim(4);
     2541            // put all non-basic variables to bounds
     2542            int iSequence;
     2543            for (iSequence=0;iSequence<numberRows_+numberColumns_;iSequence++) {
     2544              Status status = getStatus(iSequence);
     2545              if (status==atLowerBound||status==isFixed)
     2546                solution_[iSequence]=lower_[iSequence];
     2547              else if (status==atUpperBound)
     2548                solution_[iSequence]=upper_[iSequence];
     2549            }
     2550            nonLinearCost_->checkInfeasibilities(primalTolerance_);
     2551            gutsOfSolution(NULL,NULL);
     2552            checkComplementarity(info,rowArray_[2],rowArray_[3]);
     2553            // so will exit
     2554            infeasibilityCost_=1.0e30;
     2555            // reset infeasibilities
     2556            sumPrimalInfeasibilities_=nonLinearCost_->sumInfeasibilities();;
     2557            numberPrimalInfeasibilities_=
     2558              nonLinearCost_->numberInfeasibilities();
     2559          }
     2560          if (infeasibilityCost_<1.0e20) {
     2561            infeasibilityCost_ *= 5.0;
     2562            changeMade_++; // say change made
     2563            handler_->message(CLP_PRIMAL_WEIGHT,messages_)
     2564              <<infeasibilityCost_
     2565              <<CoinMessageEol;
     2566            // put back original costs and then check
     2567            createRim(4);
     2568            nonLinearCost_->checkInfeasibilities(0.0);
     2569            gutsOfSolution(NULL,NULL);
     2570            checkComplementarity(info,rowArray_[2],rowArray_[3]);
     2571            problemStatus_=-1; //continue
     2572          } else {
     2573            // say infeasible
     2574            problemStatus_ = 1;
     2575          }
     2576        }
    23552577      } else {
    2356         nonLinearCost_=NULL;
    2357         // scale
    2358         int i;
    2359         for (i=0;i<numberRows_+numberColumns_;i++)
    2360           cost_[i] *= 1.0e-100;
    2361         gutsOfSolution(NULL,NULL);
    2362         nonLinearCost_=nonLinear;
    2363         infeasibilityCost_=saveWeight;
    2364         if ((infeasibilityCost_>=1.0e18||
    2365              numberDualInfeasibilities_==0)&&perturbation_==101) {
     2578        // may be optimal
     2579        if (perturbation_==101) {
    23662580          unPerturb(); // stop any further perturbation
    2367           numberDualInfeasibilities_=1; // carry on
    2368         }
    2369         if (infeasibilityCost_>=1.0e20||
    2370             numberDualInfeasibilities_==0) {
    2371           // we are infeasible - use as ray
    2372           delete [] ray_;
    2373           ray_ = new double [numberRows_];
    2374           memcpy(ray_,dual_,numberRows_*sizeof(double));
    2375           // and get feasible duals
    2376           infeasibilityCost_=0.0;
    2377           createRim(4);
    2378           // put all non-basic variables to bounds
    2379           int iSequence;
    2380           for (iSequence=0;iSequence<numberRows_+numberColumns_;iSequence++) {
    2381             Status status = getStatus(iSequence);
    2382             if (status==atLowerBound||status==isFixed)
    2383               solution_[iSequence]=lower_[iSequence];
    2384             else if (status==atUpperBound)
    2385               solution_[iSequence]=upper_[iSequence];
     2581          lastCleaned=-1; // carry on
     2582        }
     2583        if ( lastCleaned!=numberIterations_||unflag()) {
     2584          handler_->message(CLP_PRIMAL_OPTIMAL,messages_)
     2585            <<primalTolerance_
     2586            <<CoinMessageEol;
     2587          if (numberTimesOptimal_<4) {
     2588            numberTimesOptimal_++;
     2589            changeMade_++; // say change made
     2590            if (numberTimesOptimal_==1) {
     2591              // better to have small tolerance even if slower
     2592              factorization_->zeroTolerance(1.0e-15);
     2593            }
     2594            lastCleaned=numberIterations_;
     2595            if (primalTolerance_!=dblParam_[ClpPrimalTolerance])
     2596              handler_->message(CLP_PRIMAL_ORIGINAL,messages_)
     2597                <<CoinMessageEol;
     2598            double oldTolerance = primalTolerance_;
     2599            primalTolerance_=dblParam_[ClpPrimalTolerance];
     2600            // put back original costs and then check
     2601            createRim(4);
     2602            // put all non-basic variables to bounds
     2603            int iSequence;
     2604            for (iSequence=0;iSequence<numberRows_+numberColumns_;iSequence++) {
     2605              Status status = getStatus(iSequence);
     2606              if (status==atLowerBound||status==isFixed)
     2607                solution_[iSequence]=lower_[iSequence];
     2608              else if (status==atUpperBound)
     2609                solution_[iSequence]=upper_[iSequence];
     2610            }
     2611            nonLinearCost_->checkInfeasibilities(oldTolerance);
     2612            gutsOfSolution(NULL,NULL);
     2613            checkComplementarity(info,rowArray_[2],rowArray_[3]);
     2614            problemStatus_ = -1;
     2615          } else {
     2616            problemStatus_=0; // optimal
     2617            if (lastCleaned<numberIterations_) {
     2618              handler_->message(CLP_SIMPLEX_GIVINGUP,messages_)
     2619                <<CoinMessageEol;
     2620            }
    23862621          }
    2387           nonLinearCost_->checkInfeasibilities(primalTolerance_);
    2388           gutsOfSolution(NULL,NULL);
    2389           checkComplementarity(info,rowArray_[2],rowArray_[3]);
    2390           // so will exit
    2391           infeasibilityCost_=1.0e30;
    2392           // reset infeasibilities
    2393           sumPrimalInfeasibilities_=nonLinearCost_->sumInfeasibilities();;
    2394           numberPrimalInfeasibilities_=
    2395             nonLinearCost_->numberInfeasibilities();
    2396         }
    2397         if (infeasibilityCost_<1.0e20) {
    2398           infeasibilityCost_ *= 5.0;
    2399           changeMade_++; // say change made
    2400           handler_->message(CLP_PRIMAL_WEIGHT,messages_)
    2401             <<infeasibilityCost_
    2402             <<CoinMessageEol;
    2403           // put back original costs and then check
    2404           createRim(4);
    2405           nonLinearCost_->checkInfeasibilities(0.0);
    2406           gutsOfSolution(NULL,NULL);
    2407           checkComplementarity(info,rowArray_[2],rowArray_[3]);
    2408           problemStatus_=-1; //continue
    2409         } else {
    2410           // say infeasible
    2411           problemStatus_ = 1;
    2412         }
    2413       }
    2414     } else {
    2415       // may be optimal
    2416       if (perturbation_==101) {
    2417         unPerturb(); // stop any further perturbation
    2418         lastCleaned=-1; // carry on
    2419       }
    2420       if ( lastCleaned!=numberIterations_||unflag()) {
    2421         handler_->message(CLP_PRIMAL_OPTIMAL,messages_)
    2422           <<primalTolerance_
    2423           <<CoinMessageEol;
    2424         if (numberTimesOptimal_<4) {
    2425           numberTimesOptimal_++;
    2426           changeMade_++; // say change made
    2427           if (numberTimesOptimal_==1) {
    2428             // better to have small tolerance even if slower
    2429             factorization_->zeroTolerance(1.0e-15);
    2430           }
    2431           lastCleaned=numberIterations_;
    2432           if (primalTolerance_!=dblParam_[ClpPrimalTolerance])
    2433             handler_->message(CLP_PRIMAL_ORIGINAL,messages_)
    2434               <<CoinMessageEol;
    2435           double oldTolerance = primalTolerance_;
    2436           primalTolerance_=dblParam_[ClpPrimalTolerance];
    2437           // put back original costs and then check
    2438           createRim(4);
    2439           // put all non-basic variables to bounds
    2440           int iSequence;
    2441           for (iSequence=0;iSequence<numberRows_+numberColumns_;iSequence++) {
    2442             Status status = getStatus(iSequence);
    2443             if (status==atLowerBound||status==isFixed)
    2444               solution_[iSequence]=lower_[iSequence];
    2445             else if (status==atUpperBound)
    2446               solution_[iSequence]=upper_[iSequence];
    2447           }
    2448           nonLinearCost_->checkInfeasibilities(oldTolerance);
    2449           gutsOfSolution(NULL,NULL);
    2450           checkComplementarity(info,rowArray_[2],rowArray_[3]);
    2451           problemStatus_ = -1;
    24522622        } else {
    24532623          problemStatus_=0; // optimal
    2454           if (lastCleaned<numberIterations_) {
    2455             handler_->message(CLP_SIMPLEX_GIVINGUP,messages_)
     2624        }
     2625      }
     2626    } else {
     2627      // see if looks unbounded
     2628      if (problemStatus_==-5) {
     2629        if (nonLinearCost_->numberInfeasibilities()) {
     2630          if (infeasibilityCost_>1.0e18&&perturbation_==101) {
     2631            // back off weight
     2632            infeasibilityCost_ = 1.0e13;
     2633            unPerturb(); // stop any further perturbation
     2634          }
     2635          //we need infeasiblity cost changed
     2636          if (infeasibilityCost_<1.0e20) {
     2637            infeasibilityCost_ *= 5.0;
     2638            changeMade_++; // say change made
     2639            handler_->message(CLP_PRIMAL_WEIGHT,messages_)
     2640              <<infeasibilityCost_
    24562641              <<CoinMessageEol;
     2642            // put back original costs and then check
     2643            createRim(4);
     2644            gutsOfSolution(NULL,NULL);
     2645            checkComplementarity(info,rowArray_[2],rowArray_[3]);
     2646            problemStatus_=-1; //continue
     2647          } else {
     2648            // say unbounded
     2649            problemStatus_ = 2;
    24572650          }
    2458         }
    2459       } else {
    2460         problemStatus_=0; // optimal
    2461       }
    2462     }
    2463   } else {
    2464     // see if looks unbounded
    2465     if (problemStatus_==-5) {
    2466       if (nonLinearCost_->numberInfeasibilities()) {
    2467         if (infeasibilityCost_>1.0e18&&perturbation_==101) {
    2468           // back off weight
    2469           infeasibilityCost_ = 1.0e13;
    2470           unPerturb(); // stop any further perturbation
    2471         }
    2472         //we need infeasiblity cost changed
    2473         if (infeasibilityCost_<1.0e20) {
    2474           infeasibilityCost_ *= 5.0;
    2475           changeMade_++; // say change made
    2476           handler_->message(CLP_PRIMAL_WEIGHT,messages_)
    2477             <<infeasibilityCost_
    2478             <<CoinMessageEol;
    2479           // put back original costs and then check
    2480           createRim(4);
    2481           gutsOfSolution(NULL,NULL);
    2482           checkComplementarity(info,rowArray_[2],rowArray_[3]);
    2483           problemStatus_=-1; //continue
    24842651        } else {
    24852652          // say unbounded
    24862653          problemStatus_ = 2;
    2487         }
     2654        } 
    24882655      } else {
    2489         // say unbounded
    2490         problemStatus_ = 2;
    2491       }
    2492     } else {
    2493       if(type==3&&problemStatus_!=-5)
    2494         unflag(); // odd
    2495       // carry on
    2496       problemStatus_ = -1;
    2497     }
     2656        if(type==3&&problemStatus_!=-5)
     2657          unflag(); // odd
     2658        // carry on
     2659        problemStatus_ = -1;
     2660      }
     2661    }
     2662  } else {
     2663    // don't update solution
     2664    problemStatus_=-1;
    24982665  }
    24992666  if (type==0||type==1) {
     
    25362703  const int * backRow = info->backRowSequence();
    25372704  int numberQuadraticRows = info->numberQuadraticRows();
    2538   //const int * back = info->backSequence();
    2539   //int numberQuadraticColumns = info->numberQuadraticColumns();
     2705  const int * back = info->backSequence();
     2706  int numberQuadraticColumns = info->numberQuadraticColumns();
    25402707  int numberXColumns = info->numberXColumns();
    25412708  int numberXRows = info->numberXRows();
     
    25582725    int * index = array1->getIndices();
    25592726    double * modifiedCost = array1->denseVector();
     2727    //#define CHECK
     2728    //#ifdef CHECK
     2729#if 0
    25602730    // zero out basic values
    25612731    int iRow;
     
    26912861      if (jSequence>=0) {
    26922862        jSequence += start;
     2863        //solution_[jSequence]=gradient[iSequence];
    26932864        value = solution_[jSequence];
    26942865      } else {
     
    27182889    // Value of sj is - value of pi
    27192890    for (iSequence=numberColumns_;iSequence<numberColumns_+numberXRows;iSequence++) {
    2720       int jSequence = whichRow[iSequence-numberColumns_];
     2891      int iRow = iSequence-numberColumns_;
     2892      int jSequence = whichRow[iRow];
    27212893      double value;
    27222894      if (jSequence>=0) {
    27232895        int iPi = jSequence+numberXColumns;
    27242896        value = solution_[iPi]+cost_[iSequence];
     2897        const int * columnLength = matrix_->getVectorLengths();
     2898        if (!columnLength[iPi])
     2899          value = dual_[iRow]+cost_[iSequence];
    27252900      } else {
    27262901        value=dj_[iSequence];
    2727         value=gradient[iSequence];
     2902        value=dual_[iRow]+cost_[iSequence];
    27282903      }
    27292904      double value2 = value*djWeight[iSequence];
     
    27382913      dj_[iSequence]=value;
    27392914    }
    2740   }
    2741   if (info->crucialSj()<0) {
     2915#if 0
     2916    if (info->crucialSj()<0) {
     2917      for (iSequence=numberXColumns+numberQuadraticRows;iSequence<numberColumns_;iSequence++)
     2918        if (fabs(solution_[iSequence])>1.0e-4)
     2919          printf("sj %d (%d) value %g\n",iSequence-numberXColumns-numberQuadraticRows,
     2920                 back[iSequence-numberXColumns-numberQuadraticRows],solution_[iSequence]);
     2921    }
     2922    if (info->crucialSj()<0) {
     2923      int i;
     2924      for (i=0;i<numberXRows;i++) {
     2925        double pi1 = dual_[i];
     2926        double pi2 = solution_[numberXColumns+i];
     2927        if (fabs(pi1-pi2)>1.0e-3)
     2928          printf("row %d, dual %g sol %g\n",i,pi1,pi2);
     2929      }
     2930    }
     2931#endif
     2932#endif
     2933#if 1
     2934    {
     2935      //if (info->crucialSj()>=0)
     2936      //return;
     2937#ifdef CHECK
     2938      double saveSol[2000];
     2939      memcpy(saveSol,solution_,(numberRows_+numberColumns_)*sizeof(double));
     2940      double saveDj[2000];
     2941      memcpy(saveDj,dj_,(numberRows_+numberColumns_)*sizeof(double));
     2942#endif
     2943      // Compute duals
     2944      info->createGradient(this);
     2945      double * gradient = info->gradient();
     2946      // fill in as if linear
     2947      // will not be correct unless complementary - but that should not matter
     2948      // Just save sj value in that case
     2949      //double saveValue=0.0;
     2950      //int crucialSj = info->crucialSj();
     2951      //if (crucialSj>=0)
     2952      int number=0;
     2953      int iRow;
     2954      for (iRow=0;iRow<numberRows_;iRow++) {
     2955        int iPivot=pivotVariable_[iRow];
     2956        if (gradient[iPivot]) {
     2957          modifiedCost[iRow] = gradient[iPivot];
     2958          index[number++]=iRow;
     2959        }
     2960      }
     2961      array1->setNumElements(number);
     2962      factorization_->updateColumnTranspose(array2,array1);
     2963      memcpy(dual_,modifiedCost,numberRows_*sizeof(double));
     2964      array1->clear();
     2965      memcpy(dj_,gradient,(numberColumns_+numberRows_)*sizeof(double));
     2966      matrix_->transposeTimes(-1.0,dual_,dj_,rowScale_,columnScale_);
     2967      memset(dj_+numberXColumns,0,(numberXRows+numberQuadraticColumns)*sizeof(double));
     2968      // And djs
     2969      for (iRow=0;iRow<numberRows_;iRow++) {
     2970        dj_[numberColumns_+iRow]= cost_[numberColumns_+iRow]+dual_[iRow];
     2971      }
     2972      double saveSj=0.0;
     2973      if (info->crucialSj()>=0)
     2974        saveSj=solution_[info->crucialSj()];
     2975      if (info->crucialSj()>=0&&0) {
     2976#ifdef CHECK
     2977        for (iSequence=0;iSequence<numberColumns_+numberRows_;iSequence++)
     2978          if (fabs(solution_[iSequence]-saveSol[iSequence])>1.0e-4)
     2979            printf("Badxsol %d value %g true %g\n",iSequence,
     2980                   solution_[iSequence],saveSol[iSequence]);
     2981        for (iSequence=0;iSequence<numberColumns_+numberRows_;iSequence++)
     2982          if (fabs(dj_[iSequence]-saveDj[iSequence])>1.0e-4)
     2983            printf("Badxdj %d value %g true %g\n",iSequence,
     2984                   dj_[iSequence],saveDj[iSequence]);
     2985        memcpy(solution_,saveSol,(numberRows_+numberColumns_)*sizeof(double));
     2986        memcpy(dj_,saveDj,(numberRows_+numberColumns_)*sizeof(double));
     2987#endif
     2988        if (numberIterations_==-1289) {
     2989          int i;
     2990          printf ("returning on sj\n");
     2991          for (i=0;i<numberRows_+numberColumns_;i++) {
     2992            char x='N';
     2993            if (getColumnStatus(i)==basic)
     2994              x='B';
     2995            printf("CH %d %c sol %g dj %g\n",i,x,solution_[i],dj_[i]);
     2996          }
     2997        }
     2998        return;
     2999      }
     3000      // Set dual solution
     3001      for (iRow=0;iRow<numberQuadraticRows;iRow++) {
     3002        int jRow = backRow[iRow];
     3003        solution_[iRow+numberXColumns]=dual_[jRow];
     3004      }
     3005      // clear sj
     3006      int start = numberXColumns+numberQuadraticRows;
     3007      memset(solution_+start,0,numberQuadraticColumns*sizeof(double));
     3008      memset(dj_+numberXColumns,0,(numberQuadraticRows+numberQuadraticColumns)*sizeof(double));
     3009      matrix_->times(-1.0,columnActivityWork_,modifiedCost,rowScale_,columnScale_);
     3010      memset(modifiedCost,0,numberXRows*sizeof(double));
     3011      // Do costs as rhs and sj solution values
     3012      for (iRow=0;iRow<numberQuadraticColumns;iRow++) {
     3013        int jSequence = back[iRow];
     3014        double value=cost_[jSequence];
     3015        if (fabs(value)>1.0e5) {
     3016          assert (getColumnStatus(jSequence)==basic);
     3017        }
     3018        double value2=-modifiedCost[iRow+numberXRows];
     3019        modifiedCost[iRow+numberXRows]=0.0;
     3020        jSequence = iRow+start;
     3021        if (getColumnStatus(jSequence)!=basic) {
     3022          if (fabs(value2-value)>1.0e-3)
     3023            //printf("bad nonbasic match %d %g %g\n",iRow,value,value2);
     3024          // should adjust value?
     3025          solution_[jSequence]=0.0;
     3026          //value=value2;
     3027        } else {
     3028          //if (fabs(value-value2-dj_[jSequence-start])>1.0e-3)
     3029          //printf("bad basic match %d %g %g\n",iRow,value,value2);
     3030          solution_[jSequence]=value-value2;
     3031        }
     3032        storedCost[iRow]=value;
     3033        storedUpper[iRow]=value;
     3034        storedValue[iRow]=value;
     3035      }
     3036      if (info->crucialSj()>=0)
     3037        solution_[info->crucialSj()]=saveSj;
     3038#ifdef CHECK
     3039      for (iSequence=numberXColumns+numberQuadraticRows;iSequence<numberColumns_;iSequence++)
     3040        if (fabs(solution_[iSequence])>1.0e-4||fabs(saveSol[iSequence])>1.0e-4)
     3041          printf("sj %d (%d) value %g true %g\n",iSequence-numberXColumns-numberQuadraticRows,
     3042                 back[iSequence-numberXColumns-numberQuadraticRows],solution_[iSequence],saveSol[iSequence]);
     3043      for (iSequence=0;iSequence<numberColumns_+numberRows_;iSequence++)
     3044        if (fabs(solution_[iSequence]-saveSol[iSequence])>1.0e-4)
     3045          printf("Badsol %d value %g true %g\n",iSequence,
     3046                 solution_[iSequence],saveSol[iSequence]);
     3047      for (iSequence=0;iSequence<numberColumns_+numberRows_;iSequence++)
     3048        if (fabs(dj_[iSequence]-saveDj[iSequence])>1.0e-4)
     3049          printf("Baddj %d value %g true %g\n",iSequence,
     3050                 dj_[iSequence],saveDj[iSequence]);
     3051      memcpy(solution_,saveSol,(numberRows_+numberColumns_)*sizeof(double));
     3052      memcpy(dj_,saveDj,(numberRows_+numberColumns_)*sizeof(double));
     3053#endif
     3054    }
     3055#endif
     3056  }
     3057  if (numberIterations_==-1289) {
    27423058    int i;
    2743     for (i=0;i<numberXRows;i++) {
    2744       double pi1 = dual_[i];
    2745       double pi2 = solution_[numberXColumns+i];
    2746       if (fabs(pi1-pi2)>1.0e-3)
    2747         printf("row %d, dual %g sol %g\n",i,pi1,pi2);
     3059    printf ("normal\n");
     3060    for (i=0;i<numberRows_+numberColumns_;i++) {
     3061      char x='N';
     3062      if (getColumnStatus(i)==basic)
     3063        x='B';
     3064      printf("CH %d %c sol %g dj %g\n",i,x,solution_[i],dj_[i]);
    27483065    }
    27493066  }
     
    27743091  double linearCost=0.0;
    27753092  int number0=0,number1=0,number2=0;
     3093#if 0
     3094  int numberNSj=0;
     3095  for (iColumn=numberXColumns+info->numberQuadraticRows();iColumn<numberColumns_;iColumn++) {
     3096    if (getColumnStatus(iColumn)!=basic)
     3097      numberNSj++;
     3098  }
     3099  printf("NumberN %d\n",numberNSj);
     3100#endif
    27763101  for (iColumn=0;iColumn<numberXColumns;iColumn++) {
    27773102    double valueI = solution_[iColumn];
     
    28183143            <<CoinMessageEol;
    28193144          number2++;
     3145          int crucialSj =info->crucialSj();
     3146          assert (crucialSj>=0);
     3147          assert (crucialSj==iColumn||crucialSj==jSequence);
    28203148        } else {
    28213149          number1++;
     
    28773205               iRow,solution_[iColumn],jSequence,solution_[jSequence]);
    28783206          number2++;
     3207          int crucialSj =info->crucialSj();
     3208          assert (crucialSj>=0);
     3209          assert (crucialSj==iColumn||crucialSj==jSequence);
    28793210        } else {
    28803211          number1++;
     
    29033234    }
    29043235  }
    2905   //printf("number 0 - %d, 1 - %d, 2 - %d\n",number0,number1,number2);
     3236  printf("number 0 - %d, 1 - %d, 2 - %d\n",number0,number1,number2);
    29063237  objectiveValue_ += linearCost+infeasCost;
    29073238  assert (infeasCost>=0.0);
     
    31693500  delete [] elements2;
    31703501  model2->allSlackBasis();
    3171   model2->scaling(false);
    3172   printf("scaling off\n");
     3502  //model2->scaling(false);
     3503  //printf("scaling off\n");
    31733504  model2->setLogLevel(this->logLevel());
    31743505  // Move solution across
  • branches/pre/Test/unitTest.cpp

    r200 r201  
    12101210    model.setLogLevel(63);
    12111211    //exit(77);
     1212    model.setFactorizationFrequency(1);
    12121213    model.quadraticPrimal(1);
    12131214    double objValue = model.getObjValue();
  • branches/pre/include/ClpSimplexPrimalQuadratic.hpp

    r200 r201  
    8080      Returns 0 - can do normal iteration
    8181      1 - losing complementarity
     82      cleanupIteration - 0 no, 1 yes, 2 restoring one of x/s in basis
    8283  */
    8384  int primalRow(CoinIndexedVector * rowArray,
     
    8687                CoinIndexedVector * spareArray2,
    8788                ClpQuadraticInfo * info,
    88                 bool cleanupIteration);
     89                int cleanupIteration);
    8990  /**  Refactorizes if necessary
    9091       Checks if finished.  Updates status.
Note: See TracChangeset for help on using the changeset viewer.