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

Last change on this file since 1552 was 1552, checked in by mjs, 11 years ago

Format examples with 'astyle -A4 -p'.

  • Property svn:eol-style set to native
  • Property svn:keywords set to Id
File size: 7.4 KB
Line 
1/* $Id: testQP.cpp 1552 2010-05-25 01:03:59Z mjs $ */
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          int status = m.readMps("../../Data/Sample/share2qp", "mps");
25          if (status) {
26               printf("errors on input\n");
27               exit(77);
28          }
29          ClpInterior model;
30          model.loadProblem(*m.getMatrixByCol(), m.getColLower(), m.getColUpper(),
31                            m.getObjCoefficients(),
32                            m.getRowLower(), m.getRowUpper());
33          // get quadratic part
34          int * start = NULL;
35          int * column = NULL;
36          double * element = NULL;
37          m.readQuadraticMps(NULL, start, column, element, 2);
38          int j;
39          for (j = 0; j < 79; j++) {
40               if (start[j] < start[j+1]) {
41                    int i;
42                    printf("Column %d ", j);
43                    for (i = start[j]; i < start[j+1]; i++) {
44                         printf("( %d, %g) ", column[i], element[i]);
45                    }
46                    printf("\n");
47               }
48          }
49          model.loadQuadraticObjective(model.numberColumns(), start, column, element);
50          // share2qp is in old style qp format - convert to new so other options can use
51          model.writeMps("quad.mps");
52          ClpCholeskyBase * cholesky = new ClpCholeskyBase();
53          cholesky->setKKT(true);
54          model.setCholesky(cholesky);
55          model.primalDual();
56          double *primal;
57          double *dual;
58          primal = model.primalColumnSolution();
59          dual = model.dualRowSolution();
60          int i;
61          int numberColumns = model.numberColumns();
62          int numberRows = model.numberRows();
63          for ( i = 0; i < numberColumns; i++) {
64               if (fabs(primal[i]) > 1.0e-8)
65                    printf("%d primal %g\n", i, primal[i]);
66          }
67          for ( i = 0; i < numberRows; i++) {
68               if (fabs(dual[i]) > 1.0e-8)
69                    printf("%d dual %g\n", i, dual[i]);
70          }
71     } else {
72          // Could read into ClpInterior
73          ClpSimplex model;
74          if (model.readMps(argv[1])) {
75               printf("errors on input\n");
76               exit(77);
77          }
78          model.writeMps("quad");
79          if (argc < 3) {
80               // simplex - just primal as dual does not work
81               // also I need to fix scaling of duals on output
82               // (Was okay in first place - can't mix and match scaling techniques)
83               // model.scaling(0);
84               model.primal();
85          } else {
86               // barrier
87               ClpInterior barrier;
88               barrier.borrowModel(model);
89               ClpCholeskyBase * cholesky = new ClpCholeskyBase();
90               cholesky->setKKT(true);
91               barrier.setCholesky(cholesky);
92               barrier.primalDual();
93               barrier.returnModel(model);
94          }
95          // Just check if share2qp (quad.mps here)
96          // this is because I am not checking if variables at ub
97          if (model.numberColumns() == 79) {
98               double *primal;
99               double *dual;
100               primal = model.primalColumnSolution();
101               dual = model.dualRowSolution();
102               // Check duals by hand
103               const ClpQuadraticObjective * quadraticObj =
104                    (dynamic_cast<const ClpQuadraticObjective*>(model.objectiveAsObject()));
105               assert (quadraticObj);
106               CoinPackedMatrix * quad = quadraticObj->quadraticObjective();
107               const int * columnQuadratic = quad->getIndices();
108               const CoinBigIndex * columnQuadraticStart = quad->getVectorStarts();
109               const int * columnQuadraticLength = quad->getVectorLengths();
110               const double * quadraticElement = quad->getElements();
111               int numberColumns = model.numberColumns();
112               int numberRows = model.numberRows();
113               double * gradient = new double [numberColumns];
114               // move linear objective
115               memcpy(gradient, quadraticObj->linearObjective(), numberColumns * sizeof(double));
116               int iColumn;
117               for ( iColumn = 0; iColumn < numberColumns; iColumn++) {
118                    double valueI = primal[iColumn];
119                    CoinBigIndex j;
120                    for (j = columnQuadraticStart[iColumn];
121                              j < columnQuadraticStart[iColumn] + columnQuadraticLength[iColumn]; j++) {
122                         int jColumn = columnQuadratic[j];
123                         double valueJ = primal[jColumn];
124                         double elementValue = quadraticElement[j];
125                         if (iColumn != jColumn) {
126                              double gradientI = valueJ * elementValue;
127                              double gradientJ = valueI * elementValue;
128                              gradient[iColumn] += gradientI;
129                              gradient[jColumn] += gradientJ;
130                         } else {
131                              double gradientI = valueI * elementValue;
132                              gradient[iColumn] += gradientI;
133                         }
134                    }
135                    if (fabs(primal[iColumn]) > 1.0e-8)
136                         printf("%d primal %g\n", iColumn, primal[iColumn]);
137               }
138               for (int i = 0; i < numberRows; i++) {
139                    if (fabs(dual[i]) > 1.0e-8)
140                         printf("%d dual %g\n", i, dual[i]);
141               }
142               // Now use duals to get reduced costs
143               // Can't use this as will try and use scaling
144               // model.transposeTimes(-1.0,dual,gradient);
145               // So ...
146               CoinPackedMatrix * matrix = model.matrix();
147               const int * row = matrix->getIndices();
148               const CoinBigIndex * columnStart = matrix->getVectorStarts();
149               const int * columnLength = matrix->getVectorLengths();
150               const double * element = matrix->getElements();
151               for ( iColumn = 0; iColumn < numberColumns; iColumn++) {
152                    double dj = gradient[iColumn];
153                    CoinBigIndex j;
154                    for (j = columnStart[iColumn];
155                              j < columnStart[iColumn] + columnLength[iColumn]; j++) {
156                         int jRow = row[j];
157                         dj -= element[j] * dual[jRow];
158                    }
159                    if (model.getColumnStatus(iColumn) == ClpSimplex::basic) {
160                         assert (fabs(dj) < 1.0e-5);
161                    } else {
162                         assert (dj > -1.0e-5);
163                    }
164               }
165               delete [] gradient;
166          }
167     }
168     return 0;
169}
Note: See TracBrowser for help on using the repository browser.