source: branches/pre/ClpPrimalQuadraticDantzig.cpp @ 1926

Last change on this file since 1926 was 196, checked in by forrest, 17 years ago

64 bit, Sbb etc

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 7.1 KB
Line 
1// Copyright (C) 2003, International Business Machines
2// Corporation and others.  All Rights Reserved.
3
4#include "CoinPragma.hpp"
5
6#include <cstdio>
7
8#include "CoinIndexedVector.hpp"
9
10#include "ClpSimplex.hpp"
11#include "ClpPrimalQuadraticDantzig.hpp"
12#include "ClpSimplexPrimalQuadratic.hpp"
13#include "ClpFactorization.hpp"
14#include "ClpPackedMatrix.hpp"
15
16//#############################################################################
17// Constructors / Destructor / Assignment
18//#############################################################################
19
20//-------------------------------------------------------------------
21// Default Constructor
22//-------------------------------------------------------------------
23ClpPrimalQuadraticDantzig::ClpPrimalQuadraticDantzig () 
24: ClpPrimalColumnPivot()
25{
26  type_=1;
27  quadraticInfo_=NULL;
28}
29
30//-------------------------------------------------------------------
31// Copy constructor
32//-------------------------------------------------------------------
33ClpPrimalQuadraticDantzig::ClpPrimalQuadraticDantzig (const ClpPrimalQuadraticDantzig & source) 
34: ClpPrimalColumnPivot(source)
35{ 
36  quadraticInfo_=source.quadraticInfo_;
37
38}
39// Constructor from model
40ClpPrimalQuadraticDantzig::ClpPrimalQuadraticDantzig(
41                  ClpSimplexPrimalQuadratic * model, 
42                  ClpQuadraticInfo * info)
43  :ClpPrimalColumnPivot()
44{
45  type_=1;
46  // Safe as has no extra data
47  model_=(ClpSimplex *) model;
48  quadraticInfo_ = info;
49}
50
51//-------------------------------------------------------------------
52// Destructor
53//-------------------------------------------------------------------
54ClpPrimalQuadraticDantzig::~ClpPrimalQuadraticDantzig ()
55{
56
57}
58
59//----------------------------------------------------------------
60// Assignment operator
61//-------------------------------------------------------------------
62ClpPrimalQuadraticDantzig &
63ClpPrimalQuadraticDantzig::operator=(const ClpPrimalQuadraticDantzig& rhs)
64{
65  if (this != &rhs) {
66    ClpPrimalColumnPivot::operator=(rhs);
67    quadraticInfo_ = rhs.quadraticInfo_;
68  }
69  return *this;
70}
71
72// Returns pivot column, -1 if none
73int 
74ClpPrimalQuadraticDantzig::pivotColumn(CoinIndexedVector * updates,
75                                    CoinIndexedVector * spareRow1,
76                                    CoinIndexedVector * spareRow2,
77                                    CoinIndexedVector * spareColumn1,
78                                    CoinIndexedVector * spareColumn2)
79{
80  assert(model_);
81  // Find out where stuff is
82  int originalNumberColumns = quadraticInfo_->numberXColumns();
83  int start = model_->numberRows();
84  double * solution = model_->solutionRegion();
85  //double * lower = model_->lowerRegion();
86  //double * upper = model_->upperRegion();
87  //double * trueLower = model_->columnLower();
88  //double * trueUpper = model_->columnUpper();
89 
90  //double tolerance=model_->currentPrimalTolerance();
91
92  double dualTolerance = model_->dualTolerance()*1000.0;
93  double bestDj = dualTolerance;
94  double bestInfeasibleDj = 10.0*dualTolerance;
95  int bestSequence=-1;
96  double dualIn=0.0;
97
98  int iSequence;
99  assert (!model_->scalingFlag());
100  const double * pi = solution+quadraticInfo_->numberXColumns();
101  // Matrix for linear stuff
102  CoinPackedMatrix * matrix = model_->matrix();
103  const int * row = matrix->getIndices();
104  const int * columnStart = matrix->getVectorStarts();
105  const double * element =  matrix->getElements();
106  const int * columnLength = matrix->getVectorLengths();
107  const int * which = quadraticInfo_->quadraticSequence();
108  const double * objective = quadraticInfo_->linearObjective();
109  const double * djWeight = quadraticInfo_->djWeight();
110  for (iSequence=0;iSequence<originalNumberColumns;iSequence++) {
111    if (model_->flagged(iSequence))
112      continue;
113    int jSequence = which[iSequence];
114    double value;
115    if (jSequence>=0) {
116      jSequence += start;
117      value = solution[jSequence];
118    } else {
119      value=objective[iSequence];
120      int j;
121      for (j=columnStart[iSequence];j<columnStart[iSequence]+columnLength[iSequence]; j++) {
122        int iRow = row[j];
123        value -= element[j]*pi[iRow];
124      }
125    }
126    double value2 = value*djWeight[iSequence];
127    if (fabs(value2)>dualTolerance)
128      value=value2;
129    else if (value<-dualTolerance)
130      value = -1.001*dualTolerance;
131    else if (value>dualTolerance)
132      value = 1.001*dualTolerance;
133    if (djWeight[iSequence]<1.0e-6)
134      value=value2;
135    ClpSimplex::Status status = model_->getStatus(iSequence);
136   
137    switch(status) {
138     
139    case ClpSimplex::basic:
140      if (fabs(value)>bestInfeasibleDj) {
141        // infeasible - high priority
142        bestDj=COIN_DBL_MAX;
143        bestInfeasibleDj = fabs(value);
144        bestSequence = jSequence;
145        dualIn = value;
146        abort();
147      }
148      break;
149    case ClpSimplex::isFixed:
150      break;
151    case ClpSimplex::isFree:
152    case ClpSimplex::superBasic:
153      abort();
154      break;
155    case ClpSimplex::atUpperBound:
156      if (value>bestDj) {
157        bestDj = value;
158        bestSequence = iSequence;
159        dualIn = value;
160      }
161      break;
162    case ClpSimplex::atLowerBound:
163      if (value<-bestDj) {
164        bestDj = -value;
165        bestSequence = iSequence;
166        dualIn = value;
167      }
168    }
169  }
170  // And slacks
171  // Value of sj is - value of pi
172  int firstSlack = model_->numberColumns();
173  double * lower = model_->lowerRegion()+model_->numberColumns();
174  double * upper = model_->upperRegion()+model_->numberColumns();
175  int jSequence;
176  int originalNumberRows = quadraticInfo_->numberXRows();
177  for (jSequence=0;jSequence<originalNumberRows;jSequence++) {
178    int iSequence  = jSequence + firstSlack;
179    if (model_->flagged(iSequence))
180      continue;
181    int iPi = jSequence+originalNumberColumns;
182    // for slacks either pi zero or row at bound
183    // for L rows pi negative okay so choose if positive
184    double value = solution[iPi];
185    double value2 = value*djWeight[iSequence];
186    if (fabs(value2)>dualTolerance)
187      value=value2;
188    else if (value<-dualTolerance)
189      value = -1.001*dualTolerance;
190    else if (value>dualTolerance)
191      value = 1.001*dualTolerance;
192    if (djWeight[iSequence]<1.0e-6)
193      value=value2;
194    ClpSimplex::Status status = model_->getStatus(iSequence);
195    switch(status) {
196     
197    case ClpSimplex::basic:
198      if (fabs(value)>bestDj&&upper[jSequence]>lower[jSequence]) {
199        bestDj = fabs(value);
200        bestSequence = iPi; //?
201        abort();
202      }
203      break;
204    case ClpSimplex::isFixed:
205      break;
206    case ClpSimplex::isFree:
207    case ClpSimplex::superBasic:
208      if (fabs(value)>bestDj) {
209        bestDj = fabs(value);
210        bestSequence = iSequence;
211        dualIn = value;
212      }
213      break;
214    case ClpSimplex::atUpperBound:
215      if (value>bestDj) {
216        bestDj = value;
217        bestSequence = iSequence;
218        dualIn = value;
219      }
220      break;
221    case ClpSimplex::atLowerBound:
222      if (-value>bestDj) {
223        bestDj = -value;
224        bestSequence = iSequence;
225        dualIn = value;
226      }
227    }
228  }
229  if (bestSequence>=0) {
230    model_->djRegion()[bestSequence] = dualIn;
231  }
232  return bestSequence;
233}
234 
235//-------------------------------------------------------------------
236// Clone
237//-------------------------------------------------------------------
238ClpPrimalColumnPivot * ClpPrimalQuadraticDantzig::clone(bool CopyData) const
239{
240  if (CopyData) {
241    return new ClpPrimalQuadraticDantzig(*this);
242  } else {
243    return new ClpPrimalQuadraticDantzig();
244  }
245}
246
Note: See TracBrowser for help on using the repository browser.