source: trunk/ClpSimplexPrimalQuadratic.cpp @ 284

Last change on this file since 284 was 277, checked in by forrest, 16 years ago

Allowing CoinBigIndex? as long

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 110.2 KB
Line 
1// Copyright (C) 2003, International Business Machines
2// Corporation and others.  All Rights Reserved.
3
4#include "CoinPragma.hpp"
5
6#include <math.h>
7
8#include "CoinHelperFunctions.hpp"
9#include "ClpSimplexPrimalQuadratic.hpp"
10#include "ClpPrimalQuadraticDantzig.hpp"
11#include "ClpPrimalColumnDantzig.hpp"
12#include "ClpQuadraticObjective.hpp"
13#include "ClpFactorization.hpp"
14#include "ClpNonLinearCost.hpp"
15#include "ClpPackedMatrix.hpp"
16#include "CoinIndexedVector.hpp"
17#include "CoinWarmStartBasis.hpp"
18#include "CoinMpsIO.hpp"
19#include "ClpPrimalColumnPivot.hpp"
20#include "ClpMessage.hpp"
21#include <cfloat>
22#include <cassert>
23#include <string>
24#include <stdio.h>
25#include <iostream>
26class tempMessage :
27   public CoinMessageHandler {
28
29public:
30  virtual int print() ;
31  tempMessage(ClpSimplex * model);
32  ClpSimplex * model_;
33};
34
35// Constructor with pointer to model
36tempMessage::tempMessage(ClpSimplex * model)
37  : CoinMessageHandler(),
38    model_(model)
39{
40}
41int
42tempMessage::print()
43{
44  static int numberFeasible=0;
45  if (currentSource()=="Clp") {
46    if (currentMessage().externalNumber()==5) { 
47      if (!numberFeasible&&!model_->nonLinearCost()->numberInfeasibilities()) {
48        numberFeasible++;
49        model_->setMaximumIterations(0);
50      }
51    } 
52  }
53  return CoinMessageHandler::print();
54}
55
56// A sequential LP method
57int 
58ClpSimplexPrimalQuadratic::primalSLP(int numberPasses, double deltaTolerance)
59{
60  // Are we minimizing or maximizing
61  double whichWay=optimizationDirection();
62  if (whichWay<0.0)
63    whichWay=-1.0;
64  else if (whichWay>0.0)
65    whichWay=1.0;
66
67  // This is as a user would see
68
69  int numberColumns = this->numberColumns();
70  int numberRows = this->numberRows();
71  double * columnLower = this->columnLower();
72  double * columnUpper = this->columnUpper();
73  double * objective = this->objective();
74  double * solution = this->primalColumnSolution();
75 
76  // Save objective
77 
78  double * saveObjective = new double [numberColumns];
79  memcpy(saveObjective,objective,numberColumns*sizeof(double));
80
81  // Get list of non linear columns
82  ClpQuadraticObjective * quadraticObj = (dynamic_cast< ClpQuadraticObjective*>(objective_));
83  CoinPackedMatrix * quadratic = NULL;
84  if (quadraticObj)
85    quadratic = quadraticObj->quadraticObjective();
86  if (!quadratic) {
87    // no quadratic part
88    return primal(0);
89  }
90  int numberNonLinearColumns = 0;
91  int iColumn;
92  int * listNonLinearColumn = new int[numberColumns];
93  memset(listNonLinearColumn,0,numberColumns*sizeof(int));
94  const int * columnQuadratic = quadratic->getIndices();
95  const CoinBigIndex * columnQuadraticStart = quadratic->getVectorStarts();
96  const int * columnQuadraticLength = quadratic->getVectorLengths();
97  const double * quadraticElement = quadratic->getElements();
98  for (iColumn=0;iColumn<numberColumns;iColumn++) {
99    int j;
100    for (j=columnQuadraticStart[iColumn];
101         j<columnQuadraticStart[iColumn]+columnQuadraticLength[iColumn];j++) {
102      int jColumn = columnQuadratic[j];
103      listNonLinearColumn[jColumn]=1;
104      listNonLinearColumn[iColumn]=1;
105    }
106  }
107  for (iColumn=0;iColumn<numberColumns;iColumn++) {
108    if(listNonLinearColumn[iColumn])
109      listNonLinearColumn[numberNonLinearColumns++]=iColumn;
110  }
111 
112  if (!numberNonLinearColumns) {
113    delete [] listNonLinearColumn;
114    // no quadratic part
115    return primal(0);
116  }
117
118  // get feasible
119  if (!this->status()||numberPrimalInfeasibilities())
120    primal(1);
121  // still infeasible
122  if (numberPrimalInfeasibilities())
123    return 0;
124
125  int jNon;
126  int * last[3];
127 
128  double * trust = new double[numberNonLinearColumns];
129  double * trueLower = new double[numberNonLinearColumns];
130  double * trueUpper = new double[numberNonLinearColumns];
131  for (jNon=0;jNon<numberNonLinearColumns;jNon++) {
132    iColumn=listNonLinearColumn[jNon];
133    trust[jNon]=0.5;
134    trueLower[jNon]=columnLower[iColumn];
135    trueUpper[jNon]=columnUpper[iColumn];
136    if (solution[iColumn]<trueLower[jNon])
137      solution[iColumn]=trueLower[jNon];
138    else if (solution[iColumn]>trueUpper[jNon])
139      solution[iColumn]=trueUpper[jNon];
140  }
141  int iPass;
142  double lastObjective=1.0e31;
143  double * saveSolution = new double [numberColumns];
144  double * savePi = new double [numberRows];
145  unsigned char * saveStatus = new unsigned char[numberRows+numberColumns];
146  double targetDrop=1.0e31;
147  double objectiveOffset;
148  getDblParam(ClpObjOffset,objectiveOffset);
149  // 1 bound up, 2 up, -1 bound down, -2 down, 0 no change
150  for (iPass=0;iPass<3;iPass++) {
151    last[iPass]=new int[numberNonLinearColumns];
152    for (jNon=0;jNon<numberNonLinearColumns;jNon++) 
153      last[iPass][jNon]=0;
154  }
155  // goodMove +1 yes, 0 no, -1 last was bad - just halve gaps, -2 do nothing
156  int goodMove=-2;
157  char * statusCheck = new char[numberColumns];
158  for (iPass=0;iPass<numberPasses;iPass++) {
159    // redo objective
160    double offset=0.0;
161    double objValue=-objectiveOffset;
162    double lambda=-1.0;
163    if (goodMove>=0) {
164      // get best interpolation
165      double coeff0=-objectiveOffset,coeff1=0.0,coeff2=0.0;
166      for (iColumn=0;iColumn<numberColumns;iColumn++) {
167        coeff0 += saveObjective[iColumn]*solution[iColumn];
168        coeff1 += saveObjective[iColumn]*(saveSolution[iColumn]-solution[iColumn]);
169      }
170      for (jNon=0;jNon<numberNonLinearColumns;jNon++) {
171        iColumn=listNonLinearColumn[jNon];
172        double valueI = solution[iColumn];
173        double valueISave = saveSolution[iColumn];
174        int j;
175        for (j=columnQuadraticStart[iColumn];
176             j<columnQuadraticStart[iColumn]+columnQuadraticLength[iColumn];j++) {
177          int jColumn = columnQuadratic[j];
178          double valueJ = solution[jColumn];
179          double valueJSave = saveSolution[jColumn];
180          double elementValue = 0.5*quadraticElement[j];
181          coeff0 += valueI*valueJ*elementValue;
182          coeff1 += (valueI*valueJSave+valueISave*valueJ-2.0*valueI*valueJ)*elementValue;
183          coeff2 += (valueISave*valueJSave+valueI*valueJ-valueISave*valueJ-valueI*valueJSave)*elementValue;
184        }
185      }
186      double lambdaValue;
187      if (fabs(coeff2)<1.0e-9) {
188        if (coeff1+coeff2>=0.0) 
189          lambda = 0.0;
190        else
191          lambda = 1.0;
192      } else {
193        lambda = -(0.5*coeff1)/coeff2;
194        if (lambda>1.0||lambda<0.0) {
195          if (coeff1+coeff2>=0.0) 
196            lambda = 0.0;
197          else
198            lambda = 1.0;
199        }
200      }
201      lambdaValue = lambda*lambda*coeff2+lambda*coeff1+coeff0;
202      printf("coeffs %g %g %g - lastobj %g\n",coeff0,coeff1,coeff2,lastObjective);
203      printf("obj at saved %g, obj now %g zero deriv at %g - value %g\n",
204             coeff0+coeff1+coeff2,coeff0,lambda,lambdaValue);
205      if (!iPass) lambda=0.0;
206      if (lambda>0.0&&lambda<=1.0) {
207        // update solution
208        for (iColumn=0;iColumn<numberColumns;iColumn++) 
209          solution[iColumn] = lambda * saveSolution[iColumn] 
210            + (1.0-lambda) * solution[iColumn];
211        if (lambda>0.999) {
212          memcpy(this->dualRowSolution(),savePi,numberRows*sizeof(double));
213          memcpy(status_,saveStatus,numberRows+numberColumns);
214        }
215        if (lambda>0.99999&&fabs(coeff1+coeff2)>1.0e-2) {
216          // tighten all
217          goodMove=-1;
218        }
219      }
220    }
221    memcpy(objective,saveObjective,numberColumns*sizeof(double));
222    for (iColumn=0;iColumn<numberColumns;iColumn++) 
223      objValue += objective[iColumn]*solution[iColumn];
224    for (jNon=0;jNon<numberNonLinearColumns;jNon++) {
225      iColumn=listNonLinearColumn[jNon];
226      if (getColumnStatus(iColumn)==basic) {
227        if (solution[iColumn]<columnLower[iColumn]+1.0e-8)
228          statusCheck[iColumn]='l';
229        else if (solution[iColumn]>columnUpper[iColumn]-1.0e-8)
230          statusCheck[iColumn]='u';
231        else
232          statusCheck[iColumn]='B';
233      } else {
234        if (solution[iColumn]<columnLower[iColumn]+1.0e-8)
235          statusCheck[iColumn]='L';
236        else
237          statusCheck[iColumn]='U';
238      }
239      double valueI = solution[iColumn];
240      int j;
241      for (j=columnQuadraticStart[iColumn];
242           j<columnQuadraticStart[iColumn]+columnQuadraticLength[iColumn];j++) {
243        int jColumn = columnQuadratic[j];
244        double valueJ = solution[jColumn];
245        double elementValue = quadraticElement[j];
246        elementValue *= 0.5;
247        objValue += valueI*valueJ*elementValue;
248        offset += valueI*valueJ*elementValue;
249        double gradientI = valueJ*elementValue;
250        double gradientJ = valueI*elementValue;
251        offset -= gradientI*valueI;
252        objective[iColumn] += gradientI;
253        offset -= gradientJ*valueJ;
254        objective[jColumn] += gradientJ;
255      }
256    }
257    printf("objective %g, objective offset %g\n",objValue,offset);
258    setDblParam(ClpObjOffset,objectiveOffset-offset);
259    objValue *= whichWay;
260    int * temp=last[2];
261    last[2]=last[1];
262    last[1]=last[0];
263    last[0]=temp;
264    for (jNon=0;jNon<numberNonLinearColumns;jNon++) {
265      iColumn=listNonLinearColumn[jNon];
266      double change = solution[iColumn]-saveSolution[iColumn];
267      if (change<-1.0e-5) {
268        if (fabs(change+trust[jNon])<1.0e-5) 
269          temp[jNon]=-1;
270        else
271          temp[jNon]=-2;
272      } else if(change>1.0e-5) {
273        if (fabs(change-trust[jNon])<1.0e-5) 
274          temp[jNon]=1;
275        else
276          temp[jNon]=2;
277      } else {
278        temp[jNon]=0;
279      }
280    } 
281    // goodMove +1 yes, 0 no, -1 last was bad - just halve gaps, -2 do nothing
282    double maxDelta=0.0;
283    if (goodMove>=0) {
284      if (objValue<=lastObjective) 
285        goodMove=1;
286      else
287        goodMove=0;
288    } else {
289      maxDelta=1.0e10;
290    }
291    double maxGap=0.0;
292    for (jNon=0;jNon<numberNonLinearColumns;jNon++) {
293      iColumn=listNonLinearColumn[jNon];
294      maxDelta = max(maxDelta,
295                     fabs(solution[iColumn]-saveSolution[iColumn]));
296      if (goodMove>0) {
297        if (last[0][jNon]*last[1][jNon]<0) {
298          // halve
299          trust[jNon] *= 0.5;
300        } else {
301          if (last[0][jNon]==last[1][jNon]&&
302              last[0][jNon]==last[2][jNon])
303            trust[jNon] *= 1.5; 
304        }
305      } else if (goodMove!=-2&&trust[jNon]>10.0*deltaTolerance) {
306        trust[jNon] *= 0.2;
307      }
308      maxGap = max(maxGap,trust[jNon]);
309    }
310    std::cout<<"largest gap is "<<maxGap<<std::endl;
311    if (iPass>10000) {
312      for (jNon=0;jNon<numberNonLinearColumns;jNon++) 
313        trust[jNon] *=0.0001;
314    }
315    if (goodMove>0) {
316      double drop = lastObjective-objValue;
317      std::cout<<"True drop was "<<drop<<std::endl;
318      std::cout<<"largest delta is "<<maxDelta<<std::endl;
319      if (maxDelta<deltaTolerance&&drop<1.0e-4&&goodMove&&lambda<0.99999) {
320        std::cout<<"Exiting"<<std::endl;
321        break;
322      }
323    }
324    if (!iPass)
325      goodMove=1;
326    targetDrop=0.0;
327    double * r = this->dualColumnSolution();
328    for (jNon=0;jNon<numberNonLinearColumns;jNon++) {
329      iColumn=listNonLinearColumn[jNon];
330      columnLower[iColumn]=max(solution[iColumn]
331                               -trust[jNon],
332                               trueLower[jNon]);
333      columnUpper[iColumn]=min(solution[iColumn]
334                               +trust[jNon],
335                               trueUpper[jNon]);
336    }
337    if (iPass) {
338      // get reduced costs
339      this->matrix()->transposeTimes(savePi,
340                                     this->dualColumnSolution());
341      for (jNon=0;jNon<numberNonLinearColumns;jNon++) {
342        iColumn=listNonLinearColumn[jNon];
343        double dj = objective[iColumn]-r[iColumn];
344        r[iColumn]=dj;
345        if (dj<0.0) 
346          targetDrop -= dj*(columnUpper[iColumn]-solution[iColumn]);
347        else
348          targetDrop -= dj*(columnLower[iColumn]-solution[iColumn]);
349      }
350    } else {
351      memset(r,0,numberColumns*sizeof(double));
352    }
353#if 0
354    for (jNon=0;jNon<numberNonLinearColumns;jNon++) {
355      iColumn=listNonLinearColumn[jNon];
356      if (statusCheck[iColumn]=='L'&&r[iColumn]<-1.0e-4) {
357        columnLower[iColumn]=max(solution[iColumn],
358                                 trueLower[jNon]);
359        columnUpper[iColumn]=min(solution[iColumn]
360                                 +trust[jNon],
361                                 trueUpper[jNon]);
362      } else if (statusCheck[iColumn]=='U'&&r[iColumn]>1.0e-4) {
363        columnLower[iColumn]=max(solution[iColumn]
364                                 -trust[jNon],
365                                 trueLower[jNon]);
366        columnUpper[iColumn]=min(solution[iColumn],
367                                 trueUpper[jNon]);
368      } else {
369        columnLower[iColumn]=max(solution[iColumn]
370                                 -trust[jNon],
371                                 trueLower[jNon]);
372        columnUpper[iColumn]=min(solution[iColumn]
373                                 +trust[jNon],
374                                 trueUpper[jNon]);
375      }
376    }
377#endif
378    if (goodMove) {
379      memcpy(saveSolution,solution,numberColumns*sizeof(double));
380      memcpy(savePi,this->dualRowSolution(),numberRows*sizeof(double));
381      memcpy(saveStatus,status_,numberRows+numberColumns);
382     
383      std::cout<<"Pass - "<<iPass
384               <<", target drop is "<<targetDrop
385               <<std::endl;
386      lastObjective = objValue;
387      if (targetDrop<1.0e-5&&goodMove&&iPass) {
388        printf("Exiting on target drop %g\n",targetDrop);
389        break;
390      }
391      {
392        double * r = this->dualColumnSolution();
393        for (jNon=0;jNon<numberNonLinearColumns;jNon++) {
394          iColumn=listNonLinearColumn[jNon];
395          printf("Trust %d %g - solution %d %g obj %g dj %g state %c - bounds %g %g\n",
396                 jNon,trust[jNon],iColumn,solution[iColumn],objective[iColumn],
397                 r[iColumn],statusCheck[iColumn],columnLower[iColumn],
398                 columnUpper[iColumn]);
399        }
400      }
401      setLogLevel(63);
402      this->scaling(false);
403      this->primal(1);
404      if (this->status()) {
405        CoinMpsIO writer;
406        writer.setMpsData(*matrix(), COIN_DBL_MAX,
407                          getColLower(), getColUpper(),
408                          getObjCoefficients(),
409                          (const char*) 0 /*integrality*/,
410                          getRowLower(), getRowUpper(),
411                          NULL,NULL);
412        writer.writeMps("xx.mps");
413      }
414      assert (!this->status()); // must be feasible
415      goodMove=1;
416    } else {
417      // bad pass - restore solution
418      printf("Backtracking\n");
419      memcpy(solution,saveSolution,numberColumns*sizeof(double));
420      memcpy(this->dualRowSolution(),savePi,numberRows*sizeof(double));
421      memcpy(status_,saveStatus,numberRows+numberColumns);
422      iPass--;
423      goodMove=-1;
424    }
425  }
426  // restore solution
427  memcpy(solution,saveSolution,numberColumns*sizeof(double));
428  for (jNon=0;jNon<numberNonLinearColumns;jNon++) {
429    iColumn=listNonLinearColumn[jNon];
430    columnLower[iColumn]=max(solution[iColumn],
431                             trueLower[jNon]);
432    columnUpper[iColumn]=min(solution[iColumn],
433                             trueUpper[jNon]);
434  }
435  delete [] statusCheck;
436  delete [] savePi;
437  delete [] saveStatus;
438  // redo objective
439  double offset=0.0;
440  double objValue=-objectiveOffset;
441  memcpy(objective,saveObjective,numberColumns*sizeof(double));
442  for (iColumn=0;iColumn<numberColumns;iColumn++) 
443    objValue += objective[iColumn]*solution[iColumn];
444  for (jNon=0;jNon<numberNonLinearColumns;jNon++) {
445    iColumn=listNonLinearColumn[jNon];
446    double valueI = solution[iColumn];
447    int j;
448    for (j=columnQuadraticStart[iColumn];
449         j<columnQuadraticStart[iColumn]+columnQuadraticLength[iColumn];j++) {
450      int jColumn = columnQuadratic[j];
451      double valueJ = solution[jColumn];
452      double elementValue = quadraticElement[j];
453      objValue += 0.5*valueI*valueJ*elementValue;
454      offset += 0.5*valueI*valueJ*elementValue;
455      double gradientI = valueJ*elementValue;
456      double gradientJ = valueI*elementValue;
457      offset -= gradientI*valueI;
458      objective[iColumn] += gradientI;
459      offset -= gradientJ*valueJ;
460      objective[jColumn] += gradientJ;
461    }
462  }
463  printf("objective %g, objective offset %g\n",objValue,offset);
464  setDblParam(ClpObjOffset,objectiveOffset-offset);
465  this->primal(1);
466  // redo values
467  setDblParam(ClpObjOffset,objectiveOffset);
468  objectiveValue_ += whichWay*offset;
469  for (jNon=0;jNon<numberNonLinearColumns;jNon++) {
470    iColumn=listNonLinearColumn[jNon];
471    columnLower[iColumn]= trueLower[jNon];
472    columnUpper[iColumn]= trueUpper[jNon];
473  }
474  memcpy(objective,saveObjective,numberColumns*sizeof(double));
475  delete [] saveSolution;
476  for (iPass=0;iPass<3;iPass++) 
477    delete [] last[iPass];
478  delete [] trust;
479  delete [] trueUpper;
480  delete [] trueLower;
481  delete [] saveObjective;
482  delete [] listNonLinearColumn;
483  return 0;
484}
485// Dantzig's method
486int 
487ClpSimplexPrimalQuadratic::primalQuadratic(int phase)
488{
489  // Get a feasible solution
490  if (numberPrimalInfeasibilities())
491    primal(1);
492  // still infeasible
493  if (numberPrimalInfeasibilities())
494    return 1;
495  ClpQuadraticInfo info;
496  ClpSimplexPrimalQuadratic * model2 = makeQuadratic(info);
497  if (!model2) {
498    printf("** Not quadratic\n");
499    primal(1);
500    return 0;
501  }
502#if 0
503  CoinMpsIO writer;
504  writer.setMpsData(*model2->matrix(), COIN_DBL_MAX,
505                    model2->getColLower(), model2->getColUpper(),
506                    model2->getObjCoefficients(),
507                    (const char*) 0 /*integrality*/,
508                    model2->getRowLower(), model2->getRowUpper(),
509                    NULL,NULL);
510  writer.writeMps("xx.mps");
511#endif 
512  // Now do quadratic
513  //ClpPrimalQuadraticDantzig dantzigP(model2,&info);
514  ClpPrimalColumnDantzig dantzigP;
515  model2->setPrimalColumnPivotAlgorithm(dantzigP);
516  model2->messageHandler()->setLogLevel(63);
517  model2->primalQuadratic2(&info,phase);
518  endQuadratic(model2,info);
519  return 0;
520}
521int ClpSimplexPrimalQuadratic::primalQuadratic2 (ClpQuadraticInfo * info,
522                                                 int phase)
523{
524
525  algorithm_ = +2;
526
527  // save data
528  ClpDataSave data = saveData();
529  // Only ClpPackedMatrix at present
530  assert ((dynamic_cast< ClpPackedMatrix*>(matrix_)));
531 
532  // Assume problem is feasible
533  // Stuff below will be moved into a class
534  int numberXColumns = info->numberXColumns();
535  int numberXRows = info->numberXRows();
536  // initialize - values pass coding and algorithm_ is +1
537  ClpObjective * saveObj = objectiveAsObject();
538  setObjectivePointer(info->originalObjective());
539  factorization_->setBiasLU(0);
540  if (!startup(1)) {
541
542    //setObjectivePointer(saveObj);
543    // Setup useful stuff
544    info->setCurrentPhase(phase);
545    int lastCleaned=0; // last time objective or bounds cleaned up
546    info->setSequenceIn(-1);
547    info->setCrucialSj(-1);
548    bool deleteCosts=false;
549    if (scalingFlag_>0) {
550      // scale
551      CoinPackedMatrix * quadratic = info->quadraticObjective();
552      double * objective = info->linearObjective();
553      // replace column by column
554      double * newElement = new double[numberXColumns];
555      // scale column copy
556      // get matrix data pointers
557      const int * row = quadratic->getIndices();
558      const CoinBigIndex * columnStart = quadratic->getVectorStarts();
559      const int * columnLength = quadratic->getVectorLengths(); 
560      const double * elementByColumn = quadratic->getElements();
561      double direction = optimizationDirection_;
562      // direction is actually scale out not scale in
563      if (direction)
564        direction = 1.0/direction;
565      int iColumn;
566      for (iColumn=0;iColumn<numberXColumns;iColumn++) {
567        int j;
568        double scale = columnScale_[iColumn]*direction;
569        const double * elementsInThisColumn = elementByColumn + columnStart[iColumn];
570        const int * columnsInThisColumn = row + columnStart[iColumn];
571        int number = columnLength[iColumn];
572        assert (number<=numberXColumns);
573        for (j=0;j<number;j++) {
574          int iColumn2 = columnsInThisColumn[j];
575          newElement[j] = elementsInThisColumn[j]*scale*rowScale_[iColumn2+numberXRows];
576        }
577        quadratic->replaceVector(iColumn,number,newElement);
578        objective[iColumn] *= direction*columnScale_[iColumn];
579      }
580      delete [] newElement;
581      deleteCosts=true;
582    } else if (optimizationDirection_!=1.0) {
583      CoinPackedMatrix * quadratic = info->quadraticObjective();
584      double * objective = info->linearObjective();
585      // replace column by column
586      double * newElement = new double[numberXColumns];
587      // get matrix data pointers
588      const CoinBigIndex * columnStart = quadratic->getVectorStarts();
589      const int * columnLength = quadratic->getVectorLengths(); 
590      const double * elementByColumn = quadratic->getElements();
591      double direction = optimizationDirection_;
592      // direction is actually scale out not scale in
593      if (direction)
594        direction = 1.0/direction;
595      int iColumn;
596      for (iColumn=0;iColumn<numberXColumns;iColumn++) {
597        int j;
598        const double * elementsInThisColumn = elementByColumn + columnStart[iColumn];
599        int number = columnLength[iColumn];
600        assert (number<=numberXColumns);
601        for (j=0;j<number;j++) {
602          newElement[j] = elementsInThisColumn[j]*direction;
603        }
604        quadratic->replaceVector(iColumn,number,newElement);
605        objective[iColumn] *= direction;
606      }
607      delete [] newElement;
608      deleteCosts=true;
609    }
610   
611    // Say no pivot has occurred (for steepest edge and updates)
612    pivotRow_=-2;
613   
614    // This says whether to restore things etc
615    int factorType=0;
616    /*
617      Status of problem:
618      0 - optimal
619      1 - infeasible
620      2 - unbounded
621      -1 - iterating
622      -2 - factorization wanted
623      -3 - redo checking without factorization
624      -4 - looks infeasible
625      -5 - looks unbounded
626    */
627    while (problemStatus_<0) {
628      int iRow,iColumn;
629      // clear
630      for (iRow=0;iRow<4;iRow++) {
631        rowArray_[iRow]->clear();
632      }   
633     
634      for (iColumn=0;iColumn<2;iColumn++) {
635        columnArray_[iColumn]->clear();
636      }   
637     
638      // give matrix (and model costs and bounds a chance to be
639      // refreshed (normally null)
640      matrix_->refresh(this);
641      // If we have done no iterations - special
642      if (lastGoodIteration_==numberIterations_)
643        factorType=3;
644      // may factorize, checks if problem finished
645      statusOfProblemInPrimal(lastCleaned,factorType,progress_,info);
646      if (problemStatus_>=0)
647        break; // declare victory
648     
649      checkComplementarity (info,rowArray_[0],rowArray_[1]);
650
651      // Say good factorization
652      factorType=1;
653     
654      // Say no pivot has occurred (for steepest edge and updates)
655      pivotRow_=-2;
656      // Check problem phase
657      // We assume all X are feasible
658      if (info->currentSolution()&&info->sequenceIn()<0) {
659        phase=0;
660        int nextSequenceIn=-1;
661        int numberQuadraticRows = info->numberQuadraticRows();
662        for (iColumn=0;iColumn<numberXColumns;iColumn++) {
663          double value = solution_[iColumn];
664          double lower = lower_[iColumn];
665          double upper = upper_[iColumn];
666          if ((value>lower+primalTolerance_&&value<upper-primalTolerance_)
667              ||getColumnStatus(iColumn)==superBasic) {
668            if (getColumnStatus(iColumn)!=basic&&!flagged(iColumn)) {
669              if (fabs(dj_[iColumn])>10.0*dualTolerance_||
670                  getColumnStatus(iColumn+numberXColumns+numberQuadraticRows)==basic) {
671                if (phase!=2) {
672                  phase=2;
673                  nextSequenceIn=iColumn;
674                }
675              }
676            }
677          }
678        }
679        for (iColumn=numberColumns_;iColumn<numberColumns_+numberXRows;iColumn++) {
680          double value = solution_[iColumn];
681          double lower = lower_[iColumn];
682          double upper = upper_[iColumn];
683          if ((value>lower+primalTolerance_&&value<upper-primalTolerance_)
684            ||getColumnStatus(iColumn)==superBasic) {
685            if (getColumnStatus(iColumn)!=basic&&!flagged(iColumn)&&fabs(dj_[iColumn])>10.0*dualTolerance_) {
686              if (phase!=2) {
687                phase=2;
688                nextSequenceIn=iColumn;
689              }
690            }
691          }
692        }
693        info->setSequenceIn(nextSequenceIn);
694        info->setCurrentPhase(phase);
695      }
696
697      // exit if victory declared
698      dualIn_=0.0; // so no updates
699      if (!info->currentPhase()&&info->sequenceIn()<0&&primalColumnPivot_->pivotColumn(rowArray_[1],
700                                          rowArray_[2],rowArray_[3],
701                                          columnArray_[0],
702                                          columnArray_[1]) < 0) {
703        problemStatus_=0;
704        break;
705      }
706     
707      // Iterate
708      problemStatus_=-1;
709      whileIterating(info);
710    }
711    if (deleteCosts) {
712      // delete scaled copy of objective
713      delete info->originalObjective();
714    }
715  }
716  setObjectivePointer(saveObj);
717  // clean up
718  finish();
719  restoreData(data);
720  return problemStatus_;
721}
722/*
723  Reasons to come out:
724  -1 iterations etc
725  -2 inaccuracy
726  -3 slight inaccuracy (and done iterations)
727  -4 end of values pass and done iterations
728  +0 looks optimal (might be infeasible - but we will investigate)
729  +2 looks unbounded
730  +3 max iterations
731*/
732int
733ClpSimplexPrimalQuadratic::whileIterating(
734                      ClpQuadraticInfo * info)
735{
736  checkComplementarity (info,rowArray_[0],rowArray_[1]);
737  int crucialSj = info->crucialSj();
738  if (info->crucialSj()>=0) {
739    printf("after inv %d Sj value %g\n",crucialSj,solution_[info->crucialSj()]);
740  }
741  int returnCode=-1;
742  int phase = info->currentPhase();
743  double saveObjective = objectiveValue_;
744  int numberXColumns = info->numberXColumns();
745  int numberXRows = info->numberXRows();
746  int numberQuadraticRows = info->numberQuadraticRows();
747  // Make list of implied sj
748  // And backward pointers to basic variables
749  {
750    int * impliedSj = info->impliedSj();
751    int i;
752    for (i=numberRows_;i<numberColumns_;i++) {
753      if (getColumnStatus(i)==basic)
754        impliedSj[i-numberRows_]=i;
755      else
756        impliedSj[i-numberRows_]=-1;
757    }
758    int * basicRow = info->basicRow();
759    for (i=0;i<numberRows_+numberColumns_;i++)
760      basicRow[i]=-1;
761    for (i=0;i<numberRows_;i++)
762      basicRow[pivotVariable_[i]]=i;
763  }
764  int nextSequenceIn=info->sequenceIn();
765  int oldSequenceIn=nextSequenceIn;
766  int saveSequenceIn = sequenceIn_;
767  // status stays at -1 while iterating, >=0 finished, -2 to invert
768  // status -3 to go to top without an invert
769  while (problemStatus_==-1) {
770    if (info->crucialSj()<0&&factorization_->pivots()>=10) {
771      returnCode =-2; // refactorize
772      break;
773    }
774#ifdef CLP_DEBUG
775    {
776      int i;
777      // not [1] as has information
778      for (i=0;i<4;i++) {
779        if (i!=1)
780          rowArray_[i]->checkClear();
781      }   
782      for (i=0;i<2;i++) {
783        columnArray_[i]->checkClear();
784      }   
785    }     
786#endif
787    // choose column to come in
788    // can use pivotRow_ to update weights
789    // pass in list of cost changes so can do row updates (rowArray_[1])
790    // NOTE rowArray_[0] is used by computeDuals which is a
791    // slow way of getting duals but might be used
792    // Initially Dantzig and look at s variables
793    // Only do if one not already chosen
794    int cleanupIteration;
795    if (nextSequenceIn<0) {
796      cleanupIteration=0;
797      if (phase==2) {
798        // values pass
799        // get next
800        int iColumn;
801        int iStart = oldSequenceIn+1;
802        for (iColumn=iStart;iColumn<numberXColumns;iColumn++) {
803          double value = solution_[iColumn];
804          double lower = lower_[iColumn];
805          double upper = upper_[iColumn];
806          if (value>lower+primalTolerance_&&value<upper-primalTolerance_) {
807            if (getColumnStatus(iColumn)!=basic&&!flagged(iColumn)) {
808              if (fabs(dj_[iColumn])>10.0*dualTolerance_||
809                  getColumnStatus(iColumn+numberXColumns+numberQuadraticRows)==basic) {
810                nextSequenceIn=iColumn;
811                break;
812              }
813            }
814          }
815        }
816        if (nextSequenceIn<0) {
817          iStart=max(iStart,numberColumns_);
818          int numberXRows = info->numberXRows();
819          for (iColumn=iStart;iColumn<numberColumns_+numberXRows;
820               iColumn++) {
821            double value = solution_[iColumn];
822            double lower = lower_[iColumn];
823            double upper = upper_[iColumn];
824            if (value>lower+primalTolerance_&&value<upper-primalTolerance_) {
825              if (getColumnStatus(iColumn)!=basic&&!flagged(iColumn)&&fabs(dj_[iColumn])>10.0*dualTolerance_) {
826                nextSequenceIn=iColumn;
827                break;
828              }
829            }
830          }
831        }
832        oldSequenceIn=nextSequenceIn;
833        saveSequenceIn=sequenceIn_;
834        sequenceIn_ = nextSequenceIn;
835      } else {
836        saveSequenceIn=sequenceIn_;
837        createDjs(info,rowArray_[3],rowArray_[2]);
838        dualIn_=0.0; // so no updates
839        rowArray_[1]->clear();
840        primalColumn(rowArray_[1],rowArray_[2],rowArray_[3],
841                     columnArray_[0],columnArray_[1]);
842      }
843    } else {
844      saveSequenceIn=sequenceIn_;
845      sequenceIn_ = nextSequenceIn;
846      if (phase==2) {
847        if ((sequenceIn_<numberXColumns||sequenceIn_>=numberColumns_)&&info->crucialSj()<0) 
848          cleanupIteration=0;
849        else
850          cleanupIteration=1;
851      } else {
852        cleanupIteration=1;
853      }
854    }
855    pivotRow_=-1;
856    sequenceOut_=-1;
857    rowArray_[1]->clear();
858    if (sequenceIn_>=0) {
859      nextSequenceIn=-1;
860      // we found a pivot column
861      int chosen = sequenceIn_;
862      // do second half of iteration
863      while (chosen>=0) {
864        int saveSequenceInInfo=chosen;
865        int saveCrucialSjInfo=info->crucialSj();
866        rowArray_[1]->clear();
867        checkComplementarity (info,rowArray_[3],rowArray_[1]);
868        printf("True objective is %g, infeas cost %g, sum %g\n",
869               objectiveValue_,info->infeasCost(),objectiveValue_+info->infeasCost());
870        objectiveValue_=saveObjective;
871        returnCode=-1;
872        pivotRow_=-1;
873        sequenceOut_=-1;
874        // we found a pivot column
875        // update the incoming column
876        sequenceIn_=chosen;
877        chosen=-1;
878        unpack(rowArray_[1]);
879        // compute dj in case linear
880        {
881          info->createGradient(this);
882          double * gradient = info->gradient();
883          dualIn_=gradient[sequenceIn_];
884          int j;
885          const double * element = rowArray_[1]->denseVector();
886          int numberNonZero=rowArray_[1]->getNumElements();
887          const int * whichRow = rowArray_[1]->getIndices();
888          bool packed = rowArray_[1]->packedMode();
889          if (packed) {
890            for (j=0;j<numberNonZero; j++) {
891              int iRow = whichRow[j];
892              dualIn_ -= element[j]*dual_[iRow];
893            }
894          } else {
895            for (j=0;j<numberNonZero; j++) {
896              int iRow = whichRow[j];
897              dualIn_ -= element[iRow]*dual_[iRow];
898            }
899          }
900        }
901        if ((!cleanupIteration&&sequenceIn_<numberXColumns)||
902            (cleanupIteration&&info->crucialSj()<numberColumns_&&info->crucialSj()>=numberRows_)) {
903          int iSequence;
904          if (!cleanupIteration)
905            iSequence=sequenceIn_;
906          else
907            iSequence=info->crucialSj()-numberRows_;
908          // may need to re-do row of basis
909          int * impliedSj = info->impliedSj();
910          if (impliedSj[iSequence]>=0) {
911            // mark as valid
912            impliedSj[iSequence]=-1;
913            ClpPackedMatrix* rowCopy =
914              dynamic_cast< ClpPackedMatrix*>(rowCopy_);
915            assert(rowCopy);
916            const int * column = rowCopy->getIndices();
917            const CoinBigIndex * rowStart = rowCopy->getVectorStarts();
918            const double * rowElement = rowCopy->getElements();
919            int iRow=iSequence+numberXRows;
920            int i;
921            int * index = rowArray_[2]->getIndices();
922            double * element = rowArray_[2]->denseVector();
923            int n=0;
924            int * basicRow = info->basicRow();
925            int * permute = factorization_->pivotColumn();
926            for (i=rowStart[iRow];i<rowStart[iRow+1];i++) {
927              int iColumn = column[i];
928              iColumn = basicRow[iColumn];
929              if (iColumn>=0&&iColumn!=iRow) {
930                index[n]=permute[iColumn];
931                element[n++]=rowElement[i];
932              }
933            }
934            factorization_->replaceRow(permute[iRow],n,index,element);
935            memset(element,0,n*sizeof(double));
936          }
937        }
938        // Take out elements in implied Sj rows
939        if (1) { 
940          int n=rowArray_[1]->getNumElements();
941          int * index = rowArray_[1]->getIndices();
942          double * element = rowArray_[1]->denseVector();
943          int i;
944          int * impliedSj = info->impliedSj();
945          int n2=0;
946          for (i=0;i<n;i++) {
947            int iRow=index[i]-numberXRows;
948            if (iRow<0||impliedSj[iRow]<0) {
949              index[n2++]=iRow+numberXRows;
950            } else {
951              element[iRow+numberXRows]=0.0;
952            }
953          }
954          rowArray_[1]->setNumElements(n2);
955        }
956        factorization_->updateColumnFT(rowArray_[2],rowArray_[1]);
957        if (cleanupIteration) {
958          // move back to a version of primalColumn?
959          valueIn_=solution_[sequenceIn_];
960          // should keep pivot row of crucialSj as well (speed)
961          int iSjRow=-1;
962          if (crucialSj>=0) {
963            double * work=rowArray_[1]->denseVector();
964            int number=rowArray_[1]->getNumElements();
965            int * which=rowArray_[1]->getIndices();
966            double tj = 0.0;
967            int iIndex;
968            for (iIndex=0;iIndex<number;iIndex++) {
969              int iRow = which[iIndex];
970              double alpha = work[iRow];
971              int iPivot=pivotVariable_[iRow];
972              if (iPivot==crucialSj) {
973                tj = alpha;
974                iSjRow = iRow;
975                double d2 = solution_[crucialSj]/tj;
976                if (fabs(solution_[crucialSj])<1.0e-8)
977                  printf("zero crucialSj - pivot out - get right way\n");
978                // see which way to go
979                if (d2>0)
980                  dj_[sequenceIn_]= -1.0;
981                else
982                  dj_[sequenceIn_]= 1.0;
983                break;
984              }
985            }
986            if (!tj) {
987              printf("trouble\n");
988              assert (sequenceIn_>numberXColumns&&sequenceIn_<numberColumns_);
989              dj_[sequenceIn_]=solution_[sequenceIn_];
990            //assert(tj);
991            }
992          }
993
994          dualIn_=dj_[sequenceIn_];
995          lowerIn_=lower_[sequenceIn_];
996          upperIn_=upper_[sequenceIn_];
997          if (dualIn_>0.0)
998            directionIn_ = -1;
999          else 
1000            directionIn_ = 1;
1001        } else {
1002          // not clean up
1003          lowerIn_=lower_[sequenceIn_];
1004          upperIn_=upper_[sequenceIn_];
1005          dualIn_=dj_[sequenceIn_];
1006          valueIn_=solution_[sequenceIn_];
1007          if (dualIn_>0.0)
1008            directionIn_ = -1;
1009          else 
1010            directionIn_ = 1;
1011          if (sequenceIn_<numberColumns_) {
1012            // Set dj as value of slack
1013            crucialSj = sequenceIn_+ numberQuadraticRows+numberXColumns; // sj which should go to 0.0
1014          } else {
1015            // Set dj as value of pi
1016            int iRow = sequenceIn_-numberColumns_;
1017            crucialSj = iRow+numberXColumns; // pi which should go to 0.0
1018          }
1019          if (crucialSj>=0&&getColumnStatus(crucialSj)!=basic)
1020            crucialSj=-1;
1021          info->setCrucialSj(crucialSj);
1022        }
1023        double oldSj=1.0e30;
1024        if (info->crucialSj()>=0&&cleanupIteration)
1025          oldSj= solution_[info->crucialSj()];
1026        // save reduced cost
1027        //double saveDj = dualIn_;
1028        // do ratio test and re-compute dj
1029        // Note second parameter long enough for columns
1030        int result=primalRow(rowArray_[1],rowArray_[3],
1031                             rowArray_[2],rowArray_[0],
1032                             info,
1033                             cleanupIteration);
1034        if (pivotRow_==-1&&(phase==2||cleanupIteration)&&fabs(dualIn_)<1.0e-3) {
1035          // try other way
1036          dualIn_=-dualIn_;
1037          directionIn_=-directionIn_;
1038          if (info->crucialSj()>=0)
1039            setColumnStatus(info->crucialSj(),basic);
1040          result=primalRow(rowArray_[1],rowArray_[3],
1041                           rowArray_[2],rowArray_[0],
1042                           info,
1043                           cleanupIteration);
1044        }
1045        saveObjective = objectiveValue_;
1046        if (pivotRow_>=0) {
1047          // If sj out AND not gone out since last invert then add back
1048          // if stable replace in basis
1049          int updateStatus = 0;
1050          if (result<20) {
1051            double saveCheck = factorization_->getAccuracyCheck();
1052            if (cleanupIteration)
1053              factorization_->relaxAccuracyCheck(1.0e3*saveCheck);
1054            updateStatus=factorization_->replaceColumn(rowArray_[2],
1055                                                       pivotRow_,
1056                                                       alpha_);
1057            factorization_->relaxAccuracyCheck(saveCheck);
1058          }
1059          if (result>=10) {
1060            updateStatus = max(updateStatus,result/10);
1061            result = result%10;
1062            if (updateStatus>1) {
1063              alpha_=0.0; // force re-factorization
1064              info->setSequenceIn(sequenceIn_);
1065            }
1066          }
1067          // if no pivots, bad update but reasonable alpha - take and invert
1068          if (updateStatus==2&&
1069              lastGoodIteration_==numberIterations_&&fabs(alpha_)>1.0e-5)
1070            updateStatus=4;
1071          if (updateStatus==1||updateStatus==4) {
1072            // slight error
1073            if (factorization_->pivots()>5||updateStatus==4) {
1074              returnCode=-3;
1075            }
1076          } else if (updateStatus==2) {
1077            // major error
1078            // Reset sequenceIn_
1079            // sequenceIn_=saveSequenceIn;
1080            nextSequenceIn=saveSequenceInInfo;
1081            info->setCrucialSj(saveCrucialSjInfo);
1082            if (saveCrucialSjInfo<0&&!phase)
1083              nextSequenceIn=-1;
1084            pivotRow_=-1;
1085            // better to have small tolerance even if slower
1086            factorization_->zeroTolerance(1.0e-15);
1087            int maxFactor = factorization_->maximumPivots();
1088            if (maxFactor>10) {
1089              if (forceFactorization_<0)
1090                forceFactorization_= maxFactor;
1091              forceFactorization_ = max (1,(forceFactorization_>>1));
1092            } 
1093            // later we may need to unwind more e.g. fake bounds
1094            if(lastGoodIteration_ != numberIterations_||factorization_->pivots()) {
1095              rowArray_[1]->clear();
1096              pivotRow_=-1;
1097              returnCode=-4;
1098              // retry on return
1099              if (info->crucialSj()>=0)
1100                nextSequenceIn = sequenceIn_;
1101              break;
1102            } else {
1103              // need to reject something
1104              char x = isColumn(sequenceIn_) ? 'C' :'R';
1105              handler_->message(CLP_SIMPLEX_FLAG,messages_)
1106                <<x<<sequenceWithin(sequenceIn_)
1107                <<CoinMessageEol;
1108              setFlagged(sequenceIn_);
1109              lastBadIteration_ = numberIterations_; // say be more cautious
1110              rowArray_[1]->clear();
1111              pivotRow_=-1;
1112              returnCode = -5;
1113              break;
1114             
1115            }
1116          } else if (updateStatus==3) {
1117            // out of memory
1118            // increase space if not many iterations
1119            if (factorization_->pivots()<
1120                0.5*factorization_->maximumPivots()&&
1121                factorization_->pivots()<200)
1122              factorization_->areaFactor(
1123                                         factorization_->areaFactor() * 1.1);
1124            returnCode =-2; // factorize now
1125          }
1126          // here do part of steepest - ready for next iteration
1127          primalColumnPivot_->updateWeights(rowArray_[1]);
1128        } else {
1129          if (pivotRow_==-1) {
1130            // no outgoing row is valid
1131            rowArray_[0]->clear();
1132            if (!factorization_->pivots()) {
1133              returnCode = 2; //say looks unbounded
1134              // do ray
1135              primalRay(rowArray_[1]);
1136            } else {
1137              returnCode = 4; //say looks unbounded but has iterated
1138            }
1139            break;
1140          } else {
1141            // flipping from bound to bound
1142          }
1143        }
1144
1145
1146        // update primal solution
1147
1148        double objectiveChange=0.0;
1149        // Cost on pivot row may change - may need to change dualIn
1150        double oldCost=0.0;
1151        if (pivotRow_>=0)
1152          oldCost = cost(pivotVariable_[pivotRow_]);
1153        // rowArray_[1] is not empty - used to update djs
1154        updatePrimalsInPrimal(rowArray_[1],theta_, objectiveChange,1);
1155        if (pivotRow_>=0)
1156          dualIn_ += (oldCost-cost(pivotVariable_[pivotRow_]));
1157        double oldValue = valueIn_;
1158        if (directionIn_==-1) {
1159          // as if from upper bound
1160          if (sequenceIn_!=sequenceOut_) {
1161            // variable becoming basic
1162            valueIn_ -= fabs(theta_);
1163          } else {
1164            valueIn_=lowerIn_;
1165          }
1166        } else {
1167          // as if from lower bound
1168          if (sequenceIn_!=sequenceOut_) {
1169            // variable becoming basic
1170            valueIn_ += fabs(theta_);
1171          } else {
1172            valueIn_=upperIn_;
1173          }
1174        }
1175        objectiveChange += dualIn_*(valueIn_-oldValue);
1176        // outgoing
1177        if (sequenceIn_!=sequenceOut_) {
1178          if (directionOut_>0) {
1179            valueOut_ = lowerOut_;
1180          } else {
1181            valueOut_ = upperOut_;
1182          }
1183          double lowerValue = lower_[sequenceOut_];
1184          double upperValue = upper_[sequenceOut_];
1185          if (sequenceOut_>=numberXColumns&&sequenceOut_<numberColumns_) {
1186            // really free but going to zero
1187            lowerValue=0.0;
1188            upperValue=0.0;
1189          }
1190          assert(valueOut_>=lowerValue-primalTolerance_&&
1191                 valueOut_<=upperValue+primalTolerance_);
1192          // may not be exactly at bound and bounds may have changed
1193#if 1
1194          // Make sure outgoing looks feasible
1195          double useValue=valueOut_;
1196          if (valueOut_<=lowerValue+primalTolerance_) {
1197            directionOut_=1;
1198            useValue= lower_[sequenceOut_];
1199          } else if (valueOut_>=upperValue-primalTolerance_) {
1200            directionOut_=-1;
1201            useValue= upper_[sequenceOut_];
1202          } else {
1203            printf("*** variable wandered off bound %g %g %g!\n",
1204                   lowerValue,valueOut_,upperValue);
1205            if (valueOut_-lowerValue<=upperValue-valueOut_) {
1206              useValue= lower_[sequenceOut_];
1207              valueOut_=lowerValue;
1208              directionOut_=1;
1209            } else {
1210              useValue= upper_[sequenceOut_];
1211              valueOut_=upperValue;
1212              directionOut_=-1;
1213            }
1214          }
1215          solution_[sequenceOut_]=valueOut_;
1216          nonLinearCost_->setOne(sequenceOut_,useValue);
1217#else
1218          directionOut_=nonLinearCost_->setOneOutgoing(sequenceOut_,valueOut_);
1219          solution_[sequenceOut_]=valueOut_;
1220#endif
1221        }
1222        // change cost and bounds on incoming if primal
1223        if (sequenceIn_==sequenceOut_&&(lowerOut_!=lowerIn_||upperOut_!=upperIn_)) {
1224          // linear variable superbasic
1225          setStatus(sequenceOut_,superBasic);
1226        }
1227        nonLinearCost_->setOne(sequenceIn_,valueIn_); 
1228        int whatNext=housekeeping(objectiveChange);
1229        // and again as housekeeping may have changed
1230        if (sequenceIn_==sequenceOut_&&(lowerOut_!=lowerIn_||upperOut_!=upperIn_)) {
1231          // linear variable superbasic
1232          setStatus(sequenceOut_,superBasic);
1233        }
1234        // And update backward pointers to basic variables
1235        if (sequenceIn_!=sequenceOut_) {
1236          int * basicRow = info->basicRow();
1237          basicRow[sequenceOut_]=-1;
1238          basicRow[sequenceIn_]=pivotRow_;
1239        }
1240        if (info->crucialSj()>=0) {
1241          assert(fabs(oldSj)>= fabs(solution_[info->crucialSj()]));
1242          printf("%d Sj value went from %g to %g\n",crucialSj,oldSj,solution_[info->crucialSj()]);
1243        }
1244        checkComplementarity (info,rowArray_[2],rowArray_[3]);
1245        printf("After Housekeeping True objective is %g, infeas cost %g, sum %g\n",
1246               objectiveValue_,info->infeasCost(),objectiveValue_+info->infeasCost());
1247        if (sequenceOut_>=numberXColumns&&sequenceOut_<numberColumns_) {
1248          // really free but going to zero
1249          setStatus(sequenceOut_,isFree);
1250          if (sequenceOut_==info->crucialSj())
1251            info->setCrucialSj(-1);
1252        } else if (info->crucialSj()>=0) {
1253          // Just possible crucialSj still in but original out again
1254          int crucialSj=info->crucialSj();
1255          if (crucialSj>=numberXColumns+numberQuadraticRows) {
1256            if (sequenceOut_==crucialSj-numberXColumns-numberQuadraticRows)
1257              info->setCrucialSj(-1);
1258          } else {
1259            if (sequenceOut_-numberColumns_==crucialSj-numberXColumns)
1260              info->setCrucialSj(-1);
1261          }
1262        }
1263        if (info->crucialSj()<0)
1264          result=0;
1265        if (whatNext==1) {
1266          returnCode =-2; // refactorize
1267        } else if (whatNext==2) {
1268          // maximum iterations or equivalent
1269          returnCode=3;
1270        } else if(numberIterations_ == lastGoodIteration_
1271                  + 2 * factorization_->maximumPivots()) {
1272          // done a lot of flips - be safe
1273          returnCode =-2; // refactorize
1274        }
1275        // may need to go round again
1276        cleanupIteration = 1;
1277        // may not be correct on second time
1278        if (result&&(returnCode==-1||returnCode==-2||returnCode==-3)) {
1279          int crucialSj=info->crucialSj();
1280          if (sequenceOut_<numberXColumns) {
1281            chosen =sequenceOut_ + numberQuadraticRows + numberXColumns; // sj variable
1282          } else if (sequenceOut_>=numberColumns_) {
1283            // Does this mean we can change pi
1284            int iRow = sequenceOut_-numberColumns_;
1285            if (iRow<numberXRows) {
1286              chosen=iRow+numberXColumns;
1287            } else {
1288              printf("Row %d is in column part\n",iRow);
1289              abort();
1290            }
1291          } else if (sequenceOut_<numberQuadraticRows+numberXColumns) {
1292            // pi out
1293            chosen = sequenceOut_-numberXColumns+numberColumns_;
1294          } else {
1295            // sj out
1296            chosen = sequenceOut_-numberQuadraticRows-numberXColumns;
1297          }
1298          // But double check as original incoming might have gone out
1299          if (chosen==crucialSj) {
1300            chosen=-1;
1301            break; // means original has gone out
1302          }
1303          rowArray_[0]->clear();
1304          columnArray_[0]->clear();
1305          if (returnCode==-2) {
1306            // factorization
1307            nextSequenceIn = chosen;
1308            break;
1309          }
1310        } else {
1311          break;
1312        }
1313      }
1314      if (returnCode<-1&&returnCode>-5) {
1315        problemStatus_=-2; //
1316      } else if (returnCode==-5) {
1317        // something flagged - continue;
1318      } else if (returnCode==2) {
1319        problemStatus_=-5; // looks unbounded
1320      } else if (returnCode==4) {
1321        problemStatus_=-2; // looks unbounded but has iterated
1322      } else if (returnCode!=-1) {
1323        assert(returnCode==3);
1324        problemStatus_=3;
1325      }
1326    } else {
1327      // no pivot column
1328#ifdef CLP_DEBUG
1329      if (handler_->logLevel()&32)
1330        printf("** no column pivot\n");
1331#endif
1332      if (nonLinearCost_->numberInfeasibilities())
1333        problemStatus_=-4; // might be infeasible
1334      returnCode=0;
1335      break;
1336    }
1337  }
1338  info->setSequenceIn(nextSequenceIn);
1339  return returnCode;
1340}
1341/*
1342   Row array has pivot column
1343   This chooses pivot row.
1344   For speed, we may need to go to a bucket approach when many
1345   variables go through bounds
1346   On exit rhsArray will have changes in costs of basic variables
1347*/
1348int 
1349ClpSimplexPrimalQuadratic::primalRow(CoinIndexedVector * rowArray,
1350                                     CoinIndexedVector * rhsArray,
1351                                     CoinIndexedVector * spareArray,
1352                                     CoinIndexedVector * spareArray2,
1353                                     ClpQuadraticInfo * info,
1354                                     int cleanupIteration)
1355{
1356  int result=1;
1357 
1358  // sequence stays as row number until end
1359  pivotRow_=-1;
1360  int numberSwapped=0;
1361  int numberRemaining=0;
1362  int crucialSj = info->crucialSj();
1363  if (crucialSj<0)
1364    result=0; // linear
1365  int numberThru =0; // number gone thru a barrier
1366  int lastThru =0; // number gone thru a barrier on last time
1367 
1368  double totalThru=0.0; // for when variables flip
1369  double acceptablePivot=1.0e-7;
1370  if (factorization_->pivots())
1371    acceptablePivot=1.0e-5; // if we have iterated be more strict
1372  double bestEverPivot=acceptablePivot;
1373  int lastPivotRow = -1;
1374  double lastPivot=0.0;
1375  double lastTheta=1.0e50;
1376  int lastNumberSwapped=0;
1377
1378  // use spareArrays to put ones looked at in
1379  // First one is list of candidates
1380  // We could compress if we really know we won't need any more
1381  // Second array has current set of pivot candidates
1382  // with a backup list saved in double * part of indexed vector
1383
1384  // for zeroing out arrays after
1385  int maximumSwapped=0;
1386  // pivot elements
1387  double * spare;
1388  // indices
1389  int * index, * indexSwapped;
1390  int * saveSwapped;
1391  spareArray->clear();
1392  spareArray2->clear();
1393  spare = spareArray->denseVector();
1394  index = spareArray->getIndices();
1395  saveSwapped = (int *) spareArray2->denseVector();
1396  indexSwapped = spareArray2->getIndices();
1397
1398  // we also need somewhere for effective rhs
1399  double * rhs=rhsArray->denseVector();
1400
1401  /*
1402    First we get a list of possible pivots.  We can also see if the
1403    problem looks unbounded.
1404
1405    At first we increase theta and see what happens.  We start
1406    theta at a reasonable guess.  If in right area then we do bit by bit.
1407    We save possible pivot candidates
1408
1409   */
1410
1411  // do first pass to get possibles
1412  // We can also see if unbounded
1413
1414  double * work=rowArray->denseVector();
1415  int number=rowArray->getNumElements();
1416  int * which=rowArray->getIndices();
1417
1418  // we need to swap sign if coming in from ub
1419  double way = directionIn_;
1420  double maximumMovement;
1421  if (way>0.0) 
1422    maximumMovement = min(1.0e30,upperIn_-valueIn_);
1423  else
1424    maximumMovement = min(1.0e30,valueIn_-lowerIn_);
1425
1426  int iIndex;
1427
1428  // Work out coefficients for quadratic term
1429  // This is expanded one
1430  const CoinPackedMatrix * quadratic = info->quadraticObjective();
1431  const int * columnQuadratic = quadratic->getIndices();
1432  const CoinBigIndex * columnQuadraticStart = quadratic->getVectorStarts();
1433  const int * columnQuadraticLength = quadratic->getVectorLengths();
1434  const double * quadraticElement = quadratic->getElements();
1435  //const double * originalCost = info->originalObjective();
1436  // Use rhsArray
1437  rhsArray->clear();
1438  int * index2 = rhsArray->getIndices();
1439  int numberXColumns = info->numberXColumns();
1440  int number2=0;
1441  //int numberOriginalRows = info->numberXRows();
1442  // sj
1443  int iSjRow=-1;
1444  // sj for incoming
1445  int iSjRow2=-1,crucialSj2=-1;
1446  if (sequenceIn_<numberXColumns) {
1447    crucialSj2 = sequenceIn_;
1448    crucialSj2 += numberRows_; // sj2 which should go to 0.0
1449  } else if (sequenceIn_>=numberColumns_) {
1450    int iRow=sequenceIn_-numberColumns_;
1451    crucialSj2 = iRow;
1452    crucialSj2 += numberXColumns; // sj2 which should go to 0.0
1453  }
1454  double tj = 0.0;
1455  // Change in objective will be theta*coeff1 + theta*theta*coeff2
1456  double coeff1 = 0.0;
1457  double coeff2 = 0.0;
1458  double coeffSlack=0.0;
1459  for (iIndex=0;iIndex<number;iIndex++) {
1460    int iRow = which[iIndex];
1461    double alpha = -work[iRow]*way;
1462    int iPivot=pivotVariable_[iRow];
1463    if (numberIterations_==17)
1464      printf("row %d col %d alpha %g solution %g\n",iRow,iPivot,alpha,solution_[iPivot]);
1465    if (iPivot<numberXColumns) {
1466      index2[number2++]=iPivot;
1467      rhs[iPivot]=alpha;
1468      coeff1 += alpha*cost_[iPivot];
1469      //printf("col %d alpha %g solution %g cost %g scale %g\n",iPivot,alpha,solution_[iPivot],
1470      //     cost_[iPivot],columnScale_[iPivot]);
1471    } else {
1472      if (iPivot>=numberColumns_) {
1473        // ? do we go through column of pi
1474        assert (iPivot!=crucialSj);
1475        coeffSlack += alpha*cost_[iPivot];
1476      } else if (iPivot<numberRows_) {
1477        // ? do we go through column of pi
1478        if (iPivot==crucialSj) {
1479          tj = alpha;
1480          iSjRow = iRow;
1481        }
1482      } else {
1483        if (iPivot==crucialSj) {
1484          tj = alpha;
1485          iSjRow = iRow;
1486        }
1487      }
1488    }
1489    if (iPivot==crucialSj2)
1490      iSjRow2=iRow;
1491  }
1492  // Incoming
1493  if (sequenceIn_<numberXColumns) {
1494    index2[number2++]=sequenceIn_;
1495    rhs[sequenceIn_]=way;
1496    //printf("incoming col %d alpha %g solution %g cost %g scale %g\n",sequenceIn_,way,valueIn_,
1497    //   cost_[sequenceIn_],columnScale_[sequenceIn_]);
1498    coeff1 += way*cost_[sequenceIn_];
1499  } else {
1500    //printf("incoming new col %d alpha %g solution %g\n",sequenceIn_,way,valueIn_);
1501    coeffSlack += way*cost_[sequenceIn_];
1502  }
1503  printf("coeff1 now %g\n",coeff1);
1504  rhsArray->setNumElements(number2);
1505  double largestCoeff1=1.0e-20;
1506  for (iIndex=0;iIndex<number2;iIndex++) {
1507    int iColumn=index2[iIndex];
1508    //double valueI = solution_[iColumn];
1509    double alphaI = rhs[iColumn];
1510    int j;
1511    for (j=columnQuadraticStart[iColumn];
1512         j<columnQuadraticStart[iColumn]+columnQuadraticLength[iColumn];j++) {
1513      int jColumn = columnQuadratic[j];
1514      double valueJ = solution_[jColumn];
1515      double alphaJ = rhs[jColumn];
1516      double elementValue = quadraticElement[j];
1517      double addValue = (valueJ*alphaI)*elementValue;
1518      largestCoeff1 = max(largestCoeff1,fabs(addValue));
1519      coeff1 += addValue;
1520      coeff2 += (alphaI*alphaJ)*elementValue;
1521    }
1522  }
1523  coeff2 *= 0.5;
1524  if (coeffSlack)
1525    printf("slack has non-zero cost - trouble?\n");
1526  coeff1 += coeffSlack;
1527  //coeff1 = way*dualIn_;
1528  int accuracyFlag=0;
1529  if (!cleanupIteration) {
1530    if (fabs(dualIn_)<dualTolerance_&&fabs(coeff1)<dualTolerance_) {
1531      dualIn_=0.0;
1532      coeff1=0.0;
1533    }
1534    if (fabs(way*coeff1-dualIn_)>1.0e-2*(1.0+fabs(dualIn_)))
1535      printf("primal error %g, largest %g, coeff1 %g, dualin %g\n",
1536             largestPrimalError_,largestCoeff1,way*coeff1,dualIn_);
1537    assert (fabs(way*coeff1-dualIn_)<1.0e-1*(1.0+fabs(dualIn_)));
1538    assert (way*coeff1*dualIn_>=0.0);
1539    if (way*coeff1*dualIn_<0.0) {
1540      // bad
1541      accuracyFlag=2;
1542    } else if (fabs(way*coeff1-dualIn_)>1.0e-3*(1.0+fabs(dualIn_))) {
1543      // not wonderful
1544      accuracyFlag=1;
1545    }
1546    if (crucialSj>=0) {
1547      solution_[crucialSj]=way*coeff1;
1548    }
1549  } else if (dualIn_<0.0&&way*coeff1>1.0e-2) {
1550    printf("bad pair\n");
1551    accuracyFlag=1;
1552  } else if (dualIn_>0.0&&way*coeff1<-1.0e-2) {
1553    accuracyFlag=1;
1554    printf("bad pair\n");
1555  }
1556  // interesting places are where derivative zero or sj goes to zero
1557  double d1,d2=1.0e50;
1558  if (fabs(coeff2)>1.0e-9)
1559    d1 = - 0.5*coeff1/coeff2;
1560  else if (coeff1<=1.0e-6)
1561    d1 = maximumMovement;
1562  else
1563    d1 = 0.0;
1564  if (fabs(tj)<1.0e-7||!cleanupIteration) {
1565    //if (sequenceIn_<numberXColumns)
1566      //printf("column %d is basically linear\n",sequenceIn_);
1567    //assert(!columnQuadraticLength[sequenceIn_]);
1568  } else {
1569    d2 = -solution_[crucialSj]/tj;
1570    if (d2<0.0) {
1571      printf("d2 would be negative at %g\n",d2);
1572      d2=1.0e50;
1573    }
1574  }
1575  if (d1>1.0e10&&d2>1.0e10) {
1576    // linear variable entering
1577    // maybe we should have done dual iteration to force sj to 0.0
1578    //printf("linear variable\n");
1579  }
1580  handler_->message(CLP_QUADRATIC_PRIMAL_DETAILS,messages_)
1581    <<coeff1<<coeff2<<coeffSlack<<dualIn_<<d1<<d2
1582    <<CoinMessageEol;
1583  // reality check
1584  {
1585    if (crucialSj2>=0) {
1586      // see if works out
1587      if (getColumnStatus(crucialSj2)==basic) {
1588        if (iSjRow2>=0) {
1589          //corresponding alpha nonzero
1590          double alpha=work[iSjRow2];
1591          printf("Sj alpha %g sol %g ratio %g\n",alpha,solution_[crucialSj2],solution_[crucialSj2]/alpha);
1592          if (fabs(fabs(d1)-fabs(solution_[crucialSj2]/alpha))>1.0e-3*(1.0+fabs(d1))) {
1593            printf("bad test\n");
1594            if (factorization_->pivots())
1595              accuracyFlag=2;
1596          }
1597          d1=fabs(solution_[crucialSj2]/alpha);
1598        } else {
1599          printf("Sj basic but no alpha\n");
1600        }
1601      } else {
1602        printf("Sj not basic\n");
1603      }
1604    }
1605  }
1606  //if (!cleanupIteration)
1607  //dualIn_ = way*coeff1;
1608  double mind1d2=1.0e50;
1609  if (cleanupIteration) {
1610    // we are only interested in d1 if smaller than d2
1611    mind1d2 = min(maximumMovement,d2);
1612    //if (d1>1.0e-8&&d1<0.999*d2)
1613    //mind1d2=d1;
1614  } else {
1615    // There is no d2 - d1 refers to crucialSj
1616    if (d1>1.0e-5) {
1617      mind1d2 = min(maximumMovement,d1);
1618      //assert (d1>=0.0);
1619    }
1620  }
1621  maximumMovement = min(maximumMovement,mind1d2);
1622  double trueMaximumMovement;
1623  if (way>0.0) 
1624    trueMaximumMovement = min(1.0e30,upperIn_-valueIn_);
1625  else
1626    trueMaximumMovement = min(1.0e30,valueIn_-lowerIn_);
1627  rhsArray->clear();
1628  double tentativeTheta = maximumMovement;
1629  double upperTheta = maximumMovement;
1630  bool throwAway=false;
1631  if (numberIterations_==1750) {
1632    printf("Bad iteration coming up after iteration %d\n",numberIterations_);
1633  }
1634
1635  double minimumAny=1.0e50;
1636  int whichMinimum=-1;
1637  double minimumAlpha=0.0;
1638  for (iIndex=0;iIndex<number;iIndex++) {
1639
1640    int iRow = which[iIndex];
1641    double alpha = work[iRow];
1642    int iPivot=pivotVariable_[iRow];
1643    alpha *= way;
1644    double oldValue = solution(iPivot);
1645    // get where in bound sequence
1646    if (alpha>0.0) {
1647      // basic variable going towards lower bound
1648      double bound = lower(iPivot);
1649      oldValue -= bound;
1650    } else if (alpha<0.0) {
1651      // basic variable going towards upper bound
1652      double bound = upper(iPivot);
1653      oldValue = bound-oldValue;
1654    }
1655    if (iPivot>=numberXColumns&&iPivot<numberColumns_) {
1656      double bound=0.0;
1657      double oldValue2 = solution(iPivot);
1658      if (alpha>0.0) {
1659        // basic variable going towards lower bound
1660        oldValue2 -= bound;
1661      } else if (alpha<0.0) {
1662        // basic variable going towards upper bound
1663        oldValue2 = bound-oldValue;
1664      }
1665      double value = oldValue2-minimumAny*fabs(alpha);
1666      if (value*oldValue2<0.0) {
1667        double ratio = fabs(oldValue2/alpha);
1668        if (ratio<minimumAny&&fabs(alpha)>1.0e2*acceptablePivot) {
1669          minimumAny = ratio;
1670          whichMinimum = iRow;
1671          minimumAlpha= fabs(alpha);
1672        }
1673      }
1674    }
1675    double value = oldValue-tentativeTheta*fabs(alpha);
1676    assert (oldValue>=-primalTolerance_*1.0001);
1677    if (value<-primalTolerance_) {
1678      // add to list
1679      spare[numberRemaining]=alpha;
1680      rhs[iRow]=oldValue;
1681      index[numberRemaining++]=iRow;
1682      double value=oldValue-upperTheta*fabs(alpha);
1683      if (value<-primalTolerance_&&fabs(alpha)>=acceptablePivot)
1684        upperTheta = (oldValue+primalTolerance_)/fabs(alpha);
1685    }
1686  }
1687
1688  // we need to keep where rhs non-zeros are
1689  int numberInRhs=numberRemaining;
1690  memcpy(rhsArray->getIndices(),index,numberInRhs*sizeof(int));
1691  rhsArray->setNumElements(numberInRhs);
1692
1693  theta_=maximumMovement;
1694
1695  bool goBackOne = false;
1696  double fake=1.0e-2;
1697
1698  if (numberRemaining) {
1699
1700    // looks like pivoting
1701    // now try until reasonable theta
1702    tentativeTheta = max(10.0*upperTheta,1.0e-7);
1703    tentativeTheta = min(tentativeTheta,maximumMovement);
1704   
1705    // loops increasing tentative theta until can't go through
1706   
1707    while (tentativeTheta <= maximumMovement) {
1708      double bestPivotBeforeInteresting=0.0;
1709      double thruThis = 0.0;
1710     
1711      double bestPivot=acceptablePivot;
1712      pivotRow_ = -1;
1713     
1714      numberSwapped = 0;
1715     
1716      upperTheta = maximumMovement;
1717     
1718      for (iIndex=0;iIndex<numberRemaining;iIndex++) {
1719
1720        int iRow = index[iIndex];
1721        double alpha = spare[iIndex];
1722        double oldValue = rhs[iRow];
1723        double value = oldValue-tentativeTheta*fabs(alpha);
1724
1725        if (value<-primalTolerance_) {
1726          // how much would it cost to go thru
1727          thruThis += alpha*
1728            nonLinearCost_->changeInCost(pivotVariable_[iRow],alpha);
1729          // goes on swapped list (also means candidates if too many)
1730          indexSwapped[numberSwapped++]=iRow;
1731          if (fabs(alpha)>bestPivot) {
1732            bestPivot=fabs(alpha);
1733            pivotRow_ = iRow;
1734            theta_ = oldValue/bestPivot;
1735          }
1736        } else {
1737          value = oldValue-upperTheta*fabs(alpha);
1738          if (value<-primalTolerance_ && fabs(alpha)>=acceptablePivot) 
1739            upperTheta = (oldValue+primalTolerance_)/fabs(alpha);
1740        }
1741      }
1742     
1743      maximumSwapped = max(maximumSwapped,numberSwapped);
1744      bestPivotBeforeInteresting=bestPivot;
1745
1746      //double dualCheck = - 2.0*coeff2*tentativeTheta - coeff1 - 0.1*infeasibilityCost_;
1747      double dualCheck = - 2.0*coeff2*tentativeTheta - coeff1;
1748      // but make a bit more pessimistic if pivot
1749      //if (bestPivot>acceptablePivot)
1750      //dualCheck=max(dualCheck-100.0*dualTolerance_,0.99*dualCheck);
1751      //dualCheck=0.0;
1752      if ((totalThru+thruThis>=dualCheck&&bestPivot>acceptablePivot)||fake*bestPivot>1.0e-3) {
1753        // We should be pivoting in this batch
1754        // so compress down to this lot
1755
1756        int saveNumber = numberRemaining;
1757        numberRemaining=0;
1758        for (iIndex=0;iIndex<numberSwapped;iIndex++) {
1759          int iRow = indexSwapped[iIndex];
1760          spare[numberRemaining]=way*work[iRow];
1761          index[numberRemaining++]=iRow;
1762        }
1763        memset(spare+numberRemaining,0,
1764               (saveNumber-numberRemaining)*sizeof(double));
1765        double lastThru = totalThru;
1766        int iTry;
1767#define MAXTRY 100
1768        // first get ratio with tolerance
1769        for (iTry=0;iTry<MAXTRY;iTry++) {
1770         
1771          upperTheta=maximumMovement;
1772          numberSwapped = 0;
1773         
1774          for (iIndex=0;iIndex<numberRemaining;iIndex++) {
1775           
1776            int iRow = index[iIndex];
1777            double alpha = fabs(spare[iIndex]);
1778            double oldValue = rhs[iRow];
1779            double value = oldValue-upperTheta*alpha;
1780           
1781            if (value<-primalTolerance_ && alpha>=acceptablePivot) 
1782              upperTheta = (oldValue+primalTolerance_)/alpha;
1783           
1784          }
1785         
1786          // now look at best in this lot
1787          bestPivot=acceptablePivot;
1788          pivotRow_=-1;
1789          for (iIndex=0;iIndex<numberRemaining;iIndex++) {
1790           
1791            int iRow = index[iIndex];
1792            double alpha = spare[iIndex];
1793            double oldValue = rhs[iRow];
1794            double value = oldValue-upperTheta*fabs(alpha);
1795           
1796            if (value<=0) {
1797              // how much would it cost to go thru
1798              totalThru += alpha*
1799                nonLinearCost_->changeInCost(pivotVariable_[iRow],alpha);
1800              // goes on swapped list (also means candidates if too many)
1801              indexSwapped[numberSwapped++]=iRow;
1802              if (fabs(alpha)>bestPivot) {
1803                bestPivot=fabs(alpha);
1804                theta_ = oldValue/bestPivot;
1805                pivotRow_=iRow;
1806              }
1807            } else {
1808              value = oldValue-upperTheta*fabs(alpha);
1809              if (value<-primalTolerance_ && fabs(alpha)>=acceptablePivot) 
1810                upperTheta = (oldValue+primalTolerance_)/fabs(alpha);
1811            }
1812          }
1813         
1814          maximumSwapped = max(maximumSwapped,numberSwapped);
1815          if (bestPivot<0.1*bestEverPivot&&
1816              bestEverPivot>1.0e-6&&bestPivot<1.0e-3) {
1817            // back to previous one
1818            goBackOne = true;
1819            break;
1820          } else if (pivotRow_==-1&&upperTheta>largeValue_) {
1821            if (lastPivot>acceptablePivot) {
1822              // back to previous one
1823              goBackOne = true;
1824              //break;
1825            } else {
1826              // can only get here if all pivots so far too small
1827            }
1828            break;
1829          } else {
1830            dualCheck = - 2.0*coeff2*theta_ - coeff1;
1831            if (bestPivotBeforeInteresting>1.0e-4&&bestPivot<1.0e-6)
1832              dualCheck=1.0e7;
1833            if ((totalThru>=dualCheck||fake*bestPivot>1.0e-3)
1834                &&(sequenceIn_<numberXColumns||sequenceIn_>=numberColumns_)) {
1835              if (!cleanupIteration) {
1836                // We can pivot out sj
1837                if (upperTheta==maximumMovement) {
1838                  printf("figures\n");
1839                } else if (iSjRow2>=0) {
1840                  if (lastThru>=dualCheck&&theta_<trueMaximumMovement) {
1841                    printf("totalThru %g, lastThru %g - taking sj to zero?\n",totalThru,lastThru);
1842                    // make sure will take correct path
1843                    mind1d2=1.0;
1844                    maximumMovement=1.0;
1845                    d2=0.0;
1846                    pivotRow_=-1;
1847                  }
1848                }
1849              } else {
1850                if (pivotRow_<0&&lastPivotRow<0)
1851                  throwAway=true;
1852              }
1853              break; // no point trying another loop
1854            } else if (totalThru>=dualCheck||fake*bestPivot>1.0e-3||upperTheta==maximumMovement) {
1855              break; // no point trying another loop
1856            } else {
1857              // skip this lot
1858              nonLinearCost_->goThru(numberSwapped,way,indexSwapped, work,rhs);
1859              lastPivotRow=pivotRow_;
1860              lastTheta = theta_;
1861              lastThru = numberThru;
1862              numberThru += numberSwapped;
1863              lastNumberSwapped = numberSwapped;
1864              memcpy(saveSwapped,indexSwapped,lastNumberSwapped*sizeof(int));
1865              if (bestPivot>bestEverPivot)
1866                bestEverPivot=bestPivot;
1867              lastThru = totalThru;
1868            }
1869          }
1870        }
1871        break;
1872      } else {
1873        // skip this lot
1874        nonLinearCost_->goThru(numberSwapped,way,indexSwapped, work,rhs);
1875        lastPivotRow=pivotRow_;
1876        lastTheta = theta_;
1877        lastThru = numberThru;
1878        numberThru += numberSwapped;
1879        lastNumberSwapped = numberSwapped;
1880        memcpy(saveSwapped,indexSwapped,lastNumberSwapped*sizeof(int));
1881        if (bestPivot>bestEverPivot)
1882          bestEverPivot=bestPivot;
1883        totalThru += thruThis;
1884        tentativeTheta = 2.0*upperTheta;
1885      }
1886    }
1887    // can get here without pivotRow_ set but with lastPivotRow
1888    if (goBackOne||(pivotRow_<0&&lastPivotRow>=0)) {
1889      // back to previous one
1890      pivotRow_=lastPivotRow;
1891      theta_ = lastTheta;
1892            // undo this lot
1893      nonLinearCost_->goBack(lastNumberSwapped,saveSwapped,rhs);
1894      memcpy(indexSwapped,saveSwapped,lastNumberSwapped*sizeof(int));
1895      numberSwapped = lastNumberSwapped;
1896    }
1897  }
1898  if (0&&minimumAny<theta_&&cleanupIteration&&bestEverPivot<1.0e-2
1899      &&minimumAlpha>bestEverPivot*1.0e2&&minimumAny>1.0e-1&&info->currentPhase()==0) {
1900    printf("Possible other pivot %d %g %g - alpha %g\n",
1901           whichMinimum,minimumAny,theta_,minimumAlpha);
1902    pivotRow_ = whichMinimum;
1903    alpha_ = work[pivotRow_];
1904    // translate to sequence
1905    sequenceOut_ = pivotVariable_[pivotRow_];
1906    valueOut_ = solution(sequenceOut_);
1907    lowerOut_=0.0;
1908    upperOut_=0.0;
1909    theta_= fabs(valueOut_/alpha_);
1910
1911    if (way<0.0) 
1912      theta_ = - theta_;
1913    if (alpha_*way<0.0) {
1914      directionOut_=-1;      // to upper bound
1915    } else {
1916      directionOut_=1;      // to lower bound
1917    }
1918    dualOut_ = reducedCost(sequenceOut_);
1919  } else {
1920
1921  if (pivotRow_>=0) {
1922    if (0) {
1923      double move = coeff1*theta_+coeff2*theta_*theta_;
1924      printf("Predicted change in obj is %g %g %g\n",
1925             move,objectiveValue_-move,objectiveValue_+move);
1926    }
1927#define MINIMUMTHETA 1.0e-12
1928    // will we need to increase tolerance
1929#ifdef CLP_DEBUG
1930    bool found=false;
1931#endif
1932    double largestInfeasibility = primalTolerance_;
1933    if (theta_<MINIMUMTHETA) {
1934      theta_=MINIMUMTHETA;
1935      for (iIndex=0;iIndex<numberSwapped;iIndex++) {
1936        int iRow = indexSwapped[iIndex];
1937#ifdef CLP_DEBUG
1938        if (iRow==pivotRow_)
1939          found=true;
1940#endif
1941        largestInfeasibility = max (largestInfeasibility,
1942                                    -(rhs[iRow]-fabs(work[iRow])*theta_));
1943      }
1944#ifdef CLP_DEBUG
1945      assert(found);
1946      if (largestInfeasibility>primalTolerance_&&(handler_->logLevel()&32))
1947        printf("Primal tolerance increased from %g to %g\n",
1948               primalTolerance_,largestInfeasibility);
1949#endif
1950      primalTolerance_ = max(primalTolerance_,largestInfeasibility);
1951    }
1952    alpha_ = work[pivotRow_];
1953    // translate to sequence
1954    sequenceOut_ = pivotVariable_[pivotRow_];
1955    valueOut_ = solution(sequenceOut_);
1956    lowerOut_=lower_[sequenceOut_];
1957    upperOut_=upper_[sequenceOut_];
1958
1959    if (way<0.0) 
1960      theta_ = - theta_;
1961    double newValue = valueOut_ - theta_*alpha_;
1962    if (alpha_*way<0.0) {
1963      directionOut_=-1;      // to upper bound
1964      if (fabs(theta_)>1.0e-6)
1965        upperOut_ = nonLinearCost_->nearest(sequenceOut_,newValue);
1966      else
1967        upperOut_ = newValue;
1968    } else {
1969      directionOut_=1;      // to lower bound
1970      if (fabs(theta_)>1.0e-6)
1971        lowerOut_ = nonLinearCost_->nearest(sequenceOut_,newValue);
1972      else
1973        lowerOut_ = newValue;
1974    }
1975    dualOut_ = reducedCost(sequenceOut_);
1976  } else {
1977    // If no pivot but bad move then throw away
1978    if (throwAway&&(sequenceIn_<numberXColumns||sequenceIn_>=numberColumns_)) {
1979      assert (pivotRow_<0);
1980      accuracyFlag=2;
1981    }
1982    if (maximumMovement<1.0e20&&maximumMovement==trueMaximumMovement) {
1983      // flip
1984      theta_ = maximumMovement;
1985      pivotRow_ = -2; // so we can tell its a flip
1986      result=0;
1987      sequenceOut_ = sequenceIn_;
1988      valueOut_ = valueIn_;
1989      dualOut_ = dualIn_;
1990      lowerOut_ = lowerIn_;
1991      upperOut_ = upperIn_;
1992      alpha_ = 0.0;
1993      if (accuracyFlag<2)
1994        accuracyFlag=0;
1995      if (way<0.0) {
1996        directionOut_=1;      // to lower bound
1997        theta_ = lowerOut_ - valueOut_;
1998      } else {
1999        directionOut_=-1;      // to upper bound
2000        theta_ = upperOut_ - valueOut_;
2001      }
2002      // we may still have sj to get rid of
2003    } else if (fabs(maximumMovement-mind1d2)<dualTolerance_) {
2004      // crucialSj local copy i.e. dj going to zero
2005      if (maximumMovement==d2) {
2006        // sj going to zero
2007        result=0;
2008      } else {
2009        // incoming dj going to zero
2010        if (!cleanupIteration) {
2011          result=0;
2012          iSjRow=iSjRow2;
2013          if (iSjRow>=0)
2014            crucialSj = pivotVariable_[iSjRow];
2015          else
2016            crucialSj=-1;
2017          assert (pivotRow_<0);
2018          // If crucialSj <0 then make superbasic?
2019          if (crucialSj<0) {
2020            printf("**ouch\n");
2021            theta_ = maximumMovement;
2022            pivotRow_ = -2; // so we can tell its a flip
2023            result=0;
2024            sequenceOut_ = sequenceIn_;
2025            valueOut_ = valueIn_;
2026            dualOut_ = dualIn_;
2027            lowerOut_ = lowerIn_;
2028            upperOut_ = upperIn_;
2029            alpha_ = 0.0;
2030            if (accuracyFlag<2)
2031              accuracyFlag=0;
2032            if (way<0.0) {
2033              directionOut_=1;      // to lower bound
2034              lowerIn_ = valueOut_-theta_;
2035              theta_ = lowerIn_ - valueOut_;
2036            } else {
2037              directionOut_=-1;      // to upper bound
2038              upperIn_ = valueOut_+theta_;
2039              theta_ = upperIn_ - valueOut_;
2040            }
2041          }
2042        } else {
2043          assert (fabs(dualIn_)<1.0e-3);
2044          result=1;
2045          if (minimumAlpha>0.0) {
2046            // could do this in other places if pivot looks good
2047            printf("Should take minimum\n");
2048            pivotRow_ = whichMinimum;
2049            alpha_ = work[pivotRow_];
2050            // translate to sequence
2051            sequenceOut_ = pivotVariable_[pivotRow_];
2052            valueOut_ = solution(sequenceOut_);
2053            theta_ = minimumAny;
2054            if (way<0.0) 
2055              theta_ = - theta_;
2056            lowerOut_=0.0;
2057            upperOut_=0.0;
2058            //????
2059            dualOut_ = reducedCost(sequenceOut_);
2060          }
2061        }
2062      }
2063      if (!result&&crucialSj>=0) {
2064        setColumnStatus(crucialSj,isFree);
2065        pivotRow_ = iSjRow;
2066        alpha_ = work[pivotRow_];
2067        // translate to sequence
2068        sequenceOut_ = pivotVariable_[pivotRow_];
2069        assert (sequenceOut_==crucialSj);
2070        valueOut_ = solution(sequenceOut_);
2071        theta_ = fabs(valueOut_/alpha_);
2072        if (fabs(maximumMovement-theta_)>1.0e-3*(1.0+maximumMovement))
2073          printf("maximumMovement %g mismatch with theta %g\n",maximumMovement,theta_);;
2074        if (way<0.0) 
2075          theta_ = - theta_;
2076        lowerOut_=0.0;
2077        upperOut_=0.0;
2078        //????
2079        dualOut_ = reducedCost(sequenceOut_);
2080      }
2081    } else {
2082      // need to do something
2083      accuracyFlag=2;
2084    }
2085  }
2086  }
2087
2088
2089  // clear arrays
2090
2091  memset(spare,0,numberRemaining*sizeof(double));
2092  memset(saveSwapped,0,maximumSwapped*sizeof(int));
2093
2094  // put back original bounds etc
2095  nonLinearCost_->goBackAll(rhsArray);
2096
2097  rhsArray->clear();
2098  // Change in objective will be theta*coeff1 + theta*theta*coeff2
2099  objectiveValue_ += theta_*coeff1+theta_*theta_*coeff2;
2100  {
2101    int iColumn;
2102    objectiveValue_ =0.0;
2103    CoinPackedMatrix * quadratic = 
2104      info->originalObjective()->quadraticObjective();
2105    const int * columnQuadratic = quadratic->getIndices();
2106    const CoinBigIndex * columnQuadraticStart = quadratic->getVectorStarts();
2107    const int * columnQuadraticLength = quadratic->getVectorLengths();
2108    const double * quadraticElement = quadratic->getElements();
2109    int numberColumns = info->numberXColumns();
2110    const double * objective = info->linearObjective();
2111    for (iColumn=0;iColumn<numberColumns;iColumn++) 
2112      objectiveValue_ += objective[iColumn]*solution_[iColumn];
2113    for (iColumn=0;iColumn<numberColumns;iColumn++) {
2114      double valueI = solution_[iColumn];
2115      int j;
2116      for (j=columnQuadraticStart[iColumn];
2117           j<columnQuadraticStart[iColumn]+columnQuadraticLength[iColumn];j++) {
2118        int jColumn = columnQuadratic[j];
2119        double valueJ = solution_[jColumn];
2120        double elementValue = quadraticElement[j];
2121        objectiveValue_ += 0.5*valueI*valueJ*elementValue;
2122      }
2123    }
2124  }
2125  if (accuracyFlag<2) {
2126    return result+10*accuracyFlag;
2127  } else {
2128    pivotRow_=1; // make sure correct path
2129    return 20;
2130  }
2131}
2132/* Checks if finished.  Updates status */
2133void 
2134ClpSimplexPrimalQuadratic::statusOfProblemInPrimal(int & lastCleaned,int type,
2135                                          ClpSimplexProgress * progress,
2136                                                   ClpQuadraticInfo * info)
2137{
2138  if (info->currentPhase())
2139    info->setCurrentSolution(solution_);
2140  //info->saveStatus();
2141 
2142  if (type==2) {
2143    // trouble - restore solution
2144    memcpy(status_ ,saveStatus_,(numberColumns_+numberRows_)*sizeof(char));
2145    memcpy(rowActivityWork_,savedSolution_+numberColumns_ ,
2146           numberRows_*sizeof(double));
2147    memcpy(columnActivityWork_,savedSolution_ ,
2148           numberColumns_*sizeof(double));
2149    forceFactorization_=1; // a bit drastic but ..
2150    pivotRow_=-1; // say no weights update
2151    changeMade_++; // say change made
2152    info->restoreStatus();
2153  }
2154  int saveThreshold = factorization_->sparseThreshold();
2155  int tentativeStatus = problemStatus_;
2156  if (problemStatus_>-3||problemStatus_==-4) {
2157    // factorize
2158    // later on we will need to recover from singularities
2159    // also we could skip if first time
2160    // do weights
2161    // This may save pivotRow_ for use
2162    primalColumnPivot_->saveWeights(this,1);
2163
2164    if (type) {
2165      // is factorization okay?
2166      int numberPivots = factorization_->pivots();
2167      if (internalFactorize(1)) {
2168        if (solveType_==2) {
2169          // say odd
2170          problemStatus_=5;
2171          return;
2172        }
2173        numberIterations_ = progress_->lastIterationNumber(0);
2174        // no - restore previous basis
2175        assert (info->crucialSj()<0); // need to work out how to unwind
2176        assert (type==1);
2177        memcpy(status_ ,saveStatus_,(numberColumns_+numberRows_)*sizeof(char));
2178        memcpy(rowActivityWork_,savedSolution_+numberColumns_ ,
2179               numberRows_*sizeof(double));
2180        memcpy(columnActivityWork_,savedSolution_ ,
2181               numberColumns_*sizeof(double));
2182        forceFactorization_=1; // a bit drastic but ..
2183        type = 2;
2184        assert (internalFactorize(1)==0);
2185        info->restoreStatus();
2186        // flag incoming
2187        if (numberPivots==1) {
2188          int crucialSj = info->crucialSj();
2189          if (crucialSj<0) {
2190            setFlagged(sequenceIn_);
2191          } else {
2192            int iSequence;
2193            int numberXColumns = info->numberXColumns();
2194            if (crucialSj<numberRows_) {
2195              // row
2196              iSequence = crucialSj-numberXColumns+numberColumns_;
2197            } else {
2198              // column
2199              iSequence = crucialSj-numberRows_;
2200            }
2201            if (getColumnStatus(iSequence)!=basic) {
2202              setFlagged(iSequence);
2203              info->setSequenceIn(-1);
2204            } else {
2205              printf("**** %d is basic ?\n",iSequence);
2206            }
2207          }
2208        }
2209        changeMade_++; // say change made
2210      }
2211    }
2212    if (problemStatus_!=-4)
2213      problemStatus_=-3;
2214  }
2215  if (info->crucialSj()<0) {
2216    // at this stage status is -3 or -5 if looks unbounded
2217    // get primal and dual solutions
2218    // put back original costs and then check
2219    createRim(4);
2220    if (info->currentPhase()) {
2221      memcpy(solution_,info->currentSolution(),
2222             (numberRows_+numberColumns_)*sizeof(double));
2223    }
2224    gutsOfSolution(NULL,NULL);
2225    if (info->currentPhase()) {
2226      memcpy(solution_,info->currentSolution(),
2227             (numberRows_+numberColumns_)*sizeof(double));
2228      nonLinearCost_->checkInfeasibilities(primalTolerance_);
2229    }
2230    checkComplementarity(info,rowArray_[2],rowArray_[3]);
2231    // Double check reduced costs if no action
2232    if (progress->lastIterationNumber(0)==numberIterations_) {
2233      if (primalColumnPivot_->looksOptimal()) {
2234        numberDualInfeasibilities_ = 0;
2235        sumDualInfeasibilities_ = 0.0;
2236      }
2237    }
2238    // Check if looping
2239    int loop;
2240    if (type!=2) 
2241      loop = progress->looping(); // saves iteration number
2242    else
2243      loop=-1;
2244    //loop=-1;
2245    if (loop>=0) {
2246      problemStatus_ = loop; //exit if in loop
2247      return ;
2248    } else if (loop<-1) {
2249      // something may have changed
2250      gutsOfSolution(NULL,NULL);
2251      checkComplementarity(info,rowArray_[2],rowArray_[3]);
2252    }
2253    // really for free variables in
2254    //if((progressFlag_&2)!=0)
2255    //problemStatus_=-1;;
2256    progressFlag_ = 0; //reset progress flag
2257   
2258    handler_->message(CLP_SIMPLEX_STATUS,messages_)
2259      <<numberIterations_<<objectiveValue();
2260    handler_->printing(sumPrimalInfeasibilities_>0.0)
2261      <<sumPrimalInfeasibilities_<<numberPrimalInfeasibilities_;
2262    handler_->printing(sumDualInfeasibilities_>0.0)
2263      <<sumDualInfeasibilities_<<numberDualInfeasibilities_;
2264    handler_->printing(numberDualInfeasibilitiesWithoutFree_
2265                       <numberDualInfeasibilities_)
2266                         <<numberDualInfeasibilitiesWithoutFree_;
2267    handler_->message()<<CoinMessageEol;
2268    assert (primalFeasible());
2269    // we may wish to say it is optimal even if infeasible
2270    bool alwaysOptimal = (specialOptions_&1)!=0;
2271    // give code benefit of doubt
2272    if (sumOfRelaxedDualInfeasibilities_ == 0.0&&
2273        sumOfRelaxedPrimalInfeasibilities_ == 0.0&&info->sequenceIn()<0) {
2274      // say optimal (with these bounds etc)
2275      numberDualInfeasibilities_ = 0;
2276      sumDualInfeasibilities_ = 0.0;
2277      numberPrimalInfeasibilities_ = 0;
2278      sumPrimalInfeasibilities_ = 0.0;
2279    }
2280    // had ||(type==3&&problemStatus_!=-5) -- ??? why ????
2281    if (dualFeasible()||problemStatus_==-4) {
2282      if (nonLinearCost_->numberInfeasibilities()&&!alwaysOptimal) {
2283        //may need infeasiblity cost changed
2284        // we can see if we can construct a ray
2285        // make up a new objective
2286        double saveWeight = infeasibilityCost_;
2287        // save nonlinear cost as we are going to switch off costs
2288        ClpNonLinearCost * nonLinear = nonLinearCost_;
2289        // do twice to make sure Primal solution has settled
2290        // put non-basics to bounds in case tolerance moved
2291        // put back original costs
2292        createRim(4);
2293        // put all non-basic variables to bounds
2294        int iSequence;
2295        for (iSequence=0;iSequence<numberRows_+numberColumns_;iSequence++) {
2296          Status status = getStatus(iSequence);
2297          if (status==atLowerBound||status==isFixed)
2298            solution_[iSequence]=lower_[iSequence];
2299          else if (status==atUpperBound)
2300            solution_[iSequence]=upper_[iSequence];
2301        }
2302        nonLinearCost_->checkInfeasibilities(primalTolerance_);
2303        gutsOfSolution(NULL,NULL);
2304        nonLinearCost_->checkInfeasibilities(primalTolerance_);
2305        checkComplementarity(info,rowArray_[2],rowArray_[3]);
2306       
2307        infeasibilityCost_=1.0e100;
2308        // put back original costs
2309        createRim(4);
2310        nonLinearCost_->checkInfeasibilities(primalTolerance_);
2311        // may have fixed infeasibilities - double check
2312        if (nonLinearCost_->numberInfeasibilities()==0) {
2313          // carry on
2314          problemStatus_ = -1;
2315          infeasibilityCost_=saveWeight;
2316          nonLinearCost_->checkInfeasibilities(primalTolerance_);
2317          checkComplementarity(info,rowArray_[2],rowArray_[3]);
2318        } else {
2319          nonLinearCost_=NULL;
2320          // scale
2321          int i;
2322          for (i=0;i<numberRows_+numberColumns_;i++) 
2323            cost_[i] *= 1.0e-100;
2324          gutsOfSolution(NULL,NULL);
2325          nonLinearCost_=nonLinear;
2326          infeasibilityCost_=saveWeight;
2327          if ((infeasibilityCost_>=1.0e18||
2328               numberDualInfeasibilities_==0)&&perturbation_==101) {
2329            unPerturb(); // stop any further perturbation
2330            numberDualInfeasibilities_=1; // carry on
2331          }
2332          if (infeasibilityCost_>=1.0e20||
2333              numberDualInfeasibilities_==0) {
2334            // we are infeasible - use as ray
2335            delete [] ray_;
2336            ray_ = new double [numberRows_];
2337            memcpy(ray_,dual_,numberRows_*sizeof(double));
2338            // and get feasible duals
2339            infeasibilityCost_=0.0;
2340            createRim(4);
2341            // put all non-basic variables to bounds
2342            int iSequence;
2343            for (iSequence=0;iSequence<numberRows_+numberColumns_;iSequence++) {
2344              Status status = getStatus(iSequence);
2345              if (status==atLowerBound||status==isFixed)
2346                solution_[iSequence]=lower_[iSequence];
2347              else if (status==atUpperBound)
2348                solution_[iSequence]=upper_[iSequence];
2349            }
2350            nonLinearCost_->checkInfeasibilities(primalTolerance_);
2351            gutsOfSolution(NULL,NULL);
2352            checkComplementarity(info,rowArray_[2],rowArray_[3]);
2353            // so will exit
2354            infeasibilityCost_=1.0e30;
2355            // reset infeasibilities
2356            sumPrimalInfeasibilities_=nonLinearCost_->sumInfeasibilities();;
2357            numberPrimalInfeasibilities_=
2358              nonLinearCost_->numberInfeasibilities();
2359          }
2360          if (infeasibilityCost_<1.0e20) {
2361            infeasibilityCost_ *= 5.0;
2362            changeMade_++; // say change made
2363            handler_->message(CLP_PRIMAL_WEIGHT,messages_)
2364              <<infeasibilityCost_
2365              <<CoinMessageEol;
2366            // put back original costs and then check
2367            createRim(4);
2368            nonLinearCost_->checkInfeasibilities(0.0);
2369            gutsOfSolution(NULL,NULL);
2370            checkComplementarity(info,rowArray_[2],rowArray_[3]);
2371            problemStatus_=-1; //continue
2372          } else {
2373            // say infeasible
2374            problemStatus_ = 1;
2375          }
2376        }
2377      } else {
2378        // may be optimal
2379        if (perturbation_==101) {
2380          unPerturb(); // stop any further perturbation
2381          lastCleaned=-1; // carry on
2382        }
2383        if ( lastCleaned!=numberIterations_||unflag()) {
2384          handler_->message(CLP_PRIMAL_OPTIMAL,messages_)
2385            <<primalTolerance_
2386            <<CoinMessageEol;
2387          if (numberTimesOptimal_<4) {
2388            numberTimesOptimal_++;
2389            changeMade_++; // say change made
2390            if (numberTimesOptimal_==1) {
2391              // better to have small tolerance even if slower
2392              factorization_->zeroTolerance(1.0e-15);
2393            }
2394            lastCleaned=numberIterations_;
2395            if (primalTolerance_!=dblParam_[ClpPrimalTolerance])
2396              handler_->message(CLP_PRIMAL_ORIGINAL,messages_)
2397                <<CoinMessageEol;
2398            double oldTolerance = primalTolerance_;
2399            primalTolerance_=dblParam_[ClpPrimalTolerance];
2400            // put back original costs and then check
2401            createRim(4);
2402            // put all non-basic variables to bounds
2403            int iSequence;
2404            for (iSequence=0;iSequence<numberRows_+numberColumns_;iSequence++) {
2405              Status status = getStatus(iSequence);
2406              if (status==atLowerBound||status==isFixed)
2407                solution_[iSequence]=lower_[iSequence];
2408              else if (status==atUpperBound)
2409                solution_[iSequence]=upper_[iSequence];
2410            }
2411            nonLinearCost_->checkInfeasibilities(oldTolerance);
2412            gutsOfSolution(NULL,NULL);
2413            checkComplementarity(info,rowArray_[2],rowArray_[3]);
2414            problemStatus_ = -1;
2415          } else {
2416            problemStatus_=0; // optimal
2417            if (lastCleaned<numberIterations_) {
2418              handler_->message(CLP_SIMPLEX_GIVINGUP,messages_)
2419                <<CoinMessageEol;
2420            }
2421          }
2422        } else {
2423          problemStatus_=0; // optimal
2424        }
2425      }
2426    } else {
2427      // see if looks unbounded
2428      if (problemStatus_==-5) {
2429        if (nonLinearCost_->numberInfeasibilities()) {
2430          if (infeasibilityCost_>1.0e18&&perturbation_==101) {
2431            // back off weight
2432            infeasibilityCost_ = 1.0e13;
2433            unPerturb(); // stop any further perturbation
2434          }
2435          //we need infeasiblity cost changed
2436          if (infeasibilityCost_<1.0e20) {
2437            infeasibilityCost_ *= 5.0;
2438            changeMade_++; // say change made
2439            handler_->message(CLP_PRIMAL_WEIGHT,messages_)
2440              <<infeasibilityCost_
2441              <<CoinMessageEol;
2442            // put back original costs and then check
2443            createRim(4);
2444            gutsOfSolution(NULL,NULL);
2445            checkComplementarity(info,rowArray_[2],rowArray_[3]);
2446            problemStatus_=-1; //continue
2447          } else {
2448            // say unbounded
2449            problemStatus_ = 2;
2450          }
2451        } else {
2452          // say unbounded
2453          problemStatus_ = 2;
2454        } 
2455      } else {
2456        if(type==3&&problemStatus_!=-5)
2457          unflag(); // odd
2458        // carry on
2459        problemStatus_ = -1;
2460      }
2461    }
2462  } else {
2463    // don't update solution
2464    problemStatus_=-1;
2465  }
2466  if (type==0||type==1) {
2467    if (type!=1||!saveStatus_) {
2468      // create save arrays
2469      delete [] saveStatus_;
2470      delete [] savedSolution_;
2471      saveStatus_ = new unsigned char [numberRows_+numberColumns_];
2472      savedSolution_ = new double [numberRows_+numberColumns_];
2473    }
2474    // save arrays
2475    memcpy(saveStatus_,status_,(numberColumns_+numberRows_)*sizeof(char));
2476    memcpy(savedSolution_+numberColumns_ ,rowActivityWork_,
2477           numberRows_*sizeof(double));
2478    memcpy(savedSolution_ ,columnActivityWork_,numberColumns_*sizeof(double));
2479    info->saveStatus();
2480  }
2481  // restore weights (if saved) - also recompute infeasibility list
2482  if (tentativeStatus>-3) 
2483    primalColumnPivot_->saveWeights(this,(type <2) ? 2 : 4);
2484  else
2485    primalColumnPivot_->saveWeights(this,3);
2486  if (problemStatus_<0&&!changeMade_) {
2487    problemStatus_=4; // unknown
2488  }
2489  if (saveThreshold) {
2490    // use default at present
2491    factorization_->sparseThreshold(0);
2492    factorization_->goSparse();
2493  }
2494  lastGoodIteration_ = numberIterations_;
2495}
2496// Fills in reduced costs
2497void 
2498ClpSimplexPrimalQuadratic::createDjs (ClpQuadraticInfo * info,
2499                            CoinIndexedVector * array1, CoinIndexedVector * array2)
2500{
2501  int numberQuadraticRows = info->numberQuadraticRows();
2502  int numberQuadraticColumns = info->numberQuadraticColumns();
2503  int numberXColumns = info->numberXColumns();
2504  int numberXRows = info->numberXRows();
2505  // Actually only need to do this after invert (use extra parameter to createDjs)
2506  // or when infeasibilities change
2507  // recode to do this later and update rather than recompute
2508  // When it is "change" then we don't need b (original rhs)
2509  {
2510    // update costs
2511    double * storedCost = lower_+numberColumns_+numberXRows;
2512    double * storedUpper = upper_ + numberColumns_+numberXRows;
2513    double * storedValue = solution_ + numberColumns_+numberXRows;
2514    int * index = array1->getIndices();
2515    double * modifiedCost = array1->denseVector();
2516    // Compute duals
2517    info->createGradient(this);
2518    double * gradient = info->gradient();
2519    // fill in as if linear
2520    // will not be correct unless complementary - but that should not matter
2521    // Just save sj value in that case
2522    //double saveValue=0.0;
2523    //int crucialSj = info->crucialSj();
2524    //if (crucialSj>=0)
2525    int number=0;
2526    int iRow;
2527    for (iRow=0;iRow<numberRows_;iRow++) {
2528      int iPivot=pivotVariable_[iRow];
2529      if (gradient[iPivot]) {
2530        modifiedCost[iRow] = gradient[iPivot];
2531        index[number++]=iRow;
2532      }
2533    }
2534    array1->setNumElements(number);
2535    factorization_->updateColumnTranspose(array2,array1);
2536    memcpy(dual_,modifiedCost,numberRows_*sizeof(double));
2537    array1->clear();
2538    memcpy(dj_,gradient,(numberColumns_+numberRows_)*sizeof(double));
2539    matrix_->transposeTimes(-1.0,dual_,dj_,rowScale_,columnScale_);
2540    memset(dj_+numberXColumns,0,(numberXRows+numberQuadraticColumns)*sizeof(double));
2541    // And djs
2542    for (iRow=0;iRow<numberRows_;iRow++) {
2543      dj_[numberColumns_+iRow]= cost_[numberColumns_+iRow]+dual_[iRow];
2544    }
2545    double saveSj=0.0;
2546    if (info->crucialSj()>=0) 
2547      saveSj=solution_[info->crucialSj()];
2548    // Set dual solution
2549    for (iRow=0;iRow<numberQuadraticRows;iRow++) {
2550      solution_[iRow+numberXColumns]=dual_[iRow];
2551    }
2552    // clear sj
2553    int start = numberXColumns+numberQuadraticRows;
2554    memset(solution_+start,0,numberQuadraticColumns*sizeof(double));
2555    memset(dj_+numberXColumns,0,(numberQuadraticRows+numberQuadraticColumns)*sizeof(double));
2556    matrix_->times(-1.0,columnActivityWork_,modifiedCost,rowScale_,columnScale_);
2557    memset(modifiedCost,0,numberXRows*sizeof(double));
2558    // Do costs as rhs and sj solution values
2559    for (iRow=0;iRow<numberQuadraticColumns;iRow++) {
2560      int jSequence = iRow;
2561      double value=cost_[jSequence];
2562      if (fabs(value)>1.0e5) {
2563        assert (getColumnStatus(jSequence)==basic);
2564      }
2565      double value2=-modifiedCost[iRow+numberXRows];
2566      modifiedCost[iRow+numberXRows]=0.0;
2567      jSequence = iRow+start;
2568      if (getColumnStatus(jSequence)!=basic) {
2569        if (fabs(value2-value)>1.0e-3)
2570          solution_[jSequence]=0.0;
2571        value=value2;
2572      } else {
2573        solution_[jSequence]=value-value2;
2574      }
2575      storedCost[iRow]=value;
2576      storedUpper[iRow]=value;
2577      storedValue[iRow]=value;
2578    }
2579    if (info->crucialSj()>=0) 
2580      solution_[info->crucialSj()]=saveSj;
2581  }
2582  if (numberIterations_==-1289) {
2583    int i;
2584    printf ("normal\n");
2585    for (i=0;i<numberRows_+numberColumns_;i++) {
2586      char x='N';
2587      if (getColumnStatus(i)==basic)
2588        x='B';
2589      printf("CH %d %c sol %g dj %g\n",i,x,solution_[i],dj_[i]);
2590    }
2591  }
2592}
2593
2594int 
2595ClpSimplexPrimalQuadratic::checkComplementarity (ClpQuadraticInfo * info,
2596                            CoinIndexedVector * array1, CoinIndexedVector * array2)
2597{
2598  createDjs(info,array1,array2);
2599  int numberXColumns = info->numberXColumns();
2600  int numberXRows = info->numberXRows();
2601  int start=numberXColumns+numberXRows;
2602  numberDualInfeasibilities_=0;
2603  sumDualInfeasibilities_=0.0;
2604  // Compute objective function from scratch
2605  const CoinPackedMatrix * quadratic = 
2606    info->quadraticObjective();
2607  const int * columnQuadratic = quadratic->getIndices();
2608  const CoinBigIndex * columnQuadraticStart = quadratic->getVectorStarts();
2609  const int * columnQuadraticLength = quadratic->getVectorLengths();
2610  const double * quadraticElement = quadratic->getElements();
2611  const double * originalCost = info->linearObjective();
2612  int iColumn;
2613  objectiveValue_=0.0;
2614  double infeasCost=0.0;
2615  double linearCost=0.0;
2616  int number0=0,number1=0,number2=0;
2617#if 0
2618  int numberNSj=0;
2619  for (iColumn=numberXColumns+info->numberQuadraticRows();iColumn<numberColumns_;iColumn++) {
2620    if (getColumnStatus(iColumn)!=basic)
2621      numberNSj++;
2622  }
2623  printf("NumberN %d\n",numberNSj);
2624#endif
2625  for (iColumn=0;iColumn<numberXColumns;iColumn++) {
2626    double valueI = solution_[iColumn];
2627    linearCost += valueI*originalCost[iColumn];
2628    double diff =cost_[iColumn]-originalCost[iColumn];
2629    // WE are always feasible so look at low not up!
2630    if (diff>0.0) {
2631      double above = valueI-lower_[iColumn];
2632      assert(above>=0.0);
2633      infeasCost += diff*above;
2634    } else if (diff<0.0) {
2635      double below = upper_[iColumn]-valueI;
2636      assert(below>=0.0);
2637      infeasCost -= diff*below;
2638    }
2639    int j;
2640    for (j=columnQuadraticStart[iColumn];
2641         j<columnQuadraticStart[iColumn]+columnQuadraticLength[iColumn];j++) {
2642      int jColumn = columnQuadratic[j];
2643      double valueJ = solution_[jColumn];
2644      double elementValue = 0.5*quadraticElement[j];
2645      objectiveValue_ += valueI*valueJ*elementValue;
2646    }
2647    int jSequence = iColumn+start;
2648    double value=dj_[iColumn];
2649    ClpSimplex::Status status = getColumnStatus(iColumn);
2650    if (status!=basic&&jSequence>=0) {
2651      if (getColumnStatus(jSequence)==basic) 
2652        number1++;
2653      else
2654        number0++;
2655    }
2656   
2657    switch(status) {
2658     
2659    case ClpSimplex::basic:
2660      if (getColumnStatus(jSequence)==basic) {
2661        handler_->message(CLP_QUADRATIC_BOTH,messages_)
2662          <<"Structural"<<iColumn
2663          <<solution_[iColumn]<<jSequence<<solution_[jSequence]
2664          <<CoinMessageEol;
2665        number2++;
2666        assert (info->crucialSj()>=0);
2667        assert (info->crucialSj()==iColumn||info->crucialSj()==jSequence);
2668      } else {
2669        number1++;
2670      }
2671      break;
2672    case ClpSimplex::isFixed:
2673      break;
2674    case ClpSimplex::isFree:
2675    case ClpSimplex::superBasic:
2676      if (fabs(value)>dualTolerance_) {
2677        sumDualInfeasibilities_ += fabs(value)-dualTolerance_;
2678        numberDualInfeasibilities_++;
2679      }
2680      break;
2681    case ClpSimplex::atUpperBound:
2682      if (value>dualTolerance_) {
2683        sumDualInfeasibilities_ += value-dualTolerance_;
2684        numberDualInfeasibilities_++;
2685      }
2686      break;
2687    case ClpSimplex::atLowerBound:
2688      if (value<-dualTolerance_) {
2689        sumDualInfeasibilities_ -= value+dualTolerance_;
2690        numberDualInfeasibilities_++;
2691      }
2692    }
2693  }
2694  int offset = numberXColumns;
2695  int iRow;
2696  for (iRow=0;iRow<numberXRows;iRow++) {
2697    int iColumn = iRow + numberColumns_;
2698    double diff =cost_[iColumn];
2699    double valueI = solution_[iColumn];
2700    if (diff>0.0) {
2701      double above = valueI-lower_[iColumn];
2702      assert(above>=0.0);
2703      infeasCost += diff*above;
2704    } else if (diff<0.0) {
2705      double below = upper_[iColumn]-valueI;
2706      assert(below>=0.0);
2707      infeasCost -= diff*below;
2708    }
2709    int jSequence = iRow+offset;
2710    double value=dj_[iRow+numberColumns_];
2711    ClpSimplex::Status status = getRowStatus(iRow);
2712    if (status!=basic) {
2713      if (getColumnStatus(jSequence)==basic) 
2714        number1++;
2715      else
2716        number0++;
2717    }
2718   
2719    switch(status) {
2720     
2721    case ClpSimplex::basic:
2722      if (getColumnStatus(jSequence)==basic) {
2723        printf("Row %d (%g) and %d (%g) both basic\n",
2724               iRow,solution_[iColumn],jSequence,solution_[jSequence]);
2725          number2++;
2726          assert (info->crucialSj()>=0);
2727          assert (info->crucialSj()==iColumn||info->crucialSj()==jSequence);
2728        } else {
2729          number1++;
2730        }
2731      break;
2732    case ClpSimplex::isFixed:
2733      break;
2734    case ClpSimplex::isFree:
2735    case ClpSimplex::superBasic:
2736      if (fabs(value)>dualTolerance_) {
2737        sumDualInfeasibilities_ += fabs(value)-dualTolerance_;
2738        numberDualInfeasibilities_++;
2739      }
2740      break;
2741    case ClpSimplex::atUpperBound:
2742      if (value>dualTolerance_) {
2743        sumDualInfeasibilities_ += value-dualTolerance_;
2744        numberDualInfeasibilities_++;
2745      }
2746      break;
2747    case ClpSimplex::atLowerBound:
2748      if (value<-dualTolerance_) {
2749        sumDualInfeasibilities_ -= value+dualTolerance_;
2750        numberDualInfeasibilities_++;
2751      }
2752    }
2753  }
2754  //printf("number 0 - %d, 1 - %d, 2 - %d\n",number0,number1,number2);
2755  objectiveValue_ += linearCost+infeasCost;
2756  assert (infeasCost>=0.0);
2757  sumOfRelaxedDualInfeasibilities_ =sumDualInfeasibilities_;
2758  // But not zero if cleanup iteration
2759  if (info->sequenceIn()>=0&&!numberDualInfeasibilities_)
2760    numberDualInfeasibilities_=1;
2761  numberDualInfeasibilitiesWithoutFree_= numberDualInfeasibilities_;
2762  info->setInfeasCost(infeasCost);
2763  return numberDualInfeasibilities_;
2764}
2765 
2766/* This creates the large version of QP and
2767      fills in quadratic information
2768*/
2769ClpSimplexPrimalQuadratic * 
2770ClpSimplexPrimalQuadratic::makeQuadratic(ClpQuadraticInfo & info)
2771{
2772
2773  // Get list of non linear columns
2774  ClpQuadraticObjective * quadraticObj = (dynamic_cast< ClpQuadraticObjective*>(objective_));
2775  assert (quadraticObj);
2776  CoinPackedMatrix * quadratic = quadraticObj->quadraticObjective();
2777  if (!quadratic||!quadratic->getNumElements()) {
2778    // no quadratic part
2779    return NULL;
2780  }
2781
2782  int numberColumns = this->numberColumns();
2783  double * columnLower = this->columnLower();
2784  double * columnUpper = this->columnUpper();
2785  double * objective = this->objective();
2786  double * solution = this->primalColumnSolution();
2787  double * dj = this->dualColumnSolution();
2788  double * pi = this->dualRowSolution();
2789
2790  int numberRows = this->numberRows();
2791  double * rowLower = this->rowLower();
2792  double * rowUpper = this->rowUpper();
2793
2794  // and elements
2795  CoinPackedMatrix * matrix = this->matrix();
2796  const int * row = matrix->getIndices();
2797  const CoinBigIndex * columnStart = matrix->getVectorStarts();
2798  const double * element =  matrix->getElements();
2799  const int * columnLength = matrix->getVectorLengths();
2800
2801  int iColumn;
2802  const int * columnQuadratic = quadratic->getIndices();
2803  const CoinBigIndex * columnQuadraticStart = quadratic->getVectorStarts();
2804  const int * columnQuadraticLength = quadratic->getVectorLengths();
2805  const double * quadraticElement = quadratic->getElements();
2806#if 0
2807  // deliberate bad solution
2808  // Change to use phase
2809  //double * saveO = new double[numberColumns];
2810  //memcpy(saveO,objective,numberColumns*sizeof(double));
2811  //memset(objective,0,numberColumns*sizeof(double));
2812  tempMessage messageHandler(this);;
2813  passInMessageHandler(&messageHandler);
2814  factorization()->maximumPivots(1);
2815  primal();
2816  CoinMessageHandler handler2;
2817  passInMessageHandler(&handler2);
2818  factorization()->maximumPivots(100);
2819  setMaximumIterations(1000);
2820#endif
2821  //memcpy(objective,saveO,numberColumns*sizeof(double));
2822  // Get a feasible solution
2823  //printf("For testing - deliberate bad solution\n");
2824  //columnUpper[0]=0.0;
2825  //columnLower[0]=0.0;
2826  //quadraticSLP(50,1.0e-4);
2827  //primal(1);
2828  //columnUpper[0]=COIN_DBL_MAX;
2829 
2830  // Create larger problem
2831  // First workout how many rows extra
2832  info=ClpQuadraticInfo(this);
2833  quadratic = info.quadraticObjective();
2834  int numberQuadratic = info.numberQuadraticColumns();
2835  int newNumberRows = numberRows+numberQuadratic;
2836  int newNumberColumns = numberColumns + numberRows + numberQuadratic;
2837  int numberElements = 2*matrix->getNumElements()
2838    +2*quadratic->getNumElements()
2839    + numberQuadratic;
2840  double * elements2 = new double[numberElements];
2841  CoinBigIndex * start2 = new CoinBigIndex [newNumberColumns+1];
2842  int * row2 = new int[numberElements];
2843  double * objective2 = new double[newNumberColumns];
2844  double * columnLower2 = new double[newNumberColumns];
2845  double * columnUpper2 = new double[newNumberColumns];
2846  double * rowLower2 = new double[newNumberRows];
2847  double * rowUpper2 = new double[newNumberRows];
2848  memcpy(rowLower2,rowLower,numberRows*sizeof(double));
2849  memcpy(rowUpper2,rowUpper,numberRows*sizeof(double));
2850  int iRow;
2851  for (iRow=0;iRow<numberQuadratic;iRow++) {
2852    double cost = objective[iRow];
2853    rowLower2[iRow+numberRows]=cost;
2854    rowUpper2[iRow+numberRows]=cost;
2855  }
2856  memset(objective2,0,newNumberColumns*sizeof(double));
2857  // Get a row copy of quadratic objective in standard format
2858  CoinPackedMatrix copyQ;
2859  copyQ.reverseOrderedCopyOf(*quadratic);
2860  const int * columnQ = copyQ.getIndices();
2861  const CoinBigIndex * rowStartQ = copyQ.getVectorStarts();
2862  const int * rowLengthQ = copyQ.getVectorLengths(); 
2863  const double * elementByRowQ = copyQ.getElements();
2864  // Move solution across
2865  double * solution2 = new double[newNumberColumns];
2866  memset(solution2,0,newNumberColumns*sizeof(double));
2867  newNumberColumns=0;
2868  numberElements=0;
2869  start2[0]=0;
2870  // x
2871  memcpy(dj,objective,numberColumns*sizeof(double));
2872  for (iColumn=0;iColumn<numberColumns;iColumn++) {
2873    // Original matrix
2874    columnLower2[iColumn]=columnLower[iColumn];
2875    columnUpper2[iColumn]=columnUpper[iColumn];
2876    solution2[iColumn]=solution[iColumn];
2877    // Put in cost so we know it
2878    objective2[iColumn]=objective[iColumn];
2879    int j;
2880    for (j=columnStart[iColumn];
2881         j<columnStart[iColumn]+columnLength[iColumn];
2882         j++) {
2883      elements2[numberElements]=element[j];
2884      row2[numberElements++]=row[j];
2885    }
2886    // Quadratic and modify djs
2887    for (j=columnQuadraticStart[iColumn];
2888         j<columnQuadraticStart[iColumn]+columnQuadraticLength[iColumn];
2889         j++) {
2890      int jColumn = columnQuadratic[j];
2891      double value = quadraticElement[j];
2892      if (iColumn!=jColumn) 
2893        value *= 0.5;
2894      dj[jColumn] += solution[iColumn]*value;
2895      elements2[numberElements]=-value;
2896      row2[numberElements++]=jColumn+numberRows;
2897    }
2898    for (j=rowStartQ[iColumn];
2899         j<rowStartQ[iColumn]+rowLengthQ[iColumn];
2900         j++) {
2901      int jColumn = columnQ[j];
2902      double value = elementByRowQ[j];
2903      if (iColumn!=jColumn) { 
2904        value *= 0.5;
2905        dj[jColumn] += solution[iColumn]*value;
2906        elements2[numberElements]=-value;
2907        row2[numberElements++]=jColumn+numberRows;
2908      }
2909    }
2910    start2[iColumn+1]=numberElements;
2911  }
2912  newNumberColumns=numberColumns;
2913  // pi
2914  int numberQuadraticRows = info.numberQuadraticRows();
2915  // Get a row copy in standard format
2916  CoinPackedMatrix copy;
2917  copy.reverseOrderedCopyOf(*(this->matrix()));
2918  // get matrix data pointers
2919  const int * column = copy.getIndices();
2920  const CoinBigIndex * rowStart = copy.getVectorStarts();
2921  const int * rowLength = copy.getVectorLengths(); 
2922  const double * elementByRow = copy.getElements();
2923  for (iRow=0;iRow<numberRows;iRow++) {
2924    solution2[newNumberColumns]=pi[iRow];
2925    double value = pi[iRow];
2926    columnLower2[newNumberColumns]=-COIN_DBL_MAX;
2927    columnUpper2[newNumberColumns]=COIN_DBL_MAX;
2928    int j;
2929    for (j=rowStart[iRow];
2930         j<rowStart[iRow]+rowLength[iRow];
2931         j++) {
2932      double elementValue=elementByRow[j];
2933      int jColumn = column[j];
2934      if (jColumn>=0) {
2935        elements2[numberElements]=elementValue;
2936        row2[numberElements++]=jColumn+numberRows;
2937      }
2938      dj[jColumn]-= value*elementValue;
2939    }
2940#ifndef CORRECT_ROW_COUNTS
2941    newNumberColumns++;
2942    start2[newNumberColumns]=numberElements;
2943#else
2944    if (numberElements>start2[newNumberColumns]) {
2945      newNumberColumns++;
2946      start2[newNumberColumns]=numberElements;
2947    }
2948#endif
2949  }
2950  // djs
2951  for (iColumn=0;iColumn<numberQuadratic;iColumn++) {
2952    columnLower2[newNumberColumns]=-COIN_DBL_MAX;
2953    columnUpper2[newNumberColumns]=COIN_DBL_MAX;
2954    solution2[newNumberColumns]=dj[iColumn];
2955    elements2[numberElements]=1.0;
2956    row2[numberElements++]=iColumn+numberRows;
2957    newNumberColumns++;
2958    start2[newNumberColumns]=numberElements;
2959  }
2960  // Create model
2961  ClpSimplex * model2 = new ClpSimplex(*this);
2962  model2->resize(0,0);
2963  model2->loadProblem(newNumberColumns,newNumberRows,
2964                     start2,row2, elements2,
2965                     columnLower2,columnUpper2,
2966                     objective2,
2967                     rowLower2,rowUpper2);
2968  delete [] rowLower2;
2969  delete [] rowUpper2;
2970  delete [] columnLower2;
2971  delete [] columnUpper2;
2972  // Now create expanded quadratic objective for use in primalRow
2973  // Later pack down in some way
2974  start2[0]=0;
2975  numberElements=0;
2976  for (iColumn=0;iColumn<numberColumns;iColumn++) {
2977    // Quadratic
2978    int j;
2979    for (j=columnQuadraticStart[iColumn];
2980         j<columnQuadraticStart[iColumn]+columnQuadraticLength[iColumn];
2981         j++) {
2982      int jColumn = columnQuadratic[j];
2983      double value = quadraticElement[j];
2984      if (iColumn!=jColumn) 
2985        value *= 0.5;
2986      elements2[numberElements]=value;
2987      row2[numberElements++]=jColumn;
2988    }
2989    for (j=rowStartQ[iColumn];
2990         j<rowStartQ[iColumn]+rowLengthQ[iColumn];
2991         j++) {
2992      int jColumn = columnQ[j];
2993      double value = elementByRowQ[j];
2994      if (iColumn!=jColumn) { 
2995        value *= 0.5;
2996        elements2[numberElements]=value;
2997        row2[numberElements++]=jColumn;
2998      }
2999    }
3000    start2[iColumn+1]=numberElements;
3001  }
3002  // and pad
3003  //for (;iColumn<newNumberColumns;iColumn++)
3004  //start2[iColumn+1]=numberElements;
3005  // Load up objective with expanded linear
3006  ClpQuadraticObjective * obj = 
3007    new ClpQuadraticObjective(objective2,numberColumns,
3008                              start2,row2,elements2,newNumberColumns);
3009  delete [] objective2;
3010  info.setOriginalObjective(obj);
3011  //model2->loadQuadraticObjective(newNumberColumns,start2,row2,elements2);
3012  delete [] start2;
3013  delete [] row2;
3014  delete [] elements2;
3015  model2->allSlackBasis();
3016  //model2->scaling(false);
3017  //printf("scaling off\n");
3018  model2->setLogLevel(this->logLevel());
3019  // Move solution across
3020  memcpy(model2->primalColumnSolution(),solution2,
3021         newNumberColumns*sizeof(double));
3022  columnLower2 = model2->columnLower();
3023  columnUpper2 = model2->columnUpper();
3024  delete [] solution2;
3025  solution2 = model2->primalColumnSolution();
3026  // Compute row activities and check feasible
3027  double * rowSolution2 = model2->primalRowSolution();
3028  memset(rowSolution2,0,newNumberRows*sizeof(double));
3029  model2->times(1.0,solution2,rowSolution2);
3030  rowLower2 = model2->rowLower();
3031  rowUpper2 = model2->rowUpper();
3032#if 0
3033  // redo as Dantzig says
3034  for (iColumn=0;iColumn<numberColumns;iColumn++) {
3035    Status xStatus = getColumnStatus(iColumn);
3036    bool isSuperBasic;
3037    int iS = iColumn+newNumberRows;
3038    model2->setColumnStatus(iColumn,xStatus);
3039    if (xStatus==basic) {
3040      model2->setColumnStatus(iS,isFree);
3041      solution2[iS]=0.0;
3042    } else {
3043      model2->setColumnStatus(iS,basic);
3044    }
3045    // take slack out on equality rows
3046    model2->setRowBasic(iColumn+numberRows,superBasic);
3047  }
3048  for (iRow=0;iRow<numberRows;iRow++) {
3049    Status rowStatus = getRowStatus(iRow);
3050    model2->setRowStatus(iRow,rowStatus);
3051    if (rowStatus!=basic) {
3052      model2->setColumnStatus(iRow+numberColumns,basic); // make dual basic
3053    }
3054    assert (rowSolution2[iRow]>=rowLower2[iRow]-primalTolerance_);
3055    assert (rowSolution2[iRow]<=rowUpper2[iRow]+primalTolerance_);
3056  }
3057  // why ?? - take duals out and adjust
3058  for (iRow=0;iRow<numberRows;iRow++) {
3059    model2->setRowStatus(iRow,basic);
3060    model2->setColumnStatus(iRow+numberColumns,superBasic);
3061    solution2[iRow+numberColumns]=0.0;
3062  }
3063#else
3064  for (iRow=0;iRow<numberRows;iRow++) {
3065    assert (rowSolution2[iRow]>=rowLower2[iRow]-primalTolerance_);
3066    assert (rowSolution2[iRow]<=rowUpper2[iRow]+primalTolerance_);
3067    model2->setRowStatus(iRow,basic);
3068    int jRow = iRow;
3069    model2->setColumnStatus(jRow+numberColumns,superBasic);
3070    solution2[jRow+numberColumns]=0.0;
3071  }
3072  for (iColumn=numberQuadraticRows+numberColumns;iColumn<newNumberColumns;iColumn++) {
3073    model2->setColumnStatus(iColumn,basic);
3074    model2->setRowStatus(iColumn-numberQuadraticRows-numberColumns+numberRows,isFixed);
3075  }
3076#endif
3077  memset(rowSolution2,0,newNumberRows*sizeof(double));
3078  model2->times(1.0,solution2,rowSolution2);
3079  for (iColumn=0;iColumn<numberQuadratic;iColumn++) {
3080    int iS = iColumn+numberColumns+numberQuadraticRows;
3081    int iRow = iColumn+numberRows;
3082    double value = rowSolution2[iRow];
3083    if (value>rowUpper2[iRow]) {
3084      rowSolution2[iRow] = rowUpper2[iRow];
3085      solution2[iS]-=value-rowUpper2[iRow];
3086    } else {
3087      rowSolution2[iRow] = rowLower2[iRow];
3088      solution2[iS]-=value-rowLower2[iRow];
3089    }
3090  }
3091
3092 
3093  // See if any s sub j have wrong sign and/or use djs from infeasibility objective
3094  double objectiveOffset;
3095  getDblParam(ClpObjOffset,objectiveOffset);
3096  double objValue = -objectiveOffset;
3097  for (iColumn=0;iColumn<numberColumns;iColumn++) 
3098    objValue += objective[iColumn]*solution2[iColumn];
3099  for (iColumn=0;iColumn<numberColumns;iColumn++) {
3100    double valueI = solution2[iColumn];
3101    int j;
3102    for (j=columnQuadraticStart[iColumn];
3103         j<columnQuadraticStart[iColumn]+columnQuadraticLength[iColumn];j++) {
3104      int jColumn = columnQuadratic[j];
3105      double valueJ = solution2[jColumn];
3106      double elementValue = quadraticElement[j];
3107      objValue += 0.5*valueI*valueJ*elementValue;
3108    }
3109  }
3110  printf("Objective value %g\n",objValue);
3111  //for (iColumn=0;iColumn<newNumberColumns;iColumn++)
3112  //printf("%d %g\n",iColumn,solution2[iColumn]);
3113  return (ClpSimplexPrimalQuadratic *) model2;
3114}
3115
3116// This moves solution back and deletes information
3117int 
3118ClpSimplexPrimalQuadratic::endQuadratic(ClpSimplexPrimalQuadratic * quadraticModel,
3119                   ClpQuadraticInfo & info)
3120{
3121  memcpy(dualRowSolution(),quadraticModel->dualRowSolution(),numberRows_*sizeof(double));
3122  const double * solution2 = quadraticModel->primalColumnSolution();
3123  memcpy(primalColumnSolution(),solution2,numberColumns_*sizeof(double));
3124  memset(quadraticModel->primalRowSolution(),0,
3125         quadraticModel->numberRows()*sizeof(double));
3126  quadraticModel->times(1.0,quadraticModel->primalColumnSolution(),
3127               quadraticModel->primalRowSolution());
3128
3129  int iColumn;
3130  double objectiveOffset;
3131  getDblParam(ClpObjOffset,objectiveOffset);
3132  double objValue = -objectiveOffset;
3133  double * objective = this->objective();
3134  const double * pi = dualRowSolution();
3135  for (iColumn=0;iColumn<numberColumns_;iColumn++) 
3136    objValue += objective[iColumn]*solution2[iColumn];
3137  ClpQuadraticObjective * quadraticObj = (dynamic_cast< ClpQuadraticObjective*>(objective_));
3138  assert (quadraticObj);
3139  CoinPackedMatrix * quadratic = quadraticObj->quadraticObjective();
3140  double * dj = dualColumnSolution();
3141  // Matrix for linear stuff
3142  const int * row = matrix_->getIndices();
3143  const CoinBigIndex * columnStart = matrix_->getVectorStarts();
3144  const double * element =  matrix_->getElements();
3145  const int * columnLength = matrix_->getVectorLengths();
3146  if (quadratic) {
3147    const int * columnQuadratic = quadratic->getIndices();
3148    const CoinBigIndex * columnQuadraticStart = quadratic->getVectorStarts();
3149    const int * columnQuadraticLength = quadratic->getVectorLengths();
3150    const double * quadraticElement = quadratic->getElements();
3151    int start = info.numberQuadraticRows()+numberColumns_;
3152    for (iColumn=0;iColumn<numberColumns_;iColumn++) {
3153      int jColumn = iColumn;
3154      double valueI=solution2[iColumn];
3155      double value;
3156      if (jColumn>=0) {
3157        jColumn += start;
3158        value = solution2[jColumn];
3159      } else {
3160        value=objective[iColumn];
3161        int j;
3162        for (j=columnStart[iColumn];j<columnStart[iColumn]+columnLength[iColumn]; j++) {
3163          int iRow = row[j];
3164          value -= element[j]*pi[iRow];
3165        }
3166      }
3167      dj[iColumn]=value;
3168      Status status = quadraticModel->getColumnStatus(iColumn);
3169      setColumnStatus(iColumn,status);
3170      if (status==basic) {
3171        assert(fabs(value)<1.0e-2);
3172      }
3173      int j;
3174      for (j=columnQuadraticStart[iColumn];
3175           j<columnQuadraticStart[iColumn]+columnQuadraticLength[iColumn];j++) {
3176        int jColumn = columnQuadratic[j];
3177        double valueJ = solution2[jColumn];
3178        double elementValue = quadraticElement[j];
3179        objValue += 0.5*valueI*valueJ*elementValue;
3180      }
3181    }
3182    objectiveValue_ = objValue + objectiveOffset;
3183  } else {
3184    // Do linear bit anyway
3185    for (iColumn=0;iColumn<numberColumns_;iColumn++) {
3186      double value=objective[iColumn];
3187      int j;
3188      for (j=columnStart[iColumn];j<columnStart[iColumn]+columnLength[iColumn]; j++) {
3189        int iRow = row[j];
3190        value -= element[j]*pi[iRow];
3191      }
3192      dj[iColumn]=value;
3193      Status status = quadraticModel->getColumnStatus(iColumn);
3194      setColumnStatus(iColumn,status);
3195    }
3196  }
3197  // and row status
3198  int iRow;
3199  for (iRow=0;iRow<numberRows_;iRow++) {
3200    Status status = quadraticModel->getRowStatus(iRow);
3201    setRowStatus(iRow,status);
3202  }
3203  printf("Objective value %g\n",objValue);
3204  return 0;
3205}
3206
3207/// Default constructor.
3208ClpQuadraticInfo::ClpQuadraticInfo()
3209  : originalObjective_(NULL),
3210    basicRow_(NULL),
3211    impliedSj_(NULL),
3212    currentSequenceIn_(-1),
3213    crucialSj_(-1),
3214    validSequenceIn_(-1),
3215    validCrucialSj_(-1),
3216    currentPhase_(-1),
3217    currentSolution_(NULL),
3218    validPhase_(-1),
3219    validSolution_(NULL),
3220    djWeight_(NULL),
3221    gradient_(NULL),
3222    numberXRows_(-1),
3223    numberXColumns_(-1),
3224    numberQuadraticColumns_(0),
3225    numberQuadraticRows_(0),
3226    infeasCost_(0.0)
3227{
3228}
3229// Constructor from original model
3230ClpQuadraticInfo::ClpQuadraticInfo(const ClpSimplex * model)
3231  : originalObjective_(NULL),
3232    basicRow_(NULL),
3233    impliedSj_(NULL),
3234    currentSequenceIn_(-1),
3235    crucialSj_(-1),
3236    validSequenceIn_(-1),
3237    validCrucialSj_(-1),
3238    currentPhase_(-1),
3239    currentSolution_(NULL),
3240    validPhase_(-1),
3241    validSolution_(NULL),
3242    djWeight_(NULL),
3243    gradient_(NULL),
3244    numberXRows_(-1),
3245    numberXColumns_(-1),
3246    numberQuadraticColumns_(0),
3247    numberQuadraticRows_(0),
3248    infeasCost_(0.0)
3249{
3250  if (model) {
3251    numberXRows_ = model->numberRows();
3252    numberXColumns_ = model->numberColumns();
3253    //ClpQuadraticObjective *originalObjective = (dynamic_cast< ClpQuadraticObjective*>(model->objectiveAsObject()));
3254    //assert (originalObjective);
3255    originalObjective_ = (dynamic_cast< ClpQuadraticObjective*>(model->objectiveAsObject()));
3256    assert (originalObjective_);
3257    impliedSj_ = new int[numberXColumns_];
3258    basicRow_ = new int[numberXColumns_];
3259    int i;
3260    numberQuadraticColumns_=numberXColumns_;
3261    numberQuadraticRows_=numberXRows_;
3262    int numberColumns = numberQuadraticRows_+numberXColumns_+numberQuadraticColumns_;
3263    int numberRows = numberXRows_+numberQuadraticColumns_;
3264    int size = numberRows+numberColumns;
3265    djWeight_ = new double [size];
3266    basicRow_ = new int[size];
3267    for (i=0;i<size;i++)
3268      djWeight_[i]=1.0;
3269  }
3270}
3271// Destructor
3272ClpQuadraticInfo:: ~ClpQuadraticInfo()
3273{
3274  delete [] impliedSj_;
3275  delete [] basicRow_;
3276  delete [] currentSolution_;
3277  delete [] validSolution_;
3278  delete [] djWeight_;
3279  delete [] gradient_;
3280}
3281// Copy
3282ClpQuadraticInfo::ClpQuadraticInfo(const ClpQuadraticInfo& rhs)
3283  : originalObjective_(rhs.originalObjective_),
3284    basicRow_(NULL),
3285    impliedSj_(NULL),
3286    currentSequenceIn_(rhs.currentSequenceIn_),
3287    crucialSj_(rhs.crucialSj_),
3288    validSequenceIn_(rhs.validSequenceIn_),
3289    validCrucialSj_(rhs.validCrucialSj_),
3290    currentPhase_(rhs.currentPhase_),
3291    currentSolution_(NULL),
3292    validPhase_(rhs.validPhase_),
3293    validSolution_(NULL),
3294    djWeight_(NULL),
3295    gradient_(NULL),
3296    numberXRows_(rhs.numberXRows_),
3297    numberXColumns_(rhs.numberXColumns_),
3298    numberQuadraticColumns_(rhs.numberQuadraticColumns_),
3299    numberQuadraticRows_(rhs.numberQuadraticRows_),
3300    infeasCost_(rhs.infeasCost_)
3301{
3302  if (numberXColumns_) {
3303    int numberColumns = numberQuadraticRows_+numberXColumns_+numberQuadraticColumns_;
3304    int numberRows = numberXRows_+numberQuadraticColumns_;
3305    int size = numberRows+numberColumns;
3306    impliedSj_ = new int[numberXColumns_];
3307    memcpy(impliedSj_,rhs.impliedSj_,numberXColumns_*sizeof(int));
3308    basicRow_ = new int [size];
3309    memcpy(basicRow_,rhs.basicRow_,size*sizeof(int));
3310    if (rhs.currentSolution_) {
3311      currentSolution_ = new double [size];
3312      memcpy(currentSolution_,rhs.currentSolution_,size*sizeof(double));
3313    } else {
3314      currentSolution_ = NULL;
3315    }
3316    if (rhs.validSolution_) {
3317      validSolution_ = new double [size];
3318      memcpy(validSolution_,rhs.validSolution_,size*sizeof(double));
3319    } else {
3320      validSolution_ = NULL;
3321    }
3322    if (rhs.djWeight_) {
3323      djWeight_ = new double [size];
3324      memcpy(djWeight_,rhs.djWeight_,size*sizeof(double));
3325    } else {
3326      djWeight_ = NULL;
3327    }
3328    if (rhs.gradient_) {
3329      gradient_ = new double [size];
3330      memcpy(gradient_,rhs.gradient_,size*sizeof(double));
3331    } else {
3332      gradient_ = NULL;
3333    }
3334  }
3335}
3336// Assignment
3337ClpQuadraticInfo & 
3338ClpQuadraticInfo::operator=(const ClpQuadraticInfo&rhs)
3339{
3340  if (this != &rhs) {
3341    originalObjective_ = rhs.originalObjective_;
3342    delete [] impliedSj_;
3343    delete [] basicRow_;
3344    delete [] currentSolution_;
3345    delete [] validSolution_;
3346    delete [] djWeight_;
3347    delete [] gradient_;
3348    currentSequenceIn_ = rhs.currentSequenceIn_;
3349    crucialSj_ = rhs.crucialSj_;
3350    validSequenceIn_ = rhs.validSequenceIn_;
3351    validCrucialSj_ = rhs.validCrucialSj_;
3352    currentPhase_ = rhs.currentPhase_;
3353    validPhase_ = rhs.validPhase_;
3354    numberXRows_ = rhs.numberXRows_;
3355    numberXColumns_ = rhs.numberXColumns_;
3356    infeasCost_=rhs.infeasCost_;
3357    numberQuadraticColumns_=rhs.numberQuadraticColumns_;
3358    numberQuadraticRows_=rhs.numberQuadraticRows_;
3359    int numberColumns = numberQuadraticRows_+numberXColumns_+numberQuadraticColumns_;
3360    int numberRows = numberXRows_+numberQuadraticColumns_;
3361    int size = numberRows+numberColumns;
3362    impliedSj_ = new int[numberXColumns_];
3363    memcpy(impliedSj_,rhs.impliedSj_,numberXColumns_*sizeof(int));
3364    basicRow_ = new int [size];
3365    memcpy(basicRow_,rhs.basicRow_,size*sizeof(int));
3366    if (rhs.currentSolution_) {
3367      currentSolution_ = new double [size];
3368      memcpy(currentSolution_,rhs.currentSolution_,size*sizeof(double));
3369    } else {
3370      currentSolution_ = NULL;
3371    }
3372    if (rhs.validSolution_) {
3373      validSolution_ = new double [size];
3374      memcpy(validSolution_,rhs.validSolution_,size*sizeof(double));
3375    } else {
3376      validSolution_ = NULL;
3377    }
3378    if (rhs.djWeight_) {
3379      djWeight_ = new double [size];
3380      memcpy(djWeight_,rhs.djWeight_,size*sizeof(double));
3381    } else {
3382      djWeight_ = NULL;
3383    }
3384    if (rhs.gradient_) {
3385      gradient_ = new double [size];
3386      memcpy(gradient_,rhs.gradient_,size*sizeof(double));
3387    } else {
3388      gradient_ = NULL;
3389    }
3390  }
3391  return *this;
3392}
3393// Save current In and Sj status
3394 void 
3395ClpQuadraticInfo::saveStatus()
3396{
3397  validSequenceIn_ = currentSequenceIn_;
3398  validCrucialSj_ = crucialSj_;
3399  validPhase_ = currentPhase_;
3400  int numberColumns = numberQuadraticRows_+numberXColumns_+numberQuadraticColumns_;
3401  int numberRows = numberXRows_+numberQuadraticColumns_;
3402  int size = numberRows+numberColumns;
3403  if (currentSolution_) {
3404    if (!validSolution_) 
3405      validSolution_ = new double [size];
3406    memcpy(validSolution_,currentSolution_,size*sizeof(double));
3407  }
3408}
3409// Restore previous
3410void 
3411ClpQuadraticInfo::restoreStatus()
3412{
3413  currentSequenceIn_ = validSequenceIn_;
3414  crucialSj_ = validCrucialSj_;
3415  currentPhase_ = validPhase_;
3416  delete [] currentSolution_;
3417  currentSolution_ = validSolution_;
3418  validSolution_=NULL;
3419}
3420void 
3421ClpQuadraticInfo::setCurrentSolution(const double * solution)
3422{
3423  if (currentPhase_) {
3424    int numberColumns = numberQuadraticRows_+numberXColumns_+numberQuadraticColumns_;
3425    int numberRows = numberXRows_+numberQuadraticColumns_;
3426    int size = numberRows+numberColumns;
3427    if (!currentSolution_) 
3428      currentSolution_ = new double [size];
3429    memcpy(currentSolution_,solution,size*sizeof(double));
3430  } else {
3431    delete [] currentSolution_;
3432    currentSolution_=NULL;
3433  }
3434}
3435// Quadratic objective
3436CoinPackedMatrix * 
3437ClpQuadraticInfo::quadraticObjective() const     
3438{ 
3439  return originalObjective_->quadraticObjective();
3440}
3441// Linear objective
3442double * 
3443ClpQuadraticInfo::linearObjective() const     
3444{ 
3445  return originalObjective_->linearObjective();
3446}
3447void 
3448ClpQuadraticInfo::createGradient(ClpSimplex * model)
3449{
3450  int numberColumns = numberQuadraticRows_+numberXColumns_+numberQuadraticColumns_;
3451  int numberRows = numberXRows_+numberQuadraticColumns_;
3452  int size = numberRows+numberColumns;
3453  if (!gradient_)
3454    gradient_= new double[size];
3455  memcpy(gradient_,model->costRegion(),size*sizeof(double));
3456  double * solution = model->solutionRegion();
3457  const CoinPackedMatrix * quadratic = quadraticObjective();
3458  const int * columnQuadratic = quadratic->getIndices();
3459  const CoinBigIndex * columnQuadraticStart = quadratic->getVectorStarts();
3460  const int * columnQuadraticLength = quadratic->getVectorLengths();
3461  const double * quadraticElement = quadratic->getElements();
3462  // get current costs
3463  int jSequence;
3464  for (jSequence=0;jSequence<numberQuadraticColumns_;jSequence++) {
3465    int iSequence = jSequence;
3466    // get current gradient
3467    double coeff1=gradient_[iSequence];
3468    int j;
3469    for (j=columnQuadraticStart[iSequence];
3470         j<columnQuadraticStart[iSequence]+columnQuadraticLength[iSequence];j++) {
3471      int jColumn = columnQuadratic[j];
3472      double valueJ = solution[jColumn];
3473      // maybe this is just if jColumn basic ??
3474      double elementValue = quadraticElement[j];
3475      double addValue = valueJ*elementValue;
3476      coeff1 += addValue;
3477    }
3478    gradient_[iSequence]=coeff1;
3479  }
3480}
Note: See TracBrowser for help on using the repository browser.