source: trunk/Clp/examples/testQP.cpp @ 1559

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

merge split branch into trunk

  • Property svn:eol-style set to native
  • Property svn:keywords set to Id
File size: 7.6 KB
Line 
1/* $Id: testQP.cpp 1559 2010-06-05 19:42:36Z stefan $ */
2// Copyright (C) 2005, International Business Machines
3// Corporation and others.  All Rights Reserved.
4
5#include "CoinMpsIO.hpp"
6#include "ClpInterior.hpp"
7#include "ClpSimplex.hpp"
8#include "ClpCholeskyBase.hpp"
9#include "ClpQuadraticObjective.hpp"
10#include <cassert>
11int main(int argc, const char *argv[])
12{
13     /* Read quadratic model in two stages to test loadQuadraticObjective.
14
15        And is also possible to just read into ClpSimplex/Interior which sets it all up in one go.
16        But this is only if it is in QUADOBJ format.
17
18        If no arguments does share2qp using ClpInterior (also creates quad.mps which is in QUADOBJ format)
19        If one argument uses simplex e.g. testit quad.mps
20        If > one uses barrier via ClpSimplex input and then ClpInterior borrow
21     */
22     if (argc < 2) {
23          CoinMpsIO  m;
24#if defined(COIN_HAS_SAMPLE) && defined(SAMPLEDIR)
25          int status = m.readMps(SAMPLEDIR "/share2qp", "mps");
26#else
27          fprintf(stderr, "Do not know where to find sample MPS files.\n");
28          exit(1);
29#endif
30          if (status) {
31               printf("errors on input\n");
32               exit(77);
33          }
34          ClpInterior model;
35          model.loadProblem(*m.getMatrixByCol(), m.getColLower(), m.getColUpper(),
36                            m.getObjCoefficients(),
37                            m.getRowLower(), m.getRowUpper());
38          // get quadratic part
39          int * start = NULL;
40          int * column = NULL;
41          double * element = NULL;
42          m.readQuadraticMps(NULL, start, column, element, 2);
43          int j;
44          for (j = 0; j < 79; j++) {
45               if (start[j] < start[j+1]) {
46                    int i;
47                    printf("Column %d ", j);
48                    for (i = start[j]; i < start[j+1]; i++) {
49                         printf("( %d, %g) ", column[i], element[i]);
50                    }
51                    printf("\n");
52               }
53          }
54          model.loadQuadraticObjective(model.numberColumns(), start, column, element);
55          // share2qp is in old style qp format - convert to new so other options can use
56          model.writeMps("quad.mps");
57          ClpCholeskyBase * cholesky = new ClpCholeskyBase();
58          cholesky->setKKT(true);
59          model.setCholesky(cholesky);
60          model.primalDual();
61          double *primal;
62          double *dual;
63          primal = model.primalColumnSolution();
64          dual = model.dualRowSolution();
65          int i;
66          int numberColumns = model.numberColumns();
67          int numberRows = model.numberRows();
68          for (i = 0; i < numberColumns; i++) {
69               if (fabs(primal[i]) > 1.0e-8)
70                    printf("%d primal %g\n", i, primal[i]);
71          }
72          for (i = 0; i < numberRows; i++) {
73               if (fabs(dual[i]) > 1.0e-8)
74                    printf("%d dual %g\n", i, dual[i]);
75          }
76     } else {
77          // Could read into ClpInterior
78          ClpSimplex model;
79          if (model.readMps(argv[1])) {
80               printf("errors on input\n");
81               exit(77);
82          }
83          model.writeMps("quad");
84          if (argc < 3) {
85               // simplex - just primal as dual does not work
86               // also I need to fix scaling of duals on output
87               // (Was okay in first place - can't mix and match scaling techniques)
88               // model.scaling(0);
89               model.primal();
90          } else {
91               // barrier
92               ClpInterior barrier;
93               barrier.borrowModel(model);
94               ClpCholeskyBase * cholesky = new ClpCholeskyBase();
95               cholesky->setKKT(true);
96               barrier.setCholesky(cholesky);
97               barrier.primalDual();
98               barrier.returnModel(model);
99          }
100          // Just check if share2qp (quad.mps here)
101          // this is because I am not checking if variables at ub
102          if (model.numberColumns() == 79) {
103               double *primal;
104               double *dual;
105               primal = model.primalColumnSolution();
106               dual = model.dualRowSolution();
107               // Check duals by hand
108               const ClpQuadraticObjective * quadraticObj =
109                    (dynamic_cast<const ClpQuadraticObjective*>(model.objectiveAsObject()));
110               assert(quadraticObj);
111               CoinPackedMatrix * quad = quadraticObj->quadraticObjective();
112               const int * columnQuadratic = quad->getIndices();
113               const CoinBigIndex * columnQuadraticStart = quad->getVectorStarts();
114               const int * columnQuadraticLength = quad->getVectorLengths();
115               const double * quadraticElement = quad->getElements();
116               int numberColumns = model.numberColumns();
117               int numberRows = model.numberRows();
118               double * gradient = new double [numberColumns];
119               // move linear objective
120               memcpy(gradient, quadraticObj->linearObjective(), numberColumns * sizeof(double));
121               int iColumn;
122               for (iColumn = 0; iColumn < numberColumns; iColumn++) {
123                    double valueI = primal[iColumn];
124                    CoinBigIndex j;
125                    for (j = columnQuadraticStart[iColumn];
126                              j < columnQuadraticStart[iColumn] + columnQuadraticLength[iColumn]; j++) {
127                         int jColumn = columnQuadratic[j];
128                         double valueJ = primal[jColumn];
129                         double elementValue = quadraticElement[j];
130                         if (iColumn != jColumn) {
131                              double gradientI = valueJ * elementValue;
132                              double gradientJ = valueI * elementValue;
133                              gradient[iColumn] += gradientI;
134                              gradient[jColumn] += gradientJ;
135                         } else {
136                              double gradientI = valueI * elementValue;
137                              gradient[iColumn] += gradientI;
138                         }
139                    }
140                    if (fabs(primal[iColumn]) > 1.0e-8)
141                         printf("%d primal %g\n", iColumn, primal[iColumn]);
142               }
143               for (int i = 0; i < numberRows; i++) {
144                    if (fabs(dual[i]) > 1.0e-8)
145                         printf("%d dual %g\n", i, dual[i]);
146               }
147               // Now use duals to get reduced costs
148               // Can't use this as will try and use scaling
149               // model.transposeTimes(-1.0,dual,gradient);
150               // So ...
151               CoinPackedMatrix * matrix = model.matrix();
152               const int * row = matrix->getIndices();
153               const CoinBigIndex * columnStart = matrix->getVectorStarts();
154               const int * columnLength = matrix->getVectorLengths();
155               const double * element = matrix->getElements();
156               for (iColumn = 0; iColumn < numberColumns; iColumn++) {
157                    double dj = gradient[iColumn];
158                    CoinBigIndex j;
159                    for (j = columnStart[iColumn];
160                              j < columnStart[iColumn] + columnLength[iColumn]; j++) {
161                         int jRow = row[j];
162                         dj -= element[j] * dual[jRow];
163                    }
164                    if (model.getColumnStatus(iColumn) == ClpSimplex::basic) {
165                         assert(fabs(dj) < 1.0e-5);
166                    } else {
167                         assert(dj > -1.0e-5);
168                    }
169               }
170               delete [] gradient;
171          }
172     }
173     return 0;
174}
Note: See TracBrowser for help on using the repository browser.