source: branches/devel-1/Samples/ekk_interface.cpp @ 38

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

New "example" linking osl to ekk

  • 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 "Presolve.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 Presolve * 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    clp = clpOriginal;
156    if (presolve==3||(presolve==2&&clp->status())) {
157      printf("Resolving from postsolved model\n");
158      clp->primal(1);
159      numberIterations += clp->numberIterations();
160    }
161  }
162
163  // put back solution
164
165  double * rowDual = (double *) ekk_rowduals(model);
166  int numberRows=ekk_getInumrows(model);
167  int numberColumns= ekk_getInumcols(model);
168  int * rowStatus =  (int *) ekk_rowstat(model);
169  double * rowSolution = (double *) ekk_rowacts(model);
170  int i;
171  int * columnStatus = (int *) ekk_colstat(model);
172  double * columnSolution =  (double *) ekk_colsol(model);
173  memcpy(rowSolution,clp->primalRowSolution(),numberRows*sizeof(double));
174  memcpy(rowDual,clp->dualRowSolution(),numberRows*sizeof(double));
175  for (i=0;i<numberRows;i++) {
176    if (clp->getRowStatus(i)==ClpSimplex::basic)
177      rowStatus[i]=0x80000000;
178    else
179      rowStatus[i]=0;
180  }
181 
182  double * columnDual = (double *) ekk_colrcosts(model);
183  memcpy(columnSolution,clp->primalColumnSolution(),
184         numberColumns*sizeof(double));
185  memcpy(columnDual,clp->dualColumnSolution(),numberColumns*sizeof(double));
186  for (i=0;i<numberColumns;i++) {
187    if (clp->getColumnStatus(i)==ClpSimplex::basic)
188      columnStatus[i]=0x80000000;
189    else
190      columnStatus[i]=0;
191  }
192  ekk_setIprobstat(model,clp->status());
193  ekk_setRobjvalue(model,clp->objectiveValue());
194  ekk_setInumpinf(model,clp->numberPrimalInfeasibilities());
195  ekk_setInumdinf(model,clp->numberDualInfeasibilities());
196  ekk_setIiternum(model,numberIterations);
197  ekk_setRsumpinf(model,clp->sumPrimalInfeasibilities());
198  ekk_setRsumdinf(model,clp->sumDualInfeasibilities());
199  delete clp;
200  if (compressInfo ) 
201    ekk_decompressModel(model,compressInfo);
202
203  if (scaled)
204      ekk_scaleRim(model,2);
205  return 0;
206}
207/* As ekk_primalSimplex + postsolve instructions:
208   presolve - 0 , no presolve, 1 presolve but no primal after postsolve,
209                  2 do primal if any infeasibilities,
210                  3 always do primal.
211*/
212extern "C" int ekk_primalClp(EKKModel * model, int  startup, int presolve)
213{
214  if (presolveInfo) {
215    if (!presolve)
216      presolve=3;
217    return solve(model,startup,1, presolve);
218  } else {
219    return solve(model,startup,1, 0);
220  }
221}
222
223/* As ekk_dualSimplex + postsolve instructions:
224   presolve - 0 , no presolve, 1 presolve but no primal after postsolve,
225                  2 do primal if any infeasibilities,
226                  3 always do primal.
227*/
228extern "C" int ekk_dualClp(EKKModel * model, int presolve)
229{
230  if (presolveInfo) {
231    if (!presolve)
232      presolve=3;
233    return solve(model,0,-1, presolve);
234  } else {
235    return solve(model,0,-1, 0);
236  }
237}
238/* rather like ekk_preSolve (3) plus:
239   keepIntegers - false to treat as if continuous
240   pass   - do this many passes (0==default(5))
241
242   returns 1 if infeasible
243*/
244extern "C" int ekk_preSolveClp(EKKModel * model, bool keepIntegers,
245                               int pass)
246{
247  delete presolveInfo;
248  presolveInfo = new Presolve();
249  bool scaled = ekk_scaling(model)==1;
250  if (scaled)
251      ekk_scaleRim(model,1);
252  if (!pass)
253    pass=5;
254  // very wasteful - create a clp copy of osl model
255  ClpSimplex * clp = clpmodel(model,0);
256  // could get round with tailored version of Presolve.cpp
257  ClpSimplex * newModel = presolveInfo->presolvedModel(*clp,
258                                                       ekk_getRtolpinf(model),
259                                                       keepIntegers,pass);
260  delete clp;
261  presolveInfo->setOriginalModel(NULL);
262  if (scaled)
263      ekk_scaleRim(model,2);
264  if (newModel) {
265    return 0;
266  } else {
267    delete presolveInfo;
268    presolveInfo=NULL;
269    return 1;
270  }
271}
Note: See TracBrowser for help on using the repository browser.