source: branches/pre/ClpSimplexPrimalQuadratic.cpp @ 182

Last change on this file since 182 was 182, checked in by forrest, 18 years ago

Finished for now

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 82.4 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 "ClpFactorization.hpp"
12#include "ClpNonLinearCost.hpp"
13#include "CoinPackedMatrix.hpp"
14#include "CoinIndexedVector.hpp"
15#include "CoinWarmStartBasis.hpp"
16#include "CoinMpsIO.hpp"
17#include "ClpPrimalColumnPivot.hpp"
18#include "ClpMessage.hpp"
19#include <cfloat>
20#include <cassert>
21#include <string>
22#include <stdio.h>
23#include <iostream>
24class tempMessage :
25   public CoinMessageHandler {
26
27public:
28  virtual int print() ;
29  tempMessage(ClpSimplex * model);
30  ClpSimplex * model_;
31};
32
33// Constructor with pointer to model
34tempMessage::tempMessage(ClpSimplex * model)
35  : CoinMessageHandler(),
36    model_(model)
37{
38}
39int
40tempMessage::print()
41{
42  static int numberFeasible=0;
43  if (currentSource()=="Clp") {
44    if (currentMessage().externalNumber()==5) { 
45      if (!numberFeasible&&!model_->nonLinearCost()->numberInfeasibilities()) {
46        numberFeasible++;
47        model_->setMaximumIterations(0);
48      }
49    } 
50  }
51  return CoinMessageHandler::print();
52}
53
54// A sequential LP method
55int 
56ClpSimplexPrimalQuadratic::primalSLP(int numberPasses, double deltaTolerance)
57{
58  // Are we minimizing or maximizing
59  double whichWay=optimizationDirection();
60  // This is as a user would see
61
62  int numberColumns = this->numberColumns();
63  int numberRows = this->numberRows();
64  double * columnLower = this->columnLower();
65  double * columnUpper = this->columnUpper();
66  double * objective = this->objective();
67  double * solution = this->primalColumnSolution();
68 
69  // Save objective
70 
71  double * saveObjective = new double [numberColumns];
72  memcpy(saveObjective,objective,numberColumns*sizeof(double));
73
74  // Get list of non linear columns
75  CoinPackedMatrix * quadratic = quadraticObjective();
76  if (!quadratic) {
77    // no quadratic part
78    return primal(0);
79  }
80  int numberNonLinearColumns = 0;
81  int iColumn;
82  int * listNonLinearColumn = new int[numberColumns];
83  memset(listNonLinearColumn,0,numberColumns*sizeof(int));
84  const int * columnQuadratic = quadratic->getIndices();
85  const int * columnQuadraticStart = quadratic->getVectorStarts();
86  const int * columnQuadraticLength = quadratic->getVectorLengths();
87  const double * quadraticElement = quadratic->getElements();
88  for (iColumn=0;iColumn<numberColumns;iColumn++) {
89    int j;
90    for (j=columnQuadraticStart[iColumn];
91         j<columnQuadraticStart[iColumn]+columnQuadraticLength[iColumn];j++) {
92      int jColumn = columnQuadratic[j];
93      listNonLinearColumn[jColumn]=1;
94      listNonLinearColumn[iColumn]=1;
95    }
96  }
97  for (iColumn=0;iColumn<numberColumns;iColumn++) {
98    if(listNonLinearColumn[iColumn])
99      listNonLinearColumn[numberNonLinearColumns++]=iColumn;
100  }
101 
102  if (!numberNonLinearColumns) {
103    delete [] listNonLinearColumn;
104    // no quadratic part
105    return primal(0);
106  }
107
108  // get feasible
109  if (!this->status()||numberPrimalInfeasibilities())
110    primal(1);
111  // still infeasible
112  if (numberPrimalInfeasibilities())
113    return 0;
114
115  int jNon;
116  int * last[3];
117 
118  double * trust = new double[numberNonLinearColumns];
119  double * trueLower = new double[numberNonLinearColumns];
120  double * trueUpper = new double[numberNonLinearColumns];
121  for (jNon=0;jNon<numberNonLinearColumns;jNon++) {
122    iColumn=listNonLinearColumn[jNon];
123    trust[jNon]=0.5;
124    trueLower[jNon]=columnLower[iColumn];
125    trueUpper[jNon]=columnUpper[iColumn];
126    if (solution[iColumn]<trueLower[jNon])
127      solution[iColumn]=trueLower[jNon];
128    else if (solution[iColumn]>trueUpper[jNon])
129      solution[iColumn]=trueUpper[jNon];
130  }
131  int iPass;
132  double lastObjective=1.0e31;
133  double * saveSolution = new double [numberColumns];
134  double * savePi = new double [numberRows];
135  unsigned char * saveStatus = new unsigned char[numberRows+numberColumns];
136  double targetDrop=1.0e31;
137  double objectiveOffset;
138  getDblParam(ClpObjOffset,objectiveOffset);
139  // 1 bound up, 2 up, -1 bound down, -2 down, 0 no change
140  for (iPass=0;iPass<3;iPass++) {
141    last[iPass]=new int[numberNonLinearColumns];
142    for (jNon=0;jNon<numberNonLinearColumns;jNon++) 
143      last[iPass][jNon]=0;
144  }
145  // goodMove +1 yes, 0 no, -1 last was bad - just halve gaps, -2 do nothing
146  int goodMove=-2;
147  char * statusCheck = new char[numberColumns];
148  for (iPass=0;iPass<numberPasses;iPass++) {
149    // redo objective
150    double offset=0.0;
151    double objValue=-objectiveOffset;
152    double lambda=-1.0;
153    if (goodMove>=0) {
154      // get best interpolation
155      double coeff0=-objectiveOffset,coeff1=0.0,coeff2=0.0;
156      for (iColumn=0;iColumn<numberColumns;iColumn++) {
157        coeff0 += saveObjective[iColumn]*solution[iColumn];
158        coeff1 += saveObjective[iColumn]*(saveSolution[iColumn]-solution[iColumn]);
159      }
160      for (jNon=0;jNon<numberNonLinearColumns;jNon++) {
161        iColumn=listNonLinearColumn[jNon];
162        double valueI = solution[iColumn];
163        double valueISave = saveSolution[iColumn];
164        int j;
165        for (j=columnQuadraticStart[iColumn];
166             j<columnQuadraticStart[iColumn]+columnQuadraticLength[iColumn];j++) {
167          int jColumn = columnQuadratic[j];
168          double valueJ = solution[jColumn];
169          double valueJSave = saveSolution[jColumn];
170          double elementValue = 0.5*quadraticElement[j];
171          coeff0 += valueI*valueJ*elementValue;
172          coeff1 += (valueI*valueJSave+valueISave*valueJ-2.0*valueI*valueJ)*elementValue;
173          coeff2 += (valueISave*valueJSave+valueI*valueJ-valueISave*valueJ-valueI*valueJSave)*elementValue;
174        }
175      }
176      double lambdaValue;
177      if (fabs(coeff2)<1.0e-9) {
178        if (coeff1+coeff2>=0.0) 
179          lambda = 0.0;
180        else
181          lambda = 1.0;
182      } else {
183        lambda = -(0.5*coeff1)/coeff2;
184        if (lambda>1.0||lambda<0.0) {
185          if (coeff1+coeff2>=0.0) 
186            lambda = 0.0;
187          else
188            lambda = 1.0;
189        }
190      }
191      lambdaValue = lambda*lambda*coeff2+lambda*coeff1+coeff0;
192      printf("coeffs %g %g %g - lastobj %g\n",coeff0,coeff1,coeff2,lastObjective);
193      printf("obj at saved %g, obj now %g zero deriv at %g - value %g\n",
194             coeff0+coeff1+coeff2,coeff0,lambda,lambdaValue);
195      if (!iPass) lambda=0.0;
196      if (lambda>0.0&&lambda<=1.0) {
197        // update solution
198        for (iColumn=0;iColumn<numberColumns;iColumn++) 
199          solution[iColumn] = lambda * saveSolution[iColumn] 
200            + (1.0-lambda) * solution[iColumn];
201        if (lambda>0.999) {
202          memcpy(this->dualRowSolution(),savePi,numberRows*sizeof(double));
203          memcpy(status_,saveStatus,numberRows+numberColumns);
204        }
205        if (lambda>0.99999&&fabs(coeff1+coeff2)>1.0e-2) {
206          // tighten all
207          goodMove=-1;
208        }
209      }
210    }
211    memcpy(objective,saveObjective,numberColumns*sizeof(double));
212    for (iColumn=0;iColumn<numberColumns;iColumn++) 
213      objValue += objective[iColumn]*solution[iColumn];
214    for (jNon=0;jNon<numberNonLinearColumns;jNon++) {
215      iColumn=listNonLinearColumn[jNon];
216      if (getColumnStatus(iColumn)==basic) {
217        if (solution[iColumn]<columnLower[iColumn]+1.0e-8)
218          statusCheck[iColumn]='l';
219        else if (solution[iColumn]>columnUpper[iColumn]-1.0e-8)
220          statusCheck[iColumn]='u';
221        else
222          statusCheck[iColumn]='B';
223      } else {
224        if (solution[iColumn]<columnLower[iColumn]+1.0e-8)
225          statusCheck[iColumn]='L';
226        else
227          statusCheck[iColumn]='U';
228      }
229      double valueI = solution[iColumn];
230      int j;
231      for (j=columnQuadraticStart[iColumn];
232           j<columnQuadraticStart[iColumn]+columnQuadraticLength[iColumn];j++) {
233        int jColumn = columnQuadratic[j];
234        double valueJ = solution[jColumn];
235        double elementValue = quadraticElement[j];
236        elementValue *= 0.5;
237        objValue += valueI*valueJ*elementValue;
238        offset += valueI*valueJ*elementValue;
239        double gradientI = valueJ*elementValue;
240        double gradientJ = valueI*elementValue;
241        offset -= gradientI*valueI;
242        objective[iColumn] += gradientI;
243        offset -= gradientJ*valueJ;
244        objective[jColumn] += gradientJ;
245      }
246    }
247    printf("objective %g, objective offset %g\n",objValue,offset);
248    setDblParam(ClpObjOffset,objectiveOffset-offset);
249    objValue *= whichWay;
250    int * temp=last[2];
251    last[2]=last[1];
252    last[1]=last[0];
253    last[0]=temp;
254    for (jNon=0;jNon<numberNonLinearColumns;jNon++) {
255      iColumn=listNonLinearColumn[jNon];
256      double change = solution[iColumn]-saveSolution[iColumn];
257      if (change<-1.0e-5) {
258        if (fabs(change+trust[jNon])<1.0e-5) 
259          temp[jNon]=-1;
260        else
261          temp[jNon]=-2;
262      } else if(change>1.0e-5) {
263        if (fabs(change-trust[jNon])<1.0e-5) 
264          temp[jNon]=1;
265        else
266          temp[jNon]=2;
267      } else {
268        temp[jNon]=0;
269      }
270    } 
271    // goodMove +1 yes, 0 no, -1 last was bad - just halve gaps, -2 do nothing
272    double maxDelta=0.0;
273    if (goodMove>=0) {
274      if (objValue<=lastObjective) 
275        goodMove=1;
276      else
277        goodMove=0;
278    } else {
279      maxDelta=1.0e10;
280    }
281    double maxGap=0.0;
282    for (jNon=0;jNon<numberNonLinearColumns;jNon++) {
283      iColumn=listNonLinearColumn[jNon];
284      maxDelta = max(maxDelta,
285                     fabs(solution[iColumn]-saveSolution[iColumn]));
286      if (goodMove>0) {
287        if (last[0][jNon]*last[1][jNon]<0) {
288          // halve
289          trust[jNon] *= 0.5;
290        } else {
291          if (last[0][jNon]==last[1][jNon]&&
292              last[0][jNon]==last[2][jNon])
293            trust[jNon] *= 1.5; 
294        }
295      } else if (goodMove!=-2&&trust[jNon]>10.0*deltaTolerance) {
296        trust[jNon] *= 0.2;
297      }
298      maxGap = max(maxGap,trust[jNon]);
299    }
300    std::cout<<"largest gap is "<<maxGap<<std::endl;
301    if (iPass>10000) {
302      for (jNon=0;jNon<numberNonLinearColumns;jNon++) 
303        trust[jNon] *=0.0001;
304    }
305    if (goodMove>0) {
306      double drop = lastObjective-objValue;
307      std::cout<<"True drop was "<<drop<<std::endl;
308      std::cout<<"largest delta is "<<maxDelta<<std::endl;
309      if (maxDelta<deltaTolerance&&drop<1.0e-4&&goodMove&&lambda<0.99999) {
310        std::cout<<"Exiting"<<std::endl;
311        break;
312      }
313    }
314    if (!iPass)
315      goodMove=1;
316    targetDrop=0.0;
317    double * r = this->dualColumnSolution();
318    for (jNon=0;jNon<numberNonLinearColumns;jNon++) {
319      iColumn=listNonLinearColumn[jNon];
320      columnLower[iColumn]=max(solution[iColumn]
321                               -trust[jNon],
322                               trueLower[jNon]);
323      columnUpper[iColumn]=min(solution[iColumn]
324                               +trust[jNon],
325                               trueUpper[jNon]);
326    }
327    if (iPass) {
328      // get reduced costs
329      this->matrix()->transposeTimes(savePi,
330                                     this->dualColumnSolution());
331      for (jNon=0;jNon<numberNonLinearColumns;jNon++) {
332        iColumn=listNonLinearColumn[jNon];
333        double dj = objective[iColumn]-r[iColumn];
334        r[iColumn]=dj;
335        if (dj<0.0) 
336          targetDrop -= dj*(columnUpper[iColumn]-solution[iColumn]);
337        else
338          targetDrop -= dj*(columnLower[iColumn]-solution[iColumn]);
339      }
340    } else {
341      memset(r,0,numberColumns*sizeof(double));
342    }
343#if 0
344    for (jNon=0;jNon<numberNonLinearColumns;jNon++) {
345      iColumn=listNonLinearColumn[jNon];
346      if (statusCheck[iColumn]=='L'&&r[iColumn]<-1.0e-4) {
347        columnLower[iColumn]=max(solution[iColumn],
348                                 trueLower[jNon]);
349        columnUpper[iColumn]=min(solution[iColumn]
350                                 +trust[jNon],
351                                 trueUpper[jNon]);
352      } else if (statusCheck[iColumn]=='U'&&r[iColumn]>1.0e-4) {
353        columnLower[iColumn]=max(solution[iColumn]
354                                 -trust[jNon],
355                                 trueLower[jNon]);
356        columnUpper[iColumn]=min(solution[iColumn],
357                                 trueUpper[jNon]);
358      } else {
359        columnLower[iColumn]=max(solution[iColumn]
360                                 -trust[jNon],
361                                 trueLower[jNon]);
362        columnUpper[iColumn]=min(solution[iColumn]
363                                 +trust[jNon],
364                                 trueUpper[jNon]);
365      }
366    }
367#endif
368    if (goodMove) {
369      memcpy(saveSolution,solution,numberColumns*sizeof(double));
370      memcpy(savePi,this->dualRowSolution(),numberRows*sizeof(double));
371      memcpy(saveStatus,status_,numberRows+numberColumns);
372     
373      std::cout<<"Pass - "<<iPass
374               <<", target drop is "<<targetDrop
375               <<std::endl;
376      lastObjective = objValue;
377      if (targetDrop<1.0e-5&&goodMove&&iPass) {
378        printf("Exiting on target drop %g\n",targetDrop);
379        break;
380      }
381      {
382        double * r = this->dualColumnSolution();
383        for (jNon=0;jNon<numberNonLinearColumns;jNon++) {
384          iColumn=listNonLinearColumn[jNon];
385          printf("Trust %d %g - solution %d %g obj %g dj %g state %c - bounds %g %g\n",
386                 jNon,trust[jNon],iColumn,solution[iColumn],objective[iColumn],
387                 r[iColumn],statusCheck[iColumn],columnLower[iColumn],
388                 columnUpper[iColumn]);
389        }
390      }
391      setLogLevel(63);
392      this->scaling(false);
393      this->primal(1);
394      if (this->status()) {
395        CoinMpsIO writer;
396        writer.setMpsData(*matrix(), COIN_DBL_MAX,
397                          getColLower(), getColUpper(),
398                          getObjCoefficients(),
399                          (const char*) 0 /*integrality*/,
400                          getRowLower(), getRowUpper(),
401                          NULL,NULL);
402        writer.writeMps("xx.mps");
403      }
404      assert (!this->status()); // must be feasible
405      goodMove=1;
406    } else {
407      // bad pass - restore solution
408      printf("Backtracking\n");
409      memcpy(solution,saveSolution,numberColumns*sizeof(double));
410      memcpy(this->dualRowSolution(),savePi,numberRows*sizeof(double));
411      memcpy(status_,saveStatus,numberRows+numberColumns);
412      iPass--;
413      goodMove=-1;
414    }
415  }
416  // restore solution
417  memcpy(solution,saveSolution,numberColumns*sizeof(double));
418  for (jNon=0;jNon<numberNonLinearColumns;jNon++) {
419    iColumn=listNonLinearColumn[jNon];
420    columnLower[iColumn]=max(solution[iColumn],
421                             trueLower[jNon]);
422    columnUpper[iColumn]=min(solution[iColumn],
423                             trueUpper[jNon]);
424  }
425  delete [] statusCheck;
426  delete [] savePi;
427  delete [] saveStatus;
428  // redo objective
429  double offset=0.0;
430  double objValue=-objectiveOffset;
431  memcpy(objective,saveObjective,numberColumns*sizeof(double));
432  for (iColumn=0;iColumn<numberColumns;iColumn++) 
433    objValue += objective[iColumn]*solution[iColumn];
434  for (jNon=0;jNon<numberNonLinearColumns;jNon++) {
435    iColumn=listNonLinearColumn[jNon];
436    double valueI = solution[iColumn];
437    int j;
438    for (j=columnQuadraticStart[iColumn];
439         j<columnQuadraticStart[iColumn]+columnQuadraticLength[iColumn];j++) {
440      int jColumn = columnQuadratic[j];
441      double valueJ = solution[jColumn];
442      double elementValue = quadraticElement[j];
443      objValue += 0.5*valueI*valueJ*elementValue;
444      offset += 0.5*valueI*valueJ*elementValue;
445      double gradientI = valueJ*elementValue;
446      double gradientJ = valueI*elementValue;
447      offset -= gradientI*valueI;
448      objective[iColumn] += gradientI;
449      offset -= gradientJ*valueJ;
450      objective[jColumn] += gradientJ;
451    }
452  }
453  printf("objective %g, objective offset %g\n",objValue,offset);
454  setDblParam(ClpObjOffset,objectiveOffset-offset);
455  this->primal(1);
456  // redo values
457  setDblParam(ClpObjOffset,objectiveOffset);
458  objectiveValue_ += optimizationDirection_*offset;
459  for (jNon=0;jNon<numberNonLinearColumns;jNon++) {
460    iColumn=listNonLinearColumn[jNon];
461    columnLower[iColumn]= trueLower[jNon];
462    columnUpper[iColumn]= trueUpper[jNon];
463  }
464  memcpy(objective,saveObjective,numberColumns*sizeof(double));
465  delete [] saveSolution;
466  for (iPass=0;iPass<3;iPass++) 
467    delete [] last[iPass];
468  delete [] trust;
469  delete [] trueUpper;
470  delete [] trueLower;
471  delete [] saveObjective;
472  delete [] listNonLinearColumn;
473  return 0;
474}
475// Dantzig's method
476int 
477ClpSimplexPrimalQuadratic::primalQuadratic(int phase)
478{
479#if 0
480  // Dantzig's method looks easier - lets try that
481
482  int numberColumns = this->numberColumns();
483  double * columnLower = this->columnLower();
484  double * columnUpper = this->columnUpper();
485  double * objective = this->objective();
486  double * solution = this->primalColumnSolution();
487  double * dj = this->dualColumnSolution();
488  double * pi = this->dualRowSolution();
489
490  int numberRows = this->numberRows();
491  double * rowLower = this->rowLower();
492  double * rowUpper = this->rowUpper();
493
494  // and elements
495  CoinPackedMatrix * matrix = this->matrix();
496  const int * row = matrix->getIndices();
497  const int * columnStart = matrix->getVectorStarts();
498  const double * element =  matrix->getElements();
499  const int * columnLength = matrix->getVectorLengths();
500
501  // Get list of non linear columns
502  CoinPackedMatrix * quadratic = quadraticObjective();
503  if (!quadratic||!quadratic->getNumElements()) {
504    // no quadratic part
505    return primal(1);
506  }
507
508  int iColumn;
509  const int * columnQuadratic = quadratic->getIndices();
510  const int * columnQuadraticStart = quadratic->getVectorStarts();
511  const int * columnQuadraticLength = quadratic->getVectorLengths();
512  const double * quadraticElement = quadratic->getElements();
513#if 0
514  // deliberate bad solution
515  // Change to use phase
516  //double * saveO = new double[numberColumns];
517  //memcpy(saveO,objective,numberColumns*sizeof(double));
518  //memset(objective,0,numberColumns*sizeof(double));
519  tempMessage messageHandler(this);;
520  passInMessageHandler(&messageHandler);
521  factorization()->maximumPivots(1);
522  primal();
523  CoinMessageHandler handler2;
524  passInMessageHandler(&handler2);
525  factorization()->maximumPivots(100);
526  setMaximumIterations(1000);
527#endif
528  //memcpy(objective,saveO,numberColumns*sizeof(double));
529  //printf("For testing - deliberate bad solution\n");
530  //columnUpper[0]=0.0;
531  //columnLower[0]=0.0;
532  //quadraticSLP(50,1.0e-4);
533  //primal(1);
534  //columnUpper[0]=COIN_DBL_MAX;
535 
536  // Create larger problem
537  // First workout how many rows extra
538  ClpQuadraticInfo info(this);
539  int numberQuadratic = info.numberQuadraticColumns();
540  int newNumberRows = numberRows+numberQuadratic;
541  int newNumberColumns = numberColumns + numberRows + numberQuadratic;
542  int numberElements = 2*matrix->getNumElements()
543    +2*quadratic->getNumElements()
544    + numberQuadratic;
545  double * elements2 = new double[numberElements];
546  int * start2 = new int[newNumberColumns+1];
547  int * row2 = new int[numberElements];
548  double * objective2 = new double[newNumberColumns];
549  double * columnLower2 = new double[newNumberColumns];
550  double * columnUpper2 = new double[newNumberColumns];
551  double * rowLower2 = new double[newNumberRows];
552  double * rowUpper2 = new double[newNumberRows];
553  const int * which = info.quadraticSequence();
554  const int * back = info.backSequence();
555  memcpy(rowLower2,rowLower,numberRows*sizeof(double));
556  memcpy(rowUpper2,rowUpper,numberRows*sizeof(double));
557  int iRow;
558  for (iRow=0;iRow<numberQuadratic;iRow++) {
559    double cost = objective[back[iRow]];
560    rowLower2[iRow+numberRows]=cost;
561    rowUpper2[iRow+numberRows]=cost;
562  }
563  memset(objective2,0,newNumberColumns*sizeof(double));
564  // Get a row copy of quadratic objective in standard format
565  CoinPackedMatrix copyQ;
566  copyQ.reverseOrderedCopyOf(*quadratic);
567  const int * columnQ = copyQ.getIndices();
568  const CoinBigIndex * rowStartQ = copyQ.getVectorStarts();
569  const int * rowLengthQ = copyQ.getVectorLengths();
570  const double * elementByRowQ = copyQ.getElements();
571  // Move solution across
572  double * solution2 = new double[newNumberColumns];
573  memset(solution2,0,newNumberColumns*sizeof(double));
574  newNumberColumns=0;
575  numberElements=0;
576  start2[0]=0;
577  // x
578  memcpy(dj,objective,numberColumns*sizeof(double));
579  for (iColumn=0;iColumn<numberColumns;iColumn++) {
580    // Original matrix
581    columnLower2[iColumn]=columnLower[iColumn];
582    columnUpper2[iColumn]=columnUpper[iColumn];
583    solution2[iColumn]=solution[iColumn];
584    int j;
585    for (j=columnStart[iColumn];
586         j<columnStart[iColumn]+columnLength[iColumn];
587         j++) {
588      elements2[numberElements]=element[j];
589      row2[numberElements++]=row[j];
590    }
591    if (which[iColumn]>=0) {
592      // Quadratic and modify djs
593      for (j=columnQuadraticStart[iColumn];
594           j<columnQuadraticStart[iColumn]+columnQuadraticLength[iColumn];
595           j++) {
596        int jColumn = columnQuadratic[j];
597        double value = quadraticElement[j];
598        if (iColumn!=jColumn)
599          value *= 0.5;
600        dj[jColumn] += solution[iColumn]*value;
601        elements2[numberElements]=-value;
602        row2[numberElements++]=which[jColumn]+numberRows;
603      }
604      for (j=rowStartQ[iColumn];
605           j<rowStartQ[iColumn]+rowLengthQ[iColumn];
606           j++) {
607        int jColumn = columnQ[j];
608        double value = elementByRowQ[j];
609        if (iColumn!=jColumn) {
610          value *= 0.5;
611          dj[jColumn] += solution[iColumn]*value;
612          elements2[numberElements]=-value;
613          row2[numberElements++]=which[jColumn]+numberRows;
614        }
615      }
616    }
617    start2[iColumn+1]=numberElements;
618  }
619  newNumberColumns=numberColumns;
620  // pi
621  // Get a row copy in standard format
622  CoinPackedMatrix copy;
623  copy.reverseOrderedCopyOf(*(this->matrix()));
624  // get matrix data pointers
625  const int * column = copy.getIndices();
626  const CoinBigIndex * rowStart = copy.getVectorStarts();
627  const int * rowLength = copy.getVectorLengths();
628  const double * elementByRow = copy.getElements();
629  for (iRow=0;iRow<numberRows;iRow++) {
630    solution2[newNumberColumns]=pi[iRow];
631    double value = pi[iRow];
632    columnLower2[newNumberColumns]=-COIN_DBL_MAX;
633    columnUpper2[newNumberColumns]=COIN_DBL_MAX;
634    int j;
635    for (j=rowStart[iRow];
636         j<rowStart[iRow]+rowLength[iRow];
637         j++) {
638      double elementValue=elementByRow[j];
639      int jColumn = column[j];
640      elements2[numberElements]=elementValue;
641      row2[numberElements++]=jColumn+numberRows;
642      dj[jColumn]-= value*elementValue;
643    }
644    newNumberColumns++;
645    start2[newNumberColumns]=numberElements;
646  }
647  // djs
648  for (iColumn=0;iColumn<numberQuadratic;iColumn++) {
649    columnLower2[newNumberColumns]=-COIN_DBL_MAX;
650    columnUpper2[newNumberColumns]=COIN_DBL_MAX;
651    solution2[newNumberColumns]=dj[iColumn];
652    elements2[numberElements]=1.0;
653    row2[numberElements++]=back[iColumn]+numberRows;
654    newNumberColumns++;
655    start2[newNumberColumns]=numberElements;
656  }
657  // Create model
658  ClpSimplex model2(*this);
659  model2.resize(0,0);
660  model2.loadProblem(newNumberColumns,newNumberRows,
661                     start2,row2, elements2,
662                     columnLower2,columnUpper2,
663                     objective2,
664                     rowLower2,rowUpper2);
665  delete [] objective2;
666  delete [] rowLower2;
667  delete [] rowUpper2;
668  delete [] columnLower2;
669  delete [] columnUpper2;
670  // Now create expanded quadratic objective for use in primalRow
671  // Later pack down in some way
672  start2[0]=0;
673  numberElements=0;
674  for (iColumn=0;iColumn<numberColumns;iColumn++) {
675    // Quadratic
676    int j;
677    for (j=columnQuadraticStart[iColumn];
678         j<columnQuadraticStart[iColumn]+columnQuadraticLength[iColumn];
679         j++) {
680      int jColumn = columnQuadratic[j];
681      double value = quadraticElement[j];
682      if (iColumn!=jColumn)
683        value *= 0.5;
684      elements2[numberElements]=value;
685      row2[numberElements++]=jColumn;
686    }
687    for (j=rowStartQ[iColumn];
688         j<rowStartQ[iColumn]+rowLengthQ[iColumn];
689         j++) {
690      int jColumn = columnQ[j];
691      double value = elementByRowQ[j];
692      if (iColumn!=jColumn) {
693        value *= 0.5;
694        elements2[numberElements]=value;
695        row2[numberElements++]=jColumn;
696      }
697    }
698    start2[iColumn+1]=numberElements;
699  }
700  // and pad
701  for (;iColumn<newNumberColumns;iColumn++)
702    start2[iColumn+1]=numberElements;
703  // Load up objective
704  model2.loadQuadraticObjective(newNumberColumns,start2,row2,elements2);
705  delete [] start2;
706  delete [] row2;
707  delete [] elements2;
708  model2.allSlackBasis();
709  model2.scaling(false);
710  model2.setLogLevel(this->logLevel());
711  // Move solution across
712  memcpy(model2.primalColumnSolution(),solution2,
713         newNumberColumns*sizeof(double));
714  columnLower2 = model2.columnLower();
715  columnUpper2 = model2.columnUpper();
716  delete [] solution2;
717  solution2 = model2.primalColumnSolution();
718  // Compute row activities and check feasible
719  double * rowSolution2 = model2.primalRowSolution();
720  memset(rowSolution2,0,newNumberRows*sizeof(double));
721  model2.times(1.0,solution2,rowSolution2);
722  rowLower2 = model2.rowLower();
723  rowUpper2 = model2.rowUpper();
724#if 0
725  for (iColumn=0;iColumn<numberColumns;iColumn++) {
726    Status xStatus = getColumnStatus(iColumn);
727    bool isSuperBasic;
728    int iS = iColumn+newNumberRows;
729    double value = solution2[iS];
730    if (fabs(value)>dualTolerance_)
731      isSuperBasic=true;
732    else
733      isSuperBasic=false;
734    // For moment take all x out of basis
735    // Does not seem right
736    isSuperBasic=true;
737    model2.setColumnStatus(iColumn,xStatus);
738    if (xStatus==basic) {
739      if (!isSuperBasic) {
740        model2.setRowStatus(numberRows+iColumn,basic);
741        model2.setColumnStatus(iS,superBasic);
742      } else {
743        model2.setRowStatus(numberRows+iColumn,isFixed);
744        model2.setColumnStatus(iS,basic);
745        model2.setColumnStatus(iColumn,superBasic);
746      }
747    } else {
748      model2.setRowStatus(numberRows+iColumn,isFixed);
749      model2.setColumnStatus(iS,basic);
750    }
751  }
752  for (iRow=0;iRow<numberRows;iRow++) {
753    Status rowStatus = getRowStatus(iRow);
754    model2.setRowStatus(iRow,rowStatus);
755    if (rowStatus!=basic) {
756      model2.setColumnStatus(iRow+numberColumns,basic); // make dual basic
757    }
758    assert (rowSolution2[iRow]>=rowLower2[iRow]-primalTolerance_);
759    assert (rowSolution2[iRow]<=rowUpper2[iRow]+primalTolerance_);
760  }
761  // why ?? - take duals out and adjust
762  for (iRow=0;iRow<numberRows;iRow++) {
763    model2.setRowStatus(iRow,basic);
764    model2.setColumnStatus(iRow+numberColumns,superBasic);
765    solution2[iRow+numberColumns]=0.0;
766  }
767#else
768  for (iRow=0;iRow<numberRows;iRow++) {
769    assert (rowSolution2[iRow]>=rowLower2[iRow]-primalTolerance_);
770    assert (rowSolution2[iRow]<=rowUpper2[iRow]+primalTolerance_);
771    model2.setRowStatus(iRow,basic);
772    model2.setColumnStatus(iRow+numberColumns,superBasic);
773    solution2[iRow+numberColumns]=0.0;
774  }
775  for (iColumn=numberRows+numberColumns;iColumn<newNumberColumns;iColumn++) {
776    model2.setColumnStatus(iColumn,basic);
777    model2.setRowStatus(iColumn-numberColumns,isFixed);
778  }
779#endif
780  memset(rowSolution2,0,newNumberRows*sizeof(double));
781  model2.times(1.0,solution2,rowSolution2);
782  for (iColumn=0;iColumn<numberQuadratic;iColumn++) {
783    int iS = back[iColumn]+newNumberRows;
784    int iRow = iColumn+numberRows;
785    double value = rowSolution2[iRow];
786    if (value>rowUpper2[iRow]) {
787      rowSolution2[iRow] = rowUpper2[iRow];
788      solution2[iS]-=value-rowUpper2[iRow];
789    } else {
790      rowSolution2[iRow] = rowLower2[iRow];
791      solution2[iS]-=value-rowLower2[iRow];
792    }
793  }
794
795 
796  // See if any s sub j have wrong sign and/or use djs from infeasibility objective
797  double objectiveOffset;
798  getDblParam(ClpObjOffset,objectiveOffset);
799  double objValue = -objectiveOffset;
800  for (iColumn=0;iColumn<numberColumns;iColumn++)
801    objValue += objective[iColumn]*solution2[iColumn];
802  for (iColumn=0;iColumn<numberColumns;iColumn++) {
803    double valueI = solution2[iColumn];
804    int j;
805    for (j=columnQuadraticStart[iColumn];
806         j<columnQuadraticStart[iColumn]+columnQuadraticLength[iColumn];j++) {
807      int jColumn = columnQuadratic[j];
808      double valueJ = solution2[jColumn];
809      double elementValue = quadraticElement[j];
810      objValue += 0.5*valueI*valueJ*elementValue;
811    }
812  }
813  printf("Objective value %g\n",objValue);
814  for (iColumn=0;iColumn<newNumberColumns;iColumn++)
815    printf("%d %g\n",iColumn,solution2[iColumn]);
816#if 1
817  CoinMpsIO writer;
818  writer.setMpsData(*model2.matrix(), COIN_DBL_MAX,
819                    model2.getColLower(), model2.getColUpper(),
820                    model2.getObjCoefficients(),
821                    (const char*) 0 /*integrality*/,
822                    model2.getRowLower(), model2.getRowUpper(),
823                    NULL,NULL);
824  writer.writeMps("xx.mps");
825#endif 
826  // Now do quadratic
827  // If we did not do Slp we should have primal feasible basic solution
828  // Do safe cast as no data
829  ClpSimplexPrimalQuadratic * modelPtr =
830    (ClpSimplexPrimalQuadratic *) (&model2);
831  ClpPrimalQuadraticDantzig dantzigP(modelPtr,numberRows);
832  modelPtr->setPrimalColumnPivotAlgorithm(dantzigP);
833  modelPtr->messageHandler()->setLogLevel(63);
834  modelPtr->primalQuadratic2(this,phase);
835  memcpy(dualRowSolution(),model2.primalColumnSolution()+numberColumns_,numberRows_*sizeof(double));
836  memcpy(primalColumnSolution(),model2.primalColumnSolution(),numberColumns_*sizeof(double));
837  memset(model2.primalRowSolution(),0,newNumberRows*sizeof(double));
838  model2.times(1.0,model2.primalColumnSolution(),
839               model2.primalRowSolution());
840  memcpy(dualColumnSolution(),model2.primalColumnSolution()+numberRows_+numberColumns_,
841         numberColumns_*sizeof(double));
842  objValue = -objectiveOffset;
843  for (iColumn=0;iColumn<numberColumns;iColumn++)
844    objValue += objective[iColumn]*solution2[iColumn];
845  for (iColumn=0;iColumn<numberColumns;iColumn++) {
846    double valueI = solution2[iColumn];
847    if (fabs(valueI)>1.0e-5) {
848      int djColumn = iColumn+numberRows+numberColumns;
849      assert(solution2[djColumn]<1.0e-7);
850    }
851    int j;
852    for (j=columnQuadraticStart[iColumn];
853         j<columnQuadraticStart[iColumn]+columnQuadraticLength[iColumn];j++) {
854      int jColumn = columnQuadratic[j];
855      double valueJ = solution2[jColumn];
856      double elementValue = quadraticElement[j];
857      objValue += 0.5*valueI*valueJ*elementValue;
858    }
859  }
860  printf("Objective value %g\n",objValue);
861  objectiveValue_ = objValue + objectiveOffset;
862  return 0;
863#else
864  // Get a feasible solution
865  if (numberPrimalInfeasibilities())
866    primal(1);
867  // still infeasible
868  if (numberPrimalInfeasibilities())
869    return 1;
870  ClpQuadraticInfo info;
871  ClpSimplexPrimalQuadratic * model2 = makeQuadratic(info);
872#if 0
873  CoinMpsIO writer;
874  writer.setMpsData(*model2->matrix(), COIN_DBL_MAX,
875                    model2->getColLower(), model2->getColUpper(),
876                    model2->getObjCoefficients(),
877                    (const char*) 0 /*integrality*/,
878                    model2->getRowLower(), model2->getRowUpper(),
879                    NULL,NULL);
880  writer.writeMps("xx.mps");
881#endif 
882  // Now do quadratic
883  ClpPrimalQuadraticDantzig dantzigP(model2,numberRows_);
884  model2->setPrimalColumnPivotAlgorithm(dantzigP);
885  model2->messageHandler()->setLogLevel(63);
886  model2->primalQuadratic2(this,phase);
887  endQuadratic(model2,info);
888  return 0;
889#endif
890}
891int ClpSimplexPrimalQuadratic::primalQuadratic2 (const ClpSimplexPrimalQuadratic * originalModel,int phase)
892{
893
894  algorithm_ = +2;
895
896  // save data
897  ClpDataSave data = saveData();
898 
899  // Assume problem is feasible
900  // Stuff below will be moved into a class
901  int numberXColumns = originalModel->numberColumns();
902  int numberXRows = originalModel->numberRows();
903  int baseS = numberXColumns+numberXRows;
904  int * pair = new int [numberColumns_];
905  int iColumn;
906  for (iColumn=0;iColumn<numberXColumns;iColumn++) {
907    int iS = iColumn+baseS;
908    pair[iColumn]=iS;
909    pair[iS]=iColumn;
910  }
911#if 0
912  // Throw out all x variables whose dj is nonzero
913  for (iColumn=0;iColumn<numberXColumns;iColumn++) {
914    if (getColumnStatus(iColumn)==basic) {
915      int jColumn = iColumn+numberRows_;
916      if (getColumnStatus(jColumn)==basic)
917        setColumnStatus(iColumn,superBasic);
918    }
919  } 
920  int originalNumberRows = originalModel->numberRows();
921  for (int i=0;i<originalNumberRows;i++) {
922    int jSequence = i+numberXColumns;
923    if (getRowStatus(i)==basic)
924      setColumnStatus(jSequence,superBasic);
925  }
926#endif
927  // Save solution
928  double * saveSolution = new double [numberRows_+numberColumns_];
929  assert (!scalingFlag_);
930  memcpy(saveSolution,columnActivity_,numberColumns_*sizeof(double));
931  memcpy(saveSolution+numberColumns_,rowActivity_,numberRows_*sizeof(double));
932  // initialize - values pass coding and algorithm_ is +1
933  if (!startup(1)) {
934
935    // Setup useful stuff
936    ClpQuadraticInfo info(originalModel);
937    if (phase)
938      memcpy(solution_,saveSolution,
939             (numberRows_+numberColumns_)*sizeof(double));
940    int lastCleaned=0; // last time objective or bounds cleaned up
941    int sequenceIn=-1;
942    int crucialSj=-1;
943   
944    // special nonlinear cost
945    delete nonLinearCost_;
946    nonLinearCost_ = new ClpNonLinearCost(this,originalModel->numberColumns());
947    // Say no pivot has occurred (for steepest edge and updates)
948    pivotRow_=-2;
949   
950    // This says whether to restore things etc
951    int factorType=0;
952    /*
953      Status of problem:
954      0 - optimal
955      1 - infeasible
956      2 - unbounded
957      -1 - iterating
958      -2 - factorization wanted
959      -3 - redo checking without factorization
960      -4 - looks infeasible
961      -5 - looks unbounded
962    */
963    while (problemStatus_<0) {
964      int iRow,iColumn;
965      // clear
966      for (iRow=0;iRow<4;iRow++) {
967        rowArray_[iRow]->clear();
968      }   
969     
970      for (iColumn=0;iColumn<2;iColumn++) {
971        columnArray_[iColumn]->clear();
972      }   
973     
974      // give matrix (and model costs and bounds a chance to be
975      // refreshed (normally null)
976      matrix_->refresh(this);
977      // If we have done no iterations - special
978      if (lastGoodIteration_==numberIterations_)
979        factorType=3;
980      if (phase)
981        memcpy(saveSolution,solution_,
982               (numberRows_+numberColumns_)*sizeof(double));
983      if (phase&&0) {
984        // Clean solution
985        for (iColumn=0;iColumn<numberColumns_;iColumn++) {
986          if (getColumnStatus(iColumn)==isFree)
987            solution_[iColumn]=0.0;
988        }
989      } 
990      // may factorize, checks if problem finished
991      statusOfProblemInPrimal(lastCleaned,factorType,progress_);
992      if (phase)
993        memcpy(solution_,saveSolution,
994               (numberRows_+numberColumns_)*sizeof(double));
995     
996      // Compute objective function from scratch
997      const CoinPackedMatrix * quadratic = originalModel->quadraticObjective();
998      const int * columnQuadratic = quadratic->getIndices();
999      const int * columnQuadraticStart = quadratic->getVectorStarts();
1000      const int * columnQuadraticLength = quadratic->getVectorLengths();
1001      const double * quadraticElement = quadratic->getElements();
1002      const double * originalCost = originalModel->objective();
1003      objectiveValue_=0.0;
1004      for (iColumn=0;iColumn<numberXColumns;iColumn++) {
1005        double valueI = solution_[iColumn];
1006        objectiveValue_ += valueI*originalCost[iColumn];
1007        int j;
1008        for (j=columnQuadraticStart[iColumn];
1009             j<columnQuadraticStart[iColumn]+columnQuadraticLength[iColumn];j++) {
1010          int jColumn = columnQuadratic[j];
1011          double valueJ = solution_[jColumn];
1012          double elementValue = 0.5*quadraticElement[j];
1013          objectiveValue_ += valueI*valueJ*elementValue;
1014        }
1015      }
1016
1017      // Say good factorization
1018      factorType=1;
1019     
1020      // Say no pivot has occurred (for steepest edge and updates)
1021      pivotRow_=-2;
1022      // Check problem phase
1023      // We assume all X are feasible
1024      phase=0;
1025      if (saveSolution) {
1026        sequenceIn=-1;
1027        for (iColumn=0;iColumn<numberXColumns;iColumn++) {
1028          double value = solution_[iColumn];
1029          double lower = lower_[iColumn];
1030          double upper = upper_[iColumn];
1031          if (value>lower+primalTolerance_&&value<upper-primalTolerance_) {
1032            if (getColumnStatus(iColumn)!=basic) {
1033              if (phase!=2) {
1034                phase=2;
1035                sequenceIn=iColumn;
1036              }
1037            }
1038          }
1039          if (getColumnStatus(iColumn)==basic) {
1040            int iS = pair[iColumn];
1041            assert (iS>=0&&getColumnStatus(iS)!=basic);
1042            if(fabs(solution_[iS])>dualTolerance_) {
1043              if (phase==0) {
1044                phase=1;
1045                sequenceIn=iS;
1046              }
1047            }
1048          }
1049        }
1050        int offset=numberXColumns-numberColumns_;
1051        for (iColumn=numberColumns_;iColumn<numberColumns_+numberXRows;iColumn++) {
1052          double value = solution_[iColumn];
1053          double lower = lower_[iColumn];
1054          double upper = upper_[iColumn];
1055          if (value>lower+primalTolerance_&&value<upper-primalTolerance_) {
1056            if (getColumnStatus(iColumn)!=basic) {
1057              if (phase!=2) {
1058                phase=2;
1059                sequenceIn=iColumn;
1060              }
1061            }
1062          }
1063          if (getColumnStatus(iColumn)==basic) {
1064            int iS = iColumn+offset;
1065            assert (getColumnStatus(iS)!=basic);
1066            if(fabs(solution_[iS])>dualTolerance_) {
1067              if (phase==0) {
1068                phase=1;
1069                sequenceIn=iS;
1070              }
1071            }
1072          }
1073        }
1074      }
1075      if (!phase) {
1076        delete [] saveSolution;
1077        saveSolution=NULL;
1078      }
1079
1080      // exit if victory declared
1081      if (!phase&&primalColumnPivot_->pivotColumn(rowArray_[1],
1082                                          rowArray_[2],rowArray_[3],
1083                                          columnArray_[0],
1084                                          columnArray_[1]) < 0) {
1085        problemStatus_=0;
1086        break;
1087      }
1088     
1089      // Iterate
1090      problemStatus_=-1;
1091      whileIterating(originalModel,sequenceIn,&info,crucialSj,phase);
1092    }
1093  }
1094  // clean up
1095  delete [] saveSolution;
1096  delete [] pair;
1097  finish();
1098  restoreData(data);
1099  return problemStatus_;
1100}
1101/*
1102  Reasons to come out:
1103  -1 iterations etc
1104  -2 inaccuracy
1105  -3 slight inaccuracy (and done iterations)
1106  -4 end of values pass and done iterations
1107  +0 looks optimal (might be infeasible - but we will investigate)
1108  +2 looks unbounded
1109  +3 max iterations
1110*/
1111int
1112ClpSimplexPrimalQuadratic::whileIterating(
1113                      const ClpSimplexPrimalQuadratic * originalModel,
1114                      int & sequenceIn,
1115                      ClpQuadraticInfo * info,
1116                      int & crucialSj,
1117                      int phase)
1118{
1119  checkComplementary (info);
1120
1121  int returnCode=-1;
1122  double saveObjective = objectiveValue_;
1123  int numberXColumns = originalModel->numberColumns();
1124  int oldSequenceIn=sequenceIn;
1125  // status stays at -1 while iterating, >=0 finished, -2 to invert
1126  // status -3 to go to top without an invert
1127  while (problemStatus_==-1) {
1128#ifdef CLP_DEBUG
1129    {
1130      int i;
1131      // not [1] as has information
1132      for (i=0;i<4;i++) {
1133        if (i!=1)
1134          rowArray_[i]->checkClear();
1135      }   
1136      for (i=0;i<2;i++) {
1137        columnArray_[i]->checkClear();
1138      }   
1139    }     
1140#endif
1141    // choose column to come in
1142    // can use pivotRow_ to update weights
1143    // pass in list of cost changes so can do row updates (rowArray_[1])
1144    // NOTE rowArray_[0] is used by computeDuals which is a
1145    // slow way of getting duals but might be used
1146    // Initially Dantzig and look at s variables
1147    // Only do if one not already chosen
1148    bool cleanupIteration;
1149    if (phase==2) {
1150      // values pass
1151      if (sequenceIn<0) {
1152        // get next
1153        int iColumn;
1154        int iStart = oldSequenceIn+1;
1155        for (iColumn=iStart;iColumn<numberXColumns;iColumn++) {
1156          double value = solution_[iColumn];
1157          double lower = lower_[iColumn];
1158          double upper = upper_[iColumn];
1159          if (value>lower+primalTolerance_&&value<upper-primalTolerance_) {
1160            if (getColumnStatus(iColumn)!=basic) {
1161              sequenceIn=iColumn;
1162              break;
1163            }
1164          }
1165        }
1166        if (sequenceIn<0) {
1167          iStart=max(iStart,numberColumns_);
1168          int numberXRows = originalModel->numberRows();
1169          for (iColumn=numberColumns_;iColumn<numberColumns_+numberXRows;
1170               iColumn++) {
1171            double value = solution_[iColumn];
1172            double lower = lower_[iColumn];
1173            double upper = upper_[iColumn];
1174            if (value>lower+primalTolerance_&&value<upper-primalTolerance_) {
1175              if (getColumnStatus(iColumn)!=basic) {
1176                sequenceIn=iColumn;
1177                break;
1178              }
1179            }
1180          }
1181        }
1182      }
1183      oldSequenceIn=sequenceIn;
1184      sequenceIn_ = sequenceIn;
1185      cleanupIteration=false;
1186      dualIn_ = solution_[sequenceIn_+numberRows_]; 
1187      valueIn_=solution_[sequenceIn_];
1188      if (dualIn_>0.0)
1189        directionIn_ = -1;
1190      else 
1191        directionIn_ = 1;
1192    } else {
1193      if (sequenceIn<0) {
1194        primalColumn(rowArray_[1],rowArray_[2],rowArray_[3],
1195                     columnArray_[0],columnArray_[1]);
1196        cleanupIteration=false;
1197      } else {
1198        sequenceIn_ = sequenceIn;
1199        cleanupIteration=true;
1200      }
1201    }
1202    pivotRow_=-1;
1203    sequenceOut_=-1;
1204    rowArray_[1]->clear();
1205    if (sequenceIn_>=0) {
1206      sequenceIn=-1;
1207      // we found a pivot column
1208      int chosen = sequenceIn_;
1209      // do second half of iteration
1210      while (chosen>=0) {
1211        objectiveValue_=saveObjective;
1212        returnCode=-1;
1213        pivotRow_=-1;
1214        sequenceOut_=-1;
1215        rowArray_[1]->clear();
1216        // we found a pivot column
1217        // update the incoming column
1218        sequenceIn_=chosen;
1219        chosen=-1;
1220        unpack(rowArray_[1]);
1221        factorization_->updateColumnFT(rowArray_[2],rowArray_[1]);
1222        if (cleanupIteration) {
1223          // move back to a version of primalColumn?
1224          valueIn_=solution_[sequenceIn_];
1225          // should keep pivot row of crucialSj as well (speed)
1226          int iSjRow=-1;
1227          {
1228            double * work=rowArray_[1]->denseVector();
1229            int number=rowArray_[1]->getNumElements();
1230            int * which=rowArray_[1]->getIndices();
1231            double tj = 0.0;
1232            int iIndex;
1233            for (iIndex=0;iIndex<number;iIndex++) {
1234              int iRow = which[iIndex];
1235              double alpha = work[iRow];
1236              int iPivot=pivotVariable_[iRow];
1237              if (iPivot==crucialSj) {
1238                tj = alpha;
1239                iSjRow = iRow;
1240                double d2 = solution_[crucialSj]/tj;
1241                // see which way to go
1242                if (d2>0)
1243                  dj_[sequenceIn_]= -1.0;
1244                else
1245                  dj_[sequenceIn_]= 1.0;
1246                break;
1247              }
1248            }
1249            if (!tj) {
1250              printf("trouble\n");
1251              assert (sequenceIn_>numberXColumns&&sequenceIn_<numberColumns_);
1252              dj_[sequenceIn_]=solution_[sequenceIn_];
1253            //assert(tj);
1254            }
1255          }
1256
1257          dualIn_=dj_[sequenceIn_];
1258          if (sequenceIn_>numberXColumns&&
1259              sequenceIn_<-numberColumns_) {
1260            // We can let flip to 0.0
1261            assert (cleanupIteration);
1262            if (solution_[sequenceIn_]>0.0) {
1263              nonLinearCost_->setBounds(sequenceIn_, 0.0,
1264                                        0.5*COIN_DBL_MAX);
1265              lower_[sequenceIn_]=0.0;
1266              upper_[sequenceIn_]=0.5*COIN_DBL_MAX;
1267            } else {
1268              nonLinearCost_->setBounds(sequenceIn_, -0.5*COIN_DBL_MAX,
1269                                        0.0);
1270              lower_[sequenceIn_]=-0.5*COIN_DBL_MAX;
1271              upper_[sequenceIn_]=0.0;
1272            }
1273          } else {
1274            // Let go through bounds
1275            nonLinearCost_->setBounds(sequenceIn_, -0.5*COIN_DBL_MAX,
1276                                      0.5*COIN_DBL_MAX);
1277            lower_[sequenceIn_]=-0.5*COIN_DBL_MAX;
1278            upper_[sequenceIn_]=0.5*COIN_DBL_MAX;
1279          }
1280          lowerIn_=lower_[sequenceIn_];
1281          upperIn_=upper_[sequenceIn_];
1282          if (dualIn_>0.0)
1283            directionIn_ = -1;
1284          else 
1285            directionIn_ = 1;
1286        } else {
1287          if (sequenceIn_<numberColumns_) {
1288            // Set dj as value of slack
1289            crucialSj = sequenceIn_+numberRows_; // sj which should go to 0.0
1290            //dualIn_=solution_[crucialSj];
1291          } else {
1292            // Set dj as value of pi
1293            crucialSj = sequenceIn_-numberColumns_+numberXColumns; // pi which should go to 0.0
1294            //dualIn_=solution_[crucialSj];
1295          }
1296        }
1297        // save reduced cost
1298        //double saveDj = dualIn_;
1299        // do ratio test and re-compute dj
1300        // Note second parameter long enough for columns
1301        int result=primalRow(rowArray_[1],rowArray_[3],rowArray_[2],rowArray_[0],
1302                             info,
1303                             originalModel,crucialSj,cleanupIteration);
1304        saveObjective = objectiveValue_;
1305        if (pivotRow_>=0) {
1306          // if stable replace in basis
1307          int updateStatus = factorization_->replaceColumn(rowArray_[2],
1308                                                           pivotRow_,
1309                                                           alpha_);
1310          // if no pivots, bad update but reasonable alpha - take and invert
1311          if (updateStatus==2&&
1312              lastGoodIteration_==numberIterations_&&fabs(alpha_)>1.0e-5)
1313            updateStatus=4;
1314          if (updateStatus==1||updateStatus==4) {
1315            // slight error
1316            if (factorization_->pivots()>5||updateStatus==4) {
1317              returnCode=-3;
1318            }
1319          } else if (updateStatus==2) {
1320            // major error
1321            // better to have small tolerance even if slower
1322            factorization_->zeroTolerance(1.0e-15);
1323            int maxFactor = factorization_->maximumPivots();
1324            if (maxFactor>10) {
1325              if (forceFactorization_<0)
1326                forceFactorization_= maxFactor;
1327              forceFactorization_ = max (1,(forceFactorization_>>1));
1328            } 
1329            // later we may need to unwind more e.g. fake bounds
1330            if(lastGoodIteration_ != numberIterations_) {
1331              rowArray_[1]->clear();
1332              pivotRow_=-1;
1333              returnCode=-4;
1334              // retry on return
1335              sequenceIn = sequenceIn_;
1336              break;
1337            } else {
1338              // need to reject something
1339              char x = isColumn(sequenceIn_) ? 'C' :'R';
1340              handler_->message(CLP_SIMPLEX_FLAG,messages_)
1341                <<x<<sequenceWithin(sequenceIn_)
1342                <<CoinMessageEol;
1343              setFlagged(sequenceIn_);
1344              lastBadIteration_ = numberIterations_; // say be more cautious
1345              rowArray_[1]->clear();
1346              pivotRow_=-1;
1347              returnCode = -5;
1348              break;
1349             
1350            }
1351          } else if (updateStatus==3) {
1352            // out of memory
1353            // increase space if not many iterations
1354            if (factorization_->pivots()<
1355                0.5*factorization_->maximumPivots()&&
1356                factorization_->pivots()<200)
1357              factorization_->areaFactor(
1358                                         factorization_->areaFactor() * 1.1);
1359            returnCode =-2; // factorize now
1360          }
1361          // here do part of steepest - ready for next iteration
1362          primalColumnPivot_->updateWeights(rowArray_[1]);
1363        } else {
1364          if (pivotRow_==-1) {
1365            // no outgoing row is valid
1366            rowArray_[0]->clear();
1367            if (!factorization_->pivots()) {
1368              returnCode = 2; //say looks unbounded
1369              // do ray
1370              primalRay(rowArray_[1]);
1371            } else {
1372              returnCode = 4; //say looks unbounded but has iterated
1373            }
1374            break;
1375          } else {
1376            // flipping from bound to bound
1377          }
1378        }
1379
1380
1381        // update primal solution
1382
1383        double objectiveChange=0.0;
1384        // Cost on pivot row may change - may need to change dualIn
1385        double oldCost=0.0;
1386        if (pivotRow_>=0)
1387          oldCost = cost(pivotVariable_[pivotRow_]);
1388        // rowArray_[1] is not empty - used to update djs
1389        updatePrimalsInPrimal(rowArray_[1],theta_, objectiveChange);
1390        if (pivotRow_>=0)
1391          dualIn_ += (oldCost-cost(pivotVariable_[pivotRow_]));
1392        double oldValue = valueIn_;
1393        if (directionIn_==-1) {
1394          // as if from upper bound
1395          if (sequenceIn_!=sequenceOut_) {
1396            // variable becoming basic
1397            valueIn_ -= fabs(theta_);
1398          } else {
1399            valueIn_=lowerIn_;
1400          }
1401        } else {
1402          // as if from lower bound
1403          if (sequenceIn_!=sequenceOut_) {
1404            // variable becoming basic
1405            valueIn_ += fabs(theta_);
1406          } else {
1407            valueIn_=upperIn_;
1408          }
1409        }
1410        objectiveChange += dualIn_*(valueIn_-oldValue);
1411        // outgoing
1412        if (sequenceIn_!=sequenceOut_) {
1413          if (directionOut_>0) {
1414            valueOut_ = lowerOut_;
1415          } else {
1416            valueOut_ = upperOut_;
1417          }
1418          double lowerValue = lower_[sequenceOut_];
1419          double upperValue = upper_[sequenceOut_];
1420          assert(valueOut_>=lowerValue-primalTolerance_&&
1421                 valueOut_<=upperValue+primalTolerance_);
1422          // may not be exactly at bound and bounds may have changed
1423          if (valueOut_<=lowerValue+primalTolerance_) {
1424            directionOut_=1;
1425          } else if (valueOut_>=upperValue-primalTolerance_) {
1426            directionOut_=-1;
1427          } else {
1428            printf("*** variable wandered off bound %g %g %g!\n",
1429                   lowerValue,valueOut_,upperValue);
1430            if (valueOut_-lowerValue<=upperValue-valueOut_) {
1431              valueOut_=lowerValue;
1432              directionOut_=1;
1433            } else {
1434              valueOut_=upperValue;
1435              directionOut_=-1;
1436            }
1437          }
1438          solution_[sequenceOut_]=valueOut_;
1439          nonLinearCost_->setOne(sequenceOut_,valueOut_);
1440        }
1441        // change cost and bounds on incoming if primal
1442        nonLinearCost_->setOne(sequenceIn_,valueIn_); 
1443        int whatNext=housekeeping(objectiveChange);
1444        checkComplementary (info);
1445
1446        if (whatNext==1) {
1447          returnCode =-2; // refactorize
1448        } else if (whatNext==2) {
1449          // maximum iterations or equivalent
1450          returnCode=3;
1451        } else if(numberIterations_ == lastGoodIteration_
1452                  + 2 * factorization_->maximumPivots()) {
1453          // done a lot of flips - be safe
1454          returnCode =-2; // refactorize
1455        }
1456        // may need to go round again
1457        cleanupIteration = true;
1458        double newDj;
1459        // may not be correct on second time
1460        if (sequenceIn_<numberXColumns)
1461          newDj = solution_[sequenceIn_+numberRows_];
1462        else
1463          newDj = solution_[sequenceIn_-numberColumns_+numberXColumns];
1464        newDj=1.0; // force
1465        if (result&&fabs(newDj)>dualTolerance_) {
1466          assert(sequenceOut_<numberXColumns||
1467                 sequenceOut_>=numberColumns_);
1468          if (sequenceOut_<numberXColumns) {
1469            chosen =sequenceOut_ + numberRows_; // sj variable
1470          } else {
1471            // Does this mean we can change pi
1472            int iRow = sequenceOut_-numberColumns_;
1473            if (iRow<numberRows_-numberXColumns) {
1474              int iPi = iRow+numberXColumns;
1475              printf("pi for row %d is %g\n",
1476                     iRow,solution_[iPi]);
1477              chosen=iPi;
1478            } else {
1479              printf("Row %d is in column part\n",iRow);
1480              abort();
1481            }
1482          }
1483        } else {
1484          break;
1485        }
1486      }
1487      if (returnCode<-1&&returnCode>-5) {
1488        problemStatus_=-2; //
1489      } else if (returnCode==-5) {
1490        // something flagged - continue;
1491      } else if (returnCode==2) {
1492        problemStatus_=-5; // looks unbounded
1493      } else if (returnCode==4) {
1494        problemStatus_=-2; // looks unbounded but has iterated
1495      } else if (returnCode!=-1) {
1496        assert(returnCode==3);
1497        problemStatus_=3;
1498      }
1499    } else {
1500      // no pivot column
1501#ifdef CLP_DEBUG
1502      if (handler_->logLevel()&32)
1503        printf("** no column pivot\n");
1504#endif
1505      if (nonLinearCost_->numberInfeasibilities())
1506        problemStatus_=-4; // might be infeasible
1507      returnCode=0;
1508      break;
1509    }
1510  }
1511  return returnCode;
1512}
1513/*
1514   Row array has pivot column
1515   This chooses pivot row.
1516   For speed, we may need to go to a bucket approach when many
1517   variables go through bounds
1518   On exit rhsArray will have changes in costs of basic variables
1519*/
1520int 
1521ClpSimplexPrimalQuadratic::primalRow(CoinIndexedVector * rowArray,
1522                                     CoinIndexedVector * rhsArray,
1523                                     CoinIndexedVector * spareArray,
1524                                     CoinIndexedVector * spareArray2,
1525                                     ClpQuadraticInfo * info,
1526                                     const ClpSimplexPrimalQuadratic * originalModel,
1527                                     int crucialSj,
1528                                     bool cleanupIteration)
1529{
1530  int result=1;
1531  // sequence stays as row number until end
1532  pivotRow_=-1;
1533  int numberSwapped=0;
1534  int numberRemaining=0;
1535
1536  int numberThru =0; // number gone thru a barrier
1537  int lastThru =0; // number gone thru a barrier on last time
1538 
1539  double totalThru=0.0; // for when variables flip
1540  double acceptablePivot=1.0e-7;
1541  if (factorization_->pivots())
1542    acceptablePivot=1.0e-5; // if we have iterated be more strict
1543  double bestEverPivot=acceptablePivot;
1544  int lastPivotRow = -1;
1545  double lastPivot=0.0;
1546  double lastTheta=1.0e50;
1547  int lastNumberSwapped=0;
1548
1549  // use spareArrays to put ones looked at in
1550  // First one is list of candidates
1551  // We could compress if we really know we won't need any more
1552  // Second array has current set of pivot candidates
1553  // with a backup list saved in double * part of indexed vector
1554
1555  // for zeroing out arrays after
1556  int maximumSwapped=0;
1557  // pivot elements
1558  double * spare;
1559  // indices
1560  int * index, * indexSwapped;
1561  int * saveSwapped;
1562  spareArray->clear();
1563  spareArray2->clear();
1564  spare = spareArray->denseVector();
1565  index = spareArray->getIndices();
1566  saveSwapped = (int *) spareArray2->denseVector();
1567  indexSwapped = spareArray2->getIndices();
1568
1569  // we also need somewhere for effective rhs
1570  double * rhs=rhsArray->denseVector();
1571
1572  /*
1573    First we get a list of possible pivots.  We can also see if the
1574    problem looks unbounded.
1575
1576    At first we increase theta and see what happens.  We start
1577    theta at a reasonable guess.  If in right area then we do bit by bit.
1578    We save possible pivot candidates
1579
1580   */
1581
1582  // do first pass to get possibles
1583  // We can also see if unbounded
1584
1585  double * work=rowArray->denseVector();
1586  int number=rowArray->getNumElements();
1587  int * which=rowArray->getIndices();
1588
1589  // we need to swap sign if coming in from ub
1590  double way = directionIn_;
1591  double maximumMovement;
1592  if (way>0.0) 
1593    maximumMovement = min(1.0e30,upperIn_-valueIn_);
1594  else
1595    maximumMovement = min(1.0e30,valueIn_-lowerIn_);
1596
1597  int iIndex;
1598
1599  // Work out coefficients for quadratic term
1600  // This is expanded one
1601  const CoinPackedMatrix * quadratic = quadraticObjective();
1602  const int * columnQuadratic = quadratic->getIndices();
1603  const int * columnQuadraticStart = quadratic->getVectorStarts();
1604  const int * columnQuadraticLength = quadratic->getVectorLengths();
1605  const double * quadraticElement = quadratic->getElements();
1606  const double * originalCost = originalModel->objective();
1607  // Use rhsArray
1608  rhsArray->clear();
1609  int * index2 = rhsArray->getIndices();
1610  int numberXColumns = originalModel->numberColumns();
1611  int number2=0;
1612  //int numberOriginalRows = originalModel->numberRows();
1613  // sj
1614  int iSjRow=-1;
1615  double tj = 0.0;
1616  for (iIndex=0;iIndex<number;iIndex++) {
1617    int iRow = which[iIndex];
1618    double alpha = -work[iRow]*way;
1619    int iPivot=pivotVariable_[iRow];
1620    if (iPivot<numberXColumns) {
1621      index2[number2++]=iPivot;
1622      rhs[iPivot]=alpha;
1623      //printf("col %d alpha %g solution %g\n",iPivot,alpha,solution_[iPivot]);
1624    } else {
1625      //printf("new col %d alpha %g solution %g\n",iPivot,alpha,solution_[iPivot]);
1626      if (iPivot==crucialSj) {
1627        tj = alpha;
1628        iSjRow = iRow;
1629      }
1630    }
1631  }
1632  // Incoming
1633  if (sequenceIn_<numberXColumns) {
1634    index2[number2++]=sequenceIn_;
1635    rhs[sequenceIn_]=way;
1636    printf("incoming col %d alpha %g solution %g\n",sequenceIn_,way,valueIn_);
1637  } else {
1638    printf("incoming new col %d alpha %g solution %g\n",sequenceIn_,way,valueIn_);
1639  }
1640  rhsArray->setNumElements(number2);
1641  // Change in objective will be theta*coeff1 + theta*theta*coeff2
1642  double coeff1 = 0.0;
1643  double coeff2 = 0.0;
1644  if (numberIterations_>=0||cleanupIteration) {
1645    for (iIndex=0;iIndex<number2;iIndex++) {
1646      int iColumn=index2[iIndex];
1647      //double valueI = solution_[iColumn];
1648      double alphaI = rhs[iColumn];
1649      coeff1 += alphaI*originalCost[iColumn];
1650      int j;
1651      for (j=columnQuadraticStart[iColumn];
1652           j<columnQuadraticStart[iColumn]+columnQuadraticLength[iColumn];j++) {
1653        int jColumn = columnQuadratic[j];
1654        double valueJ = solution_[jColumn];
1655        double alphaJ = rhs[jColumn];
1656        double elementValue = quadraticElement[j];
1657        coeff1 += (valueJ*alphaI)*elementValue;
1658        coeff2 += (alphaI*alphaJ)*elementValue;
1659      }
1660    }
1661  } else {
1662    const int * row = matrix_->getIndices();
1663    const int * columnStart = matrix_->getVectorStarts();
1664    const int * columnLength = matrix_->getVectorLengths();
1665    const double * element = matrix_->getElements();
1666    int j;
1667    int jRow = sequenceIn_+originalModel->numberRows();
1668    printf("sequence in %d, cost %g\n",sequenceIn_,
1669           upper_[jRow+numberColumns_]);
1670    double xx=0.0;
1671    for (j=0;j<numberColumns_;j++) {
1672      int k;
1673      for (k=columnStart[j];k<columnStart[j]+columnLength[j];k++) {
1674        if (row[k]==jRow) {
1675          printf ("col %d, el %g, sol %g, contr %g\n",
1676                  j,element[k],solution_[j],element[k]*solution_[j]);
1677          xx+= element[k]*solution_[j];
1678        }
1679      }
1680    }
1681    printf("sum %g\n",xx);
1682    for (iIndex=0;iIndex<number2;iIndex++) {
1683      int iColumn=index2[iIndex];
1684      double valueI = solution_[iColumn];
1685      double alphaI = rhs[iColumn];
1686      //rhs[iColumn]=0.0;
1687      printf("Column %d, alpha %g sol %g cost %g, contr %g\n",
1688             iColumn,alphaI,valueI,originalCost[iColumn],
1689             alphaI*originalCost[iColumn]);
1690      coeff1 += alphaI*originalCost[iColumn];
1691      int j;
1692      for (j=columnQuadraticStart[iColumn];
1693           j<columnQuadraticStart[iColumn]+columnQuadraticLength[iColumn];j++) {
1694        int jColumn = columnQuadratic[j];
1695        double valueJ = solution_[jColumn];
1696        double alphaJ = rhs[jColumn];
1697        double elementValue = quadraticElement[j];
1698        printf("j %d alphaJ %g solJ %g el %g, contr %g\n",
1699               jColumn,alphaJ,valueJ,elementValue,
1700               (valueJ*alphaI)*elementValue);
1701        coeff1 += (valueJ*alphaI)*elementValue;
1702        coeff2 += (alphaI*alphaJ)*elementValue;
1703      }
1704    }
1705  }
1706  coeff2 *= 0.5;
1707  printf("coefficients %g %g - dualIn %g\n",coeff1,coeff2,dualIn_);
1708  if (!cleanupIteration) 
1709    assert (fabs(way*coeff1-dualIn_)<1.0e-4);
1710  // interesting places are where derivative zero or sj goes to zero
1711  double d1,d2=1.0e50;
1712  if (fabs(coeff2)>1.0e-9)
1713    d1 = - 0.5*coeff1/coeff2;
1714  else if (coeff1<=1.0e-9)
1715    d1 = maximumMovement;
1716  else
1717    d1 = 0.0;
1718  if (fabs(tj)<1.0e-7) {
1719    if (sequenceIn_<numberXColumns)
1720      printf("column %d is basically linear\n",sequenceIn_);
1721    //assert(!columnQuadraticLength[sequenceIn_]);
1722  } else {
1723    d2 = -solution_[crucialSj]/tj;
1724    if (d2<0.0) {
1725      printf("d2 would be negative at %g\n",d2);
1726      d2=1.0e50;
1727    }
1728  }
1729  printf("derivative zero at %g, sj zero at %g\n",d1,d2);
1730  if (d1>1.0e10&&d2>1.0e10) {
1731    // linear variable entering
1732    // maybe we should have done dual iteration to force sj to 0.0
1733    printf("linear variable\n");
1734  }
1735  maximumMovement = min(maximumMovement,d1);
1736  maximumMovement = min(maximumMovement,d2);
1737  d2 = min(d1,d2);
1738   
1739  rhsArray->clear();
1740  double tentativeTheta = maximumMovement;
1741  double upperTheta = maximumMovement;
1742
1743
1744  for (iIndex=0;iIndex<number;iIndex++) {
1745
1746    int iRow = which[iIndex];
1747    double alpha = work[iRow];
1748    int iPivot=pivotVariable_[iRow];
1749    alpha *= way;
1750    double oldValue = solution(iPivot);
1751    // get where in bound sequence
1752    if (alpha>0.0) {
1753      // basic variable going towards lower bound
1754      double bound = lower(iPivot);
1755      oldValue -= bound;
1756    } else if (alpha<0.0) {
1757      // basic variable going towards upper bound
1758      double bound = upper(iPivot);
1759      oldValue = bound-oldValue;
1760    }
1761    double value = oldValue-tentativeTheta*fabs(alpha);
1762    assert (oldValue>=-primalTolerance_*1.0001);
1763    if (value<-primalTolerance_) {
1764      // add to list
1765      spare[numberRemaining]=alpha;
1766      rhs[iRow]=oldValue;
1767      index[numberRemaining++]=iRow;
1768      double value=oldValue-upperTheta*fabs(alpha);
1769      if (value<-primalTolerance_&&fabs(alpha)>=acceptablePivot)
1770        upperTheta = (oldValue+primalTolerance_)/fabs(alpha);
1771    }
1772  }
1773
1774  // we need to keep where rhs non-zeros are
1775  int numberInRhs=numberRemaining;
1776  memcpy(rhsArray->getIndices(),index,numberInRhs*sizeof(int));
1777  rhsArray->setNumElements(numberInRhs);
1778
1779  theta_=maximumMovement;
1780
1781  bool goBackOne = false;
1782
1783  if (numberRemaining) {
1784
1785    // looks like pivoting
1786    // now try until reasonable theta
1787    tentativeTheta = max(10.0*upperTheta,1.0e-7);
1788    tentativeTheta = min(tentativeTheta,maximumMovement);
1789   
1790    // loops increasing tentative theta until can't go through
1791   
1792    while (tentativeTheta <= maximumMovement) {
1793      double thruThis = 0.0;
1794     
1795      double bestPivot=acceptablePivot;
1796      pivotRow_ = -1;
1797     
1798      numberSwapped = 0;
1799     
1800      upperTheta = maximumMovement;
1801     
1802      for (iIndex=0;iIndex<numberRemaining;iIndex++) {
1803
1804        int iRow = index[iIndex];
1805        double alpha = spare[iIndex];
1806        double oldValue = rhs[iRow];
1807        double value = oldValue-tentativeTheta*fabs(alpha);
1808
1809        if (value<-primalTolerance_) {
1810          // how much would it cost to go thru
1811          thruThis += alpha*
1812            nonLinearCost_->changeInCost(pivotVariable_[iRow],alpha);
1813          // goes on swapped list (also means candidates if too many)
1814          indexSwapped[numberSwapped++]=iRow;
1815          if (fabs(alpha)>bestPivot) {
1816            bestPivot=fabs(alpha);
1817            pivotRow_ = iRow;
1818            theta_ = oldValue/bestPivot;
1819          }
1820        } else {
1821          value = oldValue-upperTheta*fabs(alpha);
1822          if (value<-primalTolerance_ && fabs(alpha)>=acceptablePivot) 
1823            upperTheta = (oldValue+primalTolerance_)/fabs(alpha);
1824        }
1825      }
1826     
1827      maximumSwapped = max(maximumSwapped,numberSwapped);
1828
1829      double dualCheck = - 2.0*coeff2*tentativeTheta - coeff1 - 99999999;
1830      // but make a bit more pessimistic
1831      dualCheck=max(dualCheck-100.0*dualTolerance_,0.99*dualCheck);
1832      if (totalThru+thruThis>=dualCheck) {
1833        // We should be pivoting in this batch
1834        // so compress down to this lot
1835
1836        int saveNumber = numberRemaining;
1837        numberRemaining=0;
1838        for (iIndex=0;iIndex<numberSwapped;iIndex++) {
1839          int iRow = indexSwapped[iIndex];
1840          spare[numberRemaining]=way*work[iRow];
1841          index[numberRemaining++]=iRow;
1842        }
1843        memset(spare+numberRemaining,0,
1844               (saveNumber-numberRemaining)*sizeof(double));
1845        int iTry;
1846#define MAXTRY 100
1847        // first get ratio with tolerance
1848        for (iTry=0;iTry<MAXTRY;iTry++) {
1849         
1850          upperTheta=maximumMovement;
1851          numberSwapped = 0;
1852         
1853          for (iIndex=0;iIndex<numberRemaining;iIndex++) {
1854           
1855            int iRow = index[iIndex];
1856            double alpha = fabs(spare[iIndex]);
1857            double oldValue = rhs[iRow];
1858            double value = oldValue-upperTheta*alpha;
1859           
1860            if (value<-primalTolerance_ && alpha>=acceptablePivot) 
1861              upperTheta = (oldValue+primalTolerance_)/alpha;
1862           
1863          }
1864         
1865          // now look at best in this lot
1866          bestPivot=acceptablePivot;
1867          pivotRow_=-1;
1868          for (iIndex=0;iIndex<numberRemaining;iIndex++) {
1869           
1870            int iRow = index[iIndex];
1871            double alpha = spare[iIndex];
1872            double oldValue = rhs[iRow];
1873            double value = oldValue-upperTheta*fabs(alpha);
1874           
1875            if (value<=0) {
1876              // how much would it cost to go thru
1877              totalThru += alpha*
1878                nonLinearCost_->changeInCost(pivotVariable_[iRow],alpha);
1879              // goes on swapped list (also means candidates if too many)
1880              indexSwapped[numberSwapped++]=iRow;
1881              if (fabs(alpha)>bestPivot) {
1882                bestPivot=fabs(alpha);
1883                theta_ = oldValue/bestPivot;
1884                pivotRow_=iRow;
1885              }
1886            } else {
1887              value = oldValue-upperTheta*fabs(alpha);
1888              if (value<-primalTolerance_ && fabs(alpha)>=acceptablePivot) 
1889                upperTheta = (oldValue+primalTolerance_)/fabs(alpha);
1890            }
1891          }
1892         
1893          maximumSwapped = max(maximumSwapped,numberSwapped);
1894          if (bestPivot<0.1*bestEverPivot&&
1895              bestEverPivot>1.0e-6&&bestPivot<1.0e-3) {
1896            // back to previous one
1897            goBackOne = true;
1898            break;
1899          } else if (pivotRow_==-1&&upperTheta>largeValue_) {
1900            if (lastPivot>acceptablePivot) {
1901              // back to previous one
1902              goBackOne = true;
1903              //break;
1904            } else {
1905              // can only get here if all pivots so far too small
1906            }
1907            break;
1908          } else {
1909            dualCheck = - 2.0*coeff2*theta_ - coeff1-9999999;
1910            if (totalThru>=dualCheck) {
1911              break; // no point trying another loop
1912            } else {
1913              // skip this lot
1914              nonLinearCost_->goThru(numberSwapped,way,indexSwapped, work,rhs);
1915              lastPivotRow=pivotRow_;
1916              lastTheta = theta_;
1917              lastThru = numberThru;
1918              numberThru += numberSwapped;
1919              lastNumberSwapped = numberSwapped;
1920              memcpy(saveSwapped,indexSwapped,lastNumberSwapped*sizeof(int));
1921              if (bestPivot>bestEverPivot)
1922                bestEverPivot=bestPivot;
1923            }
1924          }
1925        }
1926        break;
1927      } else {
1928        // skip this lot
1929        nonLinearCost_->goThru(numberSwapped,way,indexSwapped, work,rhs);
1930        lastPivotRow=pivotRow_;
1931        lastTheta = theta_;
1932        lastThru = numberThru;
1933        numberThru += numberSwapped;
1934        lastNumberSwapped = numberSwapped;
1935        memcpy(saveSwapped,indexSwapped,lastNumberSwapped*sizeof(int));
1936        if (bestPivot>bestEverPivot)
1937          bestEverPivot=bestPivot;
1938        totalThru += thruThis;
1939        tentativeTheta = 2.0*upperTheta;
1940      }
1941    }
1942    // can get here without pivotRow_ set but with lastPivotRow
1943    if (goBackOne||(pivotRow_<0&&lastPivotRow>=0)) {
1944      // back to previous one
1945      pivotRow_=lastPivotRow;
1946      theta_ = lastTheta;
1947            // undo this lot
1948      nonLinearCost_->goBack(lastNumberSwapped,saveSwapped,rhs);
1949      memcpy(indexSwapped,saveSwapped,lastNumberSwapped*sizeof(int));
1950      numberSwapped = lastNumberSwapped;
1951    }
1952  }
1953
1954  if (pivotRow_>=0) {
1955   
1956#define MINIMUMTHETA 1.0e-12
1957    // will we need to increase tolerance
1958#ifdef CLP_DEBUG
1959    bool found=false;
1960#endif
1961    double largestInfeasibility = primalTolerance_;
1962    if (theta_<MINIMUMTHETA) {
1963      theta_=MINIMUMTHETA;
1964      for (iIndex=0;iIndex<numberSwapped;iIndex++) {
1965        int iRow = indexSwapped[iIndex];
1966#ifdef CLP_DEBUG
1967        if (iRow==pivotRow_)
1968          found=true;
1969#endif
1970        largestInfeasibility = max (largestInfeasibility,
1971                                    -(rhs[iRow]-fabs(work[iRow])*theta_));
1972      }
1973#ifdef CLP_DEBUG
1974      assert(found);
1975      if (largestInfeasibility>primalTolerance_&&(handler_->logLevel()&32))
1976        printf("Primal tolerance increased from %g to %g\n",
1977               primalTolerance_,largestInfeasibility);
1978#endif
1979      primalTolerance_ = max(primalTolerance_,largestInfeasibility);
1980    }
1981    alpha_ = work[pivotRow_];
1982    // translate to sequence
1983    sequenceOut_ = pivotVariable_[pivotRow_];
1984    valueOut_ = solution(sequenceOut_);
1985    lowerOut_=lower_[sequenceOut_];
1986    upperOut_=upper_[sequenceOut_];
1987
1988    if (way<0.0) 
1989      theta_ = - theta_;
1990    double newValue = valueOut_ - theta_*alpha_;
1991    if (alpha_*way<0.0) {
1992      directionOut_=-1;      // to upper bound
1993      if (fabs(theta_)>1.0e-6)
1994        upperOut_ = nonLinearCost_->nearest(sequenceOut_,newValue);
1995      else
1996        upperOut_ = newValue;
1997    } else {
1998      directionOut_=1;      // to lower bound
1999      if (fabs(theta_)>1.0e-6)
2000        lowerOut_ = nonLinearCost_->nearest(sequenceOut_,newValue);
2001      else
2002        lowerOut_ = newValue;
2003    }
2004    dualOut_ = reducedCost(sequenceOut_);
2005  } else {
2006    double trueMaximumMovement;
2007    if (way>0.0) 
2008      trueMaximumMovement = min(1.0e30,upperIn_-valueIn_);
2009    else
2010      trueMaximumMovement = min(1.0e30,valueIn_-lowerIn_);
2011    if (maximumMovement<1.0e20&&maximumMovement==trueMaximumMovement) {
2012      // flip
2013      theta_ = maximumMovement;
2014      pivotRow_ = -2; // so we can tell its a flip
2015      result=0;
2016      sequenceOut_ = sequenceIn_;
2017      valueOut_ = valueIn_;
2018      dualOut_ = dualIn_;
2019      lowerOut_ = lowerIn_;
2020      upperOut_ = upperIn_;
2021      alpha_ = 0.0;
2022      if (way<0.0) {
2023        directionOut_=1;      // to lower bound
2024        theta_ = lowerOut_ - valueOut_;
2025      } else {
2026        directionOut_=-1;      // to upper bound
2027        theta_ = upperOut_ - valueOut_;
2028      }
2029      // we may still have sj to get rid of
2030    } else if (fabs(maximumMovement-d2)<dualTolerance_) {
2031      // sj going to zero
2032      result=0;
2033      assert (pivotRow_<0);
2034      nonLinearCost_->setBounds(crucialSj, 0.0,0.0);
2035      lower_[crucialSj]=0.0;
2036      upper_[crucialSj]=0.0;
2037      setColumnStatus(crucialSj,isFixed);
2038      pivotRow_ = iSjRow;
2039      alpha_ = work[pivotRow_];
2040      // translate to sequence
2041      sequenceOut_ = pivotVariable_[pivotRow_];
2042      valueOut_ = solution(sequenceOut_);
2043      lowerOut_=lower_[sequenceOut_];
2044      upperOut_=upper_[sequenceOut_];
2045      theta_ = d2;
2046      if (way<0.0) 
2047        theta_ = - theta_;
2048      double newValue = valueOut_ - theta_*alpha_;
2049      if (alpha_*way<0.0) {
2050        directionOut_=-1;      // to upper bound
2051        if (fabs(theta_)>1.0e-6)
2052          upperOut_ = nonLinearCost_->nearest(sequenceOut_,newValue);
2053        else
2054          upperOut_ = newValue;
2055      } else {
2056        directionOut_=1;      // to lower bound
2057        if (fabs(theta_)>1.0e-6)
2058          lowerOut_ = nonLinearCost_->nearest(sequenceOut_,newValue);
2059        else
2060          lowerOut_ = newValue;
2061      }
2062      //????
2063      dualOut_ = reducedCost(sequenceOut_);
2064    } else {
2065      // need to do something
2066      abort();
2067    }
2068  }
2069
2070  // clear arrays
2071
2072  memset(spare,0,numberRemaining*sizeof(double));
2073  memset(saveSwapped,0,maximumSwapped*sizeof(int));
2074
2075  // put back original bounds etc
2076  nonLinearCost_->goBackAll(rhsArray);
2077
2078  rhsArray->clear();
2079  // Change in objective will be theta*coeff1 + theta*theta*coeff2
2080  objectiveValue_ += theta_*coeff1+theta_*theta_*coeff2;
2081  printf("New objective value %g\n",objectiveValue_);
2082  {
2083    int iColumn;
2084    objectiveValue_ =0.0;
2085    CoinPackedMatrix * quadratic = originalModel->quadraticObjective();
2086    const int * columnQuadratic = quadratic->getIndices();
2087    const int * columnQuadraticStart = quadratic->getVectorStarts();
2088    const int * columnQuadraticLength = quadratic->getVectorLengths();
2089    const double * quadraticElement = quadratic->getElements();
2090    int numberColumns = originalModel->numberColumns();
2091    const double * objective = originalModel->objective();
2092    for (iColumn=0;iColumn<numberColumns;iColumn++) 
2093      objectiveValue_ += objective[iColumn]*solution_[iColumn];
2094    for (iColumn=0;iColumn<numberColumns;iColumn++) {
2095      double valueI = solution_[iColumn];
2096      int j;
2097      for (j=columnQuadraticStart[iColumn];
2098           j<columnQuadraticStart[iColumn]+columnQuadraticLength[iColumn];j++) {
2099        int jColumn = columnQuadratic[j];
2100        double valueJ = solution_[jColumn];
2101        double elementValue = quadraticElement[j];
2102        objectiveValue_ += 0.5*valueI*valueJ*elementValue;
2103      }
2104    }
2105    printf("Objective value %g\n",objectiveValue());
2106  }
2107  return result;
2108}
2109// Just for debug
2110int 
2111ClpSimplexPrimalQuadratic::checkComplementary (const 
2112                                               ClpQuadraticInfo * info)
2113{
2114  const ClpSimplex * originalModel= info->originalModel();
2115  int numberXColumns = originalModel->numberColumns();
2116  int i;
2117  for (i=0;i<numberXColumns;i++) {
2118    int jSequence = i+ numberRows_;
2119    if (getColumnStatus(i)==basic) {
2120      if (getColumnStatus(jSequence)==basic)
2121        printf("Struct %d (%g) and %d (%g) both basic\n",
2122               i,solution_[i],jSequence,solution_[jSequence]);
2123    }
2124  }
2125  int originalNumberRows = originalModel->numberRows();
2126  int offset = numberXColumns;
2127  for (i=0;i<originalNumberRows;i++) {
2128    int jSequence = i+offset;
2129    if (getRowStatus(i)==basic) {
2130      if (getColumnStatus(jSequence)==basic)
2131        printf("Row %d (%g) and %d (%g) both basic\n",
2132               i,solution_[i],jSequence,solution_[jSequence]);
2133    }
2134  }
2135  return 0;
2136}
2137 
2138/* This creates the large version of QP and
2139      fills in quadratic information
2140*/
2141ClpSimplexPrimalQuadratic * 
2142ClpSimplexPrimalQuadratic::makeQuadratic(ClpQuadraticInfo & info)
2143{
2144
2145  // Get list of non linear columns
2146  CoinPackedMatrix * quadratic = quadraticObjective();
2147  if (!quadratic||!quadratic->getNumElements()) {
2148    // no quadratic part
2149    return NULL;
2150  }
2151
2152  int numberColumns = this->numberColumns();
2153  double * columnLower = this->columnLower();
2154  double * columnUpper = this->columnUpper();
2155  double * objective = this->objective();
2156  double * solution = this->primalColumnSolution();
2157  double * dj = this->dualColumnSolution();
2158  double * pi = this->dualRowSolution();
2159
2160  int numberRows = this->numberRows();
2161  double * rowLower = this->rowLower();
2162  double * rowUpper = this->rowUpper();
2163
2164  // and elements
2165  CoinPackedMatrix * matrix = this->matrix();
2166  const int * row = matrix->getIndices();
2167  const int * columnStart = matrix->getVectorStarts();
2168  const double * element =  matrix->getElements();
2169  const int * columnLength = matrix->getVectorLengths();
2170
2171  int iColumn;
2172  const int * columnQuadratic = quadratic->getIndices();
2173  const int * columnQuadraticStart = quadratic->getVectorStarts();
2174  const int * columnQuadraticLength = quadratic->getVectorLengths();
2175  const double * quadraticElement = quadratic->getElements();
2176#if 0
2177  // deliberate bad solution
2178  // Change to use phase
2179  //double * saveO = new double[numberColumns];
2180  //memcpy(saveO,objective,numberColumns*sizeof(double));
2181  //memset(objective,0,numberColumns*sizeof(double));
2182  tempMessage messageHandler(this);;
2183  passInMessageHandler(&messageHandler);
2184  factorization()->maximumPivots(1);
2185  primal();
2186  CoinMessageHandler handler2;
2187  passInMessageHandler(&handler2);
2188  factorization()->maximumPivots(100);
2189  setMaximumIterations(1000);
2190#endif
2191  //memcpy(objective,saveO,numberColumns*sizeof(double));
2192  // Get a feasible solution
2193  //printf("For testing - deliberate bad solution\n");
2194  //columnUpper[0]=0.0;
2195  //columnLower[0]=0.0;
2196  //quadraticSLP(50,1.0e-4);
2197  //primal(1);
2198  //columnUpper[0]=COIN_DBL_MAX;
2199 
2200  // Create larger problem
2201  // First workout how many rows extra
2202  info=ClpQuadraticInfo(this);
2203  int numberQuadratic = info.numberQuadraticColumns();
2204  int newNumberRows = numberRows+numberQuadratic;
2205  int newNumberColumns = numberColumns + numberRows + numberQuadratic;
2206  int numberElements = 2*matrix->getNumElements()
2207    +2*quadratic->getNumElements()
2208    + numberQuadratic;
2209  double * elements2 = new double[numberElements];
2210  int * start2 = new int[newNumberColumns+1];
2211  int * row2 = new int[numberElements];
2212  double * objective2 = new double[newNumberColumns];
2213  double * columnLower2 = new double[newNumberColumns];
2214  double * columnUpper2 = new double[newNumberColumns];
2215  double * rowLower2 = new double[newNumberRows];
2216  double * rowUpper2 = new double[newNumberRows];
2217  const int * which = info.quadraticSequence();
2218  const int * back = info.backSequence();
2219  memcpy(rowLower2,rowLower,numberRows*sizeof(double));
2220  memcpy(rowUpper2,rowUpper,numberRows*sizeof(double));
2221  int iRow;
2222  for (iRow=0;iRow<numberQuadratic;iRow++) {
2223    double cost = objective[back[iRow]];
2224    rowLower2[iRow+numberRows]=cost;
2225    rowUpper2[iRow+numberRows]=cost;
2226  }
2227  memset(objective2,0,newNumberColumns*sizeof(double));
2228  // Get a row copy of quadratic objective in standard format
2229  CoinPackedMatrix copyQ;
2230  copyQ.reverseOrderedCopyOf(*quadratic);
2231  const int * columnQ = copyQ.getIndices();
2232  const CoinBigIndex * rowStartQ = copyQ.getVectorStarts();
2233  const int * rowLengthQ = copyQ.getVectorLengths(); 
2234  const double * elementByRowQ = copyQ.getElements();
2235  // Move solution across
2236  double * solution2 = new double[newNumberColumns];
2237  memset(solution2,0,newNumberColumns*sizeof(double));
2238  newNumberColumns=0;
2239  numberElements=0;
2240  start2[0]=0;
2241  // x
2242  memcpy(dj,objective,numberColumns*sizeof(double));
2243  for (iColumn=0;iColumn<numberColumns;iColumn++) {
2244    // Original matrix
2245    columnLower2[iColumn]=columnLower[iColumn];
2246    columnUpper2[iColumn]=columnUpper[iColumn];
2247    solution2[iColumn]=solution[iColumn];
2248    int j;
2249    for (j=columnStart[iColumn];
2250         j<columnStart[iColumn]+columnLength[iColumn];
2251         j++) {
2252      elements2[numberElements]=element[j];
2253      row2[numberElements++]=row[j];
2254    }
2255    if (which[iColumn]>=0) {
2256      // Quadratic and modify djs
2257      for (j=columnQuadraticStart[iColumn];
2258           j<columnQuadraticStart[iColumn]+columnQuadraticLength[iColumn];
2259           j++) {
2260        int jColumn = columnQuadratic[j];
2261        double value = quadraticElement[j];
2262        if (iColumn!=jColumn) 
2263          value *= 0.5;
2264        dj[jColumn] += solution[iColumn]*value;
2265        elements2[numberElements]=-value;
2266        row2[numberElements++]=which[jColumn]+numberRows;
2267      }
2268      for (j=rowStartQ[iColumn];
2269           j<rowStartQ[iColumn]+rowLengthQ[iColumn];
2270           j++) {
2271        int jColumn = columnQ[j];
2272        double value = elementByRowQ[j];
2273        if (iColumn!=jColumn) { 
2274          value *= 0.5;
2275          dj[jColumn] += solution[iColumn]*value;
2276          elements2[numberElements]=-value;
2277          row2[numberElements++]=which[jColumn]+numberRows;
2278        }
2279      }
2280    }
2281    start2[iColumn+1]=numberElements;
2282  }
2283  newNumberColumns=numberColumns;
2284  // pi
2285  // Get a row copy in standard format
2286  CoinPackedMatrix copy;
2287  copy.reverseOrderedCopyOf(*(this->matrix()));
2288  // get matrix data pointers
2289  const int * column = copy.getIndices();
2290  const CoinBigIndex * rowStart = copy.getVectorStarts();
2291  const int * rowLength = copy.getVectorLengths(); 
2292  const double * elementByRow = copy.getElements();
2293  for (iRow=0;iRow<numberRows;iRow++) {
2294    solution2[newNumberColumns]=pi[iRow];
2295    double value = pi[iRow];
2296    columnLower2[newNumberColumns]=-COIN_DBL_MAX;
2297    columnUpper2[newNumberColumns]=COIN_DBL_MAX;
2298    int j;
2299    for (j=rowStart[iRow];
2300         j<rowStart[iRow]+rowLength[iRow];
2301         j++) {
2302      double elementValue=elementByRow[j];
2303      int jColumn = column[j];
2304      elements2[numberElements]=elementValue;
2305      row2[numberElements++]=jColumn+numberRows;
2306      dj[jColumn]-= value*elementValue;
2307    }
2308    newNumberColumns++;
2309    start2[newNumberColumns]=numberElements;
2310  }
2311  // djs
2312  for (iColumn=0;iColumn<numberQuadratic;iColumn++) {
2313    columnLower2[newNumberColumns]=-COIN_DBL_MAX;
2314    columnUpper2[newNumberColumns]=COIN_DBL_MAX;
2315    solution2[newNumberColumns]=dj[iColumn];
2316    elements2[numberElements]=1.0;
2317    row2[numberElements++]=back[iColumn]+numberRows;
2318    newNumberColumns++;
2319    start2[newNumberColumns]=numberElements;
2320  }
2321  // Create model
2322  ClpSimplex * model2 = new ClpSimplex(*this);
2323  model2->resize(0,0);
2324  model2->loadProblem(newNumberColumns,newNumberRows,
2325                     start2,row2, elements2,
2326                     columnLower2,columnUpper2,
2327                     objective2,
2328                     rowLower2,rowUpper2);
2329  delete [] objective2;
2330  delete [] rowLower2;
2331  delete [] rowUpper2;
2332  delete [] columnLower2;
2333  delete [] columnUpper2;
2334  // Now create expanded quadratic objective for use in primalRow
2335  // Later pack down in some way
2336  start2[0]=0;
2337  numberElements=0;
2338  for (iColumn=0;iColumn<numberColumns;iColumn++) {
2339    // Quadratic
2340    int j;
2341    for (j=columnQuadraticStart[iColumn];
2342         j<columnQuadraticStart[iColumn]+columnQuadraticLength[iColumn];
2343         j++) {
2344      int jColumn = columnQuadratic[j];
2345      double value = quadraticElement[j];
2346      if (iColumn!=jColumn) 
2347        value *= 0.5;
2348      elements2[numberElements]=value;
2349      row2[numberElements++]=jColumn;
2350    }
2351    for (j=rowStartQ[iColumn];
2352         j<rowStartQ[iColumn]+rowLengthQ[iColumn];
2353         j++) {
2354      int jColumn = columnQ[j];
2355      double value = elementByRowQ[j];
2356      if (iColumn!=jColumn) { 
2357        value *= 0.5;
2358        elements2[numberElements]=value;
2359        row2[numberElements++]=jColumn;
2360      }
2361    }
2362    start2[iColumn+1]=numberElements;
2363  }
2364  // and pad
2365  for (;iColumn<newNumberColumns;iColumn++)
2366    start2[iColumn+1]=numberElements;
2367  // Load up objective
2368  model2->loadQuadraticObjective(newNumberColumns,start2,row2,elements2);
2369  delete [] start2;
2370  delete [] row2;
2371  delete [] elements2;
2372  model2->allSlackBasis();
2373  model2->scaling(false);
2374  model2->setLogLevel(this->logLevel());
2375  // Move solution across
2376  memcpy(model2->primalColumnSolution(),solution2,
2377         newNumberColumns*sizeof(double));
2378  columnLower2 = model2->columnLower();
2379  columnUpper2 = model2->columnUpper();
2380  delete [] solution2;
2381  solution2 = model2->primalColumnSolution();
2382  // Compute row activities and check feasible
2383  double * rowSolution2 = model2->primalRowSolution();
2384  memset(rowSolution2,0,newNumberRows*sizeof(double));
2385  model2->times(1.0,solution2,rowSolution2);
2386  rowLower2 = model2->rowLower();
2387  rowUpper2 = model2->rowUpper();
2388#if 0
2389  for (iColumn=0;iColumn<numberColumns;iColumn++) {
2390    Status xStatus = getColumnStatus(iColumn);
2391    bool isSuperBasic;
2392    int iS = iColumn+newNumberRows;
2393    double value = solution2[iS];
2394    if (fabs(value)>dualTolerance_)
2395      isSuperBasic=true;
2396    else
2397      isSuperBasic=false;
2398    // For moment take all x out of basis
2399    // Does not seem right
2400    isSuperBasic=true;
2401    model2->setColumnStatus(iColumn,xStatus);
2402    if (xStatus==basic) {
2403      if (!isSuperBasic) {
2404        model2->setRowStatus(numberRows+iColumn,basic);
2405        model2->setColumnStatus(iS,superBasic);
2406      } else {
2407        model2->setRowStatus(numberRows+iColumn,isFixed);
2408        model2->setColumnStatus(iS,basic);
2409        model2->setColumnStatus(iColumn,superBasic);
2410      }
2411    } else {
2412      model2->setRowStatus(numberRows+iColumn,isFixed);
2413      model2->setColumnStatus(iS,basic);
2414    }
2415  }
2416  for (iRow=0;iRow<numberRows;iRow++) {
2417    Status rowStatus = getRowStatus(iRow);
2418    model2->setRowStatus(iRow,rowStatus);
2419    if (rowStatus!=basic) {
2420      model2->setColumnStatus(iRow+numberColumns,basic); // make dual basic
2421    }
2422    assert (rowSolution2[iRow]>=rowLower2[iRow]-primalTolerance_);
2423    assert (rowSolution2[iRow]<=rowUpper2[iRow]+primalTolerance_);
2424  }
2425  // why ?? - take duals out and adjust
2426  for (iRow=0;iRow<numberRows;iRow++) {
2427    model2->setRowStatus(iRow,basic);
2428    model2->setColumnStatus(iRow+numberColumns,superBasic);
2429    solution2[iRow+numberColumns]=0.0;
2430  }
2431#else
2432  for (iRow=0;iRow<numberRows;iRow++) {
2433    assert (rowSolution2[iRow]>=rowLower2[iRow]-primalTolerance_);
2434    assert (rowSolution2[iRow]<=rowUpper2[iRow]+primalTolerance_);
2435    model2->setRowStatus(iRow,basic);
2436    model2->setColumnStatus(iRow+numberColumns,superBasic);
2437    solution2[iRow+numberColumns]=0.0;
2438  }
2439  for (iColumn=numberRows+numberColumns;iColumn<newNumberColumns;iColumn++) {
2440    model2->setColumnStatus(iColumn,basic);
2441    model2->setRowStatus(iColumn-numberColumns,isFixed);
2442  }
2443#endif
2444  memset(rowSolution2,0,newNumberRows*sizeof(double));
2445  model2->times(1.0,solution2,rowSolution2);
2446  for (iColumn=0;iColumn<numberQuadratic;iColumn++) {
2447    int iS = back[iColumn]+newNumberRows;
2448    int iRow = iColumn+numberRows;
2449    double value = rowSolution2[iRow];
2450    if (value>rowUpper2[iRow]) {
2451      rowSolution2[iRow] = rowUpper2[iRow];
2452      solution2[iS]-=value-rowUpper2[iRow];
2453    } else {
2454      rowSolution2[iRow] = rowLower2[iRow];
2455      solution2[iS]-=value-rowLower2[iRow];
2456    }
2457  }
2458
2459 
2460  // See if any s sub j have wrong sign and/or use djs from infeasibility objective
2461  double objectiveOffset;
2462  getDblParam(ClpObjOffset,objectiveOffset);
2463  double objValue = -objectiveOffset;
2464  for (iColumn=0;iColumn<numberColumns;iColumn++) 
2465    objValue += objective[iColumn]*solution2[iColumn];
2466  for (iColumn=0;iColumn<numberColumns;iColumn++) {
2467    double valueI = solution2[iColumn];
2468    int j;
2469    for (j=columnQuadraticStart[iColumn];
2470         j<columnQuadraticStart[iColumn]+columnQuadraticLength[iColumn];j++) {
2471      int jColumn = columnQuadratic[j];
2472      double valueJ = solution2[jColumn];
2473      double elementValue = quadraticElement[j];
2474      objValue += 0.5*valueI*valueJ*elementValue;
2475    }
2476  }
2477  printf("Objective value %g\n",objValue);
2478  for (iColumn=0;iColumn<newNumberColumns;iColumn++) 
2479    printf("%d %g\n",iColumn,solution2[iColumn]);
2480  return (ClpSimplexPrimalQuadratic *) model2;
2481}
2482
2483// This moves solution back and deletes information
2484int 
2485ClpSimplexPrimalQuadratic::endQuadratic(ClpSimplexPrimalQuadratic * quadraticModel,
2486                   ClpQuadraticInfo & info)
2487{
2488  memcpy(dualRowSolution(),quadraticModel->primalColumnSolution()+numberColumns_,numberRows_*sizeof(double));
2489  const double * solution2 = quadraticModel->primalColumnSolution();
2490  memcpy(primalColumnSolution(),solution2,numberColumns_*sizeof(double));
2491  memset(quadraticModel->primalRowSolution(),0,
2492         quadraticModel->numberRows()*sizeof(double));
2493  quadraticModel->times(1.0,quadraticModel->primalColumnSolution(),
2494               quadraticModel->primalRowSolution());
2495  memcpy(dualColumnSolution(),
2496         quadraticModel->primalColumnSolution()+numberRows_+numberColumns_,
2497         numberColumns_*sizeof(double));
2498
2499  int iColumn;
2500  double objectiveOffset;
2501  getDblParam(ClpObjOffset,objectiveOffset);
2502  double objValue = -objectiveOffset;
2503  double * objective = this->objective();
2504  for (iColumn=0;iColumn<numberColumns_;iColumn++) 
2505    objValue += objective[iColumn]*solution2[iColumn];
2506  CoinPackedMatrix * quadratic = quadraticObjective();
2507  if (quadratic) {
2508    const int * columnQuadratic = quadratic->getIndices();
2509    const int * columnQuadraticStart = quadratic->getVectorStarts();
2510    const int * columnQuadraticLength = quadratic->getVectorLengths();
2511    const double * quadraticElement = quadratic->getElements();
2512    for (iColumn=0;iColumn<numberColumns_;iColumn++) {
2513      double valueI = solution2[iColumn];
2514      if (fabs(valueI)>1.0e-5) {
2515        int djColumn = iColumn+numberRows_+numberColumns_;
2516        assert(solution2[djColumn]<1.0e-7);
2517      }
2518      int j;
2519      for (j=columnQuadraticStart[iColumn];
2520           j<columnQuadraticStart[iColumn]+columnQuadraticLength[iColumn];j++) {
2521        int jColumn = columnQuadratic[j];
2522        double valueJ = solution2[jColumn];
2523        double elementValue = quadraticElement[j];
2524        objValue += 0.5*valueI*valueJ*elementValue;
2525      }
2526    }
2527    objectiveValue_ = objValue + objectiveOffset;
2528  }
2529  printf("Objective value %g\n",objValue);
2530  return 0;
2531}
2532
2533/// Default constructor.
2534ClpQuadraticInfo::ClpQuadraticInfo()
2535  : originalModel_(NULL),
2536    quadraticSequence_(NULL),
2537    backSequence_(NULL),
2538    crucialSj_(-1),
2539    numberXRows_(-1),
2540    numberXColumns_(-1),
2541    numberQuadraticColumns_(0)
2542{
2543}
2544// Constructor from original model
2545ClpQuadraticInfo::ClpQuadraticInfo(const ClpSimplex * model)
2546  : originalModel_(model),
2547    quadraticSequence_(NULL),
2548    backSequence_(NULL),
2549    crucialSj_(-1),
2550    numberXRows_(-1),
2551    numberXColumns_(-1),
2552    numberQuadraticColumns_(0)
2553{
2554  if (originalModel_) {
2555    numberXRows_ = originalModel_->numberRows();
2556    numberXColumns_ = originalModel_->numberColumns();
2557    quadraticSequence_ = new int[numberXColumns_];
2558    backSequence_ = new int[numberXColumns_];
2559    int i;
2560    numberQuadraticColumns_=numberXColumns_;
2561    for (i=0;i<numberXColumns_;i++) {
2562      quadraticSequence_[i]=i;
2563      backSequence_[i]=i;
2564    }
2565  }
2566}
2567// Destructor
2568ClpQuadraticInfo:: ~ClpQuadraticInfo()
2569{
2570  delete [] quadraticSequence_;
2571  delete [] backSequence_;
2572}
2573// Copy
2574ClpQuadraticInfo::ClpQuadraticInfo(const ClpQuadraticInfo& rhs)
2575  : originalModel_(rhs.originalModel_),
2576    quadraticSequence_(NULL),
2577    backSequence_(NULL),
2578    crucialSj_(rhs.crucialSj_),
2579    numberXRows_(rhs.numberXRows_),
2580    numberXColumns_(rhs.numberXColumns_),
2581    numberQuadraticColumns_(rhs.numberQuadraticColumns_)
2582{
2583  if (numberXColumns_) {
2584    quadraticSequence_ = new int[numberXColumns_];
2585    memcpy(quadraticSequence_,rhs.quadraticSequence_,
2586           numberXColumns_*sizeof(int));
2587    backSequence_ = new int[numberXColumns_];
2588    memcpy(backSequence_,rhs.backSequence_,
2589           numberXColumns_*sizeof(int));
2590  }
2591}
2592// Assignment
2593ClpQuadraticInfo & 
2594ClpQuadraticInfo::operator=(const ClpQuadraticInfo&rhs)
2595{
2596  if (this != &rhs) {
2597    originalModel_ = rhs.originalModel_;
2598    delete [] quadraticSequence_;
2599    quadraticSequence_ = NULL;
2600    delete [] backSequence_;
2601    backSequence_ = NULL;
2602    crucialSj_ = rhs.crucialSj_;
2603    numberXRows_ = rhs.numberXRows_;
2604    numberXColumns_ = rhs.numberXColumns_;
2605    numberQuadraticColumns_=rhs.numberQuadraticColumns_;
2606    if (numberXColumns_) {
2607      quadraticSequence_ = new int[numberXColumns_];
2608      memcpy(quadraticSequence_,rhs.quadraticSequence_,
2609             numberXColumns_*sizeof(int));
2610      backSequence_ = new int[numberXColumns_];
2611      memcpy(backSequence_,rhs.backSequence_,
2612             numberXColumns_*sizeof(int));
2613    }
2614  }
2615  return *this;
2616}
Note: See TracBrowser for help on using the repository browser.