source: stable/BSPsplit/Clp/examples/testGub2.cpp @ 1465

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

use new pkg-config based macros to search for required and optional packages

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 12.3 KB
Line 
1// Copyright (C) 2003, International Business Machines
2// Corporation and others.  All Rights Reserved.
3
4#include "ClpSimplex.hpp"
5#include "ClpGubDynamicMatrix.hpp"
6#include "ClpPrimalColumnSteepest.hpp"
7#include "CoinSort.hpp"
8#include "CoinHelperFunctions.hpp"
9#include "CoinTime.hpp"
10#include "CoinMpsIO.hpp"
11int main (int argc, const char *argv[])
12{
13  ClpSimplex  model;
14  int status;
15  int maxIts=0;
16  int maxFactor=100;
17  if (argc<2) {
18#if defined(COIN_HAS_SAMPLE) && defined(SAMPLEDIR)
19    status=model.readMps(SAMPLEDIR "/p0033.mps",true);
20#else
21    fprintf(stderr,"Do not know where to find sample MPS files.\n");
22    exit(1);
23#endif
24  } else
25    status=model.readMps(argv[1]);
26  if (status) {
27    printf("errors on input\n");
28    exit(77);
29  }
30  if (argc>2) {
31    maxFactor = atoi(argv[2]);
32    printf("max factor %d\n",maxFactor);
33  }
34  if (argc>3) {
35    maxIts = atoi(argv[3]);
36    printf("max its %d\n",maxIts);
37  }
38  // For now scaling off
39  model.scaling(0);
40  if (maxIts) {
41    // Do partial dantzig
42    ClpPrimalColumnSteepest dantzig(5);
43    model.setPrimalColumnPivotAlgorithm(dantzig);
44    //model.messageHandler()->setLogLevel(63);
45    model.setFactorizationFrequency(maxFactor);
46    model.setMaximumIterations(maxIts);
47    model.primal();
48    if (!model.status())
49      exit(1);
50  }
51  // find gub
52  int numberRows = model.numberRows();
53  int * gubStart = new int[numberRows+1];
54  int * gubEnd = new int[numberRows];
55  int * which = new int[numberRows];
56  int * whichGub = new int[numberRows];
57  int numberColumns = model.numberColumns();
58  int * mark = new int[numberColumns];
59  int iRow,iColumn;
60  // delete variables fixed to zero
61  const double * columnLower = model.columnLower();
62  const double * columnUpper = model.columnUpper();
63  int numberDelete=0;
64  for (iColumn=0;iColumn<numberColumns;iColumn++) {
65    if (columnUpper[iColumn]==0.0&&columnLower[iColumn]==0.0)
66      mark[numberDelete++]=iColumn;
67  }
68  if (numberDelete) {
69    model.deleteColumns(numberDelete,mark);
70    numberColumns -= numberDelete;
71    columnLower = model.columnLower();
72    columnUpper = model.columnUpper();
73#if 0
74    CoinMpsIO writer;
75    writer.setMpsData(*model.matrix(), COIN_DBL_MAX,
76                      model.getColLower(), model.getColUpper(),
77                      model.getObjCoefficients(),
78                      (const char*) 0 /*integrality*/,
79                      model.getRowLower(), model.getRowUpper(),
80                      NULL, NULL);
81    writer.writeMps("cza.mps",0,0,1);
82#endif
83  }
84  double * lower = new double[numberRows];
85  double * upper = new double[numberRows];
86  const double * rowLower = model.rowLower();
87  const double * rowUpper = model.rowUpper();
88  for (iColumn=0;iColumn<numberColumns;iColumn++)
89    mark[iColumn]=-1;
90  CoinPackedMatrix * matrix = model.matrix();
91  // get row copy
92  CoinPackedMatrix rowCopy = *matrix;
93  rowCopy.reverseOrdering();
94  const int * column = rowCopy.getIndices();
95  const int * rowLength = rowCopy.getVectorLengths();
96  const CoinBigIndex * rowStart = rowCopy.getVectorStarts();
97  const double * element = rowCopy.getElements();
98  int putGub=numberRows;
99  int putNonGub=numberRows;
100  int * rowIsGub = new int [numberRows];
101  for (iRow=numberRows-1;iRow>=0;iRow--) {
102    bool gubRow=true;
103    int first=numberColumns+1;
104    int last=-1;
105    for (int j=rowStart[iRow];j<rowStart[iRow]+rowLength[iRow];j++) {
106      if (element[j]!=1.0) {
107        gubRow=false;
108        break;
109      } else {
110        int iColumn = column[j];
111        if (mark[iColumn]>=0) {
112          gubRow=false;
113          break;
114        } else {
115          last=CoinMax(last,iColumn);
116          first = CoinMin(first,iColumn);
117        }
118      }
119    }
120    if (last-first+1!=rowLength[iRow]||!gubRow) {
121      which[--putNonGub]=iRow;
122      rowIsGub[iRow]=0;
123    } else {
124      for (int j=rowStart[iRow];j<rowStart[iRow]+rowLength[iRow];j++) {
125        int iColumn = column[j];
126        mark[iColumn]=iRow;
127      }
128      rowIsGub[iRow]=-1;
129      putGub--;
130      gubStart[putGub]=first;
131      gubEnd[putGub]=last+1;
132      lower[putGub]=rowLower[iRow];
133      upper[putGub]=rowUpper[iRow];
134      whichGub[putGub]=iRow;
135    }
136  }
137  int numberNonGub=numberRows-putNonGub;
138  int numberGub = numberRows-putGub;
139  if (numberGub>0) {
140    printf("** %d gub rows\n",numberGub);
141    int numberNormal = 0;
142    const int * row = matrix->getIndices();
143    const int * columnLength = matrix->getVectorLengths();
144    const CoinBigIndex * columnStart = matrix->getVectorStarts();
145    const double * elementByColumn = matrix->getElements();
146    int numberElements = 0;
147    bool doLower=false;
148    bool doUpper=false;
149    for (iColumn=0;iColumn<numberColumns;iColumn++) {
150      if (mark[iColumn]<0) {
151        mark[numberNormal++]=iColumn;
152      } else {
153        numberElements += columnLength[iColumn];
154        if (columnLower[iColumn]!=0.0)
155          doLower=true;
156        if (columnUpper[iColumn]<1.0e20)
157          doUpper=true;
158      }
159    }
160    if (!numberNormal) {
161      printf("Putting back one gub row to make non-empty\n");
162      for (iColumn=gubStart[putGub];iColumn< gubEnd[putGub];iColumn++)
163        mark[numberNormal++]=iColumn;
164      putGub++;
165      numberGub--;
166    }
167    ClpSimplex model2(&model,numberNonGub,which+putNonGub,numberNormal,mark);
168    int numberGubColumns=numberColumns-numberNormal;
169    // sort gubs so monotonic
170    int * which = new int[numberGub];
171    int i;
172    for (i=0;i<numberGub;i++)
173      which[i]=i;
174    CoinSort_2(gubStart+putGub,gubStart+putGub+numberGub,which);
175    int * temp1 = new int [numberGub];
176    for (i=0;i<numberGub;i++) {
177      int k=which[i];
178      temp1[i]=gubEnd[putGub+k];
179    }
180    memcpy(gubEnd+putGub,temp1,numberGub*sizeof(int));
181    delete [] temp1;
182    double * temp2 = new double [numberGub];
183    for (i=0;i<numberGub;i++) {
184      int k=which[i];
185      temp2[i]=lower[putGub+k];
186    }
187    memcpy(lower+putGub,temp2,numberGub*sizeof(double));
188    for (i=0;i<numberGub;i++) {
189      int k=which[i];
190      temp2[i]=upper[putGub+k];
191    }
192    memcpy(upper+putGub,temp2,numberGub*sizeof(double));
193    delete [] temp2;
194    delete [] which;
195    numberElements -= numberGubColumns;
196    int * start2 = new int[numberGubColumns+1];
197    int * row2 = new int[numberElements];
198    double * element2 = new double[numberElements];
199    double * cost2 = new double [numberGubColumns];
200    double * lowerColumn2=NULL;
201    if (doLower) {
202      lowerColumn2 = new double [numberGubColumns];
203      CoinFillN(lowerColumn2,numberGubColumns,0.0);
204    }
205    double * upperColumn2=NULL;
206    if (doUpper) {
207      upperColumn2 = new double [numberGubColumns];
208      CoinFillN(upperColumn2,numberGubColumns,COIN_DBL_MAX);
209    }
210    numberElements = 0;
211    int numberNonGubRows=0;
212    for (iRow=0;iRow<numberRows;iRow++) {
213      if (!rowIsGub[iRow])
214        rowIsGub[iRow]= numberNonGubRows++;
215    }
216    numberColumns=0;
217    gubStart[0]=0;
218    start2[0]=0;
219    const double * cost = model.objective();
220    for (int iSet=0;iSet<numberGub;iSet++) {
221      int iStart = gubStart[iSet+putGub];
222      int iEnd = gubEnd[iSet+putGub];
223      for (int k=iStart;k<iEnd;k++) {
224        cost2[numberColumns]=cost[k];
225        if (columnLower[k])
226          lowerColumn2[numberColumns]=columnLower[k];
227        if (columnUpper[k]<1.0e20)
228          upperColumn2[numberColumns]=columnUpper[k];
229        for (int j=columnStart[k];j<columnStart[k]+columnLength[k];j++) {
230          int iRow = rowIsGub[row[j]];
231          if (iRow>=0) {
232            row2[numberElements]=iRow;
233            element2[numberElements++]=elementByColumn[j];
234          }
235        }
236        start2[++numberColumns]=numberElements;
237      }
238      gubStart[iSet+1]=numberColumns;
239    }
240    model2.replaceMatrix(new ClpGubDynamicMatrix(&model2,numberGub,
241                                                 numberColumns,gubStart,
242                                                 lower+putGub,upper+putGub,
243                                                 start2,row2,element2,cost2,
244                                                 lowerColumn2,upperColumn2));
245    delete [] rowIsGub;
246    delete [] start2;
247    delete [] row2;
248    delete [] element2;
249    delete [] cost2;
250    delete [] lowerColumn2;
251    delete [] upperColumn2;
252    // For now scaling off
253    model2.scaling(0);
254    // Do partial dantzig
255    ClpPrimalColumnSteepest dantzig(5);
256    model2.setPrimalColumnPivotAlgorithm(dantzig);
257    //model2.messageHandler()->setLogLevel(63);
258    model2.setFactorizationFrequency(maxFactor);
259    model2.setMaximumIterations(4000000);
260    double time1 = CoinCpuTime();
261    model2.primal();
262    {
263      ClpGubDynamicMatrix * gubMatrix =
264        dynamic_cast< ClpGubDynamicMatrix*>(model2.clpMatrix());
265      assert (gubMatrix);
266      const double * solution = model2.primalColumnSolution();
267      int numberGubColumns = gubMatrix->numberGubColumns();
268      int firstOdd = gubMatrix->firstDynamic();
269      int lastOdd =gubMatrix->firstAvailable();
270      int numberTotalColumns = firstOdd+numberGubColumns;
271      int numberRows = model2.numberRows();
272      char * status = new char [numberTotalColumns];
273      double * gubSolution = new double [numberTotalColumns];
274      int numberSets = gubMatrix->numberSets();
275      const int * id = gubMatrix->id();
276      int i;
277      const double * lowerColumn = gubMatrix->lowerColumn();
278      const double * upperColumn = gubMatrix->upperColumn();
279      for (i=0;i<numberGubColumns;i++) {
280        if (gubMatrix->getDynamicStatus(i)==ClpGubDynamicMatrix::atUpperBound) {
281          gubSolution[i+firstOdd]=upperColumn[i];
282          status[i+firstOdd]=2;
283        } else if (gubMatrix->getDynamicStatus(i)==ClpGubDynamicMatrix::atLowerBound&&lowerColumn) {
284          gubSolution[i+firstOdd]=lowerColumn[i];
285          status[i+firstOdd]=1;
286        } else {
287          gubSolution[i+firstOdd]=0.0;
288          status[i+firstOdd]=1;
289        }
290      }
291      for (i=0;i<firstOdd;i++) {
292        ClpSimplex::Status thisStatus = model2.getStatus(i);
293        if (thisStatus==ClpSimplex::basic)
294          status[i]=0;
295        else if (thisStatus==ClpSimplex::atLowerBound)
296          status[i]=1;
297        else if (thisStatus==ClpSimplex::atUpperBound)
298          status[i]=2;
299        else if (thisStatus==ClpSimplex::isFixed)
300          status[i]=3;
301        else
302          abort();
303        gubSolution[i]=solution[i];
304      }
305      for (i=firstOdd;i<lastOdd;i++) {
306        int iBig = id[i-firstOdd]+firstOdd;
307        ClpSimplex::Status thisStatus = model2.getStatus(i);
308        if (thisStatus==ClpSimplex::basic)
309          status[iBig]=0;
310        else if (thisStatus==ClpSimplex::atLowerBound)
311          status[iBig]=1;
312        else if (thisStatus==ClpSimplex::atUpperBound)
313          status[iBig]=2;
314        else if (thisStatus==ClpSimplex::isFixed)
315          status[iBig]=3;
316        else
317          abort();
318        gubSolution[iBig]=solution[i];
319      }
320      char * rowStatus = new char[numberRows];
321      for (i=0;i<numberRows;i++) {
322        ClpSimplex::Status thisStatus = model2.getRowStatus(i);
323        if (thisStatus==ClpSimplex::basic)
324          rowStatus[i]=0;
325        else if (thisStatus==ClpSimplex::atLowerBound)
326          rowStatus[i]=1;
327        else if (thisStatus==ClpSimplex::atUpperBound)
328          rowStatus[i]=2;
329        else if (thisStatus==ClpSimplex::isFixed)
330          rowStatus[i]=3;
331        else
332          abort();
333      }
334      char * setStatus = new char[numberSets];
335      int * keyVariable = new int[numberSets];
336      memcpy(keyVariable,gubMatrix->keyVariable(),numberSets*sizeof(int));
337      for (i=0;i<numberSets;i++) {
338        int iKey = keyVariable[i];
339        if (iKey>lastOdd)
340          iKey=numberTotalColumns+i;
341        else
342          iKey = id[iKey-firstOdd]+firstOdd;
343        keyVariable[i]=iKey;
344        ClpSimplex::Status thisStatus = gubMatrix->getStatus(i);
345        if (thisStatus==ClpSimplex::basic)
346          setStatus[i]=0;
347        else if (thisStatus==ClpSimplex::atLowerBound)
348          setStatus[i]=1;
349        else if (thisStatus==ClpSimplex::atUpperBound)
350          setStatus[i]=2;
351        else if (thisStatus==ClpSimplex::isFixed)
352          setStatus[i]=3;
353        else
354          abort();
355      }
356      FILE * fp=fopen ("xx.sol","w");
357      fwrite(gubSolution,sizeof(double),numberTotalColumns,fp);
358      fwrite(status,sizeof(char),numberTotalColumns,fp);
359      const double * rowsol = model2.primalRowSolution();
360      int originalNumberRows = model.numberRows();
361      double * rowsol2 = new double[originalNumberRows];
362      memset(rowsol2,0,originalNumberRows*sizeof(double));
363      model.times(1.0,gubSolution,rowsol2);
364      for (i=0;i<numberRows;i++) 
365        assert (fabs(rowsol[i]-rowsol2[i])<1.0e-3);
366      //for (;i<originalNumberRows;i++)
367      //printf("%d %g\n",i,rowsol2[i]);
368      delete [] rowsol2;
369      fwrite(rowsol,sizeof(double),numberRows,fp);
370      fwrite(rowStatus,sizeof(char),numberRows,fp);
371      fwrite(setStatus,sizeof(char),numberSets,fp);
372      fwrite(keyVariable,sizeof(int),numberSets,fp);
373      fclose(fp);
374      delete [] status;
375      delete [] gubSolution;
376      delete [] setStatus;
377      delete [] keyVariable;
378      // ** if going to rstart as dynamic need id_
379      // also copy coding in useEf.. from ClpGubMatrix (i.e. test for basis)
380    }
381    printf("obj offset is %g\n",model2.objectiveOffset());
382    printf("Primal took %g seconds\n",CoinCpuTime()-time1);
383    //model2.primal(1);
384  }
385  delete [] mark;
386  delete [] gubStart;
387  delete [] gubEnd;
388  delete [] which;
389  delete [] whichGub;
390  delete [] lower;
391  delete [] upper;
392  return 0;
393}   
Note: See TracBrowser for help on using the repository browser.