source: trunk/Clp/examples/ekk_interface.cpp

Last change on this file was 1662, checked in by lou, 7 years ago

Add EPL license notice in examples.

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