source: branches/pre/ClpSimplexPrimalQuadratic.cpp @ 180

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

new stuff

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