source: stable/1.6/Clp/examples/findNetwork.cpp @ 1609

Last change on this file since 1609 was 1120, checked in by forrest, 13 years ago

add find network code

File size: 7.7 KB
Line 
1// Copyright (C) 2007, International Business Machines
2// Corporation and others.  All Rights Reserved.
3
4#include "ClpSimplex.hpp"
5#include "ClpPresolve.hpp"
6#include "CoinHelperFunctions.hpp"
7#include "ClpNetworkMatrix.hpp"
8#include "CoinTime.hpp"
9#include <iomanip>
10int main (int argc, const char *argv[])
11{
12  ClpSimplex  model;
13  int status;
14  // Keep names when reading an mps file
15  if (argc<2)
16    status=model.readMps("../../Data/Sample/p0033.mps",true);
17  else
18    status=model.readMps(argv[1],true);
19 
20  if (status) {
21    fprintf(stderr,"Bad readMps %s\n",argv[1]);
22    fprintf(stdout,"Bad readMps %s\n",argv[1]);
23    exit(1);
24  }
25  ClpSimplex * model2;
26  ClpPresolve pinfo;
27  int numberPasses=5; // can change this
28  /* Use a tolerance of 1.0e-8 for feasibility, treat problem as
29     not being integer, do "numberpasses" passes and throw away names
30     in presolved model */
31#define PRESOLVE
32#ifdef PRESOLVE
33  model2 = pinfo.presolvedModel(model,1.0e-8,false,numberPasses,false);
34#else
35  model2 = &model;
36#endif
37  if (!model2) {
38    fprintf(stderr,"ClpPresolve says %s is infeasible with tolerance of %g\n",
39            argv[1],1.0e-8);
40    fprintf(stdout,"ClpPresolve says %s is infeasible with tolerance of %g\n",
41            argv[1],1.0e-8);
42    // model was infeasible - maybe try again with looser tolerances
43    model2 = pinfo.presolvedModel(model,1.0e-7,false,numberPasses,false);
44    if (!model2) {
45      fprintf(stderr,"ClpPresolve says %s is infeasible with tolerance of %g\n",
46              argv[1],1.0e-7);
47      fprintf(stdout,"ClpPresolve says %s is infeasible with tolerance of %g\n",
48              argv[1],1.0e-7);
49      exit(2);
50    }
51  }
52  // change factorization frequency from 200
53  model2->setFactorizationFrequency(100+model2->numberRows()/50);
54  model2->setPerturbation(50);
55  double time1 = CoinCpuTime();
56  int numberRows = model2->numberRows();
57  int numberColumns = model2->numberColumns();
58  char * type = new char[numberRows];
59  if (false) {
60    // scale rows
61    const int * row = model2->matrix()->getIndices();
62    const CoinBigIndex * columnStart = model2->matrix()->getVectorStarts();
63    const int * columnLength = model2->matrix()->getVectorLengths(); 
64    double * elementByColumn = model2->matrix()->getMutableElements();
65    double * largest = new double [numberRows];
66    double * smallest = new double [numberRows];
67    int iRow,iColumn;
68    for (iColumn=0;iColumn<numberColumns;iColumn++) {
69      for (CoinBigIndex j=columnStart[iColumn];j<columnStart[iColumn]+columnLength[iColumn];j++) {
70        iRow = row[j];
71        double value = fabs(elementByColumn[j]);
72        largest[iRow] = CoinMax(largest[iRow],value);
73        smallest[iRow] = CoinMin(smallest[iRow],value);
74      }
75    }
76    int nScale=0;
77    for (iRow=0;iRow<numberRows;iRow++) {
78      if (largest[iRow]==smallest[iRow]&&largest[iRow]!=1.0) {
79        nScale++;
80        largest[iRow]= 1.0/largest[iRow];
81      } else {
82        largest[iRow]=0.0;
83      }
84    }
85    if (nScale) {
86      printf("Scaling %d rows\n",nScale);
87      for (iColumn=0;iColumn<numberColumns;iColumn++) {
88        for (CoinBigIndex j=columnStart[iColumn];j<columnStart[iColumn]+columnLength[iColumn];j++) {
89          iRow = row[j];
90          double value = largest[iRow];
91          if (value)
92            elementByColumn[j] *= value;
93          assert (fabs(elementByColumn[j])==1.0);
94        }
95      } 
96      printf("would have to do rhs as well\n");
97      abort();
98    }
99    delete [] largest;
100    delete [] smallest;
101  }
102  int nFind=model2->findNetwork(type,0.001);
103  printf("%d network rows found in %g seconds\n",nFind,CoinCpuTime()-time1);
104  if (CoinAbs(nFind)>0.2*numberRows) {
105    const int * row = model2->matrix()->getIndices();
106    const CoinBigIndex * columnStart = model2->matrix()->getVectorStarts();
107    const int * columnLength = model2->matrix()->getVectorLengths(); 
108    const double * elementByColumn = model2->matrix()->getElements();
109    int * whichRow = new int [numberRows];
110    int * whichColumn = new int [numberColumns];
111    int iRow,iColumn;
112    double * newElement = new double[numberColumns];
113    int * newColumn = new int[numberColumns];
114    int nEl=0;
115    int nSmallRows=0;
116    int nSmallColumns=0;
117    for (iRow=0;iRow<numberRows;iRow++) {
118      if (type[iRow]>=0)
119        whichRow[nSmallRows++]=iRow;
120    }
121    for (iColumn=0;iColumn<numberColumns;iColumn++) {
122      double product =1.0;
123      int n=0;
124      for (CoinBigIndex j=columnStart[iColumn];j<columnStart[iColumn]+columnLength[iColumn];j++) {
125        iRow = row[j];
126        if (type[iRow]>=0) {
127          n++;
128          assert (fabs(elementByColumn[j])==1.0);
129          product *= elementByColumn[j];
130          if (type[iRow])
131            product = - product;
132        }
133      }
134      assert (n<=2);
135      if (n==2) {
136        assert (product==-1.0);
137      } else if (n==1) {
138        newColumn[nEl]=nSmallColumns;
139        newElement[nEl++] = -product;
140      }
141      if (n)
142        whichColumn[nSmallColumns++]=iColumn;
143    }
144    // Create small problem
145    ClpSimplex small(model2,nSmallRows,whichRow,nSmallColumns,whichColumn);
146    printf("Small model has %d columns\n",nSmallColumns);
147    // change elements and rhs
148    double * lower = small.rowLower();
149    double * upper = small.rowUpper();
150    for (iRow=0;iRow<nSmallRows;iRow++) {
151      if (type[whichRow[iRow]]) {
152        double temp = -lower[iRow];
153        lower[iRow]=-upper[iRow];
154        upper[iRow]=temp;
155      }
156    }
157    row = small.matrix()->getIndices();
158    columnStart = small.matrix()->getVectorStarts();
159    columnLength = small.matrix()->getVectorLengths(); 
160    double * element = small.matrix()->getMutableElements();
161    for (iColumn=0;iColumn<nSmallColumns;iColumn++) {
162      for (CoinBigIndex j=columnStart[iColumn];j<columnStart[iColumn]+columnLength[iColumn];j++) {
163        iRow = row[j];
164        if (type[whichRow[iRow]]) 
165          element[j] = - element[j];
166      }
167    }
168    if (nEl) {
169      // add row (probably should fix network)
170      small.addRow(nEl,newColumn,newElement);
171    }
172    delete [] newElement;
173    delete [] newColumn;
174    ClpNetworkMatrix * network = new ClpNetworkMatrix(*small.matrix());
175    small.replaceMatrix(network,true);
176    small.initialSolve();
177    double * fullSolution = model2->primalColumnSolution();
178    // move solution back
179    const double * solution = small.primalColumnSolution();
180    for (iColumn=0;iColumn<nSmallColumns;iColumn++) {
181      int kColumn = whichColumn[iColumn];
182      model2->setColumnStatus(kColumn,small.getColumnStatus(iColumn));
183      fullSolution[kColumn]=solution[iColumn];
184    }
185    for (iRow=0;iRow<nSmallRows;iRow++) {
186      ClpSimplex::Status status=small.getRowStatus(iRow);
187      if (type[whichRow[iRow]]) {
188        if (status==ClpSimplex::atLowerBound)
189          status=ClpSimplex::atUpperBound;
190        else if (status==ClpSimplex::atUpperBound)
191          status=ClpSimplex::atLowerBound;
192      }
193      model2->setRowStatus(whichRow[iRow],status);
194    }
195    double * rowSolution = model2->primalRowSolution();
196    memset (rowSolution,0,numberRows*sizeof(double));
197    model2->times(1.0,fullSolution,rowSolution);
198    delete [] whichRow;
199    delete [] whichColumn;
200    model2->dual();
201  } else {
202    model2->initialSolve();
203  }
204  delete [] type;
205  int numberIterations=model2->numberIterations();;
206#ifdef PRESOLVE
207  pinfo.postsolve(true);
208  delete model2;
209#endif
210  /* After this postsolve model should be optimal.
211     We can use checkSolution and test feasibility */
212  model.checkSolution();
213  if (model.numberDualInfeasibilities()||
214      model.numberPrimalInfeasibilities()) 
215    printf("%g dual %g(%d) Primal %g(%d)\n",
216           model.objectiveValue(),
217           model.sumDualInfeasibilities(),
218           model.numberDualInfeasibilities(),
219           model.sumPrimalInfeasibilities(),
220           model.numberPrimalInfeasibilities());
221  // But resolve for safety
222  model.primal(1);
223
224  numberIterations += model.numberIterations();;
225  // for running timing tests
226  std::cout<<argv[1]<<" Objective "<<model.objectiveValue()<<" took "<<
227    numberIterations<<" iterations and "<<
228    CoinCpuTime()-time1<<" seconds"<<std::endl;
229  return 0;
230}   
Note: See TracBrowser for help on using the repository browser.