source: trunk/Samples/testGub2.cpp @ 343

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

dualRanging

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