source: stable/1.13/Clp/examples/ekk_interface.cpp @ 1898

Last change on this file since 1898 was 1559, checked in by stefan, 10 years ago

merge split branch into trunk

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