source: trunk/Samples/testGub.cpp @ 493

Last change on this file since 493 was 493, checked in by forrest, 15 years ago

czprob

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 14.9 KB
Line 
1// Copyright (C) 2003, International Business Machines
2// Corporation and others.  All Rights Reserved.
3
4#include "ClpSimplex.hpp"
5#include "ClpDynamicExampleMatrix.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 maxFactor=100;
16  if (argc<2) {
17    status=model.readMps("../../Mps/Netlib/czprob.mps");
18    if (status) {
19      printf("Unable to read matrix - trying gzipped version\n");
20      status=model.readMps("../../Mps/Netlib/czprob.mps.gz");
21    }
22  } else {
23    status=model.readMps(argv[1]);
24  }
25  if (status) {
26    printf("errors on input\n");
27    exit(77);
28  }
29  if (argc>2) {
30    maxFactor = atoi(argv[2]);
31    printf("max factor %d\n",maxFactor);
32  }
33  if (argc>3) {
34    printf("Using ClpDynamicMatrix\n");
35  } else {
36    printf("Using ClpDynamicExampleMatrix\n");
37  }
38  // find gub
39  int numberRows = model.numberRows();
40  int * gubStart = new int[numberRows+1];
41  int * gubEnd = new int[numberRows];
42  int * which = new int[numberRows];
43  int * whichGub = new int[numberRows];
44  int numberColumns = model.numberColumns();
45  int * mark = new int[numberColumns];
46  int iRow,iColumn;
47  // delete variables fixed to zero
48  const double * columnLower = model.columnLower();
49  const double * columnUpper = model.columnUpper();
50  int numberDelete=0;
51  for (iColumn=0;iColumn<numberColumns;iColumn++) {
52    if (columnUpper[iColumn]==0.0&&columnLower[iColumn]==0.0)
53      mark[numberDelete++]=iColumn;
54  }
55  if (numberDelete) {
56    model.deleteColumns(numberDelete,mark);
57    numberColumns -= numberDelete;
58    columnLower = model.columnLower();
59    columnUpper = model.columnUpper();
60  }
61  double * lower = new double[numberRows];
62  double * upper = new double[numberRows];
63  const double * rowLower = model.rowLower();
64  const double * rowUpper = model.rowUpper();
65  for (iColumn=0;iColumn<numberColumns;iColumn++)
66    mark[iColumn]=-1;
67  CoinPackedMatrix * matrix = model.matrix();
68  // get row copy
69  CoinPackedMatrix rowCopy = *matrix;
70  rowCopy.reverseOrdering();
71  const int * column = rowCopy.getIndices();
72  const int * rowLength = rowCopy.getVectorLengths();
73  const CoinBigIndex * rowStart = rowCopy.getVectorStarts();
74  const double * element = rowCopy.getElements();
75  int putGub=numberRows;
76  int putNonGub=numberRows;
77  int * rowIsGub = new int [numberRows];
78  for (iRow=numberRows-1;iRow>=0;iRow--) {
79    bool gubRow=true;
80    int first=numberColumns+1;
81    int last=-1;
82    for (int j=rowStart[iRow];j<rowStart[iRow]+rowLength[iRow];j++) {
83      if (element[j]!=1.0) {
84        gubRow=false;
85        break;
86      } else {
87        int iColumn = column[j];
88        if (mark[iColumn]>=0) {
89          gubRow=false;
90          break;
91        } else {
92          last=max(last,iColumn);
93          first = min(first,iColumn);
94        }
95      }
96    }
97    if (last-first+1!=rowLength[iRow]||!gubRow) {
98      which[--putNonGub]=iRow;
99      rowIsGub[iRow]=0;
100    } else {
101      for (int j=rowStart[iRow];j<rowStart[iRow]+rowLength[iRow];j++) {
102        int iColumn = column[j];
103        mark[iColumn]=iRow;
104      }
105      rowIsGub[iRow]=-1;
106      putGub--;
107      gubStart[putGub]=first;
108      gubEnd[putGub]=last+1;
109      lower[putGub]=rowLower[iRow];
110      upper[putGub]=rowUpper[iRow];
111      whichGub[putGub]=iRow;
112    }
113  }
114  int numberNonGub=numberRows-putNonGub;
115  int numberGub = numberRows-putGub;
116  if (numberGub>0) {
117    printf("** %d gub rows\n",numberGub);
118    int numberNormal = 0;
119    const int * row = matrix->getIndices();
120    const int * columnLength = matrix->getVectorLengths();
121    const CoinBigIndex * columnStart = matrix->getVectorStarts();
122    const double * elementByColumn = matrix->getElements();
123    int numberElements = 0;
124    bool doLower=false;
125    bool doUpper=false;
126    for (iColumn=0;iColumn<numberColumns;iColumn++) {
127      if (mark[iColumn]<0) {
128        mark[numberNormal++]=iColumn;
129      } else {
130        numberElements += columnLength[iColumn];
131        if (columnLower[iColumn]!=0.0)
132          doLower=true;
133        if (columnUpper[iColumn]<1.0e20)
134          doUpper=true;
135      }
136    }
137    if (!numberNormal) {
138      printf("Putting back one gub row to make non-empty\n");
139      for (iColumn=gubStart[putGub];iColumn< gubEnd[putGub];iColumn++)
140        mark[numberNormal++]=iColumn;
141      putGub++;
142      numberGub--;
143    }
144    ClpSimplex model2(&model,numberNonGub,which+putNonGub,numberNormal,mark);
145    // and copy for restart test
146    ClpSimplex model3 = model2;
147    int numberGubColumns=numberColumns-numberNormal;
148    // sort gubs so monotonic
149    int * which = new int[numberGub];
150    int i;
151    for (i=0;i<numberGub;i++)
152      which[i]=i;
153    CoinSort_2(gubStart+putGub,gubStart+putGub+numberGub,which);
154    // move to bottom if we want to use later
155    memmove(gubStart,gubStart+putGub,numberGub*sizeof(int));
156    int * temp1 = new int [numberGub];
157    for (i=0;i<numberGub;i++) {
158      int k=which[i];
159      temp1[i]=gubEnd[putGub+k];
160    }
161    memcpy(gubEnd,temp1,numberGub*sizeof(int));
162    delete [] temp1;
163    double * temp2 = new double [numberGub];
164    for (i=0;i<numberGub;i++) {
165      int k=which[i];
166      temp2[i]=lower[putGub+k];
167    }
168    memcpy(lower,temp2,numberGub*sizeof(double));
169    for (i=0;i<numberGub;i++) {
170      int k=which[i];
171      temp2[i]=upper[putGub+k];
172    }
173    memcpy(upper,temp2,numberGub*sizeof(double));
174    delete [] temp2;
175    delete [] which;
176    numberElements -= numberGubColumns;
177    int * start2 = new int[numberGubColumns+1];
178    int * row2 = new int[numberElements];
179    double * element2 = new double[numberElements];
180    double * cost2 = new double [numberGubColumns];
181    double * lowerColumn2=NULL;
182    if (doLower) {
183      lowerColumn2 = new double [numberGubColumns];
184      CoinFillN(lowerColumn2,numberGubColumns,0.0);
185    }
186    double * upperColumn2=NULL;
187    if (doUpper) {
188      upperColumn2 = new double [numberGubColumns];
189      CoinFillN(upperColumn2,numberGubColumns,COIN_DBL_MAX);
190    }
191    numberElements = 0;
192    int numberNonGubRows=0;
193    for (iRow=0;iRow<numberRows;iRow++) {
194      if (!rowIsGub[iRow])
195        rowIsGub[iRow]= numberNonGubRows++;
196    }
197    numberColumns=0;
198    int iStart = gubStart[0];
199    gubStart[0]=0;
200    start2[0]=0;
201    const double * cost = model.objective();
202    for (int iSet=0;iSet<numberGub;iSet++) {
203      int iEnd = gubEnd[iSet];
204      for (int k=iStart;k<iEnd;k++) {
205        cost2[numberColumns]=cost[k];
206        if (columnLower[k])
207          lowerColumn2[numberColumns]=columnLower[k];
208        if (columnUpper[k]<1.0e20)
209          upperColumn2[numberColumns]=columnUpper[k];
210        for (int j=columnStart[k];j<columnStart[k]+columnLength[k];j++) {
211          int iRow = rowIsGub[row[j]];
212          if (iRow>=0) {
213            row2[numberElements]=iRow;
214            element2[numberElements++]=elementByColumn[j];
215          }
216        }
217        start2[++numberColumns]=numberElements;
218      }
219      if (iSet<numberGub-1)
220        iStart = gubStart[iSet+1];
221      gubStart[iSet+1]=numberColumns;
222    }
223    printf("** Before adding matrix there are %d rows and %d columns\n",
224           model2.numberRows(),model2.numberColumns());
225    if (argc>3) {
226      ClpDynamicMatrix * newMatrix = new ClpDynamicMatrix(&model2,numberGub,
227                                                                        numberColumns,gubStart,
228                                                                        lower,upper,
229                                                                        start2,row2,element2,cost2,
230                                                                        lowerColumn2,upperColumn2);
231      model2.replaceMatrix(newMatrix);
232      newMatrix->switchOffCheck();
233      newMatrix->setRefreshFrequency(1000);
234    } else {
235      ClpDynamicExampleMatrix * newMatrix = new ClpDynamicExampleMatrix(&model2,numberGub,
236                                                                        numberColumns,gubStart,
237                                                                        lower,upper,
238                                                                        start2,row2,element2,cost2,
239                                                                        lowerColumn2,upperColumn2);
240      model2.replaceMatrix(newMatrix);
241      newMatrix->switchOffCheck();
242      newMatrix->setRefreshFrequency(1000);
243    }
244    printf("** While after adding matrix there are %d rows and %d columns\n",
245           model2.numberRows(),model2.numberColumns());
246    model2.setSpecialOptions(4); // exactly to bound
247    // For now scaling off
248    model2.scaling(0);
249    ClpPrimalColumnSteepest steepest(5);
250    model2.setPrimalColumnPivotAlgorithm(steepest);
251    //model2.messageHandler()->setLogLevel(63);
252    model2.setFactorizationFrequency(maxFactor);
253    model2.setMaximumIterations(4000000);
254    double time1 = CoinCpuTime();
255    model2.primal();
256    // can't use values pass
257    model2.primal(0);
258    // test proper restart
259    if (argc>3) {
260      ClpDynamicMatrix * oldMatrix =
261        dynamic_cast< ClpDynamicMatrix*>(model2.clpMatrix());
262      assert (oldMatrix);
263      ClpDynamicMatrix * newMatrix = new 
264        ClpDynamicMatrix(&model3,numberGub,
265                                numberColumns,gubStart,
266                                lower,upper,
267                                start2,row2,element2,cost2,
268                                lowerColumn2,upperColumn2,
269                                oldMatrix->gubRowStatus(),oldMatrix->dynamicStatus());
270      model3.replaceMatrix(newMatrix);
271      // and ordinary status (but only NON gub rows)
272      memcpy(model3.statusArray(),model2.statusArray(),
273             (newMatrix->numberStaticRows()+model3.numberColumns())*sizeof(unsigned char));
274      newMatrix->switchOffCheck();
275      newMatrix->setRefreshFrequency(1000);
276    } else {
277      ClpDynamicExampleMatrix * oldMatrix =
278        dynamic_cast< ClpDynamicExampleMatrix*>(model2.clpMatrix());
279      assert (oldMatrix);
280      ClpDynamicExampleMatrix * newMatrix = new 
281        ClpDynamicExampleMatrix(&model3,numberGub,
282                                numberColumns,gubStart,
283                                lower,upper,
284                                start2,row2,element2,cost2,
285                                lowerColumn2,upperColumn2,
286                                oldMatrix->gubRowStatus(),oldMatrix->dynamicStatus(),
287                                oldMatrix->numberGubColumns(),oldMatrix->idGen());
288      model3.replaceMatrix(newMatrix);
289      // and ordinary status (but only NON gub rows)
290      memcpy(model3.statusArray(),model2.statusArray(),
291             (newMatrix->numberStaticRows()+model3.numberColumns())*sizeof(unsigned char));
292      newMatrix->switchOffCheck();
293      newMatrix->setRefreshFrequency(1000);
294    }
295    model3.setSpecialOptions(4); // exactly to bound
296    // For now scaling off
297    model3.scaling(0);
298    model3.setPrimalColumnPivotAlgorithm(steepest);
299    model3.messageHandler()->setLogLevel(63);
300    model3.setFactorizationFrequency(maxFactor);
301    model3.setMaximumIterations(4000000);
302    delete [] rowIsGub;
303    delete [] start2;
304    delete [] row2;
305    delete [] element2;
306    delete [] cost2;
307    delete [] lowerColumn2;
308    delete [] upperColumn2;
309    model3.primal();
310    // this code expects non gub first in original matrix
311    // and only works at present for ClpDynamicMatrix
312    ClpDynamicMatrix * gubMatrix =
313      dynamic_cast< ClpDynamicMatrix*>(model2.clpMatrix());
314    assert (gubMatrix);
315    ClpDynamicExampleMatrix * gubMatrix2 =
316      dynamic_cast< ClpDynamicExampleMatrix*>(model2.clpMatrix());
317    if (!gubMatrix2) {
318      const double * solution = model2.primalColumnSolution();
319      const double * cost =model.objective();
320      int numberGubColumns = gubMatrix->numberGubColumns();
321      int firstOdd = gubMatrix->firstDynamic();
322      int lastOdd =gubMatrix->firstAvailable();
323      int numberTotalColumns = firstOdd+numberGubColumns;
324      int originalNumberRows = model.numberRows();
325      int numberStaticRows = gubMatrix->numberStaticRows();
326      char * status = new char [numberTotalColumns];
327      double * gubSolution = new double [numberTotalColumns];
328      int numberSets = gubMatrix->numberSets();
329      const int * id = gubMatrix->id();
330      int i;
331      const float * columnLower = gubMatrix->columnLower();
332      const float * columnUpper = gubMatrix->columnUpper();
333      for (i=0;i<numberGubColumns;i++) {
334        if (gubMatrix->getDynamicStatus(i)==ClpDynamicMatrix::atUpperBound) {
335          gubSolution[i+firstOdd]=columnUpper[i];
336          status[i+firstOdd]=2;
337        } else if (gubMatrix->getDynamicStatus(i)==ClpDynamicMatrix::atLowerBound&&columnLower) {
338          gubSolution[i+firstOdd]=columnLower[i];
339          status[i+firstOdd]=1;
340        } else if (gubMatrix->getDynamicStatus(i)==ClpDynamicMatrix::soloKey) {
341          int iSet = gubMatrix->whichSet(i);
342          gubSolution[i+firstOdd]=gubMatrix->keyValue(iSet);
343          status[i+firstOdd]=0;
344        } else {
345          gubSolution[i+firstOdd]=0.0;
346          status[i+firstOdd]=1;
347        }
348      }
349      for (i=0;i<firstOdd;i++) {
350        ClpSimplex::Status thisStatus = model2.getStatus(i);
351        if (thisStatus==ClpSimplex::basic)
352          status[i]=0;
353        else if (thisStatus==ClpSimplex::atLowerBound)
354          status[i]=1;
355        else if (thisStatus==ClpSimplex::atUpperBound)
356          status[i]=2;
357        else if (thisStatus==ClpSimplex::isFixed)
358          status[i]=3;
359        else
360          abort();
361        gubSolution[i]=solution[i];
362      }
363      for (i=firstOdd;i<lastOdd;i++) {
364        int iBig = id[i-firstOdd]+firstOdd;
365        ClpSimplex::Status thisStatus = model2.getStatus(i);
366        if (thisStatus==ClpSimplex::basic)
367          status[iBig]=0;
368        else if (thisStatus==ClpSimplex::atLowerBound)
369          status[iBig]=1;
370        else if (thisStatus==ClpSimplex::atUpperBound)
371          status[iBig]=2;
372        else if (thisStatus==ClpSimplex::isFixed)
373          status[iBig]=3;
374        else
375          abort();
376        gubSolution[iBig]=solution[i];
377      }
378      char * rowStatus = new char[originalNumberRows];
379      for (i=0;i<numberStaticRows;i++) {
380        ClpSimplex::Status thisStatus = model2.getRowStatus(i);
381        if (thisStatus==ClpSimplex::basic)
382          rowStatus[i]=0;
383        else if (thisStatus==ClpSimplex::atLowerBound)
384          rowStatus[i]=1;
385        else if (thisStatus==ClpSimplex::atUpperBound)
386          rowStatus[i]=2;
387        else if (thisStatus==ClpSimplex::isFixed)
388          rowStatus[i]=3;
389        else
390          abort();
391      }
392      double objValue=0.0;
393      for (i=0;i<numberTotalColumns;i++) 
394        objValue += cost[i]*gubSolution[i];
395      printf("objective value is %g\n",objValue);
396      for (i=0;i<numberSets;i++) {
397        ClpSimplex::Status thisStatus = gubMatrix->getStatus(i);
398        if (thisStatus==ClpSimplex::basic)
399          rowStatus[numberStaticRows+i]=0;
400        else if (thisStatus==ClpSimplex::atLowerBound)
401          rowStatus[numberStaticRows+i]=1;
402        else if (thisStatus==ClpSimplex::atUpperBound)
403          rowStatus[numberStaticRows+i]=2;
404        else if (thisStatus==ClpSimplex::isFixed)
405          rowStatus[numberStaticRows+i]=3;
406        else
407          abort();
408      }
409      // Coding below may not work if gub rows not at end
410      FILE * fp=fopen ("xx.sol","w");
411      fwrite(gubSolution,sizeof(double),numberTotalColumns,fp);
412      fwrite(status,sizeof(char),numberTotalColumns,fp);
413      const double * rowsol = model2.primalRowSolution();
414      double * rowsol2 = new double[originalNumberRows];
415      memset(rowsol2,0,originalNumberRows*sizeof(double));
416      model.times(1.0,gubSolution,rowsol2);
417      for (i=0;i<numberStaticRows;i++) 
418        assert (fabs(rowsol[i]-rowsol2[i])<1.0e-3);
419      for (;i<originalNumberRows;i++) 
420        assert (rowsol2[i]>lower[i-numberStaticRows]-1.0e-3&&
421                rowsol2[i]<upper[i-numberStaticRows]+1.0e-3);
422      //for (;i<originalNumberRows;i++)
423      //printf("%d %g\n",i,rowsol2[i]);
424      fwrite(rowsol2,sizeof(double),originalNumberRows,fp);
425      delete [] rowsol2;
426      fwrite(rowStatus,sizeof(char),originalNumberRows,fp);
427      fclose(fp);
428      delete [] status;
429      delete [] gubSolution;
430      // ** if going to rstart as dynamic need id_
431      // also copy coding in useEf.. from ClpGubMatrix (i.e. test for basis)
432    }
433    printf("obj offset is %g\n",model2.objectiveOffset());
434    printf("Primal took %g seconds\n",CoinCpuTime()-time1);
435  }
436  delete [] mark;
437  delete [] gubStart;
438  delete [] gubEnd;
439  delete [] which;
440  delete [] whichGub;
441  delete [] lower;
442  delete [] upper;
443  return 0;
444}   
Note: See TracBrowser for help on using the repository browser.