source: trunk/Clp/examples/ekk_interface.cpp @ 800

Last change on this file since 800 was 225, checked in by forrest, 16 years ago

This should break everything

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 7.9 KB
Line 
1#include <cassert>
2
3#include "ClpSimplexPrimal.hpp"
4#include "ClpFactorization.hpp"
5#include "ClpPresolve.hpp"
6#include "ekk_c_api.h"
7//#include "ekk_c_api_undoc.h"
8 extern "C"{
9 OSLLIBAPI void * OSLLINKAGE ekk_compressModel(EKKModel * model);
10  OSLLIBAPI void OSLLINKAGE ekk_decompressModel(EKKModel * model,
11                                                void * compressInfo);
12 }
13
14static ClpPresolve * presolveInfo=NULL;
15
16static ClpSimplex * clpmodel(EKKModel * model,int startup)
17{
18  ClpSimplex * clp = new ClpSimplex();;
19  int numberRows=ekk_getInumrows(model);
20  int numberColumns= ekk_getInumcols(model);
21  clp->loadProblem(numberColumns,numberRows,ekk_blockColumn(model,0),
22                  ekk_blockRow(model,0),ekk_blockElement(model,0),
23                  ekk_collower(model),ekk_colupper(model),
24                  ekk_objective(model),
25                  ekk_rowlower(model),ekk_rowupper(model));
26  clp->setOptimizationDirection((int) ekk_getRmaxmin(model));
27  clp->setPrimalTolerance(ekk_getRtolpinf(model));
28  if (ekk_getRpweight(model)!=0.1)
29    clp->setInfeasibilityCost(1.0/ekk_getRpweight(model));
30  clp->setDualTolerance(ekk_getRtoldinf(model));
31  if (ekk_getRdweight(model)!=0.1)
32    clp->setDualBound(1.0/ekk_getRdweight(model));
33  clp->setDblParam(ClpObjOffset,ekk_getRobjectiveOffset(model));
34  const int * rowStatus =  ekk_rowstat(model);
35  const double * rowSolution = ekk_rowacts(model);
36  int i;
37  clp->createStatus();
38  double * clpSolution;
39  clpSolution = clp->primalRowSolution();
40  memcpy(clpSolution,rowSolution,numberRows*sizeof(double));
41  const double * rowLower = ekk_rowlower(model);
42  const double * rowUpper = ekk_rowupper(model);
43  for (i=0;i<numberRows;i++) {
44    ClpSimplex::Status status;
45    if ((rowStatus[i]&0x80000000)!=0) {
46      status = ClpSimplex::basic;
47    } else {
48      if (!startup) {
49        // believe bits
50        int ikey = rowStatus[i]&0x60000000;
51        if (ikey==0x40000000) {
52          // at ub
53          status = ClpSimplex::atUpperBound;
54          clpSolution[i]=rowUpper[i];
55        } else if (ikey==0x20000000) {
56          // at lb
57          status = ClpSimplex::atLowerBound;
58          clpSolution[i]=rowLower[i];
59        } else if (ikey==0x60000000) {
60          // free
61          status = ClpSimplex::isFree;
62          clpSolution[i]=0.0;
63        } else {
64          // fixed
65          status = ClpSimplex::atLowerBound;
66          clpSolution[i]=rowLower[i];
67        }
68      } else {
69        status = ClpSimplex::superBasic;
70      }
71    }
72    clp->setRowStatus(i,status);
73  }
74 
75  const int * columnStatus = ekk_colstat(model);
76  const double * columnSolution =  ekk_colsol(model);
77  clpSolution = clp->primalColumnSolution();
78  memcpy(clpSolution,columnSolution,numberColumns*sizeof(double));
79  const double * columnLower = ekk_collower(model);
80  const double * columnUpper = ekk_colupper(model);
81  for (i=0;i<numberColumns;i++) {
82    ClpSimplex::Status status;
83    if ((columnStatus[i]&0x80000000)!=0) {
84      status = ClpSimplex::basic;
85    } else {
86      if (!startup) {
87        // believe bits
88        int ikey = columnStatus[i]&0x60000000;
89        if (ikey==0x40000000) {
90          // at ub
91          status = ClpSimplex::atUpperBound;
92          clpSolution[i]=columnUpper[i];
93        } else if (ikey==0x20000000) {
94          // at lb
95          status = ClpSimplex::atLowerBound;
96          clpSolution[i]=columnLower[i];
97        } else if (ikey==0x60000000) {
98          // free
99          status = ClpSimplex::isFree;
100          clpSolution[i]=0.0;
101        } else {
102          // fixed
103          status = ClpSimplex::atLowerBound;
104          clpSolution[i]=columnLower[i];
105        }
106      } else {
107        status = ClpSimplex::superBasic;
108      }
109    }
110    clp->setColumnStatus(i,status);
111  }
112  return clp;
113}
114
115static int solve(EKKModel * model, int  startup,int algorithm,
116                 int presolve)
117{
118  // values pass or not
119  if (startup)
120    startup=1;
121  // if scaled then be careful
122  bool scaled = ekk_scaling(model)==1;
123  if (scaled)
124    ekk_scaleRim(model,1);
125  void * compressInfo=NULL;
126  ClpSimplex * clp;
127  if (!presolve||!presolveInfo) {
128    // no presolve or osl presolve - compact columns
129    compressInfo=ekk_compressModel(model);
130    clp = clpmodel(model,startup);;
131  } else {
132    // pick up clp model
133    clp = presolveInfo->model();
134  }
135   
136  // don't scale if alreday scaled
137  if (scaled)
138    clp->scaling(false);
139  if (clp->numberRows()>10000)
140    clp->factorization()->maximumPivots(100+clp->numberRows()/100);
141  if (algorithm>0)
142    clp->primal(startup);
143  else
144    clp->dual();
145
146  int numberIterations = clp->numberIterations();
147  if (presolve&&presolveInfo) {
148    // very wasteful - create a clp copy of osl model
149    ClpSimplex * clpOriginal = clpmodel(model,0);
150    presolveInfo->setOriginalModel(clpOriginal);
151    // do postsolve
152    presolveInfo->postsolve(true);
153    delete clp;
154    delete presolveInfo;
155    presolveInfo=NULL;
156    clp = clpOriginal;
157    if (presolve==3||(presolve==2&&clp->status())) {
158      printf("Resolving from postsolved model\n");
159      clp->primal(1);
160      numberIterations += clp->numberIterations();
161    }
162  }
163
164  // put back solution
165
166  double * rowDual = (double *) ekk_rowduals(model);
167  int numberRows=ekk_getInumrows(model);
168  int numberColumns= ekk_getInumcols(model);
169  int * rowStatus =  (int *) ekk_rowstat(model);
170  double * rowSolution = (double *) ekk_rowacts(model);
171  int i;
172  int * columnStatus = (int *) ekk_colstat(model);
173  double * columnSolution =  (double *) ekk_colsol(model);
174  memcpy(rowSolution,clp->primalRowSolution(),numberRows*sizeof(double));
175  memcpy(rowDual,clp->dualRowSolution(),numberRows*sizeof(double));
176  for (i=0;i<numberRows;i++) {
177    if (clp->getRowStatus(i)==ClpSimplex::basic)
178      rowStatus[i]=0x80000000;
179    else
180      rowStatus[i]=0;
181  }
182 
183  double * columnDual = (double *) ekk_colrcosts(model);
184  memcpy(columnSolution,clp->primalColumnSolution(),
185         numberColumns*sizeof(double));
186  memcpy(columnDual,clp->dualColumnSolution(),numberColumns*sizeof(double));
187  for (i=0;i<numberColumns;i++) {
188    if (clp->getColumnStatus(i)==ClpSimplex::basic)
189      columnStatus[i]=0x80000000;
190    else
191      columnStatus[i]=0;
192  }
193  ekk_setIprobstat(model,clp->status());
194  ekk_setRobjvalue(model,clp->objectiveValue());
195  ekk_setInumpinf(model,clp->numberPrimalInfeasibilities());
196  ekk_setInumdinf(model,clp->numberDualInfeasibilities());
197  ekk_setIiternum(model,numberIterations);
198  ekk_setRsumpinf(model,clp->sumPrimalInfeasibilities());
199  ekk_setRsumdinf(model,clp->sumDualInfeasibilities());
200  delete clp;
201  if (compressInfo ) 
202    ekk_decompressModel(model,compressInfo);
203
204  if (scaled)
205      ekk_scaleRim(model,2);
206  return 0;
207}
208/* As ekk_primalSimplex + postsolve instructions:
209   presolve - 0 , no presolve, 1 presolve but no primal after postsolve,
210                  2 do primal if any infeasibilities,
211                  3 always do primal.
212*/
213extern "C" int ekk_primalClp(EKKModel * model, int  startup, int presolve)
214{
215  if (presolveInfo) {
216    if (!presolve)
217      presolve=3;
218    return solve(model,startup,1, presolve);
219  } else {
220    return solve(model,startup,1, 0);
221  }
222}
223
224/* As ekk_dualSimplex + postsolve instructions:
225   presolve - 0 , no presolve, 1 presolve but no primal after postsolve,
226                  2 do primal if any infeasibilities,
227                  3 always do primal.
228*/
229extern "C" int ekk_dualClp(EKKModel * model, int presolve)
230{
231  if (presolveInfo) {
232    if (!presolve)
233      presolve=3;
234    return solve(model,0,-1, presolve);
235  } else {
236    return solve(model,0,-1, 0);
237  }
238}
239/* rather like ekk_preSolve (3) plus:
240   keepIntegers - false to treat as if continuous
241   pass   - do this many passes (0==default(5))
242
243   returns 1 if infeasible
244*/
245extern "C" int ekk_preSolveClp(EKKModel * model, bool keepIntegers,
246                               int pass)
247{
248  delete presolveInfo;
249  presolveInfo = new ClpPresolve();
250  bool scaled = ekk_scaling(model)==1;
251  if (scaled)
252      ekk_scaleRim(model,1);
253  if (!pass)
254    pass=5;
255  // very wasteful - create a clp copy of osl model
256  // 1 to keep solution as is
257  ClpSimplex * clp = clpmodel(model,1);
258  // could get round with tailored version of ClpPresolve.cpp
259  ClpSimplex * newModel = 
260    presolveInfo->presolvedModel(*clp,
261                                 ekk_getRtolpinf(model),
262                                 keepIntegers,pass,true);
263  delete clp;
264  presolveInfo->setOriginalModel(NULL);
265  if (scaled)
266      ekk_scaleRim(model,2);
267  if (newModel) {
268    return 0;
269  } else {
270    delete presolveInfo;
271    presolveInfo=NULL;
272    return 1;
273  }
274}
Note: See TracBrowser for help on using the repository browser.